Box rapidly expanding for rigid body simulation in NPT

[mjksho]

Hello All,
I have a system of polarizable water (SWM4-NDP) sandwiched between graphite wall to simulate a channel. The system is non-periodic in z-direction where you have the walls, but periodic in x & y. The water molecules are treated as rigid bodies. The walls are not. I have used fix spring/self to prevent the wall atoms from drifting away. The set up looks like the image below.

When I try to equilibrate using NPT at 300 K and 1 atm, the box rapidly increases, the temperature continues to shoot up and system blows up eventually. These are what I have tried:

1. Equilibrate in NVT. This actually works. System properties (PE, KE, Temp) are stable. But I need NPT for the calculations I want to perform.
2. Reduce timestep. I even went as low as 0.1 fs. Didn’t work.
3. Couple the barostat to the rigid bodies (water molecules) as opposed to the walls, since there are more water molecules than wall atoms. Didn’t work.
4. Use more aggressive damping parameter for the barostat. This slowed down the expansion. But the simulation still crashed.

I am actually out of reasonable options to try and have attached my input file here for a quick review. Please let me know what you think. Thanks.
in.equil1 (2.0 KB)

[mkanski]

Several ideas come to my mind:

1. Visualize the kinetic energy of the atoms to see which ones gain the additional energy.
2. Simplify your system - see if the simulation is stable if you have only the walls or only the water (with fix wall/reflect for z direction).
3. Change neigh_modify delay 5 to delay 1.
4. Simplify your simulation - comment out fixes to see which one introduces the observed behaviour.

[srtee]

@mjksho

But I need NPT for the calculations I want to perform.

Why?

If your system is big enough, by thermodynamic equivalence, either will give similar results (implemented properly), so you can use the NVT integrator.

If your system is too small, neither ensemble is strictly or even close to correct and you might as well just use the one which doesn’t explode.

[srtee]

See also:

Issue while hydrating the polymer membrane

Has your system equilibrated? I find that whenever I’ve tried applying force to a carbon sheet to mimic pressure, the carbon sheets’ separation fluctuates wildly – I attribute that to carbon being “too light”. A poorly compressible system will throw out big pressure waves en route to equilibration, especially if started from an improperly small or large initial volume, and a piston with too little density like carbon will fluctuate in response, often too much for a typical MD box to handle. (Or…

[mjksho]

@srtee, thanks for the suggestions, one reason I want to try NPT is that I found that in NVT, because the pressure is not controlled, I see high pressure with Pxx, Pyy, and Pzz around over 30,000 atm. I also saw in old mailing list that Pzz should go to zero because of the fixed boundary condition, but in this situation, it isn’t.

The system isn’t small. I have about 12,288 water atoms, 3600 wall atoms in a 61.32 X 61.32 X 81.32 box.
Thanks.

[srtee]

I have never tried to understand nanoscale pressure in non-periodic directions because, as the Good Book says, I do not occupy myself with things too high and lofty for me.

But Pxx and Pyy I can certainly deal with. If the average values for those are exceeding 30,000 atm (and I assume you want to run a simulation at 1 atm), then you have started from a very bad initial configuration which no sensible barostat can rescue. You are trying to confine an exploding grenade in a milk carton. Start again from a new initial configuration with sensible densities and no weird overlaps.

If the fluctuations in pressure hit 30,000 atm, then that is simply a feature of nanoscale simulations of water, which is very incompressible.

[mjksho]

Hahaha.
Thank you very much. Your response about bad initial configuration directed me to go check the system again. I found that the fix spring/self command which I used to hold the wall atoms to prevent them from falling apart is the biggest contributor to the high pressure.