Propagation of zero in step function

[Charles Brooks]

Dear Colleagues,

I am using CustomCompoundBondedForce to generate some custom restraints associated with restraining particles, or the center of mass of a collection of particles, within a sphere, cylinder or a plane. I am running into a problem for cases where a particle (or the COM of a collection of particles) are sitting along the axes that define the direction of the cylinder (or plane). I attach a simple (single particle) test case below that illustrates the issue. I note that the energy here is fine, i.e. I get zero when the position of the particles is along one of the cylinder axes, but the forces are not ok. I believe they should be zero.

Any insights into how I enforce the forces to be zero in this situation too would be appreciated.

The definition of the CostomCompoundBondedForce is:
<Force energy="mult*step(dlta)*0.5*frc*dlta*dlta;dlta=r-droff;r=sqrt(r);mult=step(abs(r)-small);r=delx*delx+dely*dely+delz*delz;delz=delz-r*zdir;dely=dely-r*ydir;delx=delx-r*xdir;r=delx*xdir+dely*ydir+delz*zdir;delz=a*z1-zref;dely=a*y1-yref;delx=a*x1-xref;" forceGroup="6" particles="1" type="CustomCompoundBondForce" version="1">
<PerBondParameters>
<Parameter name="a"/>
</PerBondParameters>
<GlobalParameters>
<Parameter default="4184.000015258789" name="frc"/>
<Parameter default=".52" name="droff"/>
<Parameter default="0" name="xref"/>
<Parameter default="0" name="yref"/>
<Parameter default="0" name="zref"/>
<Parameter default="1" name="xdir"/>
<Parameter default="0" name="ydir"/>
<Parameter default="0" name="zdir"/>
<Parameter default="9.999999747378752e-07" name="small"/>
</GlobalParameters>
<Bonds>
<Bond p1="0" param1="1"/>
</Bonds>
<Functions/>
</Force>

and the files to run this test are in the tgz file attached.

Thanks,

Charles Brooks

[Peter Eastman]

There's something wrong with your energy expression. It first defines r in terms of delx, dely, and delz. And then it defines delx, dely, and delz in terms of r! Also, delx, dely, and delz each appears in its own definition, which certainly isn't allowed. Think of the expression as a set of mathematical definitions, not as procedural code where you can assign a value to a variable and then change it.

The parser should have thrown an exception, so I'll look into why it didn't do that.

Peter

[Charles Brooks]

Something is still odd, because when I replace the expression with:

<Force energy="mult*step(dlta)*0.5*frc*dlta*dlta;dlta=r-droff;r=sqrt(r2);mult=step(abs(r2)-small);r2=delax*delax+delay*delay+delaz*delaz;delaz=delz-r1*zdir;delay=dely-r1*ydir;delax=delx-r1*xdir;r1=delx*xdir+dely*ydir+delz*zdir;delz=a*z1-zref;dely=a*y1-yref;delx=a*x1-xref;" forceGroup="6" particles="1" type="CustomCompoundBondForce" version="1">

it still produces non-zero force components.

[Charles Brooks]

Dear Peter,

I have been playing further with the problem I raised a day or so ago, and I am convinced there is an issue in the function interpreter. I append a simple model, 1 particle, only customcompoundbonded force defined, pdb coordinates and a simple python script that runs this. If you have insights into why the forces are not coming out 0, which they should be, I'd appreciate it. I note the potential is of the form:

dr_vec = r_vec - r_vec(ref)
d'r_vec = dr_vec - dr_vec*d_vec(dir)
s = sqrt(d'r_vec*d'r_vec)
energy = 1/2xKxStep(s-droff)x(s-droff)^2

where _vec denotes a 3D vector, x denotes scalar multiplication, * denotes vector dot product, K -> force constant

===========================system xml===================================
<?xml version="1.0" ?>
<System type="System" version="1">
<PeriodicBoxVectors>
<A x="2" y="0" z="0"/>
<B x="0" y="2" z="0"/>
<C x="0" y="0" z="2"/>
</PeriodicBoxVectors>
<Particles>
<Particle mass="40"/>
</Particles>
<Constraints/>
<Forces>
<Force alpha="0" cutoff="99.9" dispersionCorrection="0" ewaldTolerance=".0005" forceGroup="0" method="0" nx="0" ny="0" nz="0" recipForceGroup="-1" rfDielectric="1" switchingDistance="-1" type="NonbondedForce" useSwitchingFunction="0" version="1">
<Particles>
<Particle eps=".9977288609186554" q="0" sig=".33999999442653656"/>
</Particles>
<Exceptions/>
</Force>
<Force energy="step(r-droff)*0.5*frc*(r-droff)^2;r=sqrt(r);r=delx*delx+dely*dely+delz*delz;delz=dz-s*zdir;dely=dy-s*ydir;delx=dx-s*xdir;s=dx*xdir+dy*ydir+dz*zdir;dz=a*z1-zref;dy=a*y1-yref;dx=a*x1-xref;" forceGroup="6" particles="1" type="CustomCompoundBondForce" version="1">
<PerBondParameters>
<Parameter name="a"/>
</PerBondParameters>
<GlobalParameters>
<Parameter default="4184.000015258789" name="frc"/>
<Parameter default=".52" name="droff"/>
<Parameter default="0" name="xref"/>
<Parameter default="0" name="yref"/>
<Parameter default="0" name="zref"/>
<Parameter default="1" name="xdir"/>
<Parameter default="0" name="ydir"/>
<Parameter default="0" name="zdir"/>
<Parameter default="9.999999747378752e-07" name="small"/>
</GlobalParameters>
<Bonds>
<Bond p1="0" param1="1"/>
</Bonds>
<Functions/>
</Force>
</Forces>
</System>
===========================system pdb===================================
REMARK MMFP GEO TEST TO FIGURE OUT WHAT VARIOUS OPTIONS ARE
REMARK DATE: 12/23/14 10: 1:21 CREATED BY USER: brookscl
ATOM 1 LJ LJ 1 6.000 0.000 0.000 1.00 0.00 AR
TER 2 LJ 1
END
===========================python script=================================
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout

