GROMACS version: 2021.4
GROMACS modification: Yes/No
Here post your question : I went through the biphasic tutorial , I could easily complete it as all the files needed are already given there. But I want to use a different forcefield like charmm36, but I get stuck in the initial steps of the tutorial. Though I could get cyclohexane files from charmm GUI like pdb, itp and topol.top , but I still get a lot of errors when I am trying to energy minimise the cyclohexane box (to get equilibrated box of cyclohexane). Can anyone please help me with the steps, to create the biphasic system?
What’s the specific errors? It needs more information and error notes to determine the reason why the problem appears in simulation.
As first of all I want to prepare an equilibrated box of cyclohexane, so during the NVT equilibration step, I am getting the following error :-
Fatal error:
Group Non-Protein referenced in the .mdp file was not found in the index file.
Group names must match either [moleculetype] names or custom index group
names, in which case you must supply an index file to the ‘-n’ option
of grompp.
or
Fatal error:
Invalid T coupling input: 1 groups, 2 ref-t values and 2 tau-t values
So can you please help me know what to edit in the nvt.mdp file in the tc-grp section which by default is Protein Non-protein or Protein Ligand , two coupling groups. But here I am equilibrating a box of cyclohexane. Kindly help me with this.
- Non-Protein is a default group for the system contained protein. If your system doesn’t contain protein, you can change the Non-Protein to what you want (like System).
- Checking the
tc-grpsin the .mdp file. The number of groups intc-grpsshould be same astau_tandref_t
I changed the tc-grp section as following-
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = systm system ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
But on running the grompp command for NVT equilibration, I get the following error-
Fatal error:
Group systm referenced in the .mdp file was not found in the index file.
Group names must match either [moleculetype] names or custom index group
names, in which case you must supply an index file to the ‘-n’ option
of grompp.
Here I don’t understand which one is the index file, and how to fix this. Kindly help me with this.
The topology file for cyclohexane tht I obtained from CHARMM GUI-
;;
;; Generated by CHARMM-GUI (http://www.charmm-gui.org) v1.7
;;
;; psf2itp_mol.py
;;
;; Correspondance:
;; jul316@lehigh.edu or wonpil@lehigh.edu
;;
;; The main GROMACS topology file
;;
; Include forcefield parameters
#include “charmm36.itp”
#include “CYHE.itp”
[ system ]
; Name
Title
[ molecules ]
; Compound #mols
CYHE 207
It is no need to define the tc-grps twice time with same group. You can change the Temperature coupling method to:
tc-grps = system
tau_t = 0.1
ref_t = 300
It is a spelling error of ‘systm’ in your file, and I’m not sure if the error originated from this, or if there is an other group defined incorrectly.
Thank you very much @Seyilaxa , that was the error.
For biphasic system, first we have to create the equilibrated box of cyclohexane. So for that, after equilibration (NVT and NPT) steps, we have to do the MD production run too?
Usually the equilibrium and the production use different pressure coupling method (Berendsen and Parrinello-Rahman pressure bath). P-R barostats is stricter on theory, but you can also use a long enough NPT to equilibrate your system if you think it will not caussed an obvious issue in your system. Sometimes you can even skip the NPT and use the production parameter directly which is generally not recommanded.
I think you will have a further equilibrium when you build the biphasic system, so the NVT and NPT is enough for you. Of course, all of the simulation step depends on your purpose and the property of your system. You can read the manual to understand the different between the different pressure bath, then make a choice by yourself.
https://manual.gromacs.org/current/reference-manual/algorithms/molecular-dynamics.html#pressure-coupling
It is important to note that although the Berendsen pressure control algorithm yields a simulation with the correct average pressure, it does not yield the exact NPT ensemble, and does not compute the correct fluctuations in pressure or volume. We strongly advise against using it for new simulations. The only useful role it has had recently is to ensure fast relaxation without oscillations, e.g. at the start of a simulation for from equilibrium, but this is now provided by the stochastic cell rescaling, which should be used instead. For full anisotropic simulations you need to use the Parrinello-Rahman barostat (for now). This does have the same oscillation problems as many other correct-ensemble barostats, so if you cannot get your initial system stable you might need to use Berendsen briefly - but the warnings/errors you get are a reminder it should not be used for production runs.
When I run the box of cyclohexane through final step of MD production run, I am getting the following error-
Fatal error:
Not enough memory. Failed to allocate 18446744071606518219 aligned elements of size 4 for grid->grid
Can you please explain me what is this error and how to resolve it?
Without more information, that definitely sounds like a bug. This is probably fixed in later version of GROMACS. If you cannot switch versions, for one reason or another, the cause of it might still be possible to work-around. It looks PME related. Are you running with PME on GPU? Trying running mdrun with -pme cpu and see if it helps at all.
It could also be an ordinary system explosion that is not reported more elegantly.
As MagnusL said, have you ever paid attention to the change of the box size in the NPT process? Sometimes the system’s pressure do not be equilibrated quickly even occurs oscillations (especially when you use the Parrinello-Rahman approach), then lead to a system explosion which may cause this error. Changing to another barostat or setting a larger ref_p may accelerate the equilibrium process.
The topology file is as follows-
;;
;; Generated by CHARMM-GUI (http://www.charmm-gui.org) v1.7
;;
;; psf2itp_mol.py
;;
;; Correspondance:
;; jul316@lehigh.edu or wonpil@lehigh.edu
;;
;; The main GROMACS topology file
;;
; Include forcefield parameters
#include “charmm36.itp”
#include “CYHE.itp”
[ system ]
; Name
Title
[ molecules ]
; Compound #mols
CYHE 209
The md.mdp file-
title = charmm36 CYHE MD simulation
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 500000 ; 2 * 500000 = 1000 ps (1 ns)
dt = 0.002 ; 2 fs
; Output control
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save coordinates every 10.0 ps
; Bond parameters
continuation = yes ; continuing from NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = no ; continuing from NPT equilibration
On running the md production run on the cyclohexane box, it stops with the following comments-
starting mdrun ‘Title’
500000 steps, 1000.0 ps.
step 261800, will finish Fri Apr 19 00:27:03 2024
Step 261842, time 523.684 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.012168, max 0.484164 (between atoms 646 and 647)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
307 308 90.0 0.1089 0.1522 0.1111
646 647 90.0 0.1089 0.1649 0.1111
Step 261843, time 523.686 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 535.393616, max 18965.492188 (between atoms 322 and 323)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
307 308 107.1 0.1522 2105.7976 0.1111
322 323 52.4 0.1111 2107.1772 0.1111
322 324 30.9 0.1111 0.1105 0.1111
646 647 90.0 0.1649 0.1371 0.1111
Wrote pdb files with previous and current coordinates
Segmentation fault (core dumped)
Can you please explain this to me and help me resolve this?
It’s almost certainly a system explosion or similar crash.
Then do I need to increase the box size or change any of the parameters in the input file, to prevent this from happening again?
There are a few things I would try:
- Run the preceding NVT equilibration longer (assuming that you did that before these NPT ). Increase the number of steps by a factor 2 to 5 (depending on how long it takes to run the simulation).
- In the NPT equilibration (and production), switch to
pcoupl = c-rescaleandtau_p = 5.0ortau_p = 10.0. - Check if the compressibility is correct for the medium. You are using water compressibility for cyclohexane. It might be close enough, but I don’t know.
- Check if there are problems with the force field parameters.
Thank you very much @MagnusL , the last step worked , I changed the pcoupl=c-rescale (tau_p=2, I did not change it to 5) and the production run got completed without any error.
The density that I got after NPT equilibration is as follows-
Just to make sure if I followed all the steps correctly to build an equilibrated box of cyclohexane-
- I used charm GUI (not CGenFF website) to generate cyclohexane structure, from which i got the chx.itp, chx.pdb and chx.top files.
- chx.pdb to chx.gro using gmx editconf
- chx.gro to chx_box.gro
- energy minimization. equilibration (NVT and NPT) and finally the production run
- For the input parameter files (em.mdp, nvt.mdp, npt.mdp, md.mdp), all the .mdp files I used from the protein-ligand tutorial (as I am using charmm36 ff) with minor editing.
- In the .mdp files, I removed the position restraint (-DPOSRE)
Did I do any mistake in any of the above steps? Kindly let me know.
The above steps sound good. But looking at the density plot from the NPT equilibration, I would recommend to extend the NPT equilibration until you clearly see that the observable does not change significantly for an extended period of time, but rather fluctuates around a stable value. The alternative would be to discard the beginning of your production simulations, i.e., discarding that part of the simulation as equilibration.
In the previous density plot, NPT equilibration has already been done for 500 ps. Is much longer equilibration to build an equilibrated box of cyclohexane okay or did I do any mistake in any of the steps? Because in the biphasic system tutorial, though gromos ff is used, but equilibration-both NVT and NPT has been done for only 100 ps. So I am really worried if I did anything wrong.
Kindly see the following density plot I got after NPT equilibration of 200ps for cyclohexane box, unlike the previous plot where NPT equilibration was carried out for 500ps
Also I wanted to know if I am comparing simulation result of 3 different systems with same protein (same forcefield- charmm36 is used for all the studies), do I need to change the pcoupl= C-rescale for all the systems (as I have changed it for cyclohexane box in it’s NPT and md run .mdp files) ?
In general, you should not be too focused on equilibration times from other sources, especially if your system size or force field are different. It is always better to check that your system is in equilibrium (or in a state that you want to study). Anyhow, in the tutorial it says: " 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."
So after 10.1 ns of NPT simulations the density was reasonable (presumably stable for a long period of those 10 ns). It is of course a matter of what you call your stages and what you observe, but I would consider those additional 10 ns equilibration as well.
I would recommend using the c-rescale barostat, yes.
Ok. I got it now. I think i misunderstood, in the tutorial I thought that the cyclohexane box went through 100 ps each of NVT and NPT equilibration then 10 ns of md production run.
So for building this equilibrated box of cyclohexane, if 10.1 ns NPT equilibration was done, then what should be the appropriate time step for production run or is it okay if we skip the production run step here?