Units of force constants

[Siddharth Srinivasan]

What are the units of the various force constants for HarmonicBondForce, HarmonicAngleForce, PeriodicTorsionForce and RBTorsionForce that is expected by OpenMM?

[Siddharth Srinivasan]

I ask because Im trying some code integration, and am getting very high potential energies. I have a 33 residue protein. Psuedo code for the bonded energies looks like
{{{
HarmonicBondForce *bondForce = new HarmonicBondForce();

sys->addForce(bondForce);

numBonds = total number of bonds in system, upper triangular matrix
for (unsigned int i = 0; i < numBonds; i++) {
bondForce->addBond(atom1, atom2, bzero, kb);
// where kb is in Kcal/mol
// bzero in nanometers
}

// Set positions for OpenMM context
context->setPositions(pos);
context->setVelocities(vel);


// Ask for forces
// Ask for potential energy
}}}
Im expecting about 38 Kcal/mol bond energy, but I end up getting 31690.730312 from the OpenMM version.

[John Chodera]

Hi Siddharth,

The Doxygen documentation contains is supposed to contain detailed descriptions of the parameters expected by these force classes, and their units.

You can access the Doxygen documentation (compiled from the source headers, found in the source distribution under src/openmmapi/include/openmm/) here:

https://simtk.org/api_docs/openmm/

It looks like the developers neglected to specify the units here. I believe they are supposed to be in kJ/mol/nm^2.

The developers also neglect to specify the equation for the harmonic potential. I *think* it's

(K/2) (r - r0)^2

My guess is that you specified your force constant in kcal/mol/A^2, and were specifying the spring constant K' in

K' (r - r0)^2

This would give an energy off by a factor of 10 * 10 * 4.184 * 2, which mostly explains the difference between your observed and expected bond energies.

Good luck!

John

[Siddharth Srinivasan]

Actually, I just saw one of the examples (HelloEthane), which is much more descriptive. Based on this, the call is
{{{
k = 2 * force_constant in Kcal/mol * 4.184 * (100),
eq_distance = eq_distance_in_angstroms * 0.1f;

bondForce->addBond(atom1, atom2, k, eq_distance);
}}}

Im obviously doing something very wrong since the potential energy is 26518803.059273, but I cant make out what!

[Siddharth Srinivasan]

Whoops sorry John, I think we were both posting at the same time, I didn't see yours until after I posted! The HelloEthane example is much more descriptive, and I came to the same conclusions as you. I'm still doing something wrong though after following the example, as I mentioned in my previous post

[John Chodera]

Siddharth,

That's so odd -- it's almost as if you want to *divide* your original force constant by (10 * 10 * 4.184 * 2) to get something of the magnitude you expect, but doing that doesn't make any sense if the original units are kcal/mol/A^2...

Maybe the developers can help out here...

[Siddharth Srinivasan]

Well I have most of my code in place at least, and plan to debug each individual force component before proceeding, thats why this issue came up for me. Hopefully we can figure it out soon. I tried going through the source code, in the SimTKReferenceKernels folder (I think thats the name), but thats just wasting a bit too much time! Thanks for your help anyway

[Michael Sherman]

Hi, Siddharth. HelloEthane is probably the best documentation at the moment and your last posted conversion looked correct to me. There was also a handout at the OpenMM Developer's Workshop giving all the OpenMM functional forms, but I don't know where to find an electronic copy of that.

The energy problem you're seeing now is likely due to something other than the units -- I would suggest trying a single pair of atoms to verify that you are getting the energy you expect.

Regards, Sherm

[Siddharth Srinivasan]

Hi Sherm

I have the handout, and pretty much followed exactly the same procedure from the example, and the handout equations. I'll try narrowing down the problem further, using a single water molecule or something trivial. I've also used the createExceptionsFromBonds() method to create 1-4 exceptions for AMBER force field, does that take into account 5 sided rings having no 1-4 scaled interactions (since there is an angle term that connect them as well)?

[Michael Sherman]

- for working out the bond terms I would suggest building your trivial example with no nonbonded forces at all, just to keep it simple
- regarding 5-sided rings, all heavy atoms are 1-2 or 1-3 bonded along the shortest path, so NonbondedForce::createExceptionsFromBonds() will set all the nonbonded forces to zero. That method only works with a list of 1-2 bonds, though, so doesn't know about angles or torsions.
- you can view the exact list of exceptions that got created via getExceptionParameters().

Sherm