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.
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
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.
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:
nsttcouple x dt (ps)
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.
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.