GROMACS version: 2019.1
GROMACS modification: No
Here post your question
Simulating one peptide in vacuum in a 999 nm^3 box, rmsd shots even after recentering with trjconv. Expected structure on the top, but keeps unfolding back and forth. Any insights?
Sounds like expected behavior for a short peptide.
In some trials with same parameters, I see no unfolding of the peptide, stable rmsd within 50 ps of production run. This is what I do not understand.
I don’t think you can reasonably conclude anything from 50 ps. Even in this simulation, if you stopped it at 50 ps, you would conclude that the peptide unfolds immediately and stays that way, but clearly it re-folds quickly, and remains that way for the majority of the simulation.
@Justin, sorry I misworded it, in some trials, I see no unfolding of the peptide, stable rmsd is achieved after 50 ps till the end of the production run, which usually I run till 300 ns.
This sounds like normal inter-trajectory variability. How many replicates have you performed? Do you start from different velocities, different initial geometries, or both?
I have done atleast four runs, yes indeed so much variability between these trials, I usually use random seed generator (Set to -1) and used minimized structure from one simulation cycle (nvt, md, em) for next trial
In addition, I do not see this unfolding when I set constraints = None.
Any insight would be really appreciated.
Have you done the QM vs. MM vibrational frequency comparison that I suggested for N-methylacetamide and a short peptide (e.g. trialanine) to confirm that the force field actually models H-X vibrations correctly? If not, start there. Again, as I’ve said many times, the force field you’re using is designed to be used with constraints and in the condensed phase. If you take it out of those conditions, it may not provide a reliable outcome. Models have limitations.
Justin,
I am using h-bond constraints now with opls forcefield in vacuum. I understand thoroughly why h-bond constraints to be applied. I have another rmsd plot here from a peptide, 3 replicates, all starting with same energy minimized structure, but different velocities (set by gen_seed = -1 during nvt run)
md.mdp (2.0 KB)
Q1: In trial 1, temperature fluctuation during production run was high (rmsd 130K). Why would that happen? That temp shot coincided with rmsd shot.
Q2: Got expected structures in trial 2 (and part of trial 3). Do I perform say x replicates till I see atleast 3 runs giving expected structures?
Q3: I tried two sets of temperature coupling parameters with Nose-Hoover thermostat, dt = 2 fs:
| tau_t (ps) | nsttcouple x dt (ps) | Temp rmsd |
|---|---|---|
| 0.5 | 0.01 | 30 K |
| 0.1 | 0.002 | 19 K |
I got similar structures from these two conditions. But why is the temperature fluctuation much lower with the parameter set 2? Is there a rule on how to chose these values?
The conformational change probably corresponds to a large change in COM velocity, which then causes the thermostat to go haywire.
As long as you don’t discount the replicate that didn’t turn out as expected, do as many replicates as you think necessary (I’d probably do 5-10 here) to have a “true” representation of the ensemble.
Regarding the thermostat, I wouldn’t use N-H with a gas-phase system. It is prone to big fluctuations, which may cause unnatural outcomes. Langevin dynamics (integrator = sd) is the way to go here. Your observations about temperature variation are a direct result of how strongly you’ve got the system coupled to the thermostat. The tighter you couple the system, the lower the variation in temperature. Again, avoid this entirely and run Langevin dynamics.
Justin,
I retried with Langevin dynamics, 5 replicates, no conformational change nor temperature variation (rmsd 20 K). I also get expected structures in all, thank you for correcting me!
Q1: There is no proper guide to implementing MD simulations in gas phase. For example, I chose N-H thermostat based on a few papers. Any good resource to use?
Q2: In this mdp file nvt.mdp (1.9 KB), is there anything alarmingly wrong to use for a gas-phase system?
GROMACS currently doesn’t support proper vacuum calculations. The right way to do it is a nonperiodic simulation with angular COM motion removal. Using a huge box with NVT and PBC is just a hack. I suppose it works well enough but it’s not really how these calculations are done in other software.
I always get so many LINCS warnings with angular comm mode, hence not using it. In simple terms, can you explain what does flying ice cube effect mean if rotational COM velocity is not removed?
Basically, whatever you simulate becomes a spinning rigid body.
