Units of force constants

[Siddharth Srinivasan]

Thanks for your earlier reply Michael. So I ran my code with the second simplest of all systems, a water molecule (2 bonds). Here is a printout of the parameters that were set.

{{{
// The bonds that were set {force_constant, eq_distance}
Bond id 0 462750.406250 0.095720
Bond id 1 462750.406250 0.095720
OpenMM Platform: Reference
OpenMM Speed: double free or corruption (fasttop)

// The output forces
[-205317, -173637, -406466]
[33314.9, 380623, 112438]
[172002, -206985, 294028]

// The output potential energy
PE 81994.174235 KE 0.000748

// The positions of the atoms
[2.717, 1.873, 1.539]
[2.637, 0.959001, 1.269]
[2.304, 2.37, 0.833]
}}}

The code that I used to create the bonds was
{{{
const float
kb = force_k_in_KcalperMol * ( 4.184 * 100.0f * 2),
bzero = eq_distance_in_A * 0.1f; //nm

// Set OpenMM bond parameters (nanometers)
bondForce->addBond(atom1, atom2, bzero, kb);

}}}

The code used to get the potential energy was
{{{
printf ("PE %lf KE %lf\n",
state.getPotentialEnergy()/4.184,
state.getKineticEnergy())/4.184; // PE in KCal/mol

}}}

Anything jump out at you? I assumed here that the potential energy is returned in Kj/Mol, so I convert it to Kcal/Mol for display. As you can see, the potntial energy should be in the region of 0.000364 Kcal/mol, not 81994.174235 as I get

[John Chodera]

This is just a wild guess, but the addBond() function of HarmonicBondForce expects double arguments (and not floats) for 'length' and 'k'. I know Mark Friedrichs had some trouble with automatic casts not working as expected for him under Linux, though they worked fine for me on my Mac. Maybe you can try making sure to explicitly cast to double before calling addBond()?

Again, just a stab in the dark here.

John

[Michael Sherman]

I checked your numbers in Matlab and agree with you -- 0.000364 kcal/mol is the right answer. And I get that whether I work in Amber units or use your conversions to kJ/nm and back.

So I think that narrows it down to the parts of the program we haven't seen. For example, the non-zero KE suggests that you took a step -- is that right? Are atom1 and atom2 set to the right indices? Did you use try/catch to see if there are any error conditions? Are you using the prebuilt binaries or have you built your own? If the latter, do the test cases pass? Did you request positions, velocities, and energies in getState()? Etc.

Sherm

[Siddharth Srinivasan]

I'll try to be as precise as possible with what I have done. I'm hoping that once I get his very basic routine clear in my mind, the resrt should be much easier, so thanks for helping out. Here is what my code looks like

== openmm_init () ==
{
HarmonicBondForce *bondForce = new HarmonicBondForce();
sys->addForce(bondForce);
for (unsigned int i = 0; i < numBonds; i++) {
const double
k = force_constant_in_KcalperMolA^2 * ( 4.184 * 100.0 * 2.0),
eq = eq_distance_in_A * 0.1; //nm
bondForce->addBond(atom1, atom2, eq, k);
printf ("Bond id %d %lf %lf\n", i, k, eq);
}

// Integrator
Integrator *integrator = new VerletIntegrator(md->dt);

// Set the coordinates and masses.
vector<Vec3> pos(myopenmm->numAtoms);
vector<Vec3> vel(myopenmm->numAtoms);

// Get masses, coordinates, vel from my program datastructures
// IMPORTANT: adding particles to OpenMM system here.
for (unsigned int i = 0; i < myopenmm->numAtoms; i++) {
pos = Vec3(coords.m[0], coords.m[1], coords.m[2]);
vel = Vec3(vels.m[0], vels.m[1], vels.m[2]);
sys->addParticle (masses.m[0]);
}

// Get context. Make sure this is done after the system is fully set up
OpenMMContext *context = new OpenMMContext(*sys, *integrator);

printf("OpenMM Platform: %s\n", context->getPlatform().getName().c_str());
printf("OpenMM Speed: %s\n", context->getPlatform().getSpeed());

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

// Assign to ZymeCAD OpenMM interface
myopenmm->openmm_sys = sys;
myopenmm->integrator = integrator;
myopenmm->context = context;

return myopenmm;
}

