GROMACS version: 2019.1
GROMACS modification: No
Here post your question
I did MD run on a peptide with no bond constraints and timestep of 1 fs. I got the following note during grompp:
“The bond in molecule-type Protein_chain_A between atoms 173 O1 and 174 HO1 has an estimated oscillational period of 9.0e-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”
I was wondering how does this affect my simulation?
It certainly does. All the force fields in GROMACS are designed to have at least h-bonds constrained.
Read about why that matters at http://manual.gromacs.org/current/reference-manual/special/remove-fast-dgf.html
Thanks for the link. But when I applied h-bonds constraint, my peptide still remained extended instead expected helical structure (which I see with no constraints). Any insight?
Is it reproducible in multiple simulations? What force field are you using? Again, none of the parameter sets in GROMACS are designed to be used in the absence of constraints, so anything you see with constraints = none is probably unreliable.
I conduct these simulations in vacuum using OPLS forcefield. Yes that extended structure is reproducible in multiple simulations. And I do get expected helical structure without applying any constraints. May be it is the absence of solvent?
OPLS is designed for use in condensed-phase simulations. The charges are over-polarized for any type of gas-phase simulation. I don’t think you can draw any kind of conclusion from these types of simulations.
If you want to get to the root of the constraint issue, you would have to do a vibrational frequency analysis for the amide group to see if the N-H bond oscillation is incorrect and therefore leading to inaccurate forces. Often, little attention is paid to the quality of parameters for these kinds of bonds because they are expected to be constrained in typical simulations.
Although charges are overpolarized, it was recently tested and validated that scaling down partial charges did not significantly change the gas-phase structure of proteins (https://pubs.acs.org/doi/10.1021/acs.jpcb.9b04014)
I am still trying to understand why I get more reliable structures with no bond constraints applied.
Why does the option for Constraints = None exist?
Tried two simulations (each 2x), one without any bond constraints and integration timestep of 1 fs and the second one with constraints on h-bonds and 2 fs timestep. Shown below is the rmsd plot, I do not understand why there was more fluctuation when h-bond constraints were applied. Any thoughts appreciated.
Any insights on the above observation? Why would there be more fluctuation when constraints are applied on h-bonds?
Another observation: My peptide has equilibrated well at 300K when no constraints are applied with rmsd of only 25 K versus 45 K with h-bond constraints.
Any insight would be appreciated!
I think Dr. Lemkul has already mentioned this in a recent reply somewhere else. Make molecule whole and recenter with trjconv to eliminate rmsd shots.
Thanks Dr Masrul. I manage to fix it.
Best Regards
Veeru
Got the following note when used all-bonds constrained with opls forcefield:
“You are using constraints on all bonds, whereas the forcefield has been parametrized only with constraints involving hydrogen atoms. We suggest using constraints = h-bonds instead, this will also improve performance.”
Does it make sense to have only h-bonds constrained? Why not then all-bonds constrained?
If the force field was designed to only have bonds to H constrained, that’s how it should be used. If it was designed for all bonds to be constrained, that’s how it should be used. This is important because if you change this setting, you’re introducing preventable errors in the simulation.
Imagine the rotation of a dihedral in something simple, like butane. The 1-4 interactions and dihedral parameters are balanced to generate the potential energy surface as the molecule rotates around the C-C-C-C bond. If the force field assumes these bonds are not rigid, this is implicitly part of the dihedral parametrization. If you then make all of the C-C bonds rigid, you negate their ability to stretch, thereby shifting the potential energy surface in unexpected ways. While the errors may be small, they may not have predictable effects and are entirely avoidable through proper use of the prescribed type of constraints.
Thank you @jalemkul for the thorough explanation, I understand now clearly!
