[gmx-users] Accessing the distance_restraints energy term
chris.neale at utoronto.ca
chris.neale at utoronto.ca
Thu Nov 6 05:14:12 CET 2008
Thanks Mark for confirming that Distance restraint energies do get
written when they are present. I agree then that the issue here is
that they are not present.
I'll try to tackle this again tomorrow, but have included full
information here for completeness. Also, I would be grateful if you
could send me your scripts again off list just so that I know that I
have the right thing.
The package that I am calling "yours" was sent to me by Mark after I
requested it based on this post.
http://www.gromacs.org/pipermail/gmx-users/2008-January/031854.html
A quote from one of Marks posts from this Feb probably explains my problems:
http://www.gromacs.org/pipermail/gmx-users/2008-February/032530.html
> Anyway, the point is that using any of Yuguang's awk or shell scripts is
> not supported with the ffcharmm*itp files produced by my conversion
> script. I've no idea what the combination would do, and no interest in
> finding out! :-) The rudimentary documentation in the headers of my
> scripts discuss the necessary issues.
Still, here is what I have done. I didn't use anybody's scripts as I
received them. I began with the scripts that I downloaded from the
gromacs users contributions page to do my own conversion of our
modified charmm parameter files and I found lots of problems. For
example, the epsilon that I find in par_all27_lipid.prm is negative,
which leads to negative C6 C12 converted values and lots of nan's.
Below is an excerpt from par_all27_lipid.prm.
!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)
!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j
!
!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
!
HOL 0.0 -0.046 0.2245
HAL1 0.0 -0.022 1.3200 ! alkane, 3/92
HAL2 0.0 -0.028 1.3400 ! alkane, yin and mackerell, 4/98
HAL3 0.0 -0.024 1.3400 ! alkane, yin and mackerell, 4/98
HCL 0.0 -0.046 0.2245 ! ethanolamine
...
...
So I tried to fix what I could by taking any relevant lines from a
package that Mark sent me off list early in 2008. This generally
worked since most of my energies match very well.
I then got stuck on the angles, so I looked at the awk -f files from
the package that I got from Mark, specifically the UB-related awk
input file, and it appears to add a distance_restraints section -- I
guess that this is not in fact Mark's file and must be Yuguang's.
Below is what I have currently to deal with guanidine. Note that the
relevant [ angletypes ] line is:
NC2 C NC2 1 120.00000 435.13600 75312.00000 0.23642
which I modified to
NC2 C NC2 1 120.00000 435.13600 10.00000 0.23642
just to see if these last two columns were being picked up by gmx
3.3.1, and noted that the angle energy did not change at all. It was
then that I tried to include distance_restraints.
I did not include any -DDISRES term, as I don't protect the
distance_restraints section by an #ifdef.
[cneale at rb19 CM]$ cat test.top
#include "ffcharmm.itp"
#include "guanidine.itp"
#include "cla.itp"
#include "tip3p-charmm.itp"
[ system ]
; Name
CNtest
[ molecules ]
; Compound #mols
Guanidine 50
Chloride 50
SOL 712
###
[cneale at rb19 CM]$ cat guanidine.itp
[ moleculetype ]
; Name nrexcl
Guanidine 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
typeB chargeB massB
1 C 1 GHC CZ 1 0.64 12.011
; qtot 0.64
2 NC2 1 GHC NH1 1 -0.8 14.007
; qtot -0.16
3 HC 1 GHC HH11 1 0.46 1.008
; qtot 0.3
4 HC 1 GHC HH12 1 0.46 1.008
; qtot 0.76
5 NC2 1 GHC NH2 1 -0.8 14.007
; qtot -0.04
6 HC 1 GHC HH21 1 0.46 1.008
; qtot 0.42
7 HC 1 GHC HH22 1 0.46 1.008
; qtot 0.88
8 NC2 1 GHC NH3 1 -0.8 14.007
; qtot 0.08
9 HC 1 GHC HH31 1 0.46 1.008
; qtot 0.54
10 HC 1 GHC HH32 1 0.46 1.008 ; qtot 1
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 5 1
1 8 1
2 3 1
2 4 1
5 6 1
5 7 1
8 9 1
8 10 1
[ pairs ]
; ai aj funct c0 c1 c2 c3
2 6 1
2 7 1
2 9 1
2 10 1
3 5 1
3 8 1
4 5 1
4 8 1
5 9 1
5 10 1
6 8 1
7 8 1
[ angles ]
; ai aj ak funct c0 c1 c2
c3
2 1 5 1
2 1 8 1
5 1 8 1
1 2 3 1
1 2 4 1
3 2 4 1
1 5 6 1
1 5 7 1
6 5 7 1
1 8 9 1
1 8 10 1
9 8 10 1
[ dihedrals ]
; ai aj ak al funct c0 c1
c2 c3 c4 c5
5 1 2 3 3
5 1 2 4 3
8 1 2 3 3
8 1 2 4 3
2 1 5 6 3
2 1 5 7 3
8 1 5 6 3
8 1 5 7 3
2 1 8 9 3
2 1 8 10 3
5 1 8 9 3
5 1 8 10 3
[ distance_restraints ]
; Actually Urey-Bradley terms (thanks to Mark Abraham's script to see
how it is done)
; ai aj type index type' low up1 up2
fac
2 5 1 0 1 0.23642 0.23642 0.73642
75.312
2 8 1 1 1 0.23642 0.23642 0.73642
75.312
5 8 1 2 1 0.23642 0.23642 0.73642
75.312
###########
I used a recursive grep -v to cull any lines that were not required
from my ffcharmm files since I wanted to test it for these molecules
and then give the files to my colleague and not be responsible for any
of the other parameters that I did not check.
[cneale at rb19 CM]$ cat ffcharmmbon.itp
[ bondtypes ]
; i j func b0 kb
C C 1 0.13350 502080.0000000000
NC2 C 1 0.13650 387438.4000000000
NC2 HC 1 0.10000 380744.0000000000
OT HT 1 0.09572 376560.0000000001
[ angletypes ] ; normal potential harmonic in theta
; i j k func th0 kth
HC NC2 C 1 120.00000 410.03200 0.00000 0.00000
HC NC2 HC 1 120.00000 209.20000 0.00000 0.00000
HT OT HT 1 104.52000 460.24000 0.00000 0.00000
NC2 C NC2 1 120.00000 435.13600 75312.00000 0.23642
[ dihedraltypes ]
; i j k l func phi0 kphi n
X C C X 3 33.47200 -0.00000 -33.47200
-0.00000 0.00000 0.00000
X C NC2 X 3 18.82800 -0.00000 -18.82800
-0.00000 0.00000 0.00000
; Now the improper dihedrals
#define improper_NC2_X_X_C_ 180.00000 83.68000 2
#########
[cneale at rb19 CM]$ cat ffcharmmnb.itp
[ atomtypes ]
;name mass charge ptype c6 c12
C 12.011 0.000 A 3.7702860800e-03 7.7215458918e-06
HC 1.008 0.000 A 3.1539699357e-09 1.2921281844e-17
HT 1.008 0.000 A 3.1539699357e-09 1.2921281844e-17
NC2 14.007 0.000 A 4.2939997181e-03 5.5086142385e-06
OT 15.999 0.000 A 2.4895403920e-03 2.4347673691e-06
CLA 35.450 0.000 A 1.0991274335e-02 4.8123052707e-05
[ nonbond_params ]
; i j func c6 c12
C C 1 3.7702860800e-03 7.7215458918e-06
C HC 1 7.2125924088e-05 4.3697425971e-09
C HT 1 7.2125924088e-05 4.3697425971e-09
C NC2 1 4.0420100498e-03 6.5816043354e-06
C OT 1 3.0987482359e-03 4.4356775553e-06
C CLA 1 6.5152471328e-03 1.9745500580e-05
HC HC 1 3.1539699357e-09 1.2921281844e-17
HC HT 1 3.1539699357e-09 1.2921281844e-17
HC NC2 1 6.3972926767e-05 2.5494521165e-09
HC OT 1 4.3824474706e-05 1.3719498873e-09
HC CLA 1 1.6747395113e-04 2.0175221921e-08
HT HT 1 3.1539699357e-09 1.2921281844e-17
HT NC2 1 6.3972926767e-05 2.5494521165e-09
HT OT 1 4.3824474706e-05 1.3719498873e-09
HT CLA 1 1.6747395113e-04 2.0175221921e-08
NC2 NC2 1 4.2939997181e-03 5.5086142385e-06
NC2 OT 1 3.2745881973e-03 3.6735174929e-06
NC2 CLA 1 7.0886836461e-03 1.7334802405e-05
OT OT 1 2.4895403920e-03 2.4347673691e-06
OT CLA 1 5.4809880447e-03 1.1883812906e-05
CLA CLA 1 1.0991274335e-02 4.8123052707e-05
[ pairtypes ]
; i j func c6 c12
C C 1 3.7702860800e-03 7.7215458918e-06
C HC 1 7.2125924088e-05 4.3697425971e-09
C NC2 1 4.0420100498e-03 6.5816043354e-06
HC HC 1 3.1539699357e-09 1.2921281844e-17
HC NC2 1 6.3972926767e-05 2.5494521165e-09
NC2 NC2 1 4.2939997181e-03 5.5086142385e-06
Thanks again,
Chris.
More information about the gromacs.org_gmx-users
mailing list