My main program looks like
{{{

// openmm_t is my local interface to openmm. It contains the OpenMM System, Integrator and Context.

openmm_t *myopenmm = openmm_init ()
// Print forces, energy and positions of all atoms as initial check.
openmm_get_state (myopenmm, OPENMM_FORCE);
openmm_get_state (myopenmm, OPENMM_ENERGY);
openmm_get_state (myopenmm, OPENMM_POSITION);
sys.exit()
}}}

So to answer your questions:
>> For example, the non-zero KE suggests that you took a step -- is that right
I have initialized velocities from my programs velocities, thats true. But I dont actually take any step, I want to check my initial forces so I exit my program immediately after printing out the initial state.

>> Are atom1 and atom2 set to the right indices? Did you use try/catch to see if there are any error conditions?
I think they are set to the right values, since the positions are printed out correctly when I query the State using *openmm_get_state (myopenmm, OPENMM_POSITION);*. I should add, that openmm_get_state is simply a function that accepts an enumeration, and calls State.getPositions, State.getvelocities, state,getEnergy or state.getForces as required.

>> Are you using the prebuilt binaries or have you built your own? If the latter, do the test cases pass?
Prebuilt binaries. I've not spent too much time in getting it to compile from source, but maybe I will to check the unit tests if you think that will help

>> Did you request positions, velocities, and energies in getState()? Etc.
As described above, I did all of the above requests after my openmm_init() function.

[Siddharth Srinivasan]

Pardon the horrible formatting..a Preview button would be infinitely useful!

[Siddharth Srinivasan]

I think I figured out what the problem is..when I set positions using
{{{
// IMPORTANT: adding particles to OpenMM system here.
for (unsigned int i = 0; i < myopenmm->numAtoms; i++) {
pos = Vec3(coords.m[0], coords.m[1], coords.m[2]);
vel = Vec3(vels.m[0], vels.m[1], vels.m[2]);
sys->addParticle (masses.m[0]);
}
}}}
the distance that is computed is in Angstroms, I need to divide each coordinate by 10.0f to get it in nm! Does that sound about right? I tried it and I do get the correct potential energy...

[Michael Sherman]

That makes a lot of sense -- I assumed the positions you displayed were Angstroms, but if those were nm that was a very stretched water molecule!

Congratulations on finding that. I told you getting the units right is the hardest part!

Regards, Sherm

[John Chodera]

> I told you getting the units right is the hardest part!

Sherm, this makes me wonder if an MMTK- or Scientific Python-like Units class might make the most sense for OpenMM as well. In the "Feature Request" I submitted that Peter closed (URL below) I outlined what we did for the MMTOOLS SimTK project in Python (which is based on MMTK's concept). This made it *extremely* easy to get the units right, by making sure you are always multiplying and dividing unit-bearing numbers by the appropriate units:

https://simtk.org/tracker/index.php?fun ... 1&atid=436

For example, the call to addBond would look like:

bonds.addBond(i, j, length * Units.Angstrom, k * Units.kcal / Units.mol / Units.Angstrom);

As long as the 'Units' class represents all units in terms of fundamental units of length, mass, and time, then the user never even needs to know what these fundamental units are, and the conversions are done correctly every time. As it is, the lack of such a feature and the incomplete documentation of units in the Doxygen documentation is just causing headaches.

Scientific Python takes this concept yet further, by making unit-bearing numbers special objects. The handling of this subject by both MMTK and Scientific Python are worth looking at, especially if the goal is to facilitate making it "easy to do things right".

John

[Michael Sherman]

Let's start another thread for this discussion -- I think we wore this one out! I'll call it "Units conversion in OpenMM" and include your last post, John.

Sherm

[Siddharth Srinivasan]

Thanks for you help, John and Michael!