LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] Young's modulus
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] Young's modulus

From: J Moon <answodbs43@...24...>
Date: Wed, 2 May 2018 14:16:26 -0700


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]/vol*1.0e-4
variable pressy equal c_fy[2]/vol*1.0e-4
variable pressz equal c_fy[3]/vol*1.0e-4
variable totalp equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)*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"