[gmx-users] Possible free energy bug?

Justin A. Lemkul jalemkul at vt.edu
Thu Mar 10 19:55:54 CET 2011

Hi All,

I've been troubleshooting a problem for some time now and I wanted to report it 
here and solicit some feedback before I submit a bug report to see if there's 
anything else I can try.

Here's the situation: I ran some free energy calculations (thermodynamic 
integration) a long time ago using version 3.3.3 to determine the hydration free 
energy of a series of small molecules.  Results were good and they ended up as 
part of a paper, so I'm trying to reproduce the methodology with 4.5.3 (using 
BAR) to see if I understand the workflow completely.  The problem is my systems 
are crashing.  The runs simply stop randomly (usually within a few hundred ps) 
with lots of LINCS warnings and step*.pdb files being written.

I know the parameters are good, and produce stable trajectories, since I spent 
months on them some years ago. The system prep is steepest descents EM to Fmax < 
100 (always achieved), NVT at 298 K for 100 ps, NPT at 298K/1 bar for 100 ps, 
then 5 ns of data collection under NPT conditions.  Here's the rundown of what 
I'm seeing:

1. All LJ transformations work fine.  The problem only comes when I have a 
molecule with full LJ interaction and I am "charging" it (i.e., introducing 
charges to the partially-interacting species).

2. Simulations at lambda=1 (full interaction) work fine.

3. Simulations with the free energy code off entirely work fine under all 

4. I cannot run in serial due to http://redmine.gromacs.org/issues/715.  The bug 
seems to affect other systems and is not specifically related to my free energy 

5. Running with DD fails because my system is relatively small (more on this in 
a moment).

6. Running with mdrun -pd 2 works, but mdrun -pd 4 crashes for any value of 
lambda != 1.

7. I created a larger system (instead of a 3x3x3-nm cube of water with my 
molecule, I used 4x4x4) and ran on 4 CPU's with DD (lambda = 0, i.e. full vdW, 
no intermolecular Coulombic interactions - .mdp file is below).  This run also 
crashed with some warnings about DD cell size:

DD  load balancing is limited by minimum cell size in dimension X
DD  step 329999  vol min/aver 0.748! load imb.: force 31.5%

...and then the actual crash:

Program mdrun_4.5.3_gcc_mpi, VERSION 4.5.3
Source code file: domdec_con.c, line: 693

Fatal error:
DD cell 0 0 0 could only obtain 14 of the 15 atoms that are connected via 
constraints from the neighboring cells. This probably means your constraint 
lengths are too long compared to the domain decomposition cell size. Decrease 
the number of domain decomposition grid cells or lincs-order or use the -rcon 
option of mdrun.
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

Watching the trajectory doesn't seem to give any useful information.  The small 
molecule of interest is at a periodic boundary when the crash happens, but there 
are several crosses prior to the crash without incident, so I don't know if the 
issue is related to PBC or not, but it appears not.

8. I initially thought the problem might be related to the barostat, but 
switching from P-R to Berendsen does not alleviate the problem, nor does 
increasing tau_p (tested 0.5, 1.0, 2.0, and 5.0 - all crash).  Longer tau_p 
simply delays the crash, but does not prevent it.

So after all that, I'm wondering if (1) anyone has seen the same, or (2) if 
there's anything else I can try (environment variables, hidden tricks, etc) that 
I can use to get to the bottom of this before I give up and file a bug report.

If you made it this far, thanks for reading my novel and hopefully someone can 
give me some ideas.  The .mdp file I'm using is below, but it is just one of 
many that I've tried.  In theory, it should work, since the parameters are the 
same as my successful 3.3.3 runs, with the exception of the new free energy 
features in 4.5.3 and obvious keyword changes related to the difference in version.


--- .mdp file ---

; Run control
integrator               = sd       ; Langevin dynamics
tinit                    = 0
dt                       = 0.002
nsteps                   = 2500000  ; 5 ns
nstcomm                  = 100
; Output control
nstxout                  = 500
nstvout                  = 500
nstfout                  = 0
nstlog                   = 500
nstenergy                = 500
nstxtcout                = 0
xtc-precision            = 1000
; Neighborsearching and short-range nonbonded interactions
nstlist                  = 5
ns_type                  = grid
pbc                      = xyz
rlist                    = 0.9
; Electrostatics
coulombtype              = PME
rcoulomb                 = 0.9
; van der Waals
vdw-type                 = cutoff
rvdw                     = 1.4
; 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                = 4
ewald_rtol               = 1e-05
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                    = 298
; Pressure coupling is on for NPT
Pcoupl                   = Berendsen
tau_p                    = 2.0
compressibility          = 4.5e-05
ref_p                    = 1.0
; Free energy control stuff
free_energy              = yes
init_lambda              = 0.00
delta_lambda             = 0
foreign_lambda           = 0.05
sc-alpha                 = 0
sc-power                 = 1.0
sc-sigma                 = 0
couple-moltype           = MOR      ; name of moleculetype to couple
couple-lambda0           = vdw      ; vdW interactions
couple-lambda1           = vdw-q    ; turn on everything
couple-intramol          = no
dhdl_derivatives         = yes      ; this line (and the next two) are defaults
separate_dhdl_file       = yes      ; included only for pedantry
nstdhdl                  = 10
; Do not generate velocities
gen_vel                  = no
; options for bonds
constraints              = all-bonds
; 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              = 4


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080


More information about the gromacs.org_gmx-users mailing list