Finding elastic constants for Gold using molecular dynamics method

Continuing my previous blog here I will be showing how to calculate elastic constants of metals (illustrating gold as an example here). As stated earlier, 2 separate simulations will be performed in order to calculate the elastic constants.

Input script-Simulation 1 (for calculating C_11, C₁₂ )

                   (Tensile loading condition)

units           metal
atom_style      atomic
boundary        p p p



#FCC unit cell
variable a index 4.0782
lattice fcc ${a}
region          myreg block     0 10 &
                                0 10 &
                                0 20
#Lattice units is used, which is the  default
##in the argon example, we used box units.
create_box      1 myreg
create_atoms 1 region myreg
# Minimize  using Embedded atom method potential for Al

pair_style      eam/alloy
pair_coeff * * Au.eam.alloy Au
minimize 1e-25 1e-25 100000 100000
#=========================
#Total time is 10ps
#Total strain is 0.02
#Strain rate is 0.02 per 10 ps, which is 0.02x10^(11) strains per second
#2x10^9 strains per second!! Very high
#=========================
timestep 0.001
reset_timestep 0
compute str all pressure NULL virial
variable il equal lz
variable ezz equal (lz-${il})/${il}
fix defz all deform 1 z delta 0.0 1.63128 units box
dump dumpid all custom 100 deformtensile.lammpstrj id x y z
thermo_style custom step time c_str[1] c_str[2] c_str[3] c_str[4] c_str[5] c_str[6] v_ezz lx ly
thermo 100
run 10000
#end========================

Input script-Simulation 2 (for calculating C₄₄ )

                   (Shear loading condition)

units           metal
atom_style      atomic
boundary        p p p

#FCC unit cell
variable a equal 4.0782
#For shear deformation, we need a triclinic box
lattice custom ${a} a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.5 0.0 0.5 basis 0.0 0.5 0.5
region          myreg prism    0 10 &
                                0 10 &
                                0 10 &
0 0 0
#Lattice units is used, which is the  default
##in the argon example, we used box units.
create_box      1 myreg
create_atoms 1 region myreg basis 1 1 basis 2 1 basis 3 1 basis 4 1
# Minimize  using Embedded atom method potential for Al

pair_style      eam/alloy
pair_coeff * * Au.eam.alloy Au
minimize 1e-25 1e-25 100000 100000
#=========================
#Total time is 10ps
#5° shear strain was given.
#=========================
timestep 0.001
reset_timestep 0
compute str all pressure NULL virial
variable ily equal ly
variable ilx equal lx
variable ilz equal lz
variable exy equal xy/${ily}
variable exz equal xz/${ilz}
variable eyz equal yz/${ilz}
variable exx equal (lx-${ilx})/${ilx}
variable eyy equal (ly-${ily})/${ily}
variable ezz equal (lz-${ily})/${ilz}
fix defxy all deform 1 xy delta 3.5679 units box
dump dumpid all custom 200 deformshear.lammpstrj id x y z
thermo_style custom step time c_str[1] c_str[2] c_str[3] c_str[4] c_str[5] c_str[6] v_exx v_eyy v_ezz v_exy v_exz v_eyz
thermo 100
run 10000
#end=======================
Visit LAMMPS Manual for detailed explantion of commands used.
Potential file used: Au.eam.alloy

Results

From simulation 1:

The slope of the above plot gives C₁₁ value. So the value is 196.91 GPa.

The slope of the above plot gives C₁₂ value. So the value is 155.97 GPa.

Values:

σ₁= σ₂ (GPa)  σ₃ (GPa)       v_ezz

0.026791706 0.026791706 0

0.059960138 0.067417747 0.0002

0.093415342 0.10842185 0.0004

0.12682225 0.14939341 0.0006

0.16018103 0.19033207 0.0008

0.19349187 0.23123846 0.001

0.22675479 0.27211226 0.0012

0.25996992 0.31295347 0.0014

