[gmx-users] calculation of the desolvation energy of a charged molecule through free energy perturbation
jin.jean.zhang at gmail.com
Thu May 8 00:09:04 CEST 2014
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
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
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,
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
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
More information about the gromacs.org_gmx-users