GROMACS version: 2021.3
GROMACS modification: Yes/No
Dear users,
I’m completely noob working with MD and I performed a simulation following the Lemkul tutorial recently. Now, I would like to apply, in the same protein-ligand system, some points I read in two papers, which are:
-
“At total, each system was equilibrated using integration steps of 1 fs through 2 ns following gradually decreasing ligand restraint force constants along 4 steps of 250 ps each.”
-
“The system was heated in cycles of 100 ps simulations with gradual temperature increases (10, 50, 150, and 300 K) using the Berendsen thermostat. Then, the pressure was adjusted to 1 atm using the Berendsen barostat in a 400 ps simulation. The 100 kcal/(mol Å2) restraints were gradually removed in cycles of 500 ps simulations (50, 10, 5, and 1 kcal/(mol Å2))”
Both papers showed that restraints were removed gradually during the equilibration steps. Thus, I have some questions:
1- How can I “gradually decrease ligand restraint force constants” during the equilibration step?
2- My idea about the first question is that each step uses the previous result decreasing the force - set 1000 in the first step - that is in the file “posre_ligand.ipt” (1000, 500, 100…) to start the simulation again, is it right?
3- The “gradual temperature increases (10, 50, 150, and 300 K)” is a normal characteristic of the thermostats or should I set any function to this gradual temperature increase?
I’d like to thank you for your time and answer
Sincerely,
José Xavier
You are correct that the method as described is to do successive run of 250/500 ps, changing the parameters for each of those.
For the thermostat: customarily, for the first simulation run after minimization, the velocities for each particle will be drawn from a Boltzmann distribution at the target temperature (see gen_vel = yes), so that the first step is not effectively at 0K (and, of course, the thermostat will act as usual for the run). For the subsequent run, in the described methodology, the temperature will be raised by the thermostat at a rate given by the tau_t parameter (see description of the thermostat algorithm for exact meaning of this time constant).
Just a note: the procedure as described is very gentle. While that cannot hurt, unless there is specific concern for the equilibration of this protein system, some might say it is unnecessarily so. Also you might want to read this for more details on equilibration (and other issues).
1 Like
Dear @ebriand,
Thank you for your help.
I’ll read carefully the information from the link you sent. It looks to be very complete.
About the answers:
You are correct that the method as described is to do successive run of 250/500 ps, changing the parameters for each of those.
According to Lemkul’s tutorial, I introduced the POSRES and the POSRES_LIG in the .top and .mdp files, as well as I merged protein and ligand as the Protein_Lig group in the tc-grps during the equilibration (NVT and NPT). If I understood it correctly, the functions POSRES and POSRES_LIG restrain the protein and the ligand based on the energy threshold in the .itp files, and the tc-grps = Protein_Lig makes both molecules a single entity during the calculation.
If I want to perform a third equilibration step (NPT or NVT), “decreasing gradually the ligand restrains”,
1- should I change something else than the energy value in the posre_ligand.itp file (ex.:mdp/top files)?
2- Are there any values I shouldn’t use as an energy threshold?
Thank you again
Decreasing the ligand restraint will indeed mean changing the force constant in that file. There are no particular values to avoid.
The parameter tc-grps = Protein_Lig means that the group named Protein_Lig, that you have defined previously, will be considered together by the thermostat (tc-grps: temperature control groups). This has no influence on things other that the thermostat, and does nothing for the restraints.
A few more details on how this works mechanically. In the .mdp file, you have a line like such:
define = -DPOSRES -DPOSRES_WATER
This instructs gromacs to remember that POSRES and POSRES_WATER are defined. (Why these start with -D is because it mimics the command line argument for the C preprocessor - though this is not important)
In the topology .top, you have line like these:
#ifdef POSRES
#include "posre.itp"
#endif
The paired #ifdef and #endif directive instruct gromacs to remove whatever is in between these if POSRES is not defined. The #include directive says to look for a file named posre.itp, and paste its content at this position.
This means that, if you have define = -DPOSRES, the following line will be translated to something like this, placed at the position these directives were:
; In this topology include file, you will find position restraint
; entries for all the heavy atoms in your original pdb file.
; This means that all the protons which were added by pdb2gmx are
; not restrained.
[ position_restraints ]
; atom type fx fy fz
1 1 1000 1000 1000
5 1 1000 1000 1000
6 1 1000 1000 1000
...
But if its not defined (define = ), then it is equivalent to no text at all.
1 Like
Thank you again for your answer and explanations. I understood much better the ideas now!
I just want to ask one more question…
Following the tutorial, I performed two equilibration steps, an NVT and an NPT, using the restrains with 1000 of force constraint. Could it be considered a “gradually decrease restrain” of the ligand by performing the same two steps, and adding a third step in which I use the same NPT MDP script, but change the force constraint to 500 in the ligand posre itp file for this third step?
Yes, I would say something like an NVT equilibration (with restraint of 1000), then NPT (restraint 1000) then any number of additional NPT equilibration step with decreasing restraining force would qualify as “gradually decrease restrain” .
1 Like
Thank you very so much for your help!