[lammps-users] Questions regarding per-atom and global stress tensor calculation

Dear All,

I recently have started working with LAMMPS and having some confusions about the per-atom stress tensor and global stress tensor outuputs from LAMMPS. I am trying to run a simple case of uniaxial tensile deformation of FCC copper crystal. The input script is shown below.

My first problem is about the “compute sum” command which I used to compute the global stress tensor.
The stress output at the beginning of the simulation (dumps.stress.0) shows that the diagonal components of the per-atom stress tensor of each atom are 0. 329042. So for 50,000 atoms, the sum of the per-atom stress tensor or the diagonal components of the global stress tensor should be equal to 16452.1 (=0.329042*50,000) each. I also did the summation by writing a separate routine and got the same result.
However, the screen output for the global stress tensor shows different values for those diagonal components, 4113.0301 each.
Moreover, according to the LAMMPS documentation, the sum of the diagonal components of the per-atom stress tensor for all atoms should give the pressure value P once divided by 3V. The screen output does not agree with that. For example:

-(4113.0301+ 4113.0301+ 4113.0301)/(3*590520.42) is not equal to -0.027860375
However, if I use the value of 16452.1 for each of those diagonal components it gives the right pressure value. Am I missing something or is there anything wrong with the “compute sum” command?

********* Screen output of the global stress tensor, volume and pressure at the beginning of the simulation (0 step) ***************************
2[1] 2[2] 2[3] 2[4] 2[5] 2[6] Volume Press
4113.0301 4113.0301 4113.0301 -1.3859986e-08 -5.6610679e-08 -6.0104495e-08 590520.42 -0.027860375

When I run the fully patched current version (6 Nov 07), with the
script attached,
I get this for thermo output

       0 0 16452.12 16452.12 16452.12
590520.42 -0.027860375

which is not what you say is a problem with values of 4113. In fact
16452/V = 0.02786
as it should. So it seems per-atom stress is consitent with the
global pressure tensor.

For your 2nd question, you applied a tensile stretch in x, and a
compress in y,z. The
sum of the 3 diagonal components barely changes and the pressure stays about the
same. You are asking why xx gets big quickly and yy and zz go
negative. I think
that is normal. The changes may look large, but they are summed over
50,000 atoms
and the pressure is near 0 initially when you have a perfect lattice.
So I think it's fine.

Your 3rd Q - you can use fix npt on the y,z dimension, but only if you
don't constrain
them with the fix deform command, which your script does by "volume y volume z".
If you only fix deform on x, then see the fix npt doc for how to apply
pressure control
to only the y,z dims.

Steve

in.bug (527 Bytes)

Hi all,
I wanted to create a simulation box of alloy. For pure metal I was using create_atoms. I can't find out how to define alloy.
Please advise.

Thanks,
Sandip

Use a unit cell and 2 (or more) basis atoms in the cell. The lattice
command lets you specify that. The create_atoms command used
with the basis keyword lets you populate the unit cell, one species
at a time. Use it multiple times.

Steve

For creating a structure with more than a single atom basis, such as the 3C SiC models I have been running - I use the following commands:

lattice diamond 4.037685 # set up a diamond lattice, whose basis will be defined later
region box block 0 50 0 50 0 60 # set up a 50x50x60 cell box to fill
create_box 2 box # initialize box
create_atoms 2 box basis 1 2 basis 2 2 basis 3 2 basis 4 2 basis 5 1 basis 6 1 basis 7 1 basis 8 1 # create 3C SiC unit cell (2 atom types)

Hope that helps a bit. The tricky part was sorting out the basis command, and making sure that the basis atoms are in the right places. Starting with a 1x1x1 cell system really helps with that.

Dave