[gmx-users] Testing GROMACS: Solvation free energy tests using CGENFF

[Fahimeh Baftizadeh]

Hello everyone!

We had some discussion with JUSTIN and MARK but I dont know whether they
received my email with this details or not! So I am opening a new email for
these tests that I have done and I would like to receive your feedbacks
cause I am really lost!

I am computing solvation free energy of paracetamol in ethanol using
Cgenff. I have tested this force field for pure paracetamol and pure
ethanol and it seems good enough. (heat of sublimation, fusion and
vaporization is comparable with experimental values)! However, I have doubt
that this force field is good enough for the mixing properties of
paracetamol and ethanol. In any case I am not doing something very
complicated just look below:

TEST 1) I performed thermodynamics integration, with 20 lambda points. I
decoupled first Q and this is the one of the mdp files for production runs:
integrator               = sd       ; Langevin dynamics
pbc                      = xyz
rlist                    = 1.0
; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.0
; van der Waals
vdw-type                 = switch
rvdw-switch              = 0.8
rvdw                     = 0.9
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                  = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.12
; EWALD/PME/PPPM parameters
pme_order                = 6
ewald_rtol               = 1e-06
epsilon_surface          = 0
optimize_fft             = no
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps                  = system
tau_t                    = 1.0
ref_t                    = 300
; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman
tau_p                    = 0.5
compressibility          = 4.5e-05ref_p                    = 1.0
; Free energy control stuff
free_energy              = yes
init_lambda              = 0.00
delta_lambda             = 0
foreign_lambda           = 0.05
sc-alpha                 = 0.5
sc-power                 = 1.0
sc-sigma                 = 0.3
couple-moltype           = A  ; name of moleculetype to decouple
couple-lambda0           = vdw-q      ; only van der Waals interactions
couple-lambda1           = vdw     ; turn off everything, in this case only
vdW
couple-intramol          = no
nstdhdl                  = 10
; Do not generate velocities
gen_vel                  = no
; options for bonds
constraints              = all-bonds  ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm     = lincs
; Constrain the starting configuration
; since we are continuing from NPT
continuation             = yes
; Highest order in the expansion of the constraint coupling matrix
;lincs-order              = 12
lincs-order              = 4

The result of g_bar for this part was:
point  0.000 -  0.050,   DG  5.51 +/-  0.04
point  0.050 -  0.100,   DG  5.09 +/-  0.04
.
.
point  0.900 -  0.950,   DG  0.23 +/-  0.01
point  0.950 -  1.000,   DG  0.10 +/-  0.01

total  0.000 -  1.000,   DG 45.05 +/-  0.09

For decoupling VDW, the charges in the top file was set to zero.  and the
mdp file was the same as above, except these two lines:

couple-lambda0           = vdw      ; only van der Waals interactions
couple-lambda1           = none     ; turn off everything, in this case
only vdW

The results of the g_bar of this part is:
point  0.000 -  0.050,   DG  2.17 +/-  0.01
point  0.050 -  0.100,   DG  2.25 +/-  0.01
.
.
point  0.900 -  0.950,   DG -0.37 +/-  0.02
point  0.950 -  1.000,   DG  0.88 +/-  0.01

total  0.000 -  1.000,   DG 29.73 +/-  0.06

So I concluded that the total DG=45.05+ 29.73 = 74.78

TEST2) decoupling both VDW and Q at the same time. For this simulation, the
charges were the real charges and the mdp file had two modified lines:

couple-lambda0           = vdw-q      ; only van der Waals interactions
couple-lambda1           = none     ; turn off everything, in this case
only vdW

The results of this simulation was:
point  0.000 -  0.050,   DG 15.17 +/-  0.16
point  0.050 -  0.100,   DG  6.94 +/-  0.04
.
.
point  0.900 -  0.950,   DG  2.18 +/-  0.02
point  0.950 -  1.000,   DG  3.53 +/-  0.02

total  0.000 -  1.000,   DG 74.51 +/-  0.26

Which is more or less like before!

TEST 3: I performed another simulation, decoupling Q and VDW at the same
time by 40 lambda points:

point  0.000 -  0.025,   DG  9.75 +/-  0.02
point  0.025 -  0.050,   DG  5.66 +/-  0.01
.
.
point  0.950 -  0.975,   DG  1.60 +/-  0.00
point  0.975 -  1.000,   DG  1.93 +/-  0.00

total  0.000 -  1.000,   DG 74.68 +/-  0.06

They are all giving me same results!
Should I be worry or this is fine to receive the same results?

Thanks
Fahimeh

==============================
Fahimeh Baftizadeh, Ph.D.
MIT, Department of Chemical engineering
E19-528, 50 Ames street, Cambridge
MA, 02139

fahimeh.baftizadeh at gmail.com
fbafti at mit.edu
================================