Custom (improper) torsion, choice of function and order of atoms

[Nicholas Schafer]

Hello all -

I'm working on implementing a term that will keep the stereochemistry of the CB atom correct in my model. The term has the functional form V=k(x-x0)^2 where x = (r_{c->ca} cross r_{ca->n}) dot r_{ca->cb}. So far I have

Code: Select all

chi = CustomTorsionForce("k_chi*(theta-theta0)^2")
chi.addGlobalParameter("k_chi", k_chi)
chi.addGlobalParameter("theta0", theta0)

and, later,

Code: Select all

chi.addTorsion(ca[i], c[i], n[i], cb[i])

I have a few questions:

• 1. Is CustomTorsionForce the right function for this?
  2. If so, is the order of my atoms correct in my addTorsion call (I tried to check by tracing through the code but I'm not 100% sure that I got it right)?
  3. Finally, is there some way to get OpenMM to output the values that it calculates for "theta" so that I can check that I have a reasonable equilibrium value (by, e.g., outputting the values for my input pdb structure before doing any dynamics)?

[Peter Eastman]

The CustomTorsionForce you showed is a reasonable way of keeping four atoms in a plane. However, it isn't the same functional form as what you described:

V=k(x-x0)^2 where x = (r_{c->ca} cross r_{ca->n}) dot r_{ca->cb}

That x is equal to the cosine of the angle between the ca-cb bond and the c-ca-n plane, multiplied by the lengths of the three bonds. But in the CustomTorsionForce, theta is equal to the ca-c-n-cb dihedral angle. Those aren't the same, but either one of them can potentially be used to create the effect you want.

[Nicholas Schafer]

Thanks for the tip, Peter. Can I get OpenMM to output the theta that it calculates? If not, can you point me to the mathematical form of the ca-c-n-cb dihedral angle that CustomTorsionForce uses? Then I'll be able to calculate it myself and figure out what a good equilibrium value will be.

Thanks,
Nick

[Peter Eastman]

OpenMM doesn't have a function for reporting the angle, but MDTraj has one: http://mdtraj.org/development/analysis. ... -dihedrals. I'm not sure whether it puts zero at the same point (cis vs. trans), but other than that it should be the same.

The other thing you can do if you're just trying to debug it is create a CustomTorsionForce whose energy function is just "theta", then add only a single torsion to it. That way the energy will equal the angle of that one torsion.

[Nicholas Schafer]

MDTraj did it for me. I didn't find any evidence that there was a difference between the ways OpenMM and MDTraj are computing dihedrals. Thanks, Peter.

[Nicholas Schafer]

Upon closer inspection, it actually seems as though there is a difference between the way things are being computed in MDTraj and OpenMM. In particular, for the same atom inputs, MDTraj always yields values of the same sign (either positive or negative), while OpenMM computes dihedral angles that have both a different absolute value than those computed by MDTraj and have both positive and negative signs. I figured this out by following Peter's second tip of outputting the thetas by using only 'theta' as the expression for the energy and adding only a single torsion at a time.

[Nicholas Schafer]

I ended up going with a CustomCompoundBondForce instead and implementing the original functional form. That seems to have worked nicely. Thanks again!

[web design]

Posts were read Thanks for sharing