Hi,
I am wondering what is the unit for the stress tensor that is dumped by LAMMPS when using real units. Dump does not have units option as far as I know. units real is my command at the beginning of the input.
When I look at the website:
S = m.vi.vj + 0.5*SUM{ (ai-bi)Fij }
m in real units is given as g/cm3 in unit command, (ai-bi) is in angstroms and force is defined in kcal/mol.Angstroms in units command. v is in angstroms/fs.
When we divide by the volume later, volume is in angstroms cube.
These units are really not consistent, LAMMPS should be doing some conversions in the final output. I have found efactor, pfactor and boltz in initialize.h. But I am not sure in which units it writes in dump file for the stress tensor.
I have the similar unit problem with the velocities that are dumped. They really should be in angstroms/fs. But when I calculate temperature elsewhere using the initial velocities created by create gauss command to give 300.15K, I don’t get 300.15K but much high T.
When I calculate temperature, I use the following formula which is the same principle as for LAMMPS.
Ekin = ekfactor*(0.5mv^2)
T = 2Ekin/(3boltz)/n
I use eV for kinetic energy, A/fs for volume and amu for atom mass. K for temperature. And I have ekfactor for conversion which I calculated to be 103.64269. I calculate a much higher temperature for the same velocities but LAMMPS print 300K in the output.
I would appreciate if someone can clarify me the units in stress tensor and the velocity.
Thank you,
Burcu
Burcu,
From the dump command documentation here:
http://www.cs.sandia.gov/~sjplimp/lammps/doc/dump.html
We have:
"Note that this formula for stress does not include
virial contributions from intra-molecular interactions
(e.g. bonds, angles, torsions, etc). Also note that
this quantity is the negative of the per-atom pressure
tensor. It is also really a stress-volume formulation.
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.
Computation of stress tensor components requires a
loop thru the neighbor list and inter-processor
communication, so it can be inefficient to dump this
quantity too frequently or to have multiple dump
commands, each with stress tensor attributes."
Looking at the source code of fix_stress.cpp, we see:
// convert to pressure units (actually stress/volume =
pressure)
double nktv2p = force->nktv2p;
for (i = 0; i < nlocal; i++) {
stress[i][0] *= nktv2p;
Which means that LAMMPS is giving you units of
pressure*volume. For real units, this would be
atm*angs^3.
Mass has units of grams/mole (if units are real), not
gm/cm^3. See:
http://www.cs.sandia.gov/~sjplimp/lammps/doc/units.html
Note that nktv2p above does the units conversion.
My guess on the issue you're seeing with the
temperature and velocities not agreeing is that you
don't have the right value for n. If you're using
shake, degrees of freedom are removed and n is no
longer the number of atoms. Otherwise, the issue is
likely with the units conversion. If you're using real
units, volume would be angs^3, ke in kcal/mol, and T
in K. You'd have to use the proper conversion factors
like LAMMPS uses. For real units, you should have:
force->boltz = 0.001987191;
force->mvv2e = 48.88821 * 48.88821;
where the "boltz" is the Boltzmann factor, and mvv2e
converts from mass*velocity^2 to energy (in kcal/mol).
The conversion factors are in update.cpp.
Paul
I would appreciate if someone can clarify me the units in stress
tensor >and the velocity.
The doc page for the dump command doesn't state this explicity (I'll add it),
but all dump outputs are in the units described in the units doc pate. I.e.
if you dump velocities then they are in the same units as all velocities
that are input/output by LAMMPS, namely Ang/fmsec for real units or
Ang/psec for metal units, etc.
For the stress attributes (sxx, syy, etc), as the dump doc says, these
quantities are stress*volume, which means pressure * distance^3. So
for real units, it should be atmospheres * Ang^3.
Hope that helps,
Steve