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())
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())
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