[gmx-users] anomalous free energy dgdl.xvg values every nstlist steps while using a twin-range cutoff

[Chris Neale]

Hello,

I believe that the free-energy code dgdl contribution from the 
twin-range cutoff is being calculated incorrectly in gromacs 4.0.3 (and 
probably other versions as well).

Specifically, I notice that the dgdl values spike at nstlist intervals. 
This can be seen directly from my dgdl.xvg and also from running 
g_analyze -ac.

I suspect that while the forces from the LJ interactions that are in the 
longer range of the twin-range are added every step while the dgdl value 
is modified to have nstlist times these contributions every nstlist step.

My belief that the forces are truley added every step comes from section 
4.6.3 of the gromacs 4 manual:

"In the neighbor list all interaction pairs that fall within rlist are 
stored. Furthermore, the interactions between
pairs that do not fall within rlist but do fall within 
max(rcoulomb,rvdw) are computed during NS, and the forces
and energy are stored separately, and added to short-range forces at 
every time step between successive NS."

Perhaps the long-range LJ component forces are actually applied to the 
particles every nstlist steps as a force multipled by nstlist, as is 
done in NAMD, but I am currently under the impression that this is not 
the case.

I am happy to provide more details and files if necessary, but hopefully 
this information is sufficient for a skilled coder to take a look and 
determine if this is indeed the case.

Thank you,
Chris.

Relevant .mdp options:

integrator          =  sd
energygrps          =  SOL DPC DPN        ; annihilated group must be 
separated
gen_seed            =  -1
comm_mode           =  linear
nstcomm             =  1
comm_grps           =  System
dt                    = 0.004
nstlist             =  5
ns_type             =  grid
pbc                 =  xyz
coulombtype         =  PME
rcoulomb            =  0.9
fourierspacing      =  0.12
pme_order           =  4
vdwtype             =  cut-off
rvdw_switch         =  0
rvdw                =  1.4
rlist               =  0.9
DispCorr            =  EnerPres     
Pcoupl              =  Berendsen          ; REMOVE_FOR_EM
pcoupltype          =  isotropic          ; REMOVE_FOR_EM
compressibility     =  4.5e-5             ; REMOVE_FOR_EM
ref_p               =  1.                 ; REMOVE_FOR_EM
tau_p               =  4.0                ; REMOVE_FOR_EM
tc_grps             =  System             ; REMOVE_FOR_EM
tau_t               =  1.0                ; REMOVE_FOR_EM
ld_seed             =  -1                 ; REMOVE_FOR_EM
ref_t               =  300.               ; REMOVE_FOR_EM
gen_temp            =  300.               ; REMOVE_FOR_EM
constraints         =  all-bonds          ; REMOVE_FOR_EM
constraint_algorithm=  lincs              ; REMOVE_FOR_EM
lincs-iter          =  1                  ; REMOVE_FOR_EM
lincs-order         =  6                  ; REMOVE_FOR_EM