[gmx-users] Error in free energy calculation

[Sadaf Rani]

Dear Gromacs users
I am running a free energy calculation, my simulation crashes and generates
different PDB structures  I have well equilibrated the system before going
to the production and the system seems stable. My input file for production
run is as below:-
; Run control
integrator               = sd       ; stochastic leap-frog integrator
dt                       = 0.002
nsteps                   = 1500000   ; 3 ns
comm-mode     = Linear        ; remove center of mass translation
nstcomm     = 100           ; frequency for center of mass motion removal


; Output control
nstxout                  = 5000
nstvout                  = 5000
nstfout                  = 0
nstlog                   = 5000
nstenergy                = 5000
nstxout-compressed       = 5000


; Bond parameters
continuation            = yes        ; formerly known as
'unconstrained-start' - useful for exact continuations and reruns
constraint_algorithm    = lincs     ; holonomic constraints
constraints             = h-bonds   ; bonds to H are constrained
lincs_iter              = 1         ; accuracy of LINCS

; Neighborsearching and short-range nonbonded interactions
cutoff-scheme            = verlet
nstlist                  = 20
ns_type                  = grid
rlist                    = 1.2
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC

; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.2
pme-order       = 4        ; interpolation order for PME (default is 4)
fourierspacing   = 0.16     ; grid spacing for FFT
ewald-rtol       = 1e-6     ; relative strength of the Ewald-shifted direct
potential at rcoulomb
ewald_geometry   = 3d       ; Ewald sum is performed in all three dimensions


; van der Waals
vdwtype                  = cutoff
vdw-modifier             = Potential-shift-Verlet
verlet-buffer-tolerance = 2.5e-04
rvdw                     = 1.2
DispCorr                 = EnerPres     ; apply long range dispersion
corrections for Energy and Pressure


; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc-grps                 = system ; two coupling groups - more accurate
tau_t                   = 0.1                       ; time constant, in ps
ref_t                   = 298.15

; Pressure coupling is on for NPT
Pcoupl                  = Parrinello-Rahman
pcoupltype       = isotropic            ; uniform scaling of box vectors
tau_p           = 2.0                  ; time constant (ps)
ref_p             = 1.0                  ; reference pressure (bar)
compressibility   = 4.5e-05              ; isothermal compressibility of
water (bar^-1)
; velocities
gen_vel      = no      ; Velocity generation is on (if gen_vel is 'yes',
continuation should be 'no')
gen_seed     = -1       ; Use random seed
gen_temp     = 298.15
; Free energy control stuff
free_energy              = yes
init_lambda_state        = 10
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
 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.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
coul_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
restraint_lambdas        = 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
; Options for the decoupling
sc-alpha                 = 0.5
sc-coul                  = no       ; linear interpolation of Coulomb (none
in this case)
sc-power                 = 1
sc-sigma                 = 0.3
nstdhdl                  = 100
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 12


Any help would really be appreciated.

Thanks.
Sadaf