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

Re: [lammps-users] Young's modulus

From: "Chowdhury, Sanjib Chandra" <sanjib@...2437...>
Date: Thu, 3 May 2018 01:23:44 +0000

Experimental E that is reported in the references is usually determined from the quasi-static loading (strain rate < 0.1/s) which is quite impossible to achieve in MD simulations. That’s why MD predicted E values are usually higher than experimental E values. One approach is to extrapolate the high strain rate MD modulus data to predict the low strain rate value. Check the paper that I mentioned in the previous email.





From: J Moon <answodbs43@...24...>
Sent: Wednesday, May 2, 2018 8:53 PM
To: Chowdhury, Sanjib Chandra <sanjib@...2437...>
Subject: Re: [lammps-users] Young's modulus


Hi Sanjib, 


Thanks for pointing that out. As a person outside the field, I am quite new to calculating/measuring mechanical properties. 


If the Young's modulus depends on the strain rate, how do we define a single value for E? For example, when searching for a table of E, books/references only provide one single value. Perhaps, these values are only for linearly elastic materials. 


On the other hand, I thought c-Si was a linearly elastic material. When i changed the strain rate by 4 order of magnitude higher, I got E = 125 GPa vs 90 GPa that I used to get for lower strain rate. 


Would you happen to have some insights into this?


I'd appreciate any comments.






On Wed, May 2, 2018 at 2:51 PM, Chowdhury, Sanjib Chandra <sanjib@...2437...> wrote:

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?




From: J Moon <answodbs43@...24...>
Sent: Wednesday, May 2, 2018 5:16 PM
Subject: [lammps-users] 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]/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"