simtk.org

Code: Select all

from openff.toolkit.topology import Molecule
from openmmforcefields.generators import GAFFTemplateGenerator

from openmm import *
from openmm import XmlSerializer, LangevinMiddleIntegrator
from openmm.unit import *
from openmm.app import *    
from sys import stdout
import numpy as np

molecule = Molecule.from_smiles('C=O')
# Create the GAFF template generator
gaff = GAFFTemplateGenerator(molecules=molecule)
# Create an OpenMM ForceField object with AMBER ff14SB and TIP3P with compatible ions
forcefield = ForceField('amber/protein.ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml')
# Register the GAFF template generator
forcefield.registerTemplateGenerator(gaff.generator)

pdb_ch2o = PDBFile('ch2o.pdb')

modeller = Modeller(pdb_ch2o.topology, pdb_ch2o.positions)
system = forcefield.createSystem(modeller.topology, nonbondedCutoff=1.0*nanometer, constraints=HBonds)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.001*picoseconds)
simulation = Simulation(modeller.topology, system, integrator)
simulation.context.setPositions(modeller.positions)

# add two additional particles with zero mass, so they will be immovable
id1 = system.addParticle(0.0) 
id2 = system.addParticle(0.0) 

# try to set new positions
origpos = simulation.context.getState(getPositions=True).getPositions().copy()
origpos[4] = Vec3(10,10,10) * nanometers
simulation.context.setPositions(origpos)

# check positions, the positions are not updated in simulation.context
print(origpos[4])
print(simulation.context.getState(getPositions=True).getPositions()[4])

Code: Select all

Vec3(x=10, y=10, z=10) nm
Vec3(x=0.0, y=0.0, z=0.0) nm