I created a model containing a bilayer membrane resembling the Inner Mitochondrial Membrane (IMM) available from CHARMM-GUI. Based on how Charmm-GUI created their optimized protocol, would changing the neighbour update list (nst=1 instead of nst=10) affect the validity of using their standardized protocol?
If I use the default nstlist = 10, I get the following message:
Steepest Descents converged to machine precision in 555 steps,
but did not reach the requested Fmax < 1000.
Potential Energy = -8.0869015e+06
Maximum force = 6.2726987e+03 on atom 116466
Norm of force = 5.5581973e+01
Note: Atom 116466 is associated with one of the lipids within the membrane.
However if I set nstlist = 1, I get the following message:
Steepest Descents converged to Fmax < 1000 in 4354 steps
Potential Energy = -8.7844850e+06
Maximum force = 9.9567993e+02 on atom 43592
Norm of force = 1.3002018e+01
You shouldn’t be too worried about neighbour list updates (nstlist) during your minimization steps. The purpose of your minimization is to just make sure all the atoms are playing nice before you equilibrate them. Both messages are basically telling you this is the case.
Are all the other mdp settings identical except nstlist? From what I’ve been able to see online, steepest descent convergence without Fmax convergence is usually a by-product of your timestep being too large. Your timestep for minimization should usually be 0.001 ps for minimization runs.
Leach’s ‘Molecular Modelling Principles and Applications’ (pg 284) mentions that non-bonded neighbour lists in water shouldn’t change dramatically over 10-20 steps, (especially if your steps are just 1 fs), but that’s an older text. This paper discussed how nstlist can affect membrane topology, but that’s beyond the consideration of a minimization run.
My personal opinion is as long as your minimization has minimised it enough so that you don’t get bad contacts during your equilibration steps, it doesn’t matter too much what your energy converges to and what mdp tweaks you’ve made. Save that thinking for equilibration and production runs.
“What you are provided from CHARMM-GUI are the recommended protocols from the authors of whatever module you’re using. You’re free to depart from those inputs in whatever way(s) you see fit, but if things don’t work, it’s probably because the protocols are generally pretty robust and used for many years. See individual publications for whatever modules you’re using (Solution Builder, Membrane Builder, etc), which detail the protocol. The 2016 CHARMM-GUI paper talks about how we’ve validated the implementation across software for some general systems.”
There shouldn’t be any real need for you to change any of those mdp settings for minimization, even if it hasn’t hit Fmax by the end of the minimization step. The mdp settings are more important for equilibration and production.
If you run an equilibration and it crashes because of bad contacts or too high energy, then you would go back and force gromacs to run until a specific number of steps or until an Fmax is reached.
But as a general rule, take the time to learn what each of the inputs in your mdp file are actually doing and whether or not they are the best way to run something. Especially for bilayers. The gromacs user manual even says the default parameters for a gromacs mdp shouldn’t be used for semi-isotropic bilayer systems and it recommends some adjustments.
Also, settings that were defined 10 years ago may not be the best settings today, so it’s worth looking for papers where people have built similar systems and finding out what they used. It’s especially important for doing this when running charmm-gui built systems on gromacs. For example, the charmm small molecule library which is available online and used in gromacs actually has some incorrect molecules in it. It confuses adenine with adenosine and some of its L and D sugars are incorrect.
This influences how often you calculate the neighbour list. Also, keep in mind that usually the values you see there is minimum value, that is, GROMACS will change the value during the run depending on the hardware and the system. Take a look at the manual.
That being said, as you can see from the manual nstlist = 1 is a special case because it is forcing GROMACS to calculate the neighbour list at each step. I would say that having higher frequency for neighbour list cannot make your simulations worse, if anything they make it better, as you have less chances to get artifacts due to missing interactions. However, the list calculation is extremely expensive and will kill your simulation performance, which is why usually they are > 1 and are further fine tuned by GROMACS (usually you get values between 20 and 80 in my experience).
It is very interesting that you get such a borderline case where this value changes the outcome of the EM step. Honestly, I would just stick to the minimised structure with nstlist = 1 and go on with the rest of the simulations to relax the system with whatever protocol you have. The EM step is really to get rid of near-clashes and decrease the potential energy (and so the forces) which might explode the system if you started directly to compute forces and update positions.