I try to simulate the pressure drop dynamics of gas hydrates. To do this, I put a box of hydrate (that is simulated separately in higher pressure Phyd) between two boxes of Nitrogen gas (that are simulated separately both with the lower pressure pf Pdiss). Note that the two Nitrogen gas boxes are identical. I expect that hydrate dissociates gradually due to the pressure drop between the boxes. I use periodic boundary conditions with the pair_style"lj/cut/tip4p/long" with the cut-off radius 10 Angstrom.
I run NVT ensemble over the combined system (one box of methane hydrate between two identical boxes of N2). However, the total pressure of the new system does not drop gradually and it even becomes negative! I was wondering why this happens? Should not the total pressure at equilibrium be somewhere between PHyd and Pdiss?
I am thinking about what would happen in a similar situation in reality. If we have a box of higher pressure and just juxtapose it with another box of lower pressure, we expect that the pressure in the higher P box equilibrates gradually but not become negative, correct?
So my question is : how can I model such a behavior in LAMMPS and generally using MD simulations?
You cannot do a 1-to-1 translation from macroscopic behavior to atomic scale.
First, even in the macroscopic system, I would expect that the pressure will balance out very rapidly since the compressed system has suddenly two essentially open surfaces and can expand into it. This will also generate some kind of shockwave from the instant/rapid expansion.
Second, because you are grafting systems on top of each other that were prepared independently with PBC, they suddenly have a massive and instant change in their environment where the new interface is.
Negative pressure can easily happen. This is because of how you measure it. You cannot do the same kind of measurement macroscopically. A negative pressure is simply an indication that the system “wants” to shrink. Please note that in macroscopic environments you do not have periodic boundary conditions. Macroscopically, this would be equivalent to the force on a membrane and that can also have different directions depending on whether the system on one side wants to expand or shrink.
Microscopically, this is all the more complicated since Temperature and Pressure are cleanly defined only for very large systems (many orders of magnitude larger than what you can simulate). While it is technically possible to compute a stress tensor and compute a pressure tensor from it, its correspondence to macroscopic pressure somewhat indirect.
I’m not sure your intended simulation method can be done. But the first thing I’d try is redoing what you did in a slab geometry (boundary p p f for example) with restraining walls.
To see why, try drawing your earlier simulation system with several replicates and think about the interaction the hydrate must (self-consistently) be having with its infinite other copies.