Hi,
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.
Best,
Jack
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
###Equilibration
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]
#Variables
variable pressx equal c_fy[1]/vol1.0e-4
variable pressy equal c_fy[2]/vol1.0e-4
variable pressz equal c_fy[3]/vol1.0e-4
variable totalp equal -(c_p[1]+c_p[2]+c_p[3])/(3vol)*1.0e-4
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”