# Error met in free energy calculation **URL:** \ **Category:** User discussions\ **Tags:** mdp-parameters, mdrun\ **Created:** [April 25, 2023, 6:06am UTC](https://gromacs.bioexcel.eu/t/error-met-in-free-energy-calculation/6271 "2023-04-25T06:06:07Z")\ **Posts on this page:** 7\ **Page:** 1
**Author:** ![Wentao](https://avatars.discourse-cdn.com/v4/letter/w/a9adbd/32.png) [@Wentao](https://gromacs.bioexcel.eu/u/Wentao)\ **Post date:** [April 25, 2023, 6:06am UTC](https://gromacs.bioexcel.eu/t/error-met-in-free-energy-calculation/6271/1 "2023-04-25T06:06:08Z")
GROMACS version:2022.5 GROMACS modification: No I am trying to do free energy calculation on a protein-ligand system. I met such error when I doing the ligand part: There are perturbed non-bonded pair interactions beyond the pair-list cutoff of 1.26317 nm, which is not supported. This can happen because the system is unstable or because intra-molecular interactions at long distances are excluded. If the latter is the case, you can try to increase nstlist or rlist to avoid this.The error is likely triggered by the use of couple-intramol=no and the maximal distance in the decoupled molecule exceeding rlist Here is my mdp file for the particular lambda window: ;==================================================== ; Production simulation ;==================================================== ;---------------------------------------------------- ; RUN CONTROL ;---------------------------------------------------- integrator = sd ; stochastic leap-frog integrator nsteps = 5000000 ; 5000 ps dt = 0.001 ; 1 fs comm-mode = Linear ; remove center of mass translation nstcomm = 100 ; frequency for center of mass motion removal ;---------------------------------------------------- ; OUTPUT CONTROL ;---------------------------------------------------- nstxout = 0 ; don’t save coordinates to .trr nstvout = 0 ; don’t save velocities to .trr nstfout = 0 ; don’t save forces to .trr nstxout-compressed = 1000 ; xtc compressed trajectory output every 1000 steps (2 ps) compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file nstlog = 1000 ; update log file every 2 ps nstenergy = 1000 ; save energies every 2 ps nstcalcenergy = 100 ; calculate energies every 100 steps ;---------------------------------------------------- ; BONDS ;---------------------------------------------------- constraint\_algorithm = lincs ; holonomic constraints constraints = none ; hydrogens only are constrained lincs\_iter = 1 ; accuracy of LINCS (1 is default) lincs\_order = 4 ; also related to accuracy (4 is default) lincs-warnangle = 90 ; maximum angle that a bond can rotate before LINCS will complain (30 is default) continuation = yes ;---------------------------------------------------- ; NEIGHBOR SEARCHING ;---------------------------------------------------- cutoff-scheme = Verlet ns-type = grid ; search neighboring grid cells nstlist = 10 ; 20 fs (default is 10) rlist = 1.0 ; short-range neighborlist cutoff (in nm) pbc = xyz ; 3D PBC ;---------------------------------------------------- ; ELECTROSTATICS ;---------------------------------------------------- coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm) ewald\_geometry = 3d ; Ewald sum is performed in all three dimensions pme-order = 6 ; interpolation order for PME (default is 4) fourierspacing = 0.10 ; grid spacing for FFT ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb ;---------------------------------------------------- ; VDW ;---------------------------------------------------- vdw-type = PME rvdw = 1.0 vdw-modifier = Potential-Shift ewald-rtol-lj = 1e-3 lj-pme-comb-rule = Geometric DispCorr = EnerPres ;---------------------------------------------------- ; TEMPERATURE & PRESSURE COUPL ;---------------------------------------------------- tc\_grps = System tau\_t = 1.0 ref\_t = 300 pcoupl = Parrinello-Rahman pcoupltype = isotropic ; uniform scaling of box vectors tau\_p = 2 ; time constant (ps) ref\_p = 1.0 ; reference pressure (bar) compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1) ;---------------------------------------------------- ; VELOCITY GENERATION ;---------------------------------------------------- gen\_vel = no ; Velocity generation is off (if gen\_vel is ‘yes’, continuation should be ‘no’) gen\_seed = -1 ; Use random seed gen\_temp = 300 ;---------------------------------------------------- ; FREE ENERGY CALCULATIONS ;---------------------------------------------------- free-energy = yes couple-moltype = ligand couple-lambda0 = vdw-q couple-lambda1 = none couple-intramol = yes separate-dhdl-file = yes sc-alpha = 0.5 sc-power = 1 sc-sigma = 0.3 init-lambda-state = 12 coul-lambdas = 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 vdw-lambdas = 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 nstdhdl = 100 calc-lambda-neighbors = -1 dhdl-print-energy = potential dhdl-derivatives = yes I removed the LINCS since it always cause another error. One thing that is really confusing for me is that when I turn on the LINCS, the mdrun for lambda from 0 to 11 and 19 are all okay. The LINCS and such cut-off problem only occur in lambda 12 to 18. I don’t know how to solve such problem. Hope to get some suggestions. Thanks! ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [April 25, 2023, 6:24am UTC](https://gromacs.bioexcel.eu/t/error-met-in-free-energy-calculation/6271/2 "2023-04-25T06:24:00Z")
You need to have constraints on h-bonds. Integration with all bonds flexible is unstable with dt=0.001. And the force field is likely parameterized with constraints on h-bonds. What error do you get with constraints? ---
**Author:** ![Wentao](https://avatars.discourse-cdn.com/v4/letter/w/a9adbd/32.png) [@Wentao](https://gromacs.bioexcel.eu/u/Wentao)\ **Post date:** [April 25, 2023, 6:50am UTC](https://gromacs.bioexcel.eu/t/error-met-in-free-energy-calculation/6271/3 "2023-04-25T06:50:35Z")
Thanks for replying. The problem I met when I turn on LINCS with ds=0.002 is “bond that rotates more than 30/60/90 degree”. I build this ligand by using LigParGen and this model performed well in a typical simulation. So I think the ligand model I build is still okay. That’s the reason I tried turn off the LINCS. ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [April 25, 2023, 7:01am UTC](https://gromacs.bioexcel.eu/t/error-met-in-free-energy-calculation/6271/4 "2023-04-25T07:01:06Z")
You should never turn of constraints (except for problematic energy minimization). What do you mean with “performed well in a typical simulation”? Does the same topology run without LINCS warnings in a non free-energy simulation? ---
**Author:** ![Wentao](https://avatars.discourse-cdn.com/v4/letter/w/a9adbd/32.png) [@Wentao](https://gromacs.bioexcel.eu/u/Wentao)\ **Post date:** [April 25, 2023, 7:11am UTC](https://gromacs.bioexcel.eu/t/error-met-in-free-energy-calculation/6271/5 "2023-04-25T07:11:11Z")
The same topology run well in non free-energy simulation, no matter it is with or without LINCS. And, of course, I adopt the simulation with LINCS in that case. This cutoff problem can also appear in free-energy simulation with LINCS and ds=0.002, based on my latest simulation task. ---
**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\ **Post date:** [April 25, 2023, 7:46am UTC](https://gromacs.bioexcel.eu/t/error-met-in-free-energy-calculation/6271/6 "2023-04-25T07:46:09Z")
Some ligand force fields require constraints on all bonds. But I don’t know if that is the case for LigParGen. ---
**Author:** ![Wentao](https://avatars.discourse-cdn.com/v4/letter/w/a9adbd/32.png) [@Wentao](https://gromacs.bioexcel.eu/u/Wentao)\ **Post date:** [November 20, 2023, 12:06pm UTC](https://gromacs.bioexcel.eu/t/error-met-in-free-energy-calculation/6271/7 "2023-11-20T12:06:22Z")
Some updates: For an identical system, I successfully ran the entire decoupling process in the 2018.8 version of gmx, also using totally the same mdp settings. But I met the same problem with this: [How to solve Lincs warning in free energy calculation](https://gromacs.bioexcel.eu/t/how-to-solve-lincs-warning-in-free-energy-calculation/92) in 2023.1 version. (particularly from init-lamda-state 11 to init-lambda 18, simulation can run well in the rest of lambda states.) I successfully ran the example of gmx tutorial (the methane-water [Free Energy Calculations](http://www.mdtutorials.com/gmx/free_energy/index.html)) in gromacs version 2023.1. So I guess it’s the topology of the ligand. I tried to use LigParGen (OPLS/AA) and charmm-gui (CHARMM36) to build the topology, but they both encountered the above problem ([How to solve Lincs warning in free energy calculation](https://gromacs.bioexcel.eu/t/how-to-solve-lincs-warning-in-free-energy-calculation/92)) in such calculation.