Box changing shape after equilibration

GROMACS version: 2024.4
GROMACS modification: No
Hello!

I ran into a situation with my simulations and I was wondering if anyone would have some wisdom as to what is causing this strange phenomenon.

I create a protein-membrane system using CHARMM-GUI Membrane Builder. I used gmx editconf to transform the pdb file into a gro file. The exact command I used was gmx editconf -f step5_input.pdb -o step5_input.gro -bt cubic -c -d 1.0. Using the mdp files provided by CHARMM, I ran the minimization and equilibration steps. But when I went to check out the equilibration gro files, I had an unwelcome surprise. Somehow, my box turned from a cube into a more spheric shape. I had never had that happen to me in any of the simulations that I’ve run. It seems that something happened between the minimization and equilibration steps, as the minimization.gro is still intact.

The mdp file I have used for the first equilibration step is this one:

define = -DPOSRES -DPOSRES_FC_BB=4000.0 -DPOSRES_FC_SC=2000.0 -DPOSRES_FC_LIPID=1000.0 -DDIHRES -DDIHRES_FC=1000.0
integrator = md
dt = 0.001
nsteps = 125000
nstxout-compressed = 5000
nstxout = 0
nstvout = 0
nstfout = 0
nstcalcenergy = 100
nstenergy = 1000
nstlog = 1000
;
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
coulombtype = PME
rcoulomb = 1.2
;
tcoupl = v-rescale
tc_grps = SOLU MEMB SOLV
tau_t = 1.0 1.0 1.0
ref_t = 303.15 303.15 303.15
;
constraints = h-bonds
constraint_algorithm = LINCS
;
nstcomm = 100
comm_mode = linear
comm_grps = SOLU_MEMB SOLV
;
gen-vel = yes
gen-temp = 303.15
gen-seed = -1

Hi! I’m not familiar with the CHARMM-GUI membrane builder workflow, but I’m pretty sure this is simply a case of your box being too big. Specifically, what you get from the membrane builder is a ready-made system, waters ions and all, but then you use editconf to put it in a simulation box that’s 1 nm bigger than that system on all sides, effectively creating a vacuum. Simply leaving out the editconf step (or at least the -d 1.0 option) should do the trick.
If you think about it, what you see here is not such a strange phenomenon, but actually pretty nicely in line with what you would expect in a vacuum, with the water forming a droplet around the membrane because of surface tension.

Hi!

Thank you for such a detailed answer. I had never considered that, by adding -d 1.0, I could be inadvertently creating a vacuum. I guess I was lucky because this is the first time I ran into this problem!

Just to follow up on that, I believe there’s an option to get everything in GROMACS format straight from CHARMM-GUI. I’ve run a few systems using their setup tools and it’s worked well out-of-the-box.

You probably just stopped after doing only NVT equilibration this time. If you also did NPT, the pressure coupling would remove the vacuum.

I somehow managed to completely miss that on my first run. For some reason, I thought CHARMM-GUI only outputted pdb files and not gro files. Once I switched to using their files, everything worked out pretty well!

Unfortunately, the NPT equilibration did not solve the problem with the vacuum. Even after undergoing both equilibration steps, the vacuum didn’t budge. What did work out was to drop the editconf approach altogether and just use the gro file that was provided by CHARMM-GUI. I had thought that CHARMM-GUI only provided pdb files which is why I was converting my pdb file into a gro file. But something in that step must be messing up with the system, because even without the -d flag, it was still creating a vacuum.

As you noticed, the box provided by CHARMM-GUI is already fully prepared and ready for minimization. You should not use editconf, especially if the system is already solvated, as this will just add empty slabs on the sides of the box. And the same goes for genion etc.

Relaxing the vacuum is not necessarily so simple, also because it’s very dangerous as the NVT with vacuum borders will deform your bilayer a lot (you can see it also from the image you posted, the lipid leaflets are starting to open up). So even if the vacuum relaxed, I would still be wary of the lipid leaflets status and the deformations induced, in particular at the box edges where the vacuum will effectively meddle with the PBC image of the lipids.

I definitely learned my lesson regarding using editconf with the CHARMM-GUI files after this. At the time, I hadn’t realized how many complications would arise from adding those empty slabs.

The bilayer definitely looks very different from what I was expecting. Luckily, I ran it all again but without using editconf and now the system is looking (and behaving!) as expected.