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 -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.