0.29313723 0.35376247 0.0016

0.32625671 0.39453905 0.0018

0.3593283 0.43528247 0.002

0.39235223 0.475993        0.0022

0.42532932 0.5166719 0.0024

0.45825929 0.55731886 0.0026

0.4911424 0.5979336 0.0028

0.52397787 0.6385149 0.003

0.55676547 0.67906321 0.0032

0.5895054 0.71957877 0.0034

0.62219757 0.76006116 0.0036

0.65484255 0.80051051 0.0038

0.68744037 0.84092687 0.004

0.71999111 0.88131055 0.0042

0.75249497 0.92166156 0.0044

0.78495232 0.96197993 0.0046

0.81736289 1.0022655 0.0048

0.84972676 1.0425179 0.005

0.88204365 1.0827372 0.0052

0.91431336 1.122923         0.0054

0.94653609 1.1630757 0.0056

0.97871182 1.2031952 0.0058

1.010841         1.2432815 0.006

1.0429238 1.2833344 0.0062

1.0749602 1.3233543 0.0064

1.1069501 1.3633411 0.0066

1.1388934 1.4032945 0.0068

1.1707901 1.4432138 0.007

1.2026404 1.4830994 0.0072

1.2344443 1.5229519 0.0074

1.2662021 1.5627712 0.0076

1.2979138 1.6025569 0.0078

1.3295792 1.6423086 0.008

1.3611983 1.6820268 0.0082

1.3927711 1.7217116 0.0084

1.4242974 1.7613617 0.0086

1.4557776 1.8009775 0.0088

1.4872123 1.8405601 0.009

1.5186013 1.8801093 0.0092

1.5499447 1.9196246 0.0094

1.5812419 1.9591052 0.0096

1.6124932 1.9985515 0.0098

1.6436986 2.0379639 0.01

1.6748582 2.077342         0.0102

1.7059722 2.1166861 0.0104

1.7370406 2.1559957 0.0106

1.7680633 2.195271         0.0108

1.7990405 2.2345121 0.011

1.8299724 2.2737189 0.0112

1.8608588 2.3128913 0.0114

1.8916996 2.3520289 0.0116

1.9224948 2.3911313 0.0118

1.9532444 2.4301986 0.012

1.9839484 2.4692311 0.0122

2.0146073 2.5082293 0.0124

2.0452211 2.5471928 0.0126

2.07579         2.5861211 0.0128

2.106314         2.6250151 0.013

2.136793         2.6638744 0.0132

2.1672268 2.7026985 0.0134

2.1976152 2.7414867 0.0136

2.2279585 2.7802395 0.0138

2.2582567 2.8189567 0.014

2.2885099 2.8576384 0.0142

2.3187182 2.8962849 0.0144

2.3488817 2.9348959 0.0146

2.3790008 2.9734717 0.0148

2.4090755 3.0120125 0.015

2.4391058 3.0505183 0.0152

2.4690908 3.0889878 0.0154

2.4990304 3.1274205 0.0156

2.5289247 3.1658159 0.0158

2.5587741 3.2041754 0.016

2.5885801 3.2425004 0.0162

2.618342         3.2807897 0.0164

2.6480603 3.3190434 0.0166

2.6777338 3.3572605 0.0168

2.7073623 3.395441   0.017

2.7369459 3.4335848 0.0172

2.7664843 3.4716911 0.0174

2.7959793 3.5097615 0.0176

2.8254306 3.547796         0.0178

2.8548386 3.5857949 0.018

2.8842024 3.623757         0.0182

2.9135213 3.6616815 0.0184

2.9427955 3.6995685 0.0186

2.9720246 3.7374175 0.0188

3.0012102 3.77523         0.019

3.0303524 3.8130062 0.0192

3.0594511 3.850746         0.0194

3.0885061 3.8884487 0.0196

3.1175165 3.9261129 0.0198

3.1464828 3.9637397 0.02

From simulation 2:

The slope of the above plot gives C₄₄ value. So the value is  51.077 GPa.

Values:

σ₄                v_exy

4.27E-15         0

0.042230287 0.000874871

0.084888963 0.001749743

0.12754976 0.002624614

0.1702154 0.003499485

0.2128861 0.004374356

0.25556372 0.005249228

0.29825012 0.006124099

0.34094599 0.00699897

0.38365294 0.007873841

0.42637245 0.008748713

0.46910591 0.009623584

0.5118544 0.010498455

0.55462027 0.011373326

0.59740307 0.012248198

0.64020635 0.013123069

0.68303002 0.01399794

0.72587568 0.014872812

0.76874544 0.015747683

0.81163939 0.016622554

0.85456031 0.017497425

0.89750878 0.018372297

0.94048628 0.019247168

0.98349385 0.020122039

1.0265336 0.02099691

1.0696063 0.021871782

1.1127135 0.022746653

1.1558568 0.023621524

1.1990369 0.024496395

1.2422557 0.025371267

1.2855147 0.026246138

1.3288144 0.027121009

1.3721568 0.027995881

1.4155435 0.028870752

1.4589747 0.029745623

1.5024532 0.030620494

1.5459791 0.031495366

1.5895542 0.032370237

1.6331803 0.033245108

1.6768575 0.034119979

1.7205882 0.034994851

1.7643735 0.035869722

1.8082139 0.036744593

1.8521122 0.037619464

1.896068  0.038494336

1.9400836 0.039369207

1.9841605 0.040244078

2.0282989 0.04111895

2.0725011 0.041993821

2.1167681 0.042868692

2.1611005 0.043743563

2.2055008 0.044618435

2.2499691 0.045493306

2.2945071 0.046368177

2.3391164 0.047243048

2.3837974 0.04811792

2.4285522 0.048992791

2.4733817 0.049867662

2.5182866 0.050742533

2.563269  0.051617405

2.6083297 0.052492276

2.6534696 0.053367147

2.6986909 0.054242019

2.7439932 0.05511689

2.7893797 0.055991761

2.8348496 0.056866632

2.8804053 0.057741504

2.9260482 0.058616375

2.9717782 0.059491246

3.0175977 0.060366117

3.0635075 0.061240989

3.1095079 0.06211586

3.155602  0.062990731

3.2017892 0.063865602

3.2480708 0.064740474

3.2944493 0.065615345

3.3409242 0.066490216

3.3874976 0.067365088

3.4341711 0.068239959

3.4809439 0.06911483

3.5278193 0.069989701

3.5747973 0.070864573

3.6218788 0.071739444

3.669066   0.072614315

3.7163587 0.073489186

3.7637586 0.074364058

3.8112679 0.075238929

3.8588853 0.0761138

3.9066139 0.076988671

3.9544545 0.077863543

4.0024067 0.078738414

4.0504739 0.079613285

4.0986554 0.080488157

4.1469524 0.081363028

4.1953676 0.082237899

4.2438998 0.08311277

4.2925517 0.083987642

4.3413241 0.084862513

4.390217         0.085737384

4.4392328 0.086612255

4.4883724 0.087487127

Explanation

After calculating the C values, B, G, E are calculated using the below formula:

Bulk modulus (B) = (C₁₁ + 2 C₁₂)/ 3

Shear modulus (G) = (C₁₁ – C₁₂ + 3 C₄₄)/5

Young’s modulus (E) = (9 GB)/ (3 B + G)

Using the formula the values are calculated:

E= 108.23 GPa

G= 38.83 GPa

B= 169.61 GPa

Actual Values:

E= 78 Gpa

G= 27 GPa

B= 220 GPa

As we see, there are deviations in values, but the values are comparable. This computational study will give us a general idea of how a material behaves. This is computationally very small system on which the study was conducted (with just 8000 atoms). A larger system can give even close results. Hope this blog is useful for you. Keep visiting the page for further updates!