[gmx-users] Calculation of [constrainttypes] lengths

[Gmx QA]

Hello everyone

I know how to calculate the lengths in the [constrainttypes] section has
been brought up before on this list, but the answer from Erik Marklund was
to use bond lengths, angles and moment of inertia.

Well, I tried, but I don't quite arrive at the same answer as in the
forcefield-files, so I'd appreciate if someone could take the time to go
though my math. My aim is to implement some additional vsites, and I'd like
to be as thorough as possible.

I have tried to re-calculate the this from the oplsa.ff/ffbonded.itp -file

; constraints for 2nd type rigid CH3 (angle *-CT-HC is 110.7)
 MCH3B    CT 2    0.167031

My short python-script is as follows:

#!/usr/bin/python


import math

# opls equilibrium values

b0CH = 0.109
b0CC = 0.1522
angle = 110.7

# opls masses
mH = 1.008
mC = 12.011

# mass of the virtual site is then:
mvs = (3*mH+mC)/2

# Start to compute the moment of inertia of the CH3 group when rotating
around
# the "next" carbon atom
# Use the cosine rule to calculate the distance between each hydrogen and
the
# carbon at the center of rotation.
x2 = (b0CC * b0CC) + (b0CH * b0CH) - (2 * b0CC * b0CH *
math.cos(math.radians(angle)))

b = math.sqrt(x2)

# The final moment of inertia for three hydrogens and one carbon is thus
I = (b*b*mH)*3 + b0CC*b0CC*mC

# When using vsites, the moment of inertia must be the same, so use that to
# compute the new constraint length for two mass centers
c = math.sqrt(I/(2*mvs))
print c

When I run this script, I get c = 0.167072934788 whereas the value above is
0.167031. Ok, the difference is only in the third-fourth decimal place, but
since this is based on equilibrium values I'd have expected a better
agreement, perhaps?

Best
/PK