NPT vs. pressured applied with pull-code

GROMACS version:
GROMACS modification: Yes/No
Here post your question

What is the difference between applying a pressure using NPT equilibration versus applying a fixed force using the pull code?

I have a lipid membrane confined between two metal slabs. First I equilibrate the system under NPT at 1 atm, then I create some vacuum at top and bottom, and apply a fixed force on the metal slabs to squeeze the membrane. When that force is equivalent to 1 atm, the morphology of the membrane looks totally different from the NPT equilibrium state. Why is this happening? Thank you!

There are many things that could go wrong in such a comparison, especially unit and area conversions.

But one concrete question is what kind of pressure coupling you applied. The slabs change the volume along one dimension, a barostat in three dimensions by default, so the membrane area could be different.

The system has always been slabs+membrane.

In my NPT equilibration step, the box is a tight fit to the system, and I applied a Berendson barostat semisotropic, at 1 atm in both z and x-y plane. Then in the pull-code, where I applied fixed force onto the slabs, I expanded the box in z-direction and applied the same pressure coupling, but adjusted isothermal-compressibility in z to zero.

The force has a unit of kJ mol-1 nm-1, which is converted to force by multiplying with Avogadro’s number, and division by slab area gives the pressure.

But did you use the average area of an equilibrated NPT simulation for the simulation with a force?

Yes, the pressure is evaluated using the force in the pull-code and the area of the system, which did not change for the two methods. What’s changed is the membrane structure.

Now I think it’s because NPT gives you an equilibrium solution while replacing pressure coupling with a force gives you dynamic results. They are not equivalent as the algorithm is different. NPT equilibrium is achieved by scaling spacing and velocities of atoms, while by exerting a force the positions of atoms are changed through integration of the Langevin equation. The two methods are supposed to agree if the macroscopic property compressibility is set correctly for NPT. Is this correct?

Sorry I missed your question on the pressure coupling. It’s Berendson since I have a periodic system consists of more than one phase. Semiisotropic with tau_t=2.5 2.5 ; ref_p = 1.0 1.0; compressibility = 4.5e-5 4.5e-5; bar-1

After having done some readings on membrane compressibility, I learnt that a bulk modulus evaluated from membrane fluctuation in MD simulations is at the order of 6e-4 bar-1, which is higher than the model value. 4.5e-5, but lower than measured in our experiments 5e-3 bar-1.

I am trying to change the compressibility to the literature value 6e-4 bar-1 and run NPT, but I found it rather difficult as it prompt errors. Now I am trying to incrementally increase compressibility through ~10 simulations.

Another question is that if I change the pressure to tau_p= 1.0 2.0, will this be equivalent to adding pressure in z and be compare with experimental results?

Thank you!

Another [quote=“hess, post:2, topic:1822, full:true”]
There are many things that could go wrong in such a comparison, especially unit and area conversions.

But one concrete question is what kind of pressure coupling you applied. The slabs change the volume along one dimension, a barostat in three dimensions by default, so the membrane area could be different.
[/quote]

I don’t understand your question and what you mean with it. Do you mean changing the value of tau_p from 1 ps to 2 ps? That only changing the coupling time, it should not change the average pressure only it’s fluctuations.

Note that Gromacs 2021 now has the Bussi c-rescale barostat, which give the correct ensemble, unlike Berendsen (although the errors are usually small).