GROMACS version: 2023.3 with CUDA 12.3
GROMACS modification: No
I am currently simulating all-atom self-assembly of a lipid bilayer in a mixed water+ethanol solvent.
I was wondering if it recommended/more accurate to perform NVT and NPT equilibration prior to running the simulation?
The question arises because I wasn’t sure if they’re necessary since I am interested in the assembly process anyways during the production run. I also noticed this CG tutorial does not perform equilibration prior to running the self-assembly simulation. Moreover, it seems like I am having more issues with bilayer formation (e.g., extended boxes in z axis, holes in membrane due to permeating water) when I do NVT+NPT before the production run, so I’d prefer to avoid them if they’re not necessary in this situation. I currently randomly insert each element into a 9nm3 box and perform energy minimization, then run this simulation (see .mdp text below for reference):
Get the system into a stable state, in the ensemble you are interested in, so that simulations can run, e.g., at the temperature and with the time step you want to use.
Reach an equilibrium state, so that analyses can be made without (too much?) influence from the exact starting conditions.
In principle (or in theory), it would be possible to skip energy minimization as well and just start a production simulation. If the system is stable, i.e., the coordinates are compatible with the force field, it would work. Then you could just discard the first part of the simulation when doing analyses (e.g., using the -b option in gmx energy) - and then argue why you did it this way. There’s nothing fundamentally wrong about that. But I’d say there are good reasons why it is not common.
One could argue that the NVT equilibration stage is mostly to fulfil 1. And if you start with a simulation box that is not suitable for your system it can do more harm than good. You could do NPT equilibration directly.
If you get simulation artefacts when you equilibrate the system I would try to understand why this is happening.
I would say, if you are just interested in watching the membrane assemble, you can do it without equilibrating (if the system is stable). If you want to draw any conclusions from it I would definitely recommend equilibrating the system beforehand - otherwise you have no idea if your observations are relevant at all.
To address some problems that I did witness @MagnusL , for example, a subset of my lipid bilayer self-assemblies that I perform energy minimization + NVT + NPT + production run (with .mdp shown above) generated strange trajectories like this. (Note that I ran trjconv-whole as the last step):
As can be seen, the box that was 9x9x9nm3 continually extends in the z-direction during the 100ns long simulation. Could this be due to semi-isotropic pressure coupling? How can I prevent this from happening (either different run techniques or perhaps a different trjconv approach)? Ideally, my simulations should look like this one below everytime:
Looking at your output, I think the main problem is the orientation of your system (but it’s a bit difficult to see in the first image). When you are running with semiisotropic pressure coupling, the x and y scaling with be linked together and the z dimension scaling will be separate. You want the z dimension to be along the membrane normal. Otherwise the scaling of the system will cause artifacts.
Secondly, I would check that the simulation box is enclosing the simulation system tightly before the NVT equilibration starts, and also check that you don’t have large vaccuum volumes after NVT equilibration.
@MagnusL Thanks for your insight. Because I am interested in bilayer self-assembly, I am randomly placing the elements into a rectangular box. Thus, I understand the value of semiisotropic coupling for pre-formed membranes, but, in my case, it seems like there is a 2/3 chance that the membrane assembles with its normal not in the z-direction (i.e., x or y). Is there anyway to improve the odds of semiistropic coupling working with the normal in the z-direction? Is there a recommended precedent for bilayer self-assembly?
I noticed this paper did two rounds of NPT with isotropic coupling followed by a production run with anisotropic coupling? Would this iso+iso+aniso workflow perhaps be better for bilayer self-assembly?
Any insight you’d have on this would be super helpful. Thanks again.
Yes, I would definitely suggest that you start with isotropic pressure coupling until the membrane is formed and you know its orientation. I think semiisotropic pressure coupling should only be used (or be used with great care) when the membrane is in the x/y plane (its normal in the z direction). You could of course use aniisotropic pressure coupling, but that comes with other potential problems (similar to the first image in the post above), especially if the membrane is not yet stable or if it is very fluid.
When the membrane has been formed (after isotropic NPT equilibration) you could rotate the system to make it suitable for semiisotropic pressure coupling.
