#
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.
#

##########################################################################
# this script was generated by openmm-builder. to customize it further,
# you can save the file to disk and edit it with your favorite editor.
##########################################################################

from __future__ import print_function
from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit
from sys import stdout

print("[START] Application is now running")

pdb = app.PDBFile('input.pdb')
print("[STATUS] Loaded model")
forcefield = app.ForceField('amber03.xml', 'amber03_obc.xml')
print("[STATUS] Loaded force field")

system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.NoCutoff, 
     constraints=None, rigidWater=False)
print("[STATUS] Created system")
integrator = mm.LangevinIntegrator(300*unit.kelvin, 91/unit.picoseconds, 
    1.0*unit.femtoseconds)
print("[STATUS] Created integrator")

try:
    platform = mm.Platform.getPlatformByName('CPU')
except Exception as e:
    print("[ERROR] Could not load platform CPU. Running Reference")
    platform = mm.Platform.getPlatformByName("Reference")

simulation = app.Simulation(pdb.topology, system, integrator, platform)
print("[STATUS] Set up compute platform")
simulation.context.setPositions(pdb.positions)
print("[STATUS] Set atomic positions")

print('[STATUS] Minimizing...')
simulation.minimizeEnergy()
print('[STATUS] Equilibrating...')
simulation.step(100)

simulation.reporters.append(app.DCDReporter('trajectory.dcd', 1000))
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True, 
    potentialEnergy=True, totalEnergy=True, temperature=True, separator='\t'))
print("[STATUS] Set up reporters")

print('[STATUS] Running Production...')

increment = 1000

for i in range(0,100000,increment): 
    print("[STATUS] Step %s" % (i))
    simulation.step(increment)

print('[END] Done!')