Hi, All!
I want to dump the stress tensor in the simulation, however, lammps provide stress tensor for
every atom, and generally the stress decays very fast, so the storage is too much. I am wondering
how to record the system’s stress tensor as a whole in the log file. Thank you in advance!
How is the total stress tensor different than the system pressure which is printed
out by the thermodynamics? See the earlier mail list message about ensemble stress
for the equation relating the per-atom and total stress (pressure). The total pressure
is the averaged sum of the 3 diagonal components of the pressure tensor. If you
want to sell all 6 components, you should just put in a print statement (probably in the
pressure.cpp file) to print them to the screen when thermodynamics is computed.
Steve
Hi,
I have tried to compute the stress components of the system in dump_custom.cpp(as the atom based quantity are already available there) file by:
Sxx = (1/V) * (Sum over all atoms sxx(per atom based quantity)) ;
and print it and now if i try to find pressure by:
P = -(1/3) * (Sxx + Syy + Szz) ;
it turns out to be almost one magnitude higher than the printed pressure. I am not sure if i am missing something.
I appreciate in help with regards to this.
Thanks
Arnab
Steve Plimpton wrote:
Hi Arnab. The pressure LAMMPS gives in the log file
includes a kinetic contribution --- it looks like you
are not including that term. Here's how LAMMPS does
the pressure calculation:
double kinetic = dof * boltz * temperature / 3.0;
double volume = domain->xprd * domain->yprd *
domain->zprd;
p_component[0] = (kinetic + virial[0]) / volume *
nktv2p;
p_component[1] = (kinetic + virial[1]) / volume *
nktv2p;
p_component[2] = (kinetic + virial[2]) / volume *
nktv2p;
p_total = (p_component[0] + p_component[1] +
p_component[2]) / 3.0;
It looks like you're only including the virial term
and not the kinetic term. Also, are you doing the
units conversion correctly (see nktv2p in the LAMMPS
source code)?
Paul
According to Rob’s reply:
Hi, Zhimin. What I do to get the whole stress tensor is simply to dump
the 3 elements of p_component[] (an array defined in pressure.cpp) to
the log file. If you use thermo_style one, you can easily implement by
adding a few lines to ThermoOne::compute such as
double press = pressure->p_total; // this line is already there
double pressX = pressure->p_component[0];
double pressY = pressure->p_component[1];
double pressZ = pressure->p_component[2];
I think you could add the off-diagonal elements of the pressure tensor
by extending p_component to 6 elements and modifying pressure.cpp to
access the 3,4,5 components of virial[] as appropriate.
Rob
I changed pressure.h, pressure.cpp to add the off-diagnoal pressure tensor, and
then make some modification in thermo_one.cpp to print all the pressure tensor
into the log file. I obtained these tensors and then calculated the stress autocorrelation
function to get the shear viscosity of the lj particle system, but it seemed that
the correlation function decays very fast( about 20 timestep to decay to 0) and the
fluctuation is very large at long time. I am wondering whether something was wrong
during the modification of the code and the calculation of the shear viscosity. Any
suggestion? Thanks!
The diagonal component of pressure tensor in lammps has two part, the kinetic contribution part and the potential contribution part. It seems that every diagonal component has the same kinetic contribution. I am wondering whether it is right, and I doubt my similar extend of the off-diagonal pressure tensor.