[lammps-users] question about dump custom stress

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I have a question regarding dump custom sxx, syy, etc.

I am making uniaxial elongation of a bulk polyethylene box using
fix/uniaxial. My final strain in z is 0.5.(50 % elongation in z) and I
simulate a system of total 3000 atoms for 10 ns. I dump atomic stress values
(sxx, syy, etc.) and when I check these values, they are very large. To calculate the
normal stress along x at one time step, I add up all the x stress values for
all atoms at that time step (I calculate the other directions in the sam
way) I divide this value by total volume of the box since it's bulk
simulation. But the stress values are too high (on the order of teraPa-way
too high for polyethylene) Also the stress-strain curve does not start at
zero, there's stress even at the first time step. I know the value that
lammps gives includes a kinetic term and the inter-molecular interaction
term but not the intra-terms. I wrote down a code to calculate virial from
the intra-molecular interactions. The values from this code is much smaller
than the LAMMPS contribution. As a result, the overall stress-strain curve
is dominated by lammps results.

I was wondering what I'm doing wrong here. The units of virial is kcal/mol,
right?

Any help is appreciated, thanks so much,

Sezen

2

Sezen,

At the bottom of this page: http://lammps.sandia.gov/doc/dump.html

we read:

“It would need to be divided by a per-atom volume to have units of stress, but an individual atom’s volume is not easy to compute in a deformed solid.”

So, no, the units are not kcal/mol. The units are atm (assuming you’re using real units).

LAMMPS is essentially assuming a per-atom volume of 1 cubic angstrom (again, assuming you’re using real units), so you still have to normalize by the number of cubic angstroms per atom that you really want.

I hope this helps.

If you want to look at the source code, check out these lines in fix_stress.cpp:

// convert to pressure units (actually stress/volume = pressure)

double nktv2p = force->nktv2p;
for (i = 0; i < nlocal; i++) {
stress[i][0] *= nktv2p;
stress[i][1] *= nktv2p;
stress[i][2] *= nktv2p;
stress[i][3] *= nktv2p;
stress[i][4] *= nktv2p;
stress[i][5] *= nktv2p;
}
}

And note that nktv2p is 68589.796 for real units — see the following line in update.cpp:

force->nktv2p = 68589.796;

Note that:

kcal/(Namolangs^3) = 6.85684316 * 10^4 atm

where Na is Avagadro’s number.

Paul