Regarding addInteractionGroup() in CustomNonBodedForce

[Robin Singh]

Hello OpenMM Developers,

I am working on salt-water system and want to customize the interaction parameters of the ions of the salt using addInteractionGroups() in CustomNonbondedForce. For this, first I am looking to find energy for CustomNonBondedForce when every atom is interacting with every other atom in it and then in the second case I am defining interaction Groups. The codes and the results are given below:

The code for the first case is:

Code: Select all

from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout

pdb = PDBFile("liquid.pdb")
forcefield = ForceField('tip3p.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1.0*nanometer, constraints=HBonds)
for i, f in enumerate(system.getForces()):
    f.setForceGroup(i)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
state = simulation.context.getState(getEnergy=True)
print(state.getPotentialEnergy())
for i, f in enumerate(system.getForces()):
    state = simulation.context.getState(getEnergy=True, groups={i})
    print(f.getName(), state.getPotentialEnergy())
    

And the zero point energy (CustomNonBondedForce) I got is:

CustomNonbondedForce 6301.755701363087 kJ/mol


The code for the second case is:

Code: Select all

from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout

pdb = PDBFile("liquid.pdb")
forcefield = ForceField('tip3p.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1.0*nanometer, constraints=HBonds)
for i, f in enumerate(system.getForces()):
    f.setForceGroup(i)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)

custom_nb_force = system.getForce(2)
custom_nb_force.setUseSwitchingFunction(True)
custom_nb_force.setSwitchingDistance(2.0*nanometer)

all_indices = [int(i.index) for i in pdb.topology.atoms()]
solute = [int(i.index) for i in pdb.topology.atoms() if not (i.residue.name in ['HOH'])]

custom_nb_force.addInteractionGroup(solute, all_indices)

for i, f in enumerate(system.getForces()):
    state = simulation.context.getState(getEnergy=True, groups={i})
    print(f.getName(), state.getPotentialEnergy())
state = simulation.context.getState(getEnergy=True)
print(state.getPotentialEnergy())

And the zero point energy (CustomNonBondedForce) I got is:

CustomNonbondedForce 6301.755701363087 kJ/mol

In the second case, I am using addInteractionGroup() and giving two groups as mentioned above but still the energy due to CustomNonBondedForce is same as the default case. It should not happen like that and the energies due to CustomNonBondedForce in both of the cases should be different.

Can you please have a look at the code and let me know where is the problem?

Regards,
Robin Singh

[Peter Eastman]

Modifying a System or the Forces it contains has no effect on any Simulations that already exist. You need to make the changes before you create the Simulation. Alternatively you can call simulation.context.reinitialize().