CustomExternalForce - Odd Position behavior

[Dennis Della Corte]

Hello,

I am adding position restraints via a CustomExternalForce according to methods described here: https://onlinelibrary.wiley.com/doi/ful ... prot.25759

Unfortunately, I seem to run into problems with PBCs as discussed elsewhere on the forum/issue tracker.
The error I get suggests a blow up after ~1000 MD steps: ValueError: Particle position is NaN
Inspection of the trajectory shows very large movement in 0.5ps steps.


Adding a flat bottom potential to CA is done like this:
pdb = PDBFile('./test.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')

modeller = Modeller(pdb.topology, pdb.positions)
modeller.addSolvent(forcefield, padding=1.0*nanometers)

system = forcefield.createSystem(modeller.topology, nonbondedMethod=PME,
nonbondedCutoff=1*nanometer, constraints=HBonds)

restraint = CustomExternalForce("k0*(max(d-d0, 0.0))^2 ; d=sqrt((x-x0)^2+(y-y0)^2+(z-z0)^2)")
restraint.addGlobalParameter("d0", 5*angstroms)
restraint.addGlobalParameter("k0", 0.025 * 12 * kilocalories_per_mole/angstroms**2) #mass weighted, 12 for Carbon!!
restraint.addPerParticleParameter("x0")
restraint.addPerParticleParameter("y0")
restraint.addPerParticleParameter("z0")

names = [ i.name for i in modeller.topology.atoms() ]

for i in range(system.getNumParticles()):
if names == 'CA':
restraint.addParticle(i,modeller.positions)

system.addForce(restraint)

#SOME STUFF OMITTED HERE

integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(modeller.topology, system, integrator, platform, properties)
simulation.context.setPositions(modeller.positions)
simulation.reporters.append(PDBReporter('minimzation.pdb', 1))

###here I need some logic to check what happened to my coordinates
state = simulation.context.getState(getPositions=True)
current_positions = state.getPositions()

for i in range(10):
print(restraint.getParticleParameters(i))
idx = restraint.getParticleParameters(i)[0]
print(modeller.positions[idx])
print(current_positions[idx])

simulation.minimizeEnergy()
--------------------------------
At this stage my output reads for i == 1:
[1, (0.1458, 0.0, 0.0)]
Vec3(x=0.1458, y=0.0, z=0.0) nm
Vec3(x=0.14580016727447465, y=0.0, z=-4.2877197259372224e-07) nm

If I now look into the first .pdb frame as reported during MD I get the following position for the first CA:
ATOM 2 CA SER A 1 81.924 0.239 81.142 1.00 0.00 C

The EM finishes successfully, the error occurs about 1000 steps into the MD.

I suspect the change in position to have something to do with PBC, but not sure how to remedy it.
Is it even possible to use the above definition of CustomExternalForce with PBC?

Any help is greatly appreciated. Thank you,
Dennis

[Peter Eastman]

Yes, periodic boundary conditions will definitely cause problems. To fix it, change

Code: Select all

d=sqrt((x-x0)^2+(y-y0)^2+(z-z0)^2)

to

Code: Select all

d=periodicdistance(x, y, z, x0, y0, z0)

You can find more information about this in the docs at http://docs.openmm.org/latest/api-c++/g ... rnalForceE.