# Dt in free energy perturbation calculations **URL:** \ **Category:** User discussions\ **Tags:** mdp-parameters, free-energy\ **Created:** [March 25, 2021, 3:03pm UTC](https://gromacs.bioexcel.eu/t/dt-in-free-energy-perturbation-calculations/1824 "2021-03-25T15:03:16Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![raitis](https://avatars.discourse-cdn.com/v4/letter/r/c77e96/32.png) [@raitis](https://gromacs.bioexcel.eu/u/raitis)\ **Post date:** [March 25, 2021, 3:03pm UTC](https://gromacs.bioexcel.eu/t/dt-in-free-energy-perturbation-calculations/1824/1 "2021-03-25T15:03:16Z")
GROMACS version: 2018.2 GROMACS modification: No Dear all, I am trying to set-up FEP calculations, to calculate the relative binding energy of congeneric ligands. I used FESetup tool ([GitHub - CCPBioSim/fesetup: A tool for setting up free energy simulations.](https://github.com/CCPBioSim/fesetup)) to automate the input preparation. The prepared files seems to be fine, however, I am having some difficulties with MD simulations in Gromacs. When I use dt = 0.002 I get warning > WARNING 1 [file morph.top, line 69595]: > The bond in molecule-type LIG between atoms 45 DU45 and 47 DU47 has an > estimated oscillational period of 7.9e-03 ps, which is less than 5 times > the time step of 2.0e-03 ps. > Maybe you forgot to change the constraints mdp option. When I use dt = 0.001, it downgrades to note: > NOTE 1 [file morph.top, line 69595]: > The bond in molecule-type LIG between atoms 45 DU45 and 47 DU47 has an > estimated oscillational period of 7.9e-03 ps, which is less than 10 times > the time step of 1.0e-03 ps. > Maybe you forgot to change the constraints mdp option. The corresponding part in the pert.itp is: > … > 39 ha 1 LIG H39 39 0.137430 1.0080 ha 0.160978 1.0080 > 40 ha 1 LIG H40 40 0.132430 1.0080 du 0.000000 1.0000 > 41 ha 1 LIG H41 41 0.131930 1.0080 du 0.000000 1.0000 > 42 ha 1 LIG H42 42 0.132430 1.0080 ha 0.156978 1.0080 > 43 ha 1 LIG H43 43 0.137430 1.0080 ha 0.144978 1.0080 > 44 du 1 LIG DU44 44 0.000000 1.0000 os -0.333922 16.0000 > 45 du 1 LIG DU45 45 0.000000 1.0000 c3 0.321178 12.0100 > 46 du 1 LIG DU46 46 0.000000 1.0000 os -0.333922 16.0000 > 47 du 1 LIG DU47 47 0.000000 1.0000 h2 0.072678 1.0080 > 48 du 1 LIG DU48 48 0.000000 1.0000 h2 0.072678 1.0080 The dummy atom parameters in the pert.atp are: > du du 1.000 0.0000 A 0.000000e+00 0.000000e+00 What is the most appropriate way to solve this issue? Is there a way to run simulation using dt = 0.002 (and get accurate results)? Looking forward to hear some suggestion on what might be wrong with my FEP calculations, Raitis ---
**Author:** ![alevilla](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/alevilla/32/439_2.png) [@alevilla](https://gromacs.bioexcel.eu/u/alevilla)\ **Post date:** [March 26, 2021, 12:33pm UTC](https://gromacs.bioexcel.eu/t/dt-in-free-energy-perturbation-calculations/1824/2 "2021-03-26T12:33:39Z")
Hi, If it suitable for your case, you try to use `constraints` on the bonds (maybe only bonds involving hydrogen atoms). See the corresponding option in mdp file [Molecular dynamics parameters (.mdp options) — GROMACS 2021.1 documentation](https://manual.gromacs.org/current/user-guide/mdp-options.html#bonds). Best regards Alessandra ---
**Author:** ![raitis](https://avatars.discourse-cdn.com/v4/letter/r/c77e96/32.png) [@raitis](https://gromacs.bioexcel.eu/u/raitis)\ **Post date:** [April 1, 2021, 12:13pm UTC](https://gromacs.bioexcel.eu/t/dt-in-free-energy-perturbation-calculations/1824/3 "2021-04-01T12:13:45Z")
Dear Alevilla, Thank you for the suggestion. I checked the **constraint** options and setting it to **all-bonds** allows to run the simulation using **dt = 0.002**. However, it fails rather soon with LINCS warnings. I guess, it is not surprising. Taking this in to account, I assume there is no good way to run simulations using dt larger than 0.001 in FEP calculations (at least in my case). Since I have established that I need to use **dt = 0.001** for my system, is it better to run my production simulation with or withut any constraints? Here I need to mention that running with **dt = 0.001** and **constraints = all-bonds** , simulation fails rather soon. What are the implications of running with **constraints = h-bonds** instead of **constraints = none**? And just to be sure, does mdp file like this seems reasonable for single topology perturbation simulation (coule of atoms disapearing, and coulpe appearing)? > integrator = sd > dt = 0.001 > nsteps = 5000000 > nstcomm = 100 > comm-mode = Linear > ;---------------------------------------------------- > ; OUTPUT CONTROL > ;---------------------------------------------------- > nstxout = 0 > nstvout = 0 > nstfout = 0 > nstxout-compressed = 10000 > compressed-x-precision = 1000 > nstlog = 10000 > nstenergy = 5000 > nstcalcenergy = 100 > ;---------------------------------------------------- > ; TEMPERATURE & PRESSURE COUPL > ;---------------------------------------------------- > tcoupl = no > nsttcouple = 10 > tc\_grps = System > tau\_t = 1.0 > ref\_t = 298.0 > > pcoupl = Parrinello-Rahman > pcoupltype = isotropic > tau\_p = 10 > ref\_p = 1.0 > compressibility = 4.5e-05 > refcoord-scaling = com > ;---------------------------------------------------- > ; BONDS > ;---------------------------------------------------- > constraints = none > constraint\_algorithm = Lincs > lincs\_order = 4 > lincs\_warnangle = 30 > ;---------------------------------------------------- > ; NEIGHBOR SEARCHING > ;---------------------------------------------------- > cutoff-scheme = Verlet > nstlist = 20 > ns\_type = grid > pbc = xyz > rlist = 1.2 > ;---------------------------------------------------- > ; ELECTROSTATICS > ;---------------------------------------------------- > coulombtype = PME > ;coulomb-modifier = none > rcoulomb = 1.2 > fourierspacing = 0.12 > pme\_order = 6 > ewald\_rtol = 1.0E-6 > ;ewald\_geometry = 3d > epsilon\_surface = 0 > ;---------------------------------------------------- > ; VDW > ;---------------------------------------------------- > vdw-type = cut-off > vdw-modifier = potential-switch > rvdw-switch = 1.0 > rvdw = 1.2 > DispCorr = EnerPres > ;---------------------------------------------------- > ; FREE ENERGY CALCULATIONS > ;---------------------------------------------------- > free-energy = yes > delta-lambda = 0 > nstdhdl = 100 > calc-lambda-neighbors = -1 > sc-alpha = 0.5 > sc-coul = no > sc-power = 1 > sc-r-power = 6 > sc-sigma = 0.3 > dhdl-derivatives = yes > dhdl-print-energy = no > separate-dhdl-file = no > dh\_hist\_size = 0 > dh\_hist\_spacing = 0.1 > couple-lambda0 = vdw-q > couple-lambda1 = vdw-q > couple-intramol = no > init-lambda-state = Lambda\_here > ;---------------------------------------------------- > ; LAMBDAS > ;---------------------------------------------------- > ; lambda states 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 > fep\_lambdas = 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 Best regards, Raitis