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 stress*volume. 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 stress*volume. Here is the simple conversion from Pij* to

Pij:

Pij* = Pij * Natoms,layer / Vlayer

Here we essentially divide the peratom stress*volume 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 stress*volume.

Vlayer can be found using the box boundaries given in the header of the dumpfile:

ITEM: TIMESTEP

3390000

ITEM: NUMBER OF ATOMS

56760

ITEM: BOX BOUNDS

-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 2*22.8631*2*18*dz.

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,

Chris

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