plats = ['CUDA', 'CPU', 'Reference']

pdb = PDBFile('test1.pdb')
f = open("test1.xml", "r").read()
system = XmlSerializer.deserialize(f)

print "************TEST 3******************"
for w in plats:
pltfrm = Platform.getPlatformByName(w)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
context = Context(system,integrator,pltfrm)
context.setPositions(pdb.positions)
state = context.getState(getEnergy=True, getForces=True)
potential = state.getPotentialEnergy()
forces = state.getForces()

platform = Platform.getName(context.getPlatform())
print platform, potential
print platform, forces

del context, platform, pltfrm, integrator

del system, f, pdb
========================end of python script==============================

[Peter Eastman]

Here's what I get:

************TEST 3******************
Reference 0.0 kJ/mol
Reference [(nan, nan, nan)] kJ/(nm mol)
CPU 0.0 kJ/mol
CPU [(nan, nan, nan)] kJ/(nm mol)

The problem is coming from this part of the expression:

Code: Select all

step(r-droff)*0.5*frc*(r-droff)^2;r=sqrt(r)

(That should be r1=sqrt(r), or something like that. It should be another variable. But that's not actually the issue here.) In this case, r is 0, which is fine for computing the energy, but not the force. The derivative of sqrt(r) is 0.5/sqrt(r), which is infinite. The step function then gets evaluated to 0, and 0*inf=nan.

How about replacing sqrt(r) by sqrt(r+1e-20)? For any nonzero value of r it will make no difference, but for r=0 it will prevent the infinity, and so multiplying the result by 0 will produce 0 as intended.

Peter

[Charles Brooks]

Dear Peter,

Happy New Year. I discovered independently that I could add a "small" parameter to the sqrt and get this to work. Thanks for following through. However, I continue to be plagues using this functionality. I am attaching a gzipped tar file containing a couple of tests fail1.py and fail2.py that illustrate where I am stuck now. For the first, I am getting different answers depending on the platform I use, and the one I am getting from CUDA matches what I get from CHARMM. For the test fail1 I should get:

Energy - ( 316.2023 ) kj/mol
and forces - ( 0.00000 0.00000 0.00000 ) kj/mol/nm
(1626.64715 0.00000 0.00000 )
( 0.00000 0.00000 0.00000 )
( 0.00000 0.00000 0.00000 )

For the fail2.py test, which is a composite of the first test plus another potential
the results should be:

Energy - ( 1364.14725 ) kj/mol
and forces - ( 508.82879 90.16752-841.02259 )
( 1626.64715 0.00000 0.00000 )
( 508.82879 90.16752-841.02259 )
( 508.82879 90.16752-841.02259 )

In this case, none of the platforms seem to calculate the second restraint term defined in the xml file and I also get different answers from CUDA and the other two platforms.

I expect like the other case, this is some issue with the translation of the text to code, but I am at a bit of a loss for figuring out what it is.

Thanks,

Charlie

[Peter Eastman]

Looks like you found a bug! I just merged the fix into the master repository. The PR is here: https://github.com/pandegroup/openmm/pull/776. I also entered a bug report so it can be properly documented: https://simtk.org/tracker/index.php?fun ... 1&atid=435.

Sorry about that, and thanks for the excellent test case.

Peter

[Charles Brooks]

Dear Peter,

Following up on my post in the other forum:

The reason I compiled this was to be able to take advantage of Peter's bug fix to the issue with the CustomCompoundBondForce, however, as demonstrated with the attached tar file providing a simple test that if I either of the two CustomCoumpoundBondForces I get a result (which I've verified is correct), but if I apply both at once, only one of them is computed (it happens to be the first one). I think this too is a bug, and that the two forces should be additive. Can you let me know?

Thanks,

Charlie

[Peter Eastman]

It's getting confused because both forces define global parameters with identical names. Global parameters are stored in the Context and are shared across forces. This is sometimes useful, since it means you can have multiple forces depending on a single parameter. But it's confusing if you aren't expecting that to happen.

In this particular case, your A force defines "xdir" to have a default value of 0, and your B force defines it to have a default value of 1. With only one force or the other it gets the value you're expecting, but when you combine both forces in a single System it will get one value or the other. In this case, "xdir" is ending up set to 0, so the B force gives an energy of 0.

If you don't want global parameters to be shared between forces, it's best to give them unique prefixes: a_xdir and b_xdir, for example.

Peter