[gmx-users] maintaining stability when uncoupling ion pairs

[andrew biedermann]

Hello,

I want to calculate the free energy difference associated with decoupling a
NaCl pair from aqueous solution, but am running into a strange error. I
believe the issue is in how I define the molecule type in the topology
file, since the simulation is stable when removing individual ions, or
during regular production runs with the same mdp settings. I get the
following error soon after submitting the job (5000 steps in):

A list of missing interactions:
        LJC Pairs NB of      1 missing      1
          exclusions of  16669 missing      1

Molecule type 'tNC'
the first 10 missing interactions, except for exclusions:
        LJC Pairs NB atoms    1    2           global 11113 11114

-------------------------------------------------------
Program mdrun, VERSION 5.0.4
Source code file:
/nv/blue/amb4ht/downloads/gromacs-5.0.4/src/gromacs/mdlib/domdec_top.c,
line: 394

Fatal error:
2 of the 19448 bonded interactions could not be calculated because some
atoms involved moved further apart than the multi-body cut-off distance
(1.74968 nm) or the two-body cut-off distance (1.74968 nm), see option
-rdd, for pairs and tabulated bonds also see option -ddcheck
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------

I have modified my topology file with the following lines:

[ moleculetype ]
; molname    nrexcl
tNC          1

[ atoms ]
; id   at type     res nr   residu name  at name  cg nr  charge
1      Na          1        NA           NA        1      1.00000
2      Cl           1        CL           CL        2     -1.00000

[ bonds ]

[ angles ]

And my mdp file is:

; Run control
integrator               = sd       ; Langevin dynamics
tinit                    = 0
dt                       = 0.002
nsteps                   = 5000000   ; 10 ns
; Output control
nstxout                  = 2500
nstvout                  = 2500
nstfout                  = 0
nstlog                   = 2500
nstenergy                = 2500
nstxout-compressed       = 0
 Bond parameters
constraint_algorithm    = lincs        ; holonomic constraints
;constraints                = h-bonds  ; all bonds (even heavy atom-H
bonds) constrained
;MRS: use h-bonds
lincs_iter                  = 2        ; accuracy of LINCS
lincs_order                 = 12       ; also related to accuracy
lincs_warnangle = 30
continuation             = yes
; Neighborsearching and short-range nonbonded interactions
; Neighborsearching
cutoff-scheme   = Group    ; Group works with free energy
rlist = 1.2  ;  cut-off distance for the short-range neighbor list
ns_type             = grid              ; search neighboring grid cells
nstlist             = 10
; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.2
; van der Waals
vdw-type                 = Cutoff
vdw-modifier             = Potential-switch
rvdw-switch              = 0.90
rvdw                     = 0.95
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                  = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.10
; EWALD/PME/PPPM parameters
pme_order                = 6
ewald_rtol               = 1e-06
epsilon_surface          = 0
; Temperature coupling
tc-grps         = System
tau_t           = 2.0           ; time constant, in ps
ref_t           = 298           ; reference temperature in K
; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman
tau_p                    = 2.0
compressibility          = 4.5e-05
ref_p                    = 1.0
;############### Free energy control stuff ##############################
free_energy              = yes
init_lambda_state        = 0
delta_lambda             = 0
calc_lambda_neighbors    = 1        ; only immediate neighboring windows
; 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.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
coul_lambdas             = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80
0.90 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
; We are not transforming any bonded or restrained interactions
bonded_lambdas           = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint_lambdas        = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Masses are not changing (particle identities are the same at lambda = 0
and lambda = 1)
mass_lambdas             = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
temperature_lambdas      = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha                 = 0.5
sc-coul                  = no       ; linear interpolation of Coulomb
sc-power                 = 1.0
sc-sigma                 = 0.3
couple-moltype           = tNC ; name of moleculetype to decouple
couple-lambda0           = vdw-q      ; all interactions at lambda 0
couple-lambda1           = none     ; turn off everything at lambda 1
couple-intramol          = no
nstdhdl                  = 100

Thanks for the help,

Drew