Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- from __future__ import print_function
- from simtk.openmm import app
- import simtk.openmm as mm
- from simtk import unit
- from sys import stdout
- import mbpol
- simulation_name = "w2"
- pdb = app.PDBFile('w2.pdb')
- temperature = float(190)*unit.kelvin
- forcefield = mm.app.ForceField("/opt/conda/miniconda2/envs/openmm_env/lib/python3.6/site-packages/mbpol-1.0-py3.6-linux-x86_64.egg/mbpol.xml")
- nonbonded = mm.app.CutoffNonPeriodic
- timestep = 0.2*unit.femtoseconds
- equilib_steps = 200
- equilib_steps1 = 5
- ################################## NVT
- system = forcefield.createSystem(pdb.topology, nonbondedMethod=nonbonded,
- nonbondedCutoff=3.0*unit.nanometers)
- from openmmtools.integrators import GeodesicBAOABIntegrator
- integrator = GeodesicBAOABIntegrator(temperature, 1.0/unit.picoseconds, timestep)
- integrator.setRandomNumberSeed(1)
- platform = mm.Platform.getPlatformByName('Reference')
- simulation = app.Simulation(pdb.topology, system, integrator, platform)
- simulation.context.setPositions(pdb.positions)
- simulation.context.computeVirtualSites()
- simulation.context.setVelocitiesToTemperature(temperature)
- simulation.reporters.append(app.PDBReporter('traj_nvt.pdb', 1))
- simulation.reporters.append(app.StateDataReporter(simulation_name + "_nvt.log", 1, step=True,
- time=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True,
- progress=True, remainingTime=True, speed=True, totalSteps=equilib_steps,
- separator=' '))
- print('Equilibrating...')
- simulation.step(equilib_steps)
- simulation.saveState('nvt.xml')
- simulation.step(equilib_steps1)
- simulation.saveState('nvt1.xml')
- simulation.step(equilib_steps1)
- simulation.saveState('nvt2.xml')
- simulation.step(equilib_steps1)
- simulation.saveState('nvt3.xml')
Add Comment
Please, Sign In to add comment