Issue with lambda = 1 during free energy calculations

User discussions
,
1

GROMACS version: 2021.1
GROMACS modification: No
Hi all,

I am running a very simple test case to understand free energy calculations. I have gone through and reproduced Justin’s tutorial for methane in water with no issues.

My simpler system is 2000 LJ particles (denoted “W” for water) and one other LJ particle (denoted “P”) (see appended topology file below). I am trying to calculate the free energy of solvation of this “P” particle in the bath of “W” particles by scaling the lambda values from 0 (vdw interactions present) to 1 (vdw interactions absent). All my simulations seem to be running fine EXCEPT lambda=1 where the system weirdly solidifies.

Here is my topology file:

[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 1.0 1.0

[ atomtypes ]
; name mass charge ptype V (sig) W (eps)
W 1 0.000 A 1.000000 1.000000
P 1 0.000 A 1.000000 0.500000

[ nonbond_params ]
; i j func V(sig) W(eps)
W W 1 1.000 1.000
P P 1 1.000 0.500
W P 1 1.000 1.000

[ moleculetype ]
; name nrexcl
W 1

[ atoms ]
; nr type resnr residue atom cgnr charge
1 W 1 W W1 1 0.000

[ moleculetype ]
; name nrexcl
P 1

[ atoms ]
; nr type resnr residue atom cgnr charge
1 P 1 P P1 1 0.000

[ system ]
Binary

[ molecules ]
; name N
W 2000
P 1

And here is my mdp file for the NPT equilibration (where I see the system solidifying within 2-3 ns):
integrator = md
dt = 0.005
nsteps = 1000000
nstcomm = 1
comm-mode = Linear
comm-grps = System

nstxout = 1000 ; no. steps between printing coords to trr file
nstvout = 0 ; no. steps between printing veloc to trr file
nstfout = 0 ; no. steps between printing forces to trr file
nstlog = 1000 ; no. steps between printing energy to log file
nstcalcenergy = 1
nstenergy = 1000 ; no. steps between printing energy to edr file
nstxout-compressed = 0 ; no. steps between printing coords to xtc file

cutoff-scheme = verlet
nstcalclr = 1
pbc = xyz
rlist = 3.0

vdwtype = Cut-off
rvdw = 2.5

coulombtype = Cut-off
rcoulomb-switch = 2.0
rcoulomb = 2.5

tcoupl = v-rescale
tc-grps = System
tau-t = 0.5
ref-t = 84.19

; Pressure coupling is on
pcoupl = berendsen ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com

; FREE ENERGY CONTROLS
free-energy = yes
init-lambda-state = 20 ; Final window, lambda = 1
delta-lambda = 0
calc-lambda-neighbors = 1 ; only immediate neighboring windows
couple-moltype = P ; name of moleculetype to decouple
couple-lambda0 = vdw ; only van der Waals interactions
couple-lambda1 = none ; turn off everything, in this case only vdW
couple-intramol = no

; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for each simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
vdw_lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00

; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none in this case)
sc-power = 1.0
sc-sigma = 1.0
nstdhdl = 500

gen-vel = yes
gen-temp = 84.19
gen-seed = -1

Any help is appreciated! I want to stress again that all the windows from 0-19 seem to be working perfectly. The basic system with 2000 “W” molecules is also perfectly fine independently BUT setting lambda=1 for my system gives me strange behavior.

Thanks,
Varun

2

At lambda=1 you solute is decoupled. You can an SD integrator for correct sampling close to lambda=1. That might solve your issue.

3

Thanks for your reply, this was exactly the issue. I figured this out yesterday and was going to update my answer today.

The core issue is that the ideal particle (at lambda = 1) can move ‘through’ the box at any velocity with nothing to impede it. So to maintain the overall temperature of the system, the other component must freeze! Very interesting issue.

An alternative solution that also works is setting up two tc-grps groups with separate coupling to the thermostat. Both this solution and the sd integrator work.

Incidentally, @Justin’s methane-in-water tutorial uses the sd integrator but does not explain why. It might be useful to add a note to that tutorial stressing the importance of this (seemingly innocuous) choice.

4

I think this is now explained in the manual, but I’m not 100% sure.

Didn’t grompp give you a warning when not using SD?

5

You’re right, grompp gave me this note that I missed because I wasn’t looking at the output:

NOTE 1 [file npt_equil.mdp, line 68]:
For proper sampling of the (nearly) decoupled state, stochastic dynamics
should be used

Thanks, this is helpful!

6

Ah, I thought this was a warning. A note one can easily miss.