Units of force constants

The functionality of OpenMM will (eventually) include everything that one would need to run modern molecular simulation.
User avatar
Siddharth Srinivasan
Posts: 223
Joined: Thu Feb 12, 2009 6:49 pm

Units of force constants

Post by Siddharth Srinivasan » Fri Jul 03, 2009 5:00 pm

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

User avatar
Siddharth Srinivasan
Posts: 223
Joined: Thu Feb 12, 2009 6:49 pm

RE: Units of force constants

Post by Siddharth Srinivasan » Fri Jul 03, 2009 5:30 pm

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.

User avatar
John Chodera
Posts: 53
Joined: Wed Dec 13, 2006 6:22 pm

RE: Units of force constants

Post by John Chodera » Fri Jul 03, 2009 6:40 pm

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

User avatar
Siddharth Srinivasan
Posts: 223
Joined: Thu Feb 12, 2009 6:49 pm

RE: Units of force constants

Post by Siddharth Srinivasan » Fri Jul 03, 2009 6:43 pm

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!

User avatar
Siddharth Srinivasan
Posts: 223
Joined: Thu Feb 12, 2009 6:49 pm

RE: Units of force constants

Post by Siddharth Srinivasan » Fri Jul 03, 2009 6:46 pm

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

User avatar
John Chodera
Posts: 53
Joined: Wed Dec 13, 2006 6:22 pm

RE: Units of force constants

Post by John Chodera » Fri Jul 03, 2009 6:49 pm

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...

User avatar
Siddharth Srinivasan
Posts: 223
Joined: Thu Feb 12, 2009 6:49 pm

RE: Units of force constants

Post by Siddharth Srinivasan » Fri Jul 03, 2009 7:27 pm

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

User avatar
Michael Sherman
Posts: 807
Joined: Fri Apr 01, 2005 6:05 pm

RE: Units of force constants

Post by Michael Sherman » Sat Jul 04, 2009 10:10 am

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

User avatar
Siddharth Srinivasan
Posts: 223
Joined: Thu Feb 12, 2009 6:49 pm

RE: Units of force constants

Post by Siddharth Srinivasan » Sat Jul 04, 2009 10:59 am

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

User avatar
Michael Sherman
Posts: 807
Joined: Fri Apr 01, 2005 6:05 pm

RE: Units of force constants

Post by Michael Sherman » Sat Jul 04, 2009 11:31 am

- 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

POST REPLY