MDTraj vs openmm torsion equilibrium angles

[cwalker7]

There seems to be a mismatch between the definition of the equilibrium torsion angle in openmm for periodic torsions, and MDTraj compute_dihedrals function.

OpenMM periodic torsion defines 0 degrees as the cis conformation.

For example, if I set backbone torsion equilibrium angle to 130 degrees, and enforce no other torsions, we get a torsion distribution like this, centered at about -30 degrees (this also forms a very nice helix!). There may be sterics or nonbonded interactions preventing it reaching the set value, which I'm guessing should be -50.

torsion_dist_state10.pdf

[cwalker7]

So it definitely looks like we need to shift the torsion equilibrium angle by 180 degrees to be consistent with MDTraj. For example, if we want the torsions to be measured as +50 degrees in MDTraj or VMD, we need to input into openMM an equilibrium value of -130. I'm thinking we can have a user input the torsion as +50 degrees, and convert the angle internally in cgmodel.py.

However, this all seems to go against this convention from openmm documentation:
http://docs.openmm.org/latest/userguide/theory.html#rbtorsionforce
"For reason of convention, PeriodicTorsionForce and RBTorsonForce define the torsion angle differently. θ is zero when the first and last particles are on the same side of the bond formed by the middle two particles (the cis configuration), whereas ϕ is zero when they are on opposite sides (the trans configuration). This means that θ = ϕ - π."

The θ=0 for cis configuration is what MDTraj uses.

[mrshirts]

I'd post an issue with OpenMM to ask for clarification. I'm not sure it totally matters, as long as the documentation is clear.

[cwalker7]

Was thinking about the periodic torsion potentials again - what we are seeing for periodicity of 1 makes perfect sense actually because theta_0 is a phase offset, not an equilibrium value. At theta_0, the potential is at a maximum, and at (180-theta_0) a minimum.

Likewise, if we use periodicity of 3, the maximums are (theta_0/3, theta_0/3+120, theta_0-120),
and the minimums are at (theta_0/3+60, theta_0/3+180, theta_0/3-60). Should be able to work out a general formula for specifying the equilibrium angle rather than the phase offset.

On a related note to shirtsgroup/cg_openmm#53, it seems that using a single periodic torsion, if we want a symmetric potential about theta=0, we are limited to minima at +/-60 for periodicity=3, +/-90 for periodicity=2, and so on. That will probably work for the first paper but we will definitely want the sums of periodic terms in the future.