Does it help if you use a more correct compressibility? I think (but you should probably look it up yourself) for cyclohexane it would be closer to 1.1e-4 to 1.3e-4 per bar.
Ok Sir. I will try changing the compressibility.
Sir I have 3 systems - : protein-water, protein-cyclohexane and protein at biphasic(biphase of cyclohexane and water)interface. I have to compare these 3 systems based on different temperatures.
For protein-water system, default compressibility for water as in .mdp file is used. I can change the compressibility to cyclohexane for cyclohexane-protein system and while preparing equilibrated box of cyclohexane. But Sir, for the biphasic system of cyclohexane and water (for this only I am preparing the equilibrated box of cyclohexane), what compressibility should I use in .mdp files (default of water or compressibility of cyclohexane)?
I’m not sure exactly what compressibility that would make sense for biphasic systems. Have you looked in the literature? If you cannot find anything, I’d recommend a value closer to that of water than to cyclohexane. But you will have to decide the value of such parameters as you are the domain expert.
By the way, how many cyclohexane molecules do you have in your system? In a previous thread (Biphasic tutorial - #22 by SP06) you said you had a cubic 4.9 nm box with 218 cyclohexane molecules. In the tutorial there are 1114 molecules in a cubic 5 nm box. You might have so few molecules that it is taking unnecessarily long to reach equilibrium and the system is so small that it gets very sensitive to pressure and temperature changes. Make sure that your system is not too small (in both number of molecules and volume).
Ok Sir, I will insert more cyclohexane molecules into the box and try equilibrating it, and see if it works.
Also, in the tutorial though it’s said that gmx-insert added 1114 cyclohexane molecules, but the attached equilibrated box and topology mentions 466 molecules, this confused me a little. In the topology file given there, do we need to include water model there-
http://www.mdtutorials.com/gmx/biphasic/Files/chx.top
With this, can you please explain me the following from the biphasic tutorial-
" At this point, you should simulate the cyclohexane box such that its density stabilizes. Following steepest descent minimization, 100 ps each of NVT (298 K) and NPT equilibration (298 K, 1 bar), I was able to obtain a reasonable density with an additional 10 ns of simulation under these same NPT conditions"
In this it is said that cyclohexane box is equilibrated till density stabilizes, there is no mention of stabilizing the pressure, so is it okay only for the density to be stable and not pressure?
You don’t have to include the water model when there is no water in the system.
Remember, that information is from a tutorial. It is just an example to get people to learn what they are doing. As soon as you are doing any research you will have to decide, and defend, your criteria for when the system is in equilibrium. I would recommend checking at least temperature, pressure and potential energy. If density is important for your system, and known beforehand, it’s a good idea to monitor that as well. The density is usually more stable than the pressure, which fluctuates a lot, especially in a small system.
Roe and Brooks developed a ten-step protocol for system preparation/equilibration. Check it out, especially their density fit in step 10:
Roe, D. R.; Brooks, B. R. A Protocol for Preparing Explicitly Solvated Systems for Stable Molecular Dynamics Simulations. J. Chem. Phys. 2020, 153 (5), 054123. A protocol for preparing explicitly solvated systems for stable molecular dynamics simulations | The Journal of Chemical Physics | AIP Publishing.
I understand. Actually, I was little confused from lysozyme tutorial where average pressure of system is 7.5 bar (though reference is 1 bar) with RMSD value +/-160.5, and there it is said it’s okay for such large fluctuation in pressure. (So while preparing my system I was ignoring the pressure fluctuations at first and focused on stabilizing the density, and once I got the desired density I was proceeding to production run). So is it really okay for any system to have such fluctuations in pressure? Can you please explain and clarify my this confusion?
It is a bit unfortunate that the tutorial uses the RMSD as specification of the uncertainty of the average value. It is better to use the Err.Est. value for that (from gmx help energy):
An error estimate of the average is given based on a block averages
over 5 blocks using the full-precision averages. The error estimate
can be performed over multiple block lengths with the options
-nbmin and -nbmax.
Very high fluctuations (very high RMSDs) are expected for small systems (see the manual: Terminology - GROMACS 2024.2 documentation). Even for large systems the RMSD can be high. Pressure RMSDs between 250 and 500 bar are not uncommon at all. The pressure may or may not end up within 2 SEMs (standard errors of the mean) of the pressure set by the barostat. But it should, in most cases, be fairly close. What you want to make sure, when evaluating whether your equilibration stage is finished, is that the pressure is stable. Check the Tot-Drift from gmx energy to make sure it is not drifting very much (you will have to judge what is OK) and/or check the running average of the pressure to make sure it is flat (with large enough average windows). I would pay more attention to other observables than pressure, in most cases.
For the lysozyme tutorial, I think it’s a good example system. But keep in mind, the equilibration stage is extremely short in that example (like in most tutorials, so that they can be run in reasonable time). You would not find many NPT equilibrations, of proteins in water (even if they are simple proteins), as short as 100 ps in modern literature.
Thank you very much Sir for explaining in detail.
It’s worrying me that the pressure is still not stabilizing after 50 ns. Sir would you recommend changing the barostat ?
The pressure coupling I am using for NPT-
pcoupl = C-rescale ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 5.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in
Those barostat settings should be fine (if the compressibility is reasonably correct). When you say that the pressure is not stabilizing after 50 ns, do you mean that its average is drifting significantly during the last 10 or 25 ns? That sounds strange to me.
Yes, till 50 ns the average was drifting. But, after another 10 ns of equilibration, i.e. from 50-60 ns -
Avg. pressure = 1.20019
Err.Est. = 0.76
RMSD = 106.226
Tot-Drift = 0.123671
(bar)
Sir, I had another confusion that I wanted to clarify, does the use of barostat type affected by or has compatibility issues with GROMACS version?
Good that you reached an equilibrated state. It’s strange that it took so long.
I can’t think of any issues with the c-rescale barostat.
Sir
Though you explained me about the drift and err.est, can you please explain this in form of values like what range should be for Err.Est and Tot-Drift, for pressure, that will be much easier for me to analyse after gmx energy if my system is well-equilibrated.
It is difficult to give any specific numbers. You could focus on the density and use the criteria outlined in the article recommended above, in Longer NPT equilibration - #27 by slovas.
That said, if pressure would be an important criterion for a simulation of mine I would roughly consider it equilibrated when the average pressure (over the last 5-25 ns of the equilibration, depending on system) is no more than 2 SEM (err.est.) from the target and the tot-drift is within 1 SEM of 0. If that is difficult to fulfill (e.g. if the SEM is very low) I would visually inspect the pressure profile (plotting running averages with wide windows) to make sure that it looks stable. Then I would probably (if I think it can be motivated if challenged) accept a tot-drift of ~0.5 bar and an average pressure within 2 bar of the target, provided that the potential energy has clearly reached a plateau.
By plateau of potential energy, do you mean by the one calculated after energy minimisation step initially, or, do I need to check it after NPT equilibration too?
I think the potential energy is good to check during NPT equilibration. I think that is at least as important as temperature and density, as I’ve stated repeatedly above.
Got it. Thank you very much Sir, for patiently explaining some of the very important concepts.
Hi Sir
For one of my systems of protein-water (which I studied long before while just starting out to learn and followed the tutorial completely, and did not focus on pressure stability after NPT), I just checked density after NPT equilibration (density reached the target), did not check pressure at that point and proceeded to production run(3 runs, 100 ns each). Now, when I went back and checked the pressure of that system after NPT, I found it was -6 bar with larger values of err.est. and Tot-Drift. But when I checked the pressure during the 3 production runs, I found it was stable during each run (around 1 bar, with err.est. varying from 0.011- 0.003 and Tot-Drift of -0.034 to +0.02 ). Can you please help me with this. Do I need to go back to extend NPT again?
If it is stable during production it sounds OK.
For some systems, for which I have run an NPT equilibration way longer than required (e.g. I set NPT for 50 ns, but from plot I found after 3 ns, reached plateau, for the entire plot from 3-50 ns.) Should I or can I shorten my range (maybe for 3-15 or 20 ns) for NPT using gmx energy ( using gmx energy -b 3000, it’s calculating for the entire range of 3-50 ns)