[gmx-users] need help with free energy of solvation calculation
Moore, Jonathan (J)
JMoore2 at dow.com
Fri Jan 6 20:15:15 CET 2006
I'm trying to reproduce the free energy of solvation difference between
cellononaose and 2,3,6-trimethylcellononaose as reported in as a way to get
my feet wet for using GROMACS to calculate free energies in these types of
systems:
H. Yu, M. Amann, T. Hansson, J. Köhler, G. Wich, and W. F. van Gunsteren.
"Effect of methylation on the stability and solvation free energy of amylose
and cellulose fragments: a molecular dynamics study." Carbohydr. Res., 339,
1697-1709 (2004).
http://dx.doi.org/10.1016/j.carres.2004.05.003
The reported value in Table 3 of that paper is 393 kJ per mol. I tried to
calculate the value with Gromacs 3.3rc3 using thermodynamic integration with
soft cores (since the hydrogens that change to methyls have no sigma or
epsilon) and 20 lambda values between 0 and 1 (50 ps of equilibration
followed by 150 ps of production at each lambda value). I have reasonable
confidence that I have the force field implemented correctly because I
calculated exactly the same values for the various energy terms as one of
the paper's authors did for a configuration of an isolated cellobiose
molecule provided by the author. However, the free energy difference that I
calculated was the wrong sign and an order of magnitude too large, so I must
have done something very wrong. One thing that has me perplexed is that if
I do a gmxdump of the .tpr file, it looks like the A and B states are
exactly the same for things like atom type, charge, and mass...suggesting
that I'm not specifying the B state successfully. But if the A and B states
are exactly the same, then why am I getting anything at all for dgdl? I'm
relatively new to GROMACS and have never done a coupling-parameter-type
calculation before.
I'll copy below the content of my .mdp file and the first part of my .itp
file. Could someone with free energy experience have a look at them and let
me know if something looks wrong? Also, how would you recommend that I try
to figure out what is going wrong?
.mdp
======
;
; Input file
;
title = FE run
cpp = /usr/bin/cpp
constraints = all-bonds
constraint_algorithm= lincs
integrator = md
dt = 0.002 ; ps !
nsteps = 100000 ;
nstcomm = 500
nstxout = 10000
nstvout = 10000
nstfout = 0
nstlog = 10000
nstenergy = 250
nstxtcout = 250
nstlist = 5
ns_type = grid
rlist = 0.8
coulombtype = Reaction-Field
rcoulomb = 1.4
epsilon_r = 78.5
vdwtype = Cut-off
rvdw = 1.4
; temperature coupling is on in two groups
Tcoupl = nose-hoover
tc-grps = polymer SOL
tau_t = 0.4 0.4
ref_t = 300 300
; Energy monitoring
energygrps = polymer SOL
; Isotropic pressure coupling is now on
Pcoupl = berendsen
Pcoupltype = isotropic
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.01325
; Generate velocites is off at 300 K.
;gen_vel = no
;gen_temp = 363.0
;gen_seed = 173529
free_energy = yes
init_lambda = 0.55
sc_alpha = 1.51
.itp
===========
[ moleculetype ]
; name nrexcl
cellononaose 3
[ atoms ]
; nr type resnr residu atom cgnr charge
1 H 1 H000 HO4 1 0.410
2 OA 1 H000 O4 1 -0.642
3 CH1 1 H000 C4 1 0.232
4 CH1 1 H000 C3 2 0.232 CH1 0.180
5 OA 1 H000 O3 2 -0.642 OA -0.360
6 H 1 H000 HO3 2 0.410 CH3 0.180
7 CH1 1 H000 C2 3 0.232 CH1 0.180
8 OA 1 H000 O2 3 -0.642 OA -0.360
9 H 1 H000 HO2 3 0.410 CH3 0.180
10 CH2 1 H000 C6 4 0.232 CH2 0.180
11 OA 1 H000 O6 4 -0.642 OA -0.360
12 H 1 H000 HO6 4 0.410 CH3 0.180
13 CH1 1 H000 C5 5 0.376
14 OA 1 H000 O5 5 -0.480
15 CH1 1 H000 C1 5 0.232
16 OA 1 H000 O1 5 -0.360
17 CH1 2 000 C4 5 0.232
18 CH1 2 000 C3 6 0.232 CH1 0.180
19 OA 2 000 O3 6 -0.642 OA -0.360
20 H 2 000 HO3 6 0.410 CH3 0.180
21 CH1 2 000 C2 7 0.232 CH1 0.180
22 OA 2 000 O2 7 -0.642 OA -0.360
23 H 2 000 HO2 7 0.410 CH3 0.180
24 CH2 2 000 C6 8 0.232 CH2 0.180
25 OA 2 000 O6 8 -0.642 OA -0.360
26 H 2 000 HO6 8 0.410 CH3 0.180
27 CH1 2 000 C5 9 0.376
28 OA 2 000 O5 9 -0.480
29 CH1 2 000 C1 9 0.232
30 OA 2 000 O1 9 -0.360
31 CH1 3 000 C4 9 0.232
32 CH1 3 000 C3 10 0.232 CH1 0.180
33 OA 3 000 O3 10 -0.642 OA -0.360
34 H 3 000 HO3 10 0.410 CH3 0.180
35 CH1 3 000 C2 11 0.232 CH1 0.180
36 OA 3 000 O2 11 -0.642 OA -0.360
37 H 3 000 HO2 11 0.410 CH3 0.180
38 CH2 3 000 C6 12 0.232 CH2 0.180
39 OA 3 000 O6 12 -0.642 OA -0.360
40 H 3 000 HO6 12 0.410 CH3 0.180
41 CH1 3 000 C5 13 0.376
42 OA 3 000 O5 13 -0.480
43 CH1 3 000 C1 13 0.232
44 OA 3 000 O1 13 -0.360
45 CH1 4 000 C4 13 0.232
46 CH1 4 000 C3 14 0.232 CH1 0.180
47 OA 4 000 O3 14 -0.642 OA -0.360
48 H 4 000 HO3 14 0.410 CH3 0.180
49 CH1 4 000 C2 15 0.232 CH1 0.180
50 OA 4 000 O2 15 -0.642 OA -0.360
51 H 4 000 HO2 15 0.410 CH3 0.180
52 CH2 4 000 C6 16 0.232 CH2 0.180
53 OA 4 000 O6 16 -0.642 OA -0.360
54 H 4 000 HO6 16 0.410 CH3 0.180
55 CH1 4 000 C5 17 0.376
56 OA 4 000 O5 17 -0.480
57 CH1 4 000 C1 17 0.232
58 OA 4 000 O1 17 -0.360
59 CH1 5 000 C4 17 0.232
60 CH1 5 000 C3 18 0.232 CH1 0.180
61 OA 5 000 O3 18 -0.642 OA -0.360
62 H 5 000 HO3 18 0.410 CH3 0.180
63 CH1 5 000 C2 19 0.232 CH1 0.180
64 OA 5 000 O2 19 -0.642 OA -0.360
65 H 5 000 HO2 19 0.410 CH3 0.180
66 CH2 5 000 C6 20 0.232 CH2 0.180
67 OA 5 000 O6 20 -0.642 OA -0.360
68 H 5 000 HO6 20 0.410 CH3 0.180
69 CH1 5 000 C5 21 0.376
70 OA 5 000 O5 21 -0.480
71 CH1 5 000 C1 21 0.232
72 OA 5 000 O1 21 -0.360
73 CH1 6 000 C4 21 0.232
74 CH1 6 000 C3 22 0.232 CH1 0.180
75 OA 6 000 O3 22 -0.642 OA -0.360
76 H 6 000 HO3 22 0.410 CH3 0.180
77 CH1 6 000 C2 23 0.232 CH1 0.180
78 OA 6 000 O2 23 -0.642 OA -0.360
79 H 6 000 HO2 23 0.410 CH3 0.180
80 CH2 6 000 C6 24 0.232 CH2 0.180
81 OA 6 000 O6 24 -0.642 OA -0.360
82 H 6 000 HO6 24 0.410 CH3 0.180
83 CH1 6 000 C5 25 0.376
84 OA 6 000 O5 25 -0.480
85 CH1 6 000 C1 25 0.232
86 OA 6 000 O1 25 -0.360
87 CH1 7 000 C4 25 0.232
88 CH1 7 000 C3 26 0.232 CH1 0.180
89 OA 7 000 O3 26 -0.642 OA -0.360
90 H 7 000 HO3 26 0.410 CH3 0.180
91 CH1 7 000 C2 27 0.232 CH1 0.180
92 OA 7 000 O2 27 -0.642 OA -0.360
93 H 7 000 HO2 27 0.410 CH3 0.180
94 CH2 7 000 C6 28 0.232 CH2 0.180
95 OA 7 000 O6 28 -0.642 OA -0.360
96 H 7 000 HO6 28 0.410 CH3 0.180
97 CH1 7 000 C5 29 0.376
98 OA 7 000 O5 29 -0.480
99 CH1 7 000 C1 29 0.232
100 OA 7 000 O1 29 -0.360
101 CH1 8 000 C4 29 0.232
102 CH1 8 000 C3 30 0.232 CH1 0.180
103 OA 8 000 O3 30 -0.642 OA -0.360
104 H 8 000 HO3 30 0.410 CH3 0.180
105 CH1 8 000 C2 31 0.232 CH1 0.180
106 OA 8 000 O2 31 -0.642 OA -0.360
107 H 8 000 HO2 31 0.410 CH3 0.180
108 CH2 8 000 C6 32 0.232 CH2 0.180
109 OA 8 000 O6 32 -0.642 OA -0.360
110 H 8 000 HO6 32 0.410 CH3 0.180
111 CH1 8 000 C5 33 0.376
112 OA 8 000 O5 33 -0.480
113 CH1 8 000 C1 33 0.232
114 OA 8 000 O1 33 -0.360
115 CH1 9 T000 C4 33 0.232
116 CH1 9 T000 C3 34 0.232 CH1 0.180
117 OA 9 T000 O3 34 -0.642 OA -0.360
118 H 9 T000 HO3 34 0.410 CH3 0.180
119 CH1 9 T000 C2 35 0.232 CH1 0.180
120 OA 9 T000 O2 35 -0.642 OA -0.360
121 H 9 T000 HO2 35 0.410 CH3 0.180
122 CH2 9 T000 C6 36 0.232 CH2 0.180
123 OA 9 T000 O6 36 -0.642 OA -0.360
124 H 9 T000 HO6 36 0.410 CH3 0.180
125 CH1 9 T000 C5 37 0.376
126 OA 9 T000 O5 37 -0.480
127 CH1 9 T000 C1 37 0.232
128 OA 9 T000 O1 37 -0.538
129 H 9 T000 HO1 37 0.410
[ bonds ]
; ai aj funct type
1 2 2 gb_1
5 6 2 gb_1 gb_19
8 9 2 gb_1 gb_19
11 12 2 gb_1 gb_19
2 3 2 gb_19
3 4 2 gb_25
3 13 2 gb_25
4 5 2 gb_19
4 7 2 gb_25
7 8 2 gb_19
7 15 2 gb_25
10 11 2 gb_19
10 13 2 gb_25
13 14 2 gb_19
14 15 2 gb_19
15 16 2 gb_19
16 17 2 gb_19
19 20 2 gb_1 gb_19
22 23 2 gb_1 gb_19
25 26 2 gb_1 gb_19
17 18 2 gb_25
17 27 2 gb_25
18 19 2 gb_19
18 21 2 gb_25
21 22 2 gb_19
21 29 2 gb_25
24 25 2 gb_19
24 27 2 gb_25
27 28 2 gb_19
28 29 2 gb_19
29 30 2 gb_19
30 31 2 gb_19
33 34 2 gb_1 gb_19
36 37 2 gb_1 gb_19
39 40 2 gb_1 gb_19
31 32 2 gb_25
31 41 2 gb_25
32 33 2 gb_19
32 35 2 gb_25
35 36 2 gb_19
35 43 2 gb_25
38 39 2 gb_19
38 41 2 gb_25
41 42 2 gb_19
42 43 2 gb_19
43 44 2 gb_19
44 45 2 gb_19
47 48 2 gb_1 gb_19
50 51 2 gb_1 gb_19
53 54 2 gb_1 gb_19
45 46 2 gb_25
45 55 2 gb_25
46 47 2 gb_19
46 49 2 gb_25
49 50 2 gb_19
49 57 2 gb_25
52 53 2 gb_19
52 55 2 gb_25
55 56 2 gb_19
56 57 2 gb_19
57 58 2 gb_19
58 59 2 gb_19
61 62 2 gb_1 gb_19
64 65 2 gb_1 gb_19
67 68 2 gb_1 gb_19
59 60 2 gb_25
59 69 2 gb_25
60 61 2 gb_19
60 63 2 gb_25
63 64 2 gb_19
63 71 2 gb_25
66 67 2 gb_19
66 69 2 gb_25
69 70 2 gb_19
70 71 2 gb_19
71 72 2 gb_19
72 73 2 gb_19
75 76 2 gb_1 gb_19
78 79 2 gb_1 gb_19
81 82 2 gb_1 gb_19
73 74 2 gb_25
73 83 2 gb_25
74 75 2 gb_19
74 77 2 gb_25
77 78 2 gb_19
77 85 2 gb_25
80 81 2 gb_19
80 83 2 gb_25
83 84 2 gb_19
84 85 2 gb_19
85 86 2 gb_19
86 87 2 gb_19
89 90 2 gb_1 gb_19
92 93 2 gb_1 gb_19
95 96 2 gb_1 gb_19
87 88 2 gb_25
87 97 2 gb_25
88 89 2 gb_19
88 91 2 gb_25
91 92 2 gb_19
91 99 2 gb_25
94 95 2 gb_19
94 97 2 gb_25
97 98 2 gb_19
98 99 2 gb_19
99 100 2 gb_19
100 101 2 gb_19
103 104 2 gb_1 gb_19
106 107 2 gb_1 gb_19
109 110 2 gb_1 gb_19
101 102 2 gb_25
101 111 2 gb_25
102 103 2 gb_19
102 105 2 gb_25
105 106 2 gb_19
105 113 2 gb_25
108 109 2 gb_19
108 111 2 gb_25
111 112 2 gb_19
112 113 2 gb_19
113 114 2 gb_19
114 115 2 gb_19
117 118 2 gb_1 gb_19
120 121 2 gb_1 gb_19
123 124 2 gb_1 gb_19
128 129 2 gb_1
115 116 2 gb_25
115 125 2 gb_25
116 117 2 gb_19
116 119 2 gb_25
119 120 2 gb_19
119 127 2 gb_25
122 123 2 gb_19
122 125 2 gb_25
125 126 2 gb_19
126 127 2 gb_19
127 128 2 gb_19
Thanks,
Jonathan
____________________________
Jonathan Moore, Ph.D.
Research and Engineering Sciences - New Products
Core R&D
The Dow Chemical Company
1702 Building, Office 4E
Midland, MI 48674 USA
Phone: (989) 636-9765
Fax: (989) 636-4019
E Mail: jmoore2 at dow.com
More information about the gromacs.org_gmx-users
mailing list