Simple OpenMM system for testing code development

OpenMM
Author

Quyen Vu

Published

February 9, 2023

This post contains a simple system definition of two beads with custom nonbonded for testing


import openmm as mm
import openmm.app as app
import numpy as np

print(mm.__version__)
sys = mm.System()
# 2 GLY
sys.addParticle(57.05*mm.unit.amu)
sys.addParticle(57.05*mm.unit.amu)

# build topology
top=mm.app.topology.Topology()
top.addChain('prot')
prot_chain = list(top.chains())[-1]
top.addResidue('GLY', prot_chain)
res1 = list(top.residues())[-1]
top.addAtom('CA', app.Element.getBySymbol('C'), res1)

top.addResidue('GLY', prot_chain)
res2 = list(top.residues())[-1]
top.addAtom('CA', app.Element.getBySymbol('C'), res2)

# # first implementation
# energy_function = '((rc/r)^(4)-1)^(2)'
# wang_frenkel_force = mm.CustomNonbondedForce(energy_function)
# wang_frenkel_force.addGlobalParameter('rc', 1.408533*mm.unit.nanometer)

# second implementation
energy_function = 'step(rc-r)*((rc/r)^(2*mu)-1)^(2*nu);'
wang_frenkel_force = mm.CustomNonbondedForce(energy_function)
wang_frenkel_force.addGlobalParameter('rc', 1.408533*mm.unit.nanometer)
wang_frenkel_force.addGlobalParameter('mu', 2)
wang_frenkel_force.addGlobalParameter('nu', 1)

# # # third implementation
# energy_function = '((rc/r)^(4)-1)^(2*nu);'
# wang_frenkel_force = mm.CustomNonbondedForce(energy_function)
# wang_frenkel_force.addGlobalParameter('rc', 1.408533*mm.unit.nanometer)
# # wang_frenkel_force.addGlobalParameter('mu', 2)
# wang_frenkel_force.addGlobalParameter('nu', 1)

for _ in top.atoms():
    wang_frenkel_force.addParticle(())

sys.addForce(wang_frenkel_force)

integrator = mm.LangevinIntegrator(1*mm.unit.kelvin, 1/mm.unit.picosecond, 0.000001*mm.unit.picosecond)

# # CPU
# platform = mm.Platform.getPlatformByName('CPU')
# properties = {'Threads': str(8)}

# GPU
platform = mm.Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'single', "DeviceIndex": "0"}
sim = mm.app.Simulation(top, sys, integrator, platform, properties)

# pos=np.array([[0,0,0],[1.409,0,0]])
# pos=np.array([[0,0,0],[1.4,0,0]])
pos=np.array([[0,0,0],[1.9,0,0]])

sim.context.setPositions(pos)
state = sim.context.getState(getForces=True, getEnergy=True)
print(platform)
print(properties)
print(state.getPotentialEnergy())
print(state.getForces())