[gmx-users] Charging free energies are unaffected by couple-intramol

[Maksim Mišin]

Dear Gromacs users,


I was trying to compute charging free energies of simple solutes such as
ions and dipoles using BAR.

In the manual it says that setting mdp option couple-intramol to no makes
decoupled state of the molecule correspond to the proper vacuum state
without periodicity effects. However, two runs with this option on and off
gave practically identical results:

couple-intramol = no
LJ particle:  -59.667  +-  0.037 kcal/mol
Dipole:        -35.172  +-  0.032 kcal/mol

couple-intramol = yes
LJ particle:  -59.729  +-  0.037 kcal/mol
Dipole:        -35.177  +-  0.032 kcal/mol

If I understand everything correctly (and this has been mentioned recently
on mail list:
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2015-July/099237.html)
these energies should differ by the charging free energy in vacuum, which
for these molecules is:

LJ particle: 18.12  kcal/mol
Dipole:        0.074 kcal/mol

I have also taken look at tpr files produced by grompp using gmxdump.
Except for different initial velocities tpr files are identical with either
option, which is somewhat strange.

Is this a bug or intended behaviour?


Parameters:
LJ particle: sigma = 3.8 A, eps = 0.125 kcal/mol
Dipole: same LJ parameters for both atoms with bond length = 2 A
Water model is SPC/E; cubic box with a length of 2.5 nm was used for both
solvent and vacuum simulations. Other free energy estimation methods such
as TI and MBAR gave very similar results.

Here is typical mdp file I used for these calculation:

; MD Run control
integrator               = sd
tinit                    = 0
dt                       = 0.002
nsteps                   = 2500000
nstcomm                  = 100

 ; Output control
nstxout                  = 0
nstvout                  = 0
nstfout                  = 0
nstlog                   = 500
nstenergy                = 500
nstxout-compressed       = 0

; Neighboursearch and short-range nondbonded interactions
cutoff-scheme            = verlet
nstlist                  = 40
ns_type                  = grid
pbc                      = xyz
rlist                    = 1.2

; Electrostatics
 coulombtype              = pme
rcoulomb                 = 1.2
epsilon-rf               = 73.5
; van der Waals
vdwtype                  = cutoff
vdw-modifier             = potential-switch
rvdw-switch              = 1.0
rvdw                     = 1.2
DispCorr                 = EnerPres
fourierspacing           = 0.15
; EWALD/PME/PPPM parameters
pme_order                = 6
ewald_rtol               = 1e-06
epsilon_surface          = 0

; Temperature coupling
tc_grps                  = system
tau-t                    = 1.0
ref-t                    = 298.15

; Presure coupling
Pcoupl                   = berendsen
tau_p                    = 1.0
compressibility          = 4.5e-05
ref_p                    = 1.0

 ; Free energy control
free_energy              = yes
init_lambda_state        = 0
delta_lambda             = 0
calc_lambda_neighbors    = -1
coul_lambdas      = 0.0 0.25 0.5 0.75 1.0

; Options for the decoupling
sc-alpha                 = 0.5
sc-coul                  = no
sc-power                 = 1.0
sc-sigma                 = 0.3
couple-moltype           = MOL
couple-lambda0           = vdw-q
couple-lambda1           = vdw
couple-intramol          = no  ; (or yes)
nstdhdl                  = 100

; Velocity generation
gen-vel                 = yes
gen-temp                = 298.15

; Bond constraint
constraints              = all-bonds
constraint-algorithm     = lincs
continuation             = no
lincs-order              = 12


Thank you,
Maksim Mišin