[gmx-users] Some notes on parameterising dreiding
Laura Leay
laura.leay at postgrad.manchester.ac.uk
Tue Nov 20 18:29:25 CET 2012
It took me quite a while to get to grips with this and, since the question was asked about a year ago by someone else with no response, I thought I'd share my insight.
LJ potentials: Its important to note the slight difference from the widely used form. You need to convert Ro to sigma by dividing by 2^(1/6) (see wikipedia on the LJ interaction). I've also found that although the paper states that you can either include, exclude or scale the LJ 1-4 interactions I've found it work best to exclude them. As far as I can tell this means setting nrexcl to 4 in your topology.
Bond stretch: This is a straightforward harmonic potential (type 1 in gromacs). All you need to do is convert the length from A to nm and the force constant from kcal/mol/A^2 to kJ/mol/nm^2 by multiplying by 4.184/0.01
Angles: Although the paper states that the cosine harmonic (type 2 in gormacs) is preferred I've seen this generate some unstable systems. Both the paper and the gromacs manual tell you how to convert to the normal harmonic form (type 1 in gromacs). Just convert the force constant from kcal/mol/rad^2 to kJ/mol/rad^2. Leave the angles as degrees.
Dihedrals: There is already a thread on here about this and the conversion does work (search for dihedral form) for methanol at least. Just remember to convert from kcal to kJ. The force constant for resonant bonds seems to be too high though (the paper gives 25 kcal/mol which converts to 25*4.184*0.5=52.3 again see previous posting for the conversion) but a force constant of 40 seems to result in a stable system. Higher values may work but I've not tested them yet. Its possible that including some improper dihedrals (inversions in the dreiding paper) may fix this but its not been tested.
If anyone has any coments on why this force constant seems unacceptable I'd like to hear it!
I'm still testing this out so if the system may still blow up but this is a good start for anyone trying to parameterise dreiding. I hope this helps someone.
Below is the force field file that I've generated for my system. The atoms names follow the opls convention since I originally simulated this system with opls-aa and it would have taken extra effort to change all my topologies (I have some complicated molecules). I would advise someone starting from scratch to stick with the dreiding naming convention, edit this (if you're using it) and check for any double entries that may occur e.g. HC should be exactly the same as HA resulting in the same line twice under [ atomtypes ] so one of them can be removed.
-------------------------------------------------------------
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 3 yes 1 1
[ atomtypes ]
; full atom descriptions are available in ffoplsaa.atp
; name bond_type mass charge ptype sigma epsilon
opls_135 CT 6 12.01100 -0.180 A 3.47299e-01 3.97898e-01
opls_136 CT 6 12.01100 -0.120 A 3.47299e-01 3.97898e-01
opls_137 CT 6 12.01100 -0.060 A 3.47299e-01 3.97898e-01
opls_140 HC 1 1.00800 0.060 A 2.84642e-01 6.35968e-02
opls_145 CA 6 12.01100 -0.115 A 3.47299e-01 3.97898e-01
opls_146 HA 1 1.00800 0.115 A 2.84642e-01 6.35968e-02
opls_186 OS 8 15.99940 -0.400 A 3.03315e-01 4.00409e-01
opls_260 CA 6 12.01100 0.035 A 3.47299e-01 3.97898e-01
opls_902 NT 7 14.00670 -0.630 A 3.26256e-01 3.23842e-01
opls_154 OH 8 15.99940 -0.683 A 3.03315e-01 4.00409e-01 ; methanol
opls_155 HO 1 1.00800 0.418 A 2.84642e-01 6.35968e-02 ; methanol
opls_156 HC 1 1.00800 0.040 A 2.84642e-01 6.35968e-02 ; methanol
opls_157 CT 6 12.01100 0.145 A 3.47299e-01 3.97898e-01 ; methanol
[ bondtypes ]
; i j func b0 kb
CA CT 1 0.14600 292880.0
CT HC 1 0.10900 292880.0
CA CA 1 0.13900 439320.0
CA NT 1 0.13920 292880.0
CT NT 1 0.14620 292880.0
CA HA 1 0.10200 292880.0
CA OS 1 0.13500 292880.0
CT CT 1 0.15300 292880.0
OH HO 1 0.09800 292880.0 ; methanol
OH CT 1 0.14200 292880.0 ; methanol
[ pairtypes ]
[ constrainttypes ]
[ angletypes ]
; i j k func th0 cth
HC CT CT 1 109.471 418.400
HC CT NT 1 109.471 418.400
NT CT NT 1 109.471 418.400
CT CT CT 1 109.471 418.400
HC CT CA 1 109.471 418.400
CA CT CA 1 109.471 418.400
NT CT CA 1 109.471 418.400
CA CA CA 1 120.000 418.400
HA CA CA 1 120.000 418.400
HA CA CT 1 120.000 418.400
CA CA CT 1 120.000 418.400
CA CA NT 1 120.000 418.400
CT CA NT 1 120.000 418.400
HA CA NT 1 120.000 418.400
CA NT CT 1 106.700 418.400
CT NT CT 1 106.700 418.400
CA OS CA 1 104.510 418.400
HO OH CT 1 104.510 418.400 ;methanol
OH CT HC 1 109.471 418.400 ;methanol
HC CT HC 1 109.471 418.400 ;methanol
HC CT CT 2 109.471 470.741
HC CT NT 2 109.471 470.741
NT CT NT 2 109.471 470.741
CT CT CT 2 109.471 470.741
HC CT CA 2 109.471 470.741
CA CT CA 2 109.471 470.741
NT CT CA 2 109.471 470.741
CA CA CA 2 120.000 557.956
HA CA CA 2 120.000 557.956
HA CA CT 2 120.000 557.956
CA CA CT 2 120.000 557.956
CA CA NT 2 120.000 557.956
CT CA NT 2 120.000 557.956
HA CA NT 2 120.000 557.956
CA NT CT 2 106.700 456.093
CT NT CT 2 106.700 456.093
CA OS CA 2 104.510 446.452
HO OH CT 2 104.510 446.452 ;methanol
OH CT HC 2 109.471 470.741 ;methanol
HC CT HC 2 109.471 470.741 ;methanol
[ dihedraltypes ]
; i j k l func phi0 cp mult
CA CA CT X 1 180.0 2.092 6
X CT CA X 1 0.0 4.184 3
X CT CT X 1 0.0 4.184 3
X CT NT X 1 0.0 4.184 3
X CA NT X 1 0.0 4.184 3
CA CA NT X 1 180.0 2.092 6
X CA CA X 1 180.0 40.30 2
X OS CA X 1 180.0 4.184 2
X OH CT X 1 0.0 4.184 3 ;methanol
More information about the gromacs.org_gmx-users
mailing list