Problems with npt in one periodic dimension

Dear all,

I am trying to simulate a system, with a zeolite pore and water next to it.

To keep the water molecules from interacting with the zeolite (which is cut from the volume structure and hence doesn’t have a proper surface to study the acutual interface) I build a box around the zeolite and made the walls repulsive for water molecules.

Now for my problem: To get the right density in the water reservior I want to apply pressure only along the z axis. After the system gets compressed for the first 1000 to 10000 steps the system starts to expand into the z direction.

Since I want the system to be under ambient conditions I set the barostat to 1 atm, but I also tried higher pressures and different damping factors. The system always expands after initially getting compressed and I can’t figure out why.

Here you can find the input and the datafile.
in.zeo_h2o_clean (1.7 KB) (73.6 KB)
h2otip4p.txt (437 Bytes)

Any help is greatly appreciated and thanks in advance


That is just normal behavior of how the barostat coupling works.

You have - at the core - a fictitious harmonic oscillator where the target pressure represents the equilibrium distance of the oscillator and the current pressure the current distance. The pdamp parameter is translated into a (fictitious) mass. Now at the beginning of the simulation, the fictitious mass is at rest but will experience a force because its current pressure is not the target. Over some steps, the mass will pick up speed and the box size adjust accordingly. Depending on how far away your initial configuration is from a state where the pressure changes to reverse the process, the mass can pick up significant speed. Because of that, the box will change too much, i.e. the fictitious oscillator with “overshoot”, and thus the size of the box will move away from the target. But it should return and the amplitude of the oscillations should become smaller until there are just small oscillations around the target value.

You have the same behavior also with the thermostat, but it is not so visible, because you typically have a different fictitious mass (different tdamp). In general, a Nose-Hoover thermostat/barostat is not so effective if your system is far away from the target/equilibrium and thus you need to give it some time, i.e. increase the number of MD steps for the run. The choice of the tdamp/pdamp parameters is crucial. Generally, using larger values reduce the “overshoot”, but they also have to be chosen in such a way that the characteristic frequency of the fictitious oscillator(s) does not couple to any of the intrinsic frequencies of your system, e.g. libration motions or stretch/bend oscillations of your molecules since then you would get a resonance and the system may not stabilize at all.

In summary,

  • it helps to start from a geometry that is closer to the equilibrium. you may be able to achieve this right from the original input, but also adding a minimization with fix box/relax may help.
  • you need to be more patient and run longer to wait until the system equilibrates
  • varying the tdamp/pdamp parameter somewhat may make a difference. for box relaxation, you usually would make it larger, which means to run longer.

Thank you very much for the detailled answer!

I got it to work by modifying the system a bit to be closer to equilibrium and running more MD steps. Also the fix box/relax looks like a handy tool, that I will use in the future.

Kind regards!

1 Like

Thanks for reporting back. This will be extremely useful for people that have a similar issue and search for it in the archives and find this discussion.