I'm writing a paper on surface tension right now too using LAMMPS. I
suggest that you:
1. Use kspace style Ewald/n for long-range corrections for the LJ
fluid (or just use a long cutoff)
2. Use the mechanical definition of surface tension: gamma =
0.5*Lz*(<Pn> - <Pt>) where the normal pressure Pn=Pz and the
tangential pressure Pt=0.5*(Px+Py). You can easily calculate the
pressure tensor in LAMMPS using something like this:
compute pt all pressure thermo_temp
fix 1 all ave/time \{Nevery\} {Nrepeat} ${Nfreq} c_pt[1] c_pt[2]
c_pt[3] file ptensor.txt
Take a look at the following reference for more details on using the
mechanical definition of surface tension:
Alejandre, J.; Tildesley, D. J. & Chapela, G. A.
MOLECULAR-DYNAMICS SIMULATION OF THE ORTHOBARIC DENSITIES AND
SURFACE-TENSION OF WATER
J. Chem. Phys., 1995, 102, 4574-4583
Stan
I appreciate your suggestion.
I did this few weeks ago excluding the first step, although my cutoff was about 5*sigma for the LJ fluid. FYI, I had an elongated simulation box with a stable liquid slab in the middle and 2 vapor phases on each side. The problem was when I used fix ave/spatial to get Pz as a function of z, where z is the direction normal to the interface, I saw that Pz is not the same everywhere. At both interfaces, there was a pick in Pz plot. Physically, I expect Pz to be the same across the interface. That’s why I decided to calculate Pz (or say gamma) numerically. You can guess the very first thing that I did was waiting more with a hope to see the pick is damped. Didn’t happen though…
Also, This Pz is somehow the vapor pressure, isn’t it? But I couldn’t reproduce Marcus Martin’s (TraPPE force field) vapor pressure for n-alkanes while vapor and liquid densities were successfully reproduced.
Any help/suggestion is appreciated. Thanks.
I appreciate your suggestion.
I did this few weeks ago excluding the first step, although my cutoff was
about 5*sigma for the LJ fluid. FYI, I had an elongated simulation box with
a stable liquid slab in the middle and 2 vapor phases on each side. The
problem was when I used fix ave/spatial to get Pz as a function of z, where
z is the direction normal to the interface, I saw that Pz is not the same
everywhere. At both interfaces, there was a pick in Pz plot. Physically, I
expect Pz to be the same across the interface. That's why I decided to
calculate Pz (or say gamma) numerically. You can guess the very first thing
that I did was waiting more with a hope to see the pick is damped. Didn't
happen though...
The pick you see in your Pz plot is expected. Using fix ave/spatial
amounts to the Harasima contour, which won't give you the straight
line in Pz you expect; you would need to use the Irving-Kirkwood
contour instead (hard to do in LAMMPS). However, this doesn't matter
because total amount of surface tension is *completely independent* of
this pick, and you can still use the formula I suggested.
Also, This Pz is somehow the vapor pressure, isn't it? But I couldn't
reproduce Marcus Martin's (TraPPE force field) vapor pressure for n-alkanes
while vapor and liquid densities were successfully reproduced.
The *volume-averaged* value of Pz is also independent of the pick and
will also be equal to the vapor pressure. I'd check to make sure
that:
1. You divided by the correct volume (of the bin)
2. You properly included the kinetic part of the pressure tensor
In short, trying to obtain an accurate local (position-dependent)
pressure tensor or surface tension is troublesome, but when you
average over the entire system it is relatively easy. I have lots of
references for this subject--let me know if you have further
questions.
Stan
Yes, as you are saying, I have got to be more careful about the spatial averages. I was multiplying ave(Pz) by ave(number density) of each bin and assuming that it should be z dependent vapor pressure. To check that, I summed over number density across the interface and multiplied that by total volume. The result, I assumed, should be total number of atoms in the system which was not. That’s probably why I couldn’t reproduce vapor pressure.
Thanks for your help. For the sake of educational purposes, I will update this post with my future results.
Right, another good way to check your results is to calculate the
pressure using both "compute pressure" and "fix ave/spatial" and make
sure that they match up.
Stan
Now I can reproduce TraPPE results for both vapor pressure and surface tension of Methane and Ethane. As you mentioned the cutoff distance has to be large enough (at least 6*sig for LJ system) and one has to add the “NORM SAMPLE” keyword to the fix ave/spatial command. When averaging over bins, it is necessary to obtain results consistent with Lammps built-in function for P calculation.
Thanks for your comments.
For Lammps users: to read a description on what “Norm Sample” keyword does, one can refer to Lammps manual, FIX ave/spatial command.