I encountered an issue where shear stress components were not properly relaxed when performing lattice relaxation on a cell with a specific crystal orientation using the orient option in LAMMPS’s lattice command. This article shares the phenomenon and its solution.

In short, I was able to resolve this by first fully relaxing a minimal unit cell and then expanding it to the desired size using the replicate command. Ensuring that the initial structure is sufficiently relaxed is crucial for the reliability of subsequent calculations. Incidentally, I overlooked this confirmation and had to redo many of my calculations.

The Problem

Let’s use the lattice relaxation of BCC iron as an example. Here, instead of the typical (100), (010), (001) crystal orientations, we prepare a simulation cell where the X, Y, and Z axes are oriented along the (111), (-110), and (-1-12) directions, respectively. As you can verify, these directions are orthogonal.

Below is the LAMMPS input script that demonstrates the problem.

1
# initialization
2
units metal
3
dimension 3
4
boundary p p p
5
atom_style atomic
6

7
# create geometry
8
lattice bcc 2.85 orient x 1 1 1 orient y -1 1 0 orient z -1 -1 2
9
region box block 0 20 0 20 0 20
10
create_box 1 box
11
labelmap atom 1 Fe
12
create_atoms Fe box
13

14
# potentials
15
pair_style eam/fs
16
pair_coeff * * ./Fe_Ackland04_ZBL.eam.fs Fe
17

18
# calc
19
timestep 0.001
20

21
dump 1 all custom 1000 relax_*.dump id type x y z xs ys zs vx vy vz fx fy fz
22

23
thermo 500
24
thermo_style custom step temp pe pxx pyy pzz pxy pxz pyz
25

26
# Minimize atomic coordinates
27
neigh_modify delay 0 every 1 check yes
28
min_style cg
29
min_modify dmax 0.001
30
minimize 0.0 1.0e-12 10000 10000
31

32
# Minimize atomic coordinates + volume
33
fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 1.0e-5 nreset 100
34
min_modify line quadratic
35
minimize 0.0 1.0e-6 10000 10000

In line 8, the crystal orientation is defined using orient, and in line 9, the simulation cell size is defined simultaneously. A key point is that the thermo_style in line 24 is set to output not only the hydrostatic pressure components (Pxx, Pyy, Pzz) but also the shear stress components (Pxy, Pxz, Pyz).

Running this script (lmp_serial -in relax.cmd) produces the following log output:

1
   Step          Temp          PotEng          Pxx            Pyy            Pzz            Pxy            Pxz            Pyz
2
      1056   0             -256805.04      9934.46        7092.2933      12886.431      4.0878753e-12 -6326.5485     -3.1600352e-12
3
      1392   0             -256819.73     -1.5004335e-07 -6.7607257e-06  6.55549e-06   -5.7846837e-12 -5449.1859     -2.4248676e-11

The relaxation calculation converges in 1392 steps, but a very large residual stress of -5449.18 bar remains in Pxz. This can hardly be called a properly relaxed initial structure.

While the exact cause has not been identified, I suspect that when creating a large simulation domain directly with a specified orient, accumulated rounding errors in atomic coordinates may result in residual shear stress.

The Solution

To solve this problem, I tried a two-step approach:

1. First, create and relax a minimal unit cell (0 1 0 1 0 1) using the region command.
2. Replicate the relaxed unit cell to the desired size (20 20 20) using the replicate command.

The final input script based on this approach is as follows:

1
# initialization
2
units metal
3
dimension 3
4
boundary p p p
5
atom_style atomic
6

7
# create geometry
8
lattice bcc 2.85 orient x 1 1 1 orient y -1 1 0 orient z -1 -1 2
9
region box block 0 1 0 1 0 1
10
create_box 1 box
11
labelmap atom 1 Fe
12
create_atoms Fe box
13

14
# potentials
15
pair_style eam/fs
16
pair_coeff * * ./Fe_Ackland04_ZBL.eam.fs Fe
17

18
# calc
19
timestep 0.001
20

21
dump 1 all custom 1000 relax_*.dump id type x y z xs ys zs vx vy vz fx fy fz
22

23
thermo 500
24
thermo_style custom step temp pe pxx pyy pzz pxy pxz pyz
25

26
# Minimize atomic coordinates
27
neigh_modify delay 0 every 1 check yes
28
min_style cg
29
min_modify dmax 0.001
30
minimize 0.0 1.0e-12 10000 10000
31

32
# Minimize atomic coordinates + volume
33
fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 1.0e-5 nreset 100
34
min_modify line quadratic
35
minimize 0.0 1.0e-6 10000 10000
36

37
# Expand cell
38
replicate 20 20 20
39

40
# Minimize atomic coordinates
41
neigh_modify delay 0 every 1 check yes
42
min_style cg
43
min_modify dmax 0.001
44
minimize 0.0 1.0e-12 10000 10000
45

46
# Minimize atomic coordinates + volume
47
fix 1 all box/relax x 0.0 y 0.0 z 0.0 vmax 1.0e-5 nreset 100
48
min_modify line quadratic
49
minimize 0.0 1.0e-6 10000 10000

The result of the relaxation calculation with this script is as follows:

1
   Step          Temp          PotEng          Pxx            Pyy            Pzz            Pxy            Pxz            Pyz
2
      2135   0             -252931.74     -0.0006457296  -0.0034065106   0.0067741329  -3.3083e-11     1.3272091e-09 -1.2592623e-09
3
      2136   0             -252931.74     -0.0006457296  -0.0034065106   0.0067741329  -3.3083e-11     1.3272091e-09 -1.2592623e-09

It is confirmed that both the hydrostatic and shear stress components have been relaxed to near-zero values.

Conclusion

I learned that when dealing with cells of complex crystal orientations using the orient option in LAMMPS, a simple relaxation calculation may leave unintended residual stresses. In such cases, the issue can be resolved by first relaxing a minimal unit cell and then expanding the system based on it.

This experience reminded me that when performing simulations, it is essential not only to confirm that the relaxation calculation has converged but also to specifically check that physical quantities like the stress tensor have reached their expected values to ensure the validity of the results.