[lammps-users] surface tension values of polymers from md

Hi all,
I am surprised at the low values I am calculating for the surface tension of my polymer system. I have simple 10 bead chains comprised of LJ particles. I have used the following LJ parameters: epsilon=1.6 ; sigma=0.84. I treat the unit of energy as kT at 298 and the unit of length as 1 nm. For reference, the values for water are [ sigma=.31nm and eps=0.26 kT(298)]. The surface tension values I measure (with an air interface) are ~ 4dyne/cm or (milliNewton/meter). The surface tension of water is about 80 dyne/cm. The surface tension values are about an order of magnitude less that where I expected them to be.

So, that’s a summary of the problem. Here are some details about how I do the calculation.

The surface tension is determined by summing/integrating:

Sum_over_z of [(Pxx+Pyy)/2 - Pzz]dz

, where Pii(z) is the ii component of the stress tensor as a function of z and dz is the bin length normal to the interface. The Pii values are found using the following lines in my input file:

compute atomstress all stress/atom
fix 1 ave ave/spatial 2 500 1000 z lower .00390625 c_atomstress[1] c_atomstress[2] c_atomstress[3] units reduced file ptens.dat

The data written to my ptens.dat is as follows:

Spatial-averaged data for fix avespat and group motion

TimeStep Number-of-layers

Layer Coordinate Natoms c_atomstress c_atomstress c_atomstress

We have a polymer sim in which the peratom stress was averaged in each of 256 vertical slices.
As lammps documentation details, the peratom stresses are really in units of stressvolume. Let’s give this lammps value the symbol Pij.
To determine the surface tension, we will need the actual Pij values, not peratom and not in units of stressvolume. Here is the simple conversion from Pij to

Pij* = Pij * Natoms,layer / Vlayer

Here we essentially divide the peratom stressvolume by the volume per atom. If you prefer, you can think of it as multiplying the peratom value by number of atoms, getting a total pressure in the layer and then dividing out the volume component of the stressvolume.
Vlayer can be found using the box boundaries given in the header of the dumpfile:

-22.8631 22.8631
-18 18
-20.0837 30.9

The x,y,z box lengths are obtained from the last three lines. dz from the expression above is thus (30.9+20.0837)/256. Vlayer is 222.8631218dz.

What is the mapping to physical dimensions of our reduced surface tension values? Surface tension units = [ E ] / [L2] = kT / nm2 = 4.11 X 10-3 N/m
The common unit of surface tension is dyne/cm = 0.001 N / m. Therefore, units of surface tension in our sims will equal approximately 4 dyne/cm. So my final answer is simply scaled (multiplied) by 4.11 to get dyne/cm.

After all that, I get surface tension values about an order of magnitude less that is expected.

I have attached a plot of [(Pxx+Pyy)/2 - Pzz] * Natoms,layer/Vlayer for the z range covering the air interface. The peak is beautiful, but the magnitudes are off. Any ideas on what’s going wrong here?

Thanks and happy lammping,

(Attachment File:Tens_calc_t1.4_redblck.png is missing)

Hi, Chris:

Assuming there are no problems with unit conversions going on, the most likely culprit is the potential you used, including the cutoff used. The cutoff can have a significant effect on the value of surface tension, as the effect of the long-range correction needs to be taken into account. (These can be 30-40% of the total value.) In addition, if the model is uncharged, that could lead to significant differences compared with a more chemically detailed model.

Another issue which needs to be considered is whether you’re coarse-grained system accurately represents z polymer at 298 K, or if the effective temperature is in fact much higher. If the system is acting like it would at higher temperatures–say 350 or 400 K, the resulting surface tension will be lower than you would expect.

Also, note that even when corrections are included, the calculated surface tensions are often 15-20% lower than experiment. (There are more details about these issues in in 't Veld, Ismail, and Grest, JCP (2007).)


Rather than do the unit conversion early, you could also try to compute
the surface tension in LJ units, then do one final conversion to real
units. That would be a consistency check on what you've done.