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.) 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
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)?