Young's modulus


I am trying to calculate Young’s modulus of crystalline silicon.

I am aware that there’s the “Elastic” example provided in the package but I wanted to write my own input file using fix deform and compute stress/atom commands and learn how to compute mechanical properties.

I calculated pressure in the x-directon using the compute stress/atom command and used the fix deform command to deform my structure in x direction and recorded the corresponding strain.

Then, I calculated the E = pressure in x/ strain in x and I got ~90 GPa vs ~160 GPa (expected value).

Could you please take a look at my input file and see what could be wrong?

I’d really appreciate your help.



bulk Si via tersoff

units metal
atom_style atomic

dimension 3
boundary p p p
lattice diamond 5.431
region box block 0 8 0 8 0 8
create_box 1 box
create_atoms 1 box

pair_style tersoff
pair_coeff * * Si.tersoff Si
mass 1 28.06
group every region box
neighbor 2.0 multi
neigh_modify every 2 delay 4 check yes

##Simulation Variables
variable temp equal 1 # temperature to be set
variable dt equal 0.0005 # timestep to be used 0.5 fs
variable volume equal 82018.0367 # molecule volume
variable L equal 43.448 #when I used L to be the last lx of the equilibration run 1, it didn’t make differences.
###Do a quick initial run to see initial state of the system

dump initial all atom 1000 ini.atom
timestep ${dt}
thermo_style custom step temp press pe ke etotal vol pxx pyy pzz lx ly lz
run 0
undump initial

###Minimize Energy

min_style cg
minimize 0 1.0e-6 500 5000


velocity all create 1 87625467 rot yes mom yes dist gaussian loop geom # Create random velocity distribution using random seed
reset_timestep 0

##Equilibration Run 1
dump NPT all atom 50 NPT.atom
fix 3 all npt temp 1 1 0.05 iso 0.0 0.0 0.5
thermo 50
run 100000
unfix 3
undump NPT

compute peratom all stress/atom NULL #Units are pressure*volume
compute fy all reduce sum c_peratom[1] c_peratom[2] c_peratom[3] c_peratom[4] c_peratom[5] c_peratom[6]
compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3]


variable pressx equal c_fy[1]/vol1.0e-4
variable pressy equal c_fy[2]/vol
variable pressz equal c_fy[3]/vol1.0e-4
variable totalp equal -(c_p[1]+c_p[2]+c_p[3])/(3

variable strainx equal (lx-v_L)/v_L
variable strainy equal (ly-v_L)/v_L
variable strainz equal (lz-v_L)/v_L

variable tmp equal temp

fix 4 all npt temp 1 1 0.05 y 0 0 1000 z 0 0 1000 drag 2
fix 5 all deform 1 x erate 1e-5 units box remap x
fix def1 all print 50 “{pressx} {pressy} {pressz} {totalp} {strainx} {strainy} {strainz} {tmp}” file def1.txt screen no
dump Deform all atom 100 deform.atom
run 1000000
unfix 4
unfix 5
unfix def1
undump Deform
print “All done”

Depending on strain rate, temp and loading type (plane stress or plane strain) you can get different modulus values.

Are you using the same loading and BCs to get 160 GPa?