[gmx-users] LINCS warnings for some lambdas.

[Sana Saeed]

hii am performing protein ligand (pdb:4HBV)binding free energy calculations, i am done with complex simulation , now stuck in ligand simulation since last 2 weeks. first i was getting zero energy, and when i changed some of parameters and position restraint my ligand , now i am getting LINCS warnings sometimes with NPT  and sometimes during production md. the strange part is that some of the lambdas are perfectly running while some of them crashes due to LINCS warnings, i know how to fix LINCS warnings (extending table size, decreasing stepsize, properly minimizing structure etc) but in  this case some of lambdas are running perfectly, which means the structure and parameters are fine. i decresed step size to 1 fs and got the same result, bow i decreased it to 0.0001 ps. if someone ever got this kind of problem , please share your views. my second question is that the simple simulation for solvation free energy of small molecule and the ligand simulation in binding free energy system are same? i mean can i use the solvation free energy parameters from BEVAN lab's tutorial in this situation? in some of tutorials they use position restraints for ligand simulation in water, that is the reason i used the same. i provide npt and md.mdp parameters. if something is wrong kindly give comments. ignore the stepsize , i know its too small, i tried 0.001, 0.002, and this is the recent one, but couldnt see any difference.thanks in advance 
NPT.mdp
define              =  -DPOSRESintegrator          =  sdtinit               =  0dt                  =  0.0001   nsteps              =  500000  ; 100 psnstcomm             =  100
; Boundary conditionspbc                 =  xyz
; Outputnstxout             =  0nstvout             =  0nstfout             =  0nstlog              =  100   ; 10 psnstenergy           =  100   ; 10 psnstcalcenergy       =  20nstxtcout           =  100   ; 10 psxtc-precision       =  1000energygrps          =  System
; Neighbour searchingcutoff-scheme       = Verletnstlist             = 10      ; 20 fsns-type             = gridrlist               = 1.2               ; short-range neighborlist cutoff (in nm)rcoulomb            = 1.2               ; short-range electrostatic cutoff (in nm)rvdw                = 1.2               ; short-range van der Waals cutoff (in nm)
; Constraintsconstraints         =  noneconstraint-algorithm = lincslincs_iter          = 1         ; accuracy of LINCS (1 is default)lincs_order         = 12        ; also related to accuracy (4 is default)lincs-warnangle     = 30        ; maximum angle that a bond can rotate before LINCS will complain^M; Electrostaticscoulombtype     = PME           ; Particle Mesh Ewald for long-range electrostaticspme-order       = 6             fourierspacing  = 0.10          ; grid spacing for FFTewald-rtol      = 1e-6          ; relative strength of the Ewald-shifted direct potential at rcoulomboptimize-fft    = noewald_geometry   = 3d  
; van der Waal'svdwtype            =  Switchrvdw-switch         =  0.9DispCorr            =  EnerPres
;Temperature Coupling (SD integrator => Langevin dynamics)tc_grps             =  Systemtau_t               =  1.0ref_t               =  298.15
; Pressure couplingPcoupl              =  Berendsenpcoupltype          =  isotropictau_p               =  2compressibility     =  4.5e-05ref_p               =  1.0
; Generate velocitiesgen_vel             = yesgen_seed            = -1gen_temp            = 298.15
continuation        = no 
; Free energy control stuff = free-energy              = yessc-alpha                 = 0.5sc-power                 = 1sc-sigma                 = 0.3init-lambda-state        = 0coul-lambdas             = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0vdw-lambdas              = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0  couple-moltype          = 15E couple-intramol         = no couple-lambda0          = none couple-lambda1          = vdw-qnstdhdl                  = 100calc-lambda-neighbors    = -1refcoord-scaling         = com
PRODUCTION MD:md.mdp
; Production MD;title           = production
; Run parametersintegrator      = sd            ; stochastic leap-frog integratornsteps          = 2500000       dt              = 0.0001                nstcalcenergy   =  100energygrps      =  System
; Periodic boundary conditionspbc             = xyz           ; 3-D PBC
; Output controlnstvout         = 500        ; save velocities to .trr file every 200 psnstfout         = 500        ; save forces to .trr file every 100 psnstxout         = 500           ; save coordinates every 5 psnstxtcout       = 500           ; xtc compressed trajectory output every 2 psnstenergy       = 500           ; save energies every 2 psnstlog          = 500           ; update log file every 2 ps
; Bond parametersconstraint_algorithm = lincs    ; holonomic constraints constraints     = h-bonds       ; all bonds (even heavy atom-H bonds) constrainedlincs_iter      = 1             ; accuracy of LINCS (1 is default)lincs_order     = 6             ; also related to accuracy (4 is default)lincs-warnangle = 30            ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
; Neighbor searchingns-type         = grid          ; search neighboring grid cellsnstlist         = 5             ; 10 fsrlist           = 1.2           ; short-range neighborlist cutoff (in nm)rcoulomb        = 1.2           ; short-range electrostatic cutoff (in nm)rvdw            = 1.2           ; short-range van der Waals cutoff (in nm)
; Electrostaticscoulombtype     = PME           ; Particle Mesh Ewald for long-range electrostaticspme-order       = 6            fourierspacing  = 0.10          ; grid spacing for FFTewald-rtol      = 1e-6          ; relative strength of the Ewald-shifted direct potential at rcoulomb; Van der Waalsvdwtype        =  Switchrvdw-switch     =  0.9DispCorr        =  EnerPres     ; account for cut-off vdW scheme
;Temperature Coupling (SD integrator => Langevin dynamics)
tc_grps             =  Systemtau_t               =  1.0^Mref_t               =  298.15
; Pressure coupling is onpcoupl          = Parrinello-Rahmanpcoupltype      = isotropic     ; uniform scaling of box vectorstau_p           = 2             ; time constant, in psref_p           = 1.0           ; reference pressure, in barcompressibility = 4.5e-05       ; isothermal compressibility of water, bar^-1  
; Velocity generationgen_vel         = no            ; Velocity generation is ongen_seed        = -1            ; Use random seedgen_temp        = 298.15
continuation    = yes           ; Restarting after NPT
; Free energy control stuff free-energy              = yessc-alpha                 = 0.5sc-power                 = 1sc-sigma                 = 0.3sc-coul                  = noinit-lambda-state        = 0coul-lambdas             = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0vdw-lambdas              = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 couple-moltype          = 15E couple-intramol         = no couple-lambda0          = none couple-lambda1          = vdw-qnstdhdl                  = 100dhdl-print-energy        = yescalc-lambda-neighbors    = -1
 Sana Saeed Khan,Research AssistantChemoinformatics LabGraduate Student, MS bioinfoDepartment of BioinformaticsSoongsil University, Seoul, South Korea.