Sina,
I looked through your script, and I think the low Young’s modulus may come from a combination of several points rather than only from residual stress.
First, please note that with your current crystal orientation,
lattice diamond 5.431 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
and tensile loading in the y direction, you are calculating the modulus along the [010] direction, which is equivalent to [100] for cubic silicon. For single-crystal silicon, the Young’s modulus is anisotropic. The value along [100] is closer to about 130 GPa, not 160 GPa. A value around 160–170 GPa is more consistent with loading close to the [110] direction. So the expected reference value should first be checked according to the crystal direction.
However, 70 GPa is still too low. One important issue I noticed is your region definition:
region SiReg block 0 20 0 20 0 10 units box
Because you used units box, LAMMPS interprets the dimensions as box-distance units, not lattice-cell units. Therefore, your simulation box is much smaller than a 20 × 20 × 10 lattice-cell crystal. If your intention is to create 20 × 20 × 10 lattice cells, you should use:
region SiReg block 0 20 0 20 0 10 units lattice
This may significantly affect the result, especially because the system becomes much smaller than expected.
Another point is the deformation procedure. During tensile loading, you use fix deform together with fix npt. In this case, it is better to use a deformation-aware temperature compute, for example:
compute Tdef all temp/deform
fix loadeq2 all npt temp 1.0 1.0 10 x 0 0 10 z 0 0 10 couple none
fix_modify loadeq2 temp Tdef
fix def all deform 1 y erate 1.0e-5 units box remap x
This avoids including the affine deformation velocity in the temperature calculation.
I would also reduce the strain rate. Your current strain rate,
fix def all deform 1 y erate 1.0e-2 units box remap x
is very high for extracting an elastic modulus. For modulus calculation, it is better to use a much lower strain rate and fit the stress–strain curve only in the small elastic region, for example 0–0.5% or 0–1% strain.
Your stress conversion is generally correct, since in metal units the pressure is in bar and 1 GPa = 10000 bar:
variable stressY_GPa equal -((c_myPress[2]-v_s0)/10000.0)
The negative sign is also reasonable because LAMMPS pressure is positive in compression and negative in tension.
I also noticed a small syntax issue in the force variable:
variable forceY_N equal -c_myPress[2]v_area1.0e-15
It should be written with multiplication signs:
variable forceY_N equal -c_myPress[2]*v_area*1.0e-15
although this does not directly affect the modulus if you are using the stress variable.
In summary, I suggest the following checks:
-
Use units lattice in the region command if you want 20 × 20 × 10 lattice cells.
-
Check the reference modulus according to the loading direction. Your current loading direction corresponds to [100], not [110].
-
Use compute temp/deform during tensile loading.
-
Reduce the strain rate.
-
Fit only the initial linear elastic part of the stress–strain curve.
-
Make sure the initial stress is taken after the final equilibration step, just before deformation.
If your goal is to obtain a modulus close to 160–170 GPa for silicon, you may want to change the orientation so that the loading direction is [110], for example:
lattice diamond 5.431 orient x 1 -1 0 orient y 1 1 0 orient z 0 0 1
Then tensile loading in the y direction would correspond to the [110] direction.
Best regards,
Bahman