[gmx-users] Strange pressure coupling behaviour, rapid expanding box

[Kamps, M.]

Dear GMX users,

I'm experiencing difficulties with my simulation.
I have created a ionic bonded, rough surface of gold atoms (3276
atoms) in the OPLS FF, which is held together with Lennard-Jones
interactions. In the same box I've placed a C4H10 (14 atoms) molecule,
although the problem also occurs with proteins. This box is then
solved using SPCE water, although I've also tried TIP4P water.

The box size is originally 5.4x5.4x4.5 nm, and is first energy
minimized and NVT equilibrated. Up to this point the system behaves
normal and I can't see any red flags. When I then want to NPT
equilibrate the system, something strange happens. The periodic box
seems to ´explode´, since it rapidly (but constantly) expands. The
final dimensions are 5.4x5.4x24 nm after 50ps.

Analyzing the outputs of the NPT equilibrium shows some strange
things: The initial pressure starts at roughly 14000 bar and rapidly
decreases to around 3000 bar at t=5ps, after which it linearly
decreases to 2000 bar. Since the box is rapidly expanding, the volume
of the system is linearly growing, starting at roughly 160nm3 and
expanding to 850nm3 at 50ps.

I've added two screenshots, one before and one after the simulation.
I've also uploaded a screenshot of my mdp file, since it is easier to
read from an image than in here (i think).
https://s23.postimg.org/yceosj0q3/Before.png
https://s23.postimg.org/9kjceaobv/After.png
https://s24.postimg.org/xrzza3mt1/NPT_mdp.png

My system uses a semiisotrpic pcoupletype, since I want different
behaviour in xy and z direction. The reference pressure has been set
to 1 atmosphere (almost 1 bar), while the compressibility in the xy
direction is that of the material, while in the z direction it is that
of water. I only want periodic boundary conditions in the xy
directions, and both the z planes are LJ 12-6 walls. At z=0 the wall
is behaving like the same gold atoms, with the corresponding LJ
potential, while the upper wall is exerting an LJ potential of a dummy
atom (so no LJ interaction). Might this cause the error? I varied a
lot with the wall-r-linpot option, but putting in an positive,
negative or zero value does not seem to make a difference.

Are there any ideas as to what I am doing wrong? Thanks in advance,

Mark