I am a user of LAMMPS. Recently I do research about coarse-grained simulations for lipid bilayer membrane. I adopt Cooke model (Cooke et al. PRE. 2005). Now I want to generate a large planar membrane in LAMMPS, however I really have met some problems.
Some researchers used barostat to control the lateral tension and generate planar membranes (Reynwar et al. Nature. 2007). I don’t know how to controll the lateral tension in LAMMPS, and I just use barostat to control the pressure. I adopt the NVE ensemble with Berendsen barostat(Fu et al. Computer Physics Communications. 2017. Section 5). So I use (units lj):
fix 1 all nve
fix 2 all press/berendsen x 0.0 0.0 1.0 y 0.0 0.0 1.0 couple none
The membrane is in x-y plane at time 0 tau. However, during simulations the planar membrane always rolling over and over. Other researches can get stable planar membrane(like Reynwar and Fu). And my questions are:
Does the berendsen barostat control the lateral tension of membrane?
How can I have a stable planar membrane. Did I misunderstand the zero-lateral-tension state? Or should I apply some constraints?
Thank you and best wishes,
1. Does the berendsen barostat control the lateral tension of membrane?
The Berendsen barostat (and thermostat) does not generate the
correct distribution of volumes (or energies) for a system under NPT
conditions. If possible, avoid using this old barostat. Your
membrane fluctuations will be incorrect.
2. How can I have a stable planar membrane. Did I misunderstand the
zero-lateral-tension state? Or should I apply some constraints?
I am currently using:
fix fxnph all nph x 0 0 1000 y 0 0 1000 couple xy
This applies a 0-tension barostat in the x and y directions with no
shearing. (Note: "couple xy" preserves the aspect ratio of the
simulation rectangle, and "1000" is the Pdamp parameter, which you can
For details, see:
For temperature control, I use:
fix fxlan all langevin 1.1 1.1 10.0 123456
# Note: In this example, the temperature is 1.1*ε
# (ε = "epsilon" used by the coarse-grained lipid)
# Note: The langevin damping parameter is now "10.0" in units of time
# (which are (m/(ε*σ^2))^(1/2) because we are using "units lj"
# This is underdamped compared to the real system, but it
# leads to more efficient sampling. Feel free to play with this parameter.)