[gmx-users] Calculation of [constrainttypes] lengths
Gmx QA
gmxquestions at gmail.com
Wed Dec 2 15:46:28 CET 2015
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
More information about the gromacs.org_gmx-users
mailing list