[gmx-users] calculation of the desolvation energy of a charged molecule through free energy perturbation

[Jin Zhang]

Hi,

I'm doing fep simulations to get the desolvation energy of a small
molecules with +1 net charge. I'm new to free energy calculation, so any
suggestion and help would be appreciated. Thanks in advance.

To be more specific, I have several questions:

1.Since the molecule is charged, it will leave the system a net charge of
-1 when the coulombic interaction between the molecule and its surround is
fully decoupled. This will cause problem due to the use of PME. I
understand there's a number of discussions about this issue.
Here, I just decoupled the molecule (denotes as M) together with a clorine.
Then there're several ways to treat the interaction between M and the
clorine:

I. treat them as a single molecule (then the interaction between M and the
clorine become intramolecular interaction).
If the "couple-moltype" option is used without touching topology, to my
understanding, the M-clorine interaction would not be coupled, and the
(M-clorine)-surrounding interaction would be coupled. This is the direct
way. Is this correct?
If the "couple-moltype" option is used together with setting the to be
decoupled M-Cl charge and epsion to be 0 in topology, then the M
intramolecular, M-Cl intermolecular interactions would be decoupled
together with the (M-Cl) surrounding interaction.
Or add "couple-intramol=yes" without touching topology, then the 1-4
nonbond interaction within M would not be decoupled.

II. M and Cl can also be treated separately, I think. I did in this way:
use "couple-moltype=M" with blank type B, and edit Cl.itp to make type B
charge and epsion to be 0 (since clorine is simple, it doesn't has any
intramolecular interaction)

However, it crashes at lambda=0 during NVT simulation, which is followed by
two minimization processes.
It has been equilibrated before doing free energy calculation.
The error I got is like this:

A list of missing interactions:
        LJC Pairs NB of    428 missing      1
          exclusions of   6958 missing      1

Molecule type 'M'
the first 10 missing interactions, except for exclusions:
        LJC Pairs NB atoms   28   34           global    28    34
-------------------------------------------------------
Program mdrun, VERSION 4.6.5
Source code file: /home/apps/tarballs/gromacs-4.6.5/src/mdlib/domdec_top.c,
line: 393

Fatal error:
2 of the 9695 bonded interactions could not be calculated because some
atoms involved moved further apart than the multi-body cut-off distance
(0.715402 nm) or the two-body cut-off distance (1.2 nm), see option -rdd,
for pairs and tabulated bonds also see option -ddcheck


I found several ways can make it run, but I'm still confused.
Can anyone please explain a little bit to me?

free-energy=no                        works fine
use mdrun -pd rather than dd    works fine
use mdrun -rdd 1.25 or larger    works fine
play with cut-off settings           works fine

  the distance between atom 28 and  34 is slightly larger than 1.2nm
  the interesting thing is both increase and decrease cut-off can make it
work. I list the cut-off relate setting here, they different sets are
separated by ","
rlist                           = 1.2, 1.0, 1.6
; Electrostatics
coulombtype              = PME
rcoulomb                   = 1.2, 1.0, 1.0
; van der Waals
vdw-type                    = switch
rvdw-switch                = 1.0, 0.8, 1.4
rvdw                          = 1.1, 0.9, 1.5
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                    = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing             = 0.12, 0.10, 0.16