[lammps-users] For free boundary conditions, why stress/atom isn't zero everywhere?

Dear all,

I compute stress/atom with two free surfaces normal to z axis, x and y directions are periodic.

I think, if z is free, atoms are free to move in z direction to bring sigma_zz (stress/atom_zz) to zero. Since in continuum, for free surfaces, Cauchy stress is indeed zero everywhere, from surface to bulk.

I accept that sigma_xx and sigma_yy are not zero because surface stress, but sigma_zz should be zero everywhere.

I also realize if I sum over all blue numbers (see below) in dump file, I do get a value that is close to zero, but theory says it should be zero everywhere.

I am not quite sure what virial stress that LAMMPS calculates is, but I do get the knowledge from somewhere says virial stress is equivalent to Cauchy stress in continuum.

Could someone explain to me where I am wrong?

Best,

Weina

Input file

…………………………….

boundary p p s

lattice fcc 4.08 origin 0 0 0 orient x 1 -1 0 orient y 1 1 0 orient z 0 0 1

region mybox block 0 1 0 1 -1 1 units lattice

create_box 1 mybox

create_atoms 1 region mybox

compute 1 all stress/atom

dump 1 all custom 5000 001stress.dump.* id type x y z c_1[1] c_1[2] c_1[3]

minimize 1e-100 1e-120 10000000 100000000

…………………………….

Dump file

…………………………….

ITEM: ATOMS id type x y z c_1[1] c_1[2] c_1[3]

1 1 2.885 -7.22443e-17 -3.96282 1.5637e+06 1.5637e+06 195121

2 1 4.32749 1.4425 -2.05164 -251372 -251372 -262882

3 1 -5.85237e-17 -3.05374e-18 -3.96282 1.5637e+06 1.5637e+06 195121

4 1 2.9288e-17 2.885 -3.96282 1.5637e+06 1.5637e+06 195121

5 1 1.4425 1.4425 -2.05164 -251372 -251372 -262882

6 1 2.885 2.885 -3.96282 1.5637e+06 1.5637e+06 195121

7 1 4.32749 4.32749 -2.05164 -251372 -251372 -262882

8 1 1.4425 4.32749 -2.05164 -251372 -251372 -262882

9 1 2.885 -2.14746e-16 -2.397e-16 -9253.95 -9253.95 135522

10 1 4.32749 1.4425 2.05164 -251372 -251372 -262882

11 1 -1.97155e-16 6.30718e-17 -2.25604e-16 -9253.95 -9253.95 135522

12 1 -1.15288e-16 2.885 1.09957e-16 -9253.95 -9253.95 135522

13 1 1.4425 1.4425 2.05164 -251372 -251372 -262882

14 1 2.885 2.885 -1.79755e-16 -9253.95 -9253.95 135522

15 1 4.32749 4.32749 2.05164 -251372 -251372 -262882

16 1 1.4425 4.32749 2.05164 -251372 -251372 -262882

17 1 2.885 1.09667e-16 3.96282 1.5637e+06 1.5637e+06 195121

18 1 -1.56657e-16 3.50824e-17 3.96282 1.5637e+06 1.5637e+06 195121

19 1 -6.18136e-17 2.885 3.96282 1.5637e+06 1.5637e+06 195121

20 1 2.885 2.885 3.96282 1.5637e+06 1.5637e+06 195121

…………………………….

Try boundary p p f instead. Or you firstly dump the atomic coordination and run another simulation with ‘read_data’ command. And dont forget to enlarge the box size to avoid interaction of the two ends of z direction if you wish to use ‘boundary p p p’. This way you actually create a vacuum in between the ends.

AC

Hi, Ajing. Thank you very much for your help and advice.

I tried ppf and larger box. For atoms which are far away from surface, they do have a near-zero stress, but surface atoms still have a big sigma_zz.

My concern is shouldn’t stress_zz be zero entire system?

At least in my opinion, stress_zz for surface atoms should be zero. Analog to continuum, free surface’s boundary condition is sigma_zz=0.

Or virial stress is just another stress, we can’t compare it to Cauchy stress in continuum?

Dump file:

ITEM: ATOMS id type x y z c_1[1] c_1[2] c_1[3]

1 1 2.885 7.96951e-07 -40.6829 1.56423e+06 1.56423e+06 194740

2 1 4.32749 1.4425 -38.7716 -254159 -254159 -255979

3 1 7.96951e-07 7.96951e-07 -40.6829 1.56423e+06 1.56423e+06 194740

4 1 7.96951e-07 2.885 -40.6829 1.56423e+06 1.56423e+06 194740

5 1 1.4425 1.4425 -38.7716 -254159 -254159 -255979

6 1 2.885 2.885 -40.6829 1.56423e+06 1.56423e+06 194740

7 1 4.32749 4.32749 -38.7716 -254159 -254159 -255979

8 1 1.4425 4.32749 -38.7716 -254159 -254159 -255979

9 1 2.885 7.96951e-07 -36.7199 -4326.88 -4326.88 67947.7

10 1 4.32749 1.4425 -34.6799 3148.86 3148.86 -6952.01

.

.

.

41 1 2.885 7.96951e-07 -20.4 -0.118811 -0.118811 0.00125513

42 1 4.32749 1.4425 -18.36 -0.118464 -0.118464 -8.77682e-05

PS. Or you firstly dump the atomic coordination and run another simulation with ‘read_data’ command. I did this, still non-zero stress_zz. What’s special with this move?

Best,

Weina

Have the forces on those atoms gone to (near) zero?
If so, that may be the best the minimizer can do.

Steve

Hi Weina,

This issue is addressed in:

Model. Simul. Mater. Sci. Eng.; vol. 12, pp S319-S332 (2004)

and your observations seem consistent with what we saw. What follows should not be considered a rigorous explanation but the discrepancy arises in that the virial expression, when averaged over an appropriate atomic ensemble, does indeed correspond with the Cauchy stress. When considered on a per atom basis, the virial exhibits certain behaviors that do not correspond with continuum notions. For instance, in the reference above, we show that integrating over the per atom virial “stress” for only a few atomic layers near the free surface did indeed give zero (w/in statistical limits), even though contributions within a given plane may be finite.

You can read, in the reference above and also in:

Math. Mech. Phys. Sol.; vol. 13, pp. 221-266 (2008)

about efforts to evaluate atomic scale stress in ways that give consistent behavior with continuum notions (note also the references in the two papers above … a fairly large number of very good people have worked on this challenge). In fact, the “Hardy” implementation of evaluating stress in atomic scale simulations (and other tools that go along with this) can be accessed in LAMMPS via the AtC package. Honestly, I am only now familiarizing myself with the implementation in LAMMPS but I encourage you to explore this if you are so inclined. Caveat emptor: this is not a “standard” package (you can read about standard versus user supplied packages in the LAMMPS online documentation).

Happy computing!

Ed

Hi, Ed.

Thanks for your reply and great paper.

Good to know you’re already observed this, cause I’ve been suspecting maybe I did something wrong.

Maybe surface is a transition zone, which is composed of a few atom-layers. As Gibbs dividing surface model, we could find a surface with no volume that is equivalent to real surface zone. This is maybe the theory about.

Thanks again for your help.

Best regards,

Weina

I tried ppf and larger box. For atoms which are far away from surface, they do have a near-zero stress, but surface atoms still have a big sigma_zz.

I would say this is reasonable. In solids, surfaces are several layers with "re-arranged" atomic positions. Continum mechanics does not take care of microstructure, thus the stress around surface layers can not be predicted very well.

Also usually, the average of Virial stress (volume/time/esemble) is equivilant to Cauch stress.

AC