Dear all,
I am trying to model bulk Si and SiC in LAMMPS using the Tersoff potential under uniaxial tensile loading, then calculate the stress–strain response, and finally determine the elastic/Young’s modulus.
However, I am facing a problem: I cannot obtain the expected Young’s modulus value. For bulk silicon, the expected value is approximately 160 GPa, but my simulation gives around 70 GPa. I have already checked several possible factors, including strain rate, timestep, number of run steps, units, number of atoms, and system size. I also tried to minimize the initial stress, but I think some residual stress may still remain.
I would be very grateful if you could take a look at my LAMMPS code and help me identify and fix the problem.
Best regards,
Sina S.
#Bulk Si
units metal
dimension 3
boundary p p p
atom_style atomic
neighbor 2.5 bin
neigh_modify every 10 delay 0 check yes
geometry
lattice diamond 5.431 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region SiReg block 0 20 0 20 0 10 units box
create_box 1 SiReg
create_atoms 1 region SiReg
mass 1 28.0855
tersoff potential
pair_style tersoff
pair_coeff * * Sinew.tersoff Si
#define groups
#group 1 region SiReg
initial velocity
min_style cg
minimize 1.0e-12 1.0e-12 10000 100000
velocity all create 1.0 12345 mom yes rot no dist gaussian
Temperature and pressure computes
compute myTemp all temp
compute myPress all pressure myTemp
Per-atom outputs for dump
compute disp all displace/atom
compute st all stress/atom NULL
compute peatom all pe/atom
compute keatom all ke/atom
fixes
fix 1 all nve
run nve
timestep 0.001
thermo 100
thermo_style custom step temp pe ke etotal lx ly lz press c_myPress[1] c_myPress[2] c_myPress[3]
thermo_modify format float %12.6f
dump d_nve all custom 100 dump_nve.lammpstrj &
id type x y z &
c_disp[1] c_disp[2] c_disp[3] &
c_st[1] c_st[2] c_st[3] c_st[4] c_st[5] c_st[6] &
c_peatom c_keatom
dump_modify d_nve sort id
run 20000
unfix 1
undump d_nve
#run npt
fix EqDy1 all npt temp 1.0 1.0 10 iso 0.0 0.0 10
dump d_npt all custom 100 dump_npt.lammpstrj &
id type x y z &
c_disp[1] c_disp[2] c_disp[3] &
c_st[1] c_st[2] c_st[3] c_st[4] c_st[5] c_st[6] &
c_peatom c_keatom
dump_modify d_npt sort id
run 20000
unfix EqDy1
undump d_npt
run nvt
fix EqDy2 all nvt temp 1.0 1.0 10
dump d_nvt all custom 100 dump_nvt.lammpstrj &
id type x y z &
c_disp[1] c_disp[2] c_disp[3] &
c_st[1] c_st[2] c_st[3] c_st[4] c_st[5] c_st[6] &
c_peatom c_keatom
dump_modify d_nvt sort id
run 20000
unfix EqDy2
undump d_nvt
Reference values
variable Ly0 equal ly
run 0
variable s0 equal c_myPress[2]
run 0
variable Ly0 equal ${Ly0}
variable s0 equal ${s0}
Strain and stress variables
variable strain equal (ly-v_Ly0)/v_Ly0
variable deflection equal (ly-v_Ly0)
variable area equal lx*lz
variable stressY_bar equal c_myPress[2]
variable stressY_GPa equal -((c_myPress[2]-v_s0)/10000.0)
variable forceY_N equal -c_myPress[2]v_area1.0e-15
Tensile loading in y
thermo 100
thermo_style custom step temp lx ly lz v_strain v_deflection v_stressY_bar v_stressY_GPa v_forceY_N pe ke etotal press
run loading
fix loadeq2 all npt temp 1.0 1.0 10 x 0 0 1 z 0 0 1 couple none
fix def all deform 1 y erate 1.0e-2 units box remap x
dump d_load all custom 100 dump_load.lammpstrj &
id type x y z &
c_disp[1] c_disp[2] c_disp[3] &
c_st[1] c_st[2] c_st[3] c_st[4] c_st[5] c_st[6] &
c_peatom c_keatom
dump_modify d_load sort id
fix out all print 100 “${strain} ${stressY_GPa}” &
file stress_strain_100.txt screen no &
title “strain stress(GPa)”
run 20000
unfix loadeq2
unfix def
undump d_load