Dear developers of LAMMPS,
These days I am using a system as shown below to explore water-decane interaction. I apply OPLS-AA force field to decane, use large cut off for LJ pair potential instead of applying tail correction, use compute stress/atom to calculate local pressure, and average pressure curve along horizontal direction from 0.5 ns to 1 ns. The input file is shown at the bottom (remove some output lines for brevity), and a figure plotting the pressure curve along horizontal direction is shown in the attachment. I find that bulk water pressure is always about 85-100 atm higher than decane pressure if I use rigid water models such as SPC/E. For flexible water models such as SPC/Fw, this pressure difference is reduced but still exists (bulk water pressure is about 15-25 atm lower than decane pressure). I have tried to change simulation box size, LJ cutoff, run time, and shake tolerance, but the difference could not be apparently reduced.
The fix shake and fix rigid commands in LAMMPS compute
a virial term which is due to the constraint forces
they add to the system. That contribution is include
in the compute pressure and compute stress/atom calculations
(unless you turn it off, see the “fix” keyword option).
As far as I know, those calculations are correct. You could
test it by modeling a box of all water or all decane with and w/out
SHAKE constraints, with NPT, and see if you get densities
reported in the literature.
Note that in water like SPC (and presumably decance) models
I don’t think you are free to choose your own cutoff, if you want
to reproduce what it was parameterized for. Changing the cutoff
will change the pressure, and thus the density.
Thank you very much for your reply. But I think I haven’t expressed my idea clearly in last email. What I study is a multicomponent multiphase problem. And I use a system like this:
If you can get the pressure/density correct for pure decane or pure water,
then it is not clear there is any issue with LAMMPS. It doesn’t
do anything different for a 2 phase system. How do you know
there is anything wrong with the result you are getting?
To summarize, the observed pressure differences between the interiors of the two phases are :
Timestep (fs) SHAKE deltaP (atm)
1.0 yes 90
1.0 yes 30
0.5 no -20
0.5 no -10
You did not provide error bars, but I will assume these numbers are statistically distinct from 0. This is not too surprising. In a non-homogeneous system, especially one where local dynamics are coupled to global properties (kinetic energy, virial), you are always going to see subtle things like this, if you look carefully enough, which you did. Well done. The best thing would be to use fix nve and fix langevin, so that there are no constraints coupled to global properties. However, you might also want to take a step back and ask “What am I really trying to simulate?” In an inhomogeneous system, pressure is always a somewhat ill-defined quantity.
One more detail: your tdamp and pdamp values are very large. It is possible that your temperature profile is also not uniform.
Thank you very much. I really appreciate your help.
I want to study some surface forces’s effect on wettability in nanoscale. To make my work recognized by researchers specialized in mechanics, I’d better validate Young-Laplace equation first. I would be seriously judged if I could not reduce the unphysical pressure difference. That’s why I care about this.
Yesterday I tried fix nve and fix langevin as you suggested, but the pressure profile becomes very strange (highly fluctuated and SPC/E water pressure is 200 atm lower than decane). Maybe I should run a period of npt to relax the simulation box before nve. And I plot the temperature profile and it is truly not uniform. I reduce Tdump and Pdump, and temperature profile becomes more uniform though there is no much difference to pressure profile.
I think flexible water & small time step & small Tdump have reached good enough results. If I reduce time step until no much variation occurs, could I believe my results? Is there any subtle things could not be reduced by small time step?
Looking forward to your reply.
Another thing that might help is making the system larger.
Thank you very much for your reply. According to Young-Laplace equation, pressure difference between two phases should be zero if interface curvature is zero. That’s why I doubt my results. And I want to reduce such unphysical pressure difference because I want my work recognized by researchers specialized in mechanics.
I have tried to enlarge my simulation box before, but that pressure difference could not be reduced. Actually, I find the interface thickness is no larger than 2 nm. I think my system is large enough now. Based on this premise, maybe small time step could lead to convincing results?
I find that the method I calculated temperature curve before is not accurate. Actually temperature variation is already within 1 K for Tdump=100, Pdump=200.