Thanks for the picture. (Pictures are always helpful.)

In this case, I would not barostat the Z direction. Use something similar to my first suggestion instead:

fix 2 all nph x 0 0 1000 y 0 0 1000 couple xy
If my understanding is correct, then the width of simulation box in the Z direction should remain constant.

However people (and reviewers) may raise questions about why you did not fill the entire simulation box with water. (If you only did it to increase the speed of the simulation, this may not placate your critics. Thin water layers can behave differently than bulk water.)


The picture confirms what I suspected and I also have to fully support what Andrew said.

In summary:

  • you must not have barostatting in z-direction. It is meaningless, since you are effectively simulating a non-periodic system, so the pressure would be always zero due to the infinite volume. You only get a non-zero pressure because you arbitrarily set a boundary. But since you could do the same simulation with a large box (and thus lower pressure) for the same physics, it should be obvious that pressure in z direction is meaningless and so would be barostatting. on top of that, when using pppm with slab correction, you have to make certain, that your box does not shrink in z direction and system itself not expand or buckle significantly, so that the requirements for the slab corrections are maintained, otherwise you will need to increase the vacuum even more (the usually used slab expansion factor of 3.0 is for well behaved systems).

  • I have to strongly emphasize that the system size should be determined by the physics not computational efficiency. what you have is a membrane with a thin water coating. if that is what you want to simulate, fine, but if you want to simulate your membrane in contact with bulk water, you need much more water.

  • furthermore, you have to deal with serious finite size effects, and that is probably going to hurt you much more than the reduction of water. if you are pulling a particle this big through a membrane this small, with a box this small, you can forget getting meaningful results. due to periodic boundary conditions you are moving a large (and rather dense) array of particles through your membrane, that is going to be very different from just a single nanoparticle. for that you have to have a box so large, that the pressure/density fluctuations caused by moving the nanoparticle through the membrane will not affect the membrane particles across the periodic boundary. on top of that, you will have to move it so (very) slow, that the system has sufficient time to relax (and it still is likely going to be much faster than any equivalent experiment).

so unless you can round up orders of magnitude more computing power, and do a system with meaningful choices of (relative) dimensions, your results will be bogus and useless, or at least not representative of the kind of system that you are thinking of.
of course, it also is just common sense to not start with simulating such a large system given your limited of experience. such large amounts of CPU time require hard work to come by and thus need to be used wisely. The only way to confirm that, is to start with a much smaller nanoparticle, where you can employ more manageable system sizes (you still need to increase the amount of water significantly, though).

…and there is one more problem to contend with: when using fix rigid for your nanoparticle, you are - as far as the other atoms in the system are concerned - moving an object with effectively 0K around. that will also taint the results. the more the large the nanoparticle is. collaborators of mind have been simulating (lots) of (different size) buckyballs in coarse grained membrane systems and that has been a major concern for their results (and that was 10 years ago, so the barrier for what is meaningful and how much compute power may be used is much higher).


Thanks for your complete response , can please just restate the last part about rigid ?

Thanks for your complete response , can please just restate the last part about rigid ?

I am not sure, but I am guessing that Axel is referring to problems
with the thermostat caused by the reducing the degrees of freedom as a
result of making large objects rigid (instead of letting their
internal atoms move). If you are not using fix rigid, then you don't
have to worry about this issue. If you are, I think there are ways
around it, but you will have to read the fix rigid and compute
temperature documentation carefully.

As for how to make your simulations faster, I suggest taking a look at
the MARTINI force field. (Or use DRYMARTINI which is much faster and
does not require water). Trying to simulate an all atom system with
explicit water and long range electrostatics will be very, very slow
and costly. I'm not sure it is a doable project.

My own simulations of much simpler faster coarse grained models (such
as the Cook and Deserno 2005 model, or the Brannigan2005 model), it
can take a day on a fast computer with 24 cores to simulate a hole
opening in a lipid membrane. And this kind of coarse grained model is
about 1000 times faster than an all-atom water similar to what you are

Here are some links.

1) The MARTINI force field
(Efforts to make MARTINI available in LAMMPS are in progress but are
not working well yet for all lipid types.)

2) "Tunable generic model for fluid bilayer membranes"
Cooke IR, Kremer K, Deserno M, Physical Review E, 2005

3) G. Brannigan, P.F. Philips, and F.L.H. Brown,
Physical Review E, Vol 72, 011915 (2005)

With fix rigid the object is effectively frozen and thermostatting it won’t change that as that thermostat only applies to the 6 remaining degrees of freedom.