AMOEBA small molecule topology and energy decomposition

[Madhuranga Rathnayake]

Dear all,

I extracted DMF parameters from tinker AMOEBA09 parameter file and trying to run it in OpenMM. Initially, it looks okay but when I looked at energy decomposition, some energy terms are zero. i.e AmoebaAngleForce.

Code: Select all

from simtk.openmm import app
import simtk.openmm as mm
from simtk import unit as u
import parmed as pmd

input_pdb = 'dmf.pdb'

report_freq  = 10
sim_steps    = 500
timestep     = 1.0 * u.femtoseconds
pressure     = 1 * u.atmosphere
temperature  = 300 * u.kelvin
forcefield   = app.ForceField('amoeba2009.xml', 'dmf.xml')

# add PDB file and load topology
pdb = app.PDBFile(input_pdb)
topologyFile = 'dmf-residue.xml'
pdb.topology.loadBondDefinitions(topologyFile)
pdb = app.PDBFile(input_pdb) 

print('Number of Residues:', pdb.topology.getNumResidues())
print('Number of Atoms:', pdb.topology.getNumAtoms())
print('Number of bonds in Topology:', pdb.topology.getNumBonds())

# Number of Residues: 500
# Number of Atoms: 6000
# Number of bonds in Topology: 5500

# Set platform
platform = mm.Platform.getPlatformByName('CUDA')
properties = {'CudaPrecision': 'mixed'}
# Create system for simulation
system = forcefield.createSystem(pdb.topology, 
    polarization='mutual',
    mutualInducedTargetEpsilon=0.00001, 
    nonbondedMethod=app.PME,
    nonbondedCutoff=1.2 * u.nanometers, 
    vdwCutoff=1.2 * u.nanometer,
    useDispersionCorrection=True, 
    rigidWater=False, 
    constraints=None,
    ewaldErrorTolerance=0.0005)
    
# Integrator
integrator = mm.VerletIntegrator(timestep)
integrator.setConstraintTolerance(0.00001)

# Thermostat
thermostat = mm.AndersenThermostat(temperature, 0.1 / u.picoseconds)
system.addForce(thermostat)

# Barostat
barostat = mm.MonteCarloBarostat(pressure, temperature)
system.addForce(barostat)

# Set simulation
simulation = app.Simulation(pdb.topology, system, integrator, platform, properties)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()

simulation.step(100)

for i in range(system.getNumForces()): 
    force = system.getForce(i) 
    force.setForceGroup(i) 
    print(i, force.__class__.__name__, simulation.context.getState(getEnergy=True, groups={i}).getPotentialEnergy())
    
0 AmoebaBondForce -476.64913181999026 kJ/mol
1 AmoebaOutOfPlaneBendForce 0.0 kJ/mol
2 AmoebaAngleForce 0.0 kJ/mol
3 AmoebaInPlaneAngleForce 0.0 kJ/mol
4 AmoebaStretchBendForce 0.0 kJ/mol
5 PeriodicTorsionForce 0.0 kJ/mol
6 AmoebaPiTorsionForce 0.0 kJ/mol
7 AmoebaTorsionTorsionForce 0.0 kJ/mol
8 AmoebaVdwForce -3142.0272416332896 kJ/mol
9 AmoebaMultipoleForce -33538.827972081395 kJ/mol
10 HarmonicBondForce 0.0 kJ/mol
11 CMMotionRemover 0.0 kJ/mol
12 AndersenThermostat 0.0 kJ/mol
13 MonteCarloBarostat 0.0 kJ/mol

but if I save the current state as a PDB and do an energy decomposition in ParmEd, it gives energy components correctly. which means parameter file and topology is correct.

Code: Select all

state = simulation.context.getState(getPositions=True)
with open('current_state.pdb', 'w') as output:
  app.PDBFile.writeFile(simulation.topology, state.getPositions(), output)
  
struct=pmd.load_file('current_state.pdb')
ecomps=(pmd.openmm.energy_decomposition_system(struct,system,nrg=u.kilojoules_per_mole))

[('AmoebaBondForce', 839.722412109375),
 ('AmoebaOutOfPlaneBendForce', 27.751140594482422),
 ('AmoebaAngleForce', 1355.9627685546875),
 ('AmoebaInPlaneAngleForce', 1121.425537109375),
 ('AmoebaStretchBendForce', -15.018678665161133),
 ('PeriodicTorsionForce', -3842.664794921875),
 ('AmoebaPiTorsionForce', 37.94041442871094),
 ('AmoebaTorsionTorsionForce', 0.0),
 ('AmoebaVdwForce', -3166.4592501336283),
 ('AmoebaMultipoleForce', -33350.7421875),
 ('HarmonicBondForce', 0.0),
 ('CMMotionRemover', 0.0),
 ('AndersenThermostat', 0.0),
 ('MonteCarloBarostat', 0.0)]

Any thoughts why OpenMM doesn't pick up the rest of energy terms?

Thank you in advance!

-Ranga

[Peter Eastman]

Code: Select all

for i in range(system.getNumForces()): 
    force = system.getForce(i) 
    force.setForceGroup(i) 
    print(i, force.__class__.__name__, simulation.context.getState(getEnergy=True, groups={i}).getPotentialEnergy())

That code doesn't do what you intend. Any changes to a System or the Forces it contains, such as setting force groups, must be done before you create the Simulation. Once you create it, it won't see any later changes you make to the System. So you need to split this into two parts. First, before you create the Simulation, assign force groups:

Code: Select all

for i,f in enumerate(system.getForces()):
  f.setForceGroup(i)

And then you can loop over them later to query the energies of the groups as you're doing now.

[Madhuranga Rathnayake]

Thank you! appreciate it. solved the problem.