Skip to content

Instantly share code, notes, and snippets.

@amrhamedp
Created December 1, 2021 22:49
Show Gist options
  • Select an option

  • Save amrhamedp/187e6aceedf67d3e8be783763f20a983 to your computer and use it in GitHub Desktop.

Select an option

Save amrhamedp/187e6aceedf67d3e8be783763f20a983 to your computer and use it in GitHub Desktop.
from MDAnalysis.coordinates.DCD import DCDFile
dms = app.DesmondDMSFile('a.dms')
system = dms.createSystem(nonbondedMethod=app.PME,
nonbondedCutoff=9.0*unit.angstroms,OPLS=True)
import simtk.openmm as openmm
system.addForce(openmm.MonteCarloBarostat(1 * units.bar, 300*unit.kelvin))
integrator = mm.LangevinIntegrator(300*unit.kelvin, 1.0/unit.picoseconds,
2.0*unit.femtoseconds)
integrator.setConstraintTolerance(0.00001)
platform = mm.Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'single'}
simulation = app.Simulation(dms.topology, system, integrator, platform,
properties)
import mdtraj as md
a=md.load('a.dcd', top='a.pdb')
index1=0
with DCDFile("a.dcd") as dcd:
header = dcd.header
for frame in dcd.readframes()[0]:
simulation.context.setPositions(frame * unit.angstrom)
simulation.context.setPeriodicBoxVectors(a=a.unitcell_vectors[index1][0],b=a.unitcell_vectors[index1][1],c=a.unitcell_vectors[index1][2])
state = simulation.context.getState(getForces=True,getEnergy=True)
print(state.getPeriodicBoxVectors())
print(state.getPotentialEnergy()/unit.kilocalories_per_mole)
index1=index1+1
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment