From 313d850dd26e05f772237e7d1694c7adaff73c34 Mon Sep 17 00:00:00 2001 From: Aidan Thompson Date: Mon, 21 Feb 2022 19:38:28 -0700 Subject: [PATCH] Created a Born matrix example for silicon --- examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw | 1 + .../ELASTIC_T/BORN_MATRIX/Silicon/in.elastic | 76 ++++++++++ .../ELASTIC_T/BORN_MATRIX/Silicon/init.in | 38 +++++ .../ELASTIC_T/BORN_MATRIX/Silicon/output.in | 139 ++++++++++++++++++ .../BORN_MATRIX/Silicon/potential.in | 24 +++ 5 files changed, 278 insertions(+) create mode 120000 examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw create mode 100644 examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic create mode 100644 examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in create mode 100644 examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in create mode 100644 examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw b/examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw new file mode 120000 index 0000000000..e575921334 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw @@ -0,0 +1 @@ +../../../../potentials/Si.sw \ No newline at end of file diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic b/examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic new file mode 100644 index 0000000000..b010d9daf9 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic @@ -0,0 +1,76 @@ +# Compute elastic constant tensor for a crystal at finite temperature +# + +include init.in + +# Compute initial state + +include potential.in +run ${nequil} + +# Run dynamics + +include potential.in +include output.in + +run ${nrun} + +# Output final values + +# Average moduli for cubic crystals + +variable C11cubic equal (${C11}+${C22}+${C33})/3.0 +variable C12cubic equal (${C12}+${C13}+${C23})/3.0 +variable C44cubic equal (${C44}+${C55}+${C66})/3.0 + +variable bulkmodulus equal (${C11cubic}+2*${C12cubic})/3.0 +variable shearmodulus1 equal ${C44cubic} +variable shearmodulus2 equal (${C11cubic}-${C12cubic})/2.0 +variable poissonratio equal 1.0/(1.0+${C11cubic}/${C12cubic}) + +# For Stillinger-Weber silicon, the analytical results +# are known to be (E. R. Cowley, 1988): +# C11 = 151.4 GPa +# C12 = 76.4 GPa +# C44 = 56.4 GPa + +#print "=========================================" +#print "Components of the Elastic Constant Tensor" +#print "=========================================" + +print "Elastic Constant C11 = ${C11} ${cunits}" +print "Elastic Constant C22 = ${C22} ${cunits}" +print "Elastic Constant C33 = ${C33} ${cunits}" + +print "Elastic Constant C12 = ${C12} ${cunits}" +print "Elastic Constant C13 = ${C13} ${cunits}" +print "Elastic Constant C23 = ${C23} ${cunits}" + +print "Elastic Constant C44 = ${C44} ${cunits}" +print "Elastic Constant C55 = ${C55} ${cunits}" +print "Elastic Constant C66 = ${C66} ${cunits}" + +print "Elastic Constant C14 = ${C14} ${cunits}" +print "Elastic Constant C15 = ${C15} ${cunits}" +print "Elastic Constant C16 = ${C16} ${cunits}" + +print "Elastic Constant C24 = ${C24} ${cunits}" +print "Elastic Constant C25 = ${C25} ${cunits}" +print "Elastic Constant C26 = ${C26} ${cunits}" + +print "Elastic Constant C34 = ${C34} ${cunits}" +print "Elastic Constant C35 = ${C35} ${cunits}" +print "Elastic Constant C36 = ${C36} ${cunits}" + +print "Elastic Constant C45 = ${C45} ${cunits}" +print "Elastic Constant C46 = ${C46} ${cunits}" +print "Elastic Constant C56 = ${C56} ${cunits}" + +print "=========================================" +print "Average properties for a cubic crystal" +print "=========================================" + +print "Bulk Modulus = ${bulkmodulus} ${cunits}" +print "Shear Modulus 1 = ${shearmodulus1} ${cunits}" +print "Shear Modulus 2 = ${shearmodulus2} ${cunits}" +print "Poisson Ratio = ${poissonratio}" diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in new file mode 100644 index 0000000000..aeac99aaf0 --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in @@ -0,0 +1,38 @@ +# NOTE: This script can be modified for different atomic structures, +# units, etc. See in.elastic for more info. +# + +# Define MD parameters +variable nthermo index 1000 # interval for thermo output +variable nfreq equal ${nthermo} # intrrval for averaging output +variable nevery index 10 # sampling interval +variable nrepeat equal ${nfreq}/${nevery} # number of samples +variable neveryborn index 100 # sampling interval +variable nrepeatborn equal ${nfreq}/${neveryborn} # number of samples +variable nequil equal 10*${nthermo} # length of equilibration run +variable nrun equal 100*${nthermo} # length of equilibrated run +variable nlat index 3 # number of lattice unit cells +variable temp index 888.0 # temperature of initial sample +variable timestep index 0.002 # timestep +variable mass1 index 28.06 # mass +variable adiabatic index 0 # adiabatic (1) or isothermal (2) +variable tdamp index 0.01 # time constant for thermostat +variable seed index 123457 # seed for thermostat +variable a index 5.45 # lattice constant +variable thermostat index 1 # 0 if NVE, 1 if NVT +variable delta index 1.0e-6 # numdiff strain magnitude + +# generate the box and atom positions using a diamond lattice + +units metal + +boundary p p p + +lattice diamond $a +region box prism 0 ${nlat} 0 ${nlat} 0 ${nlat} 0.0 0.0 0.0 +create_box 1 box +create_atoms 1 box +mass 1 ${mass1} +velocity all create ${temp} 87287 + + diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in new file mode 100644 index 0000000000..29c27e6a5b --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in @@ -0,0 +1,139 @@ +# Setup output + +# Stress fluctuation term F + +compute stress all pressure thermo_temp +variable s1 equal c_stress[1] +variable s2 equal c_stress[2] +variable s3 equal c_stress[3] +variable s4 equal c_stress[6] +variable s5 equal c_stress[5] +variable s6 equal c_stress[4] + +variable s11 equal v_s1*v_s1 +variable s22 equal v_s2*v_s2 +variable s33 equal v_s3*v_s3 +variable s44 equal v_s4*v_s4 +variable s55 equal v_s5*v_s5 +variable s66 equal v_s6*v_s6 +variable s33 equal v_s3*v_s3 +variable s12 equal v_s1*v_s2 +variable s13 equal v_s1*v_s3 +variable s14 equal v_s1*v_s4 +variable s15 equal v_s1*v_s5 +variable s16 equal v_s1*v_s6 +variable s23 equal v_s2*v_s3 +variable s24 equal v_s2*v_s4 +variable s25 equal v_s2*v_s5 +variable s26 equal v_s2*v_s6 +variable s34 equal v_s3*v_s4 +variable s35 equal v_s3*v_s5 +variable s36 equal v_s3*v_s6 +variable s45 equal v_s4*v_s5 +variable s46 equal v_s4*v_s6 +variable s56 equal v_s5*v_s6 + +variable mytemp equal temp +variable mypress equal press +fix avt all ave/time ${nevery} ${nrepeat} ${nfreq} v_mytemp ave running +fix avp all ave/time ${nevery} ${nrepeat} ${nfreq} v_mypress ave running +fix avs all ave/time ${nevery} ${nrepeat} ${nfreq} v_s1 v_s2 v_s3 v_s4 v_s5 v_s6 ave running +fix avssq all ave/time ${nevery} ${nrepeat} ${nfreq} & +v_s11 v_s22 v_s33 v_s44 v_s55 v_s66 & +v_s12 v_s13 v_s14 v_s15 v_s16 & +v_s23 v_s24 v_s25 v_s26 & +v_s34 v_s35 v_s36 & +v_s45 v_s46 & +v_s56 & + ave running + +thermo ${nthermo} + +thermo_style custom step temp pe press f_avt f_avp f_avs[*] f_avssq[*] +thermo_modify norm no + +# bar to GPa +variable pconv equal 1.0e5/1.0e9 +variable cunits index GPa +# metal unit constants from LAMMPS +# force->nktv2p = 1.6021765e6; +# force->boltz = 8.617343e-5; +variable boltz equal 8.617343e-5 +variable nktv2p equal 1.6021765e6 +variable vkt equal vol/(${boltz}*${temp})/${nktv2p} +variable ffac equal ${pconv}*${vkt} + +variable F11 equal -(f_avssq[1]-f_avs[1]*f_avs[1])*${ffac} +variable F22 equal -(f_avssq[2]-f_avs[2]*f_avs[2])*${ffac} +variable F33 equal -(f_avssq[3]-f_avs[3]*f_avs[3])*${ffac} +variable F44 equal -(f_avssq[4]-f_avs[4]*f_avs[4])*${ffac} +variable F55 equal -(f_avssq[5]-f_avs[5]*f_avs[5])*${ffac} +variable F66 equal -(f_avssq[6]-f_avs[6]*f_avs[6])*${ffac} + +variable F12 equal -(f_avssq[7]-f_avs[1]*f_avs[2])*${ffac} +variable F13 equal -(f_avssq[8]-f_avs[1]*f_avs[3])*${ffac} +variable F14 equal -(f_avssq[9]-f_avs[1]*f_avs[4])*${ffac} +variable F15 equal -(f_avssq[10]-f_avs[1]*f_avs[5])*${ffac} +variable F16 equal -(f_avssq[11]-f_avs[1]*f_avs[6])*${ffac} + +variable F23 equal -(f_avssq[12]-f_avs[2]*f_avs[3])*${ffac} +variable F24 equal -(f_avssq[13]-f_avs[2]*f_avs[4])*${ffac} +variable F25 equal -(f_avssq[14]-f_avs[2]*f_avs[5])*${ffac} +variable F26 equal -(f_avssq[15]-f_avs[2]*f_avs[6])*${ffac} + +variable F34 equal -(f_avssq[16]-f_avs[3]*f_avs[4])*${ffac} +variable F35 equal -(f_avssq[17]-f_avs[3]*f_avs[5])*${ffac} +variable F36 equal -(f_avssq[18]-f_avs[3]*f_avs[6])*${ffac} + +variable F45 equal -(f_avssq[19]-f_avs[4]*f_avs[5])*${ffac} +variable F46 equal -(f_avssq[20]-f_avs[4]*f_avs[6])*${ffac} + +variable F56 equal -(f_avssq[21]-f_avs[5]*f_avs[6])*${ffac} + +# Born term + +compute virial all pressure NULL virial +compute born all born/matrix numdiff ${delta} virial +fix avborn all ave/time ${neveryborn} ${nrepeatborn} ${nfreq} c_born[*] ave running + +variable bfac equal ${pconv}*${nktv2p}/vol +variable B vector f_avborn*${bfac} + +# Kinetic term + +variable kfac equal ${pconv}*${nktv2p}*atoms*${boltz}*${temp}/vol +variable K11 equal 4.0*${kfac} +variable K22 equal 4.0*${kfac} +variable K33 equal 4.0*${kfac} +variable K44 equal 2.0*${kfac} +variable K55 equal 2.0*${kfac} +variable K66 equal 2.0*${kfac} + +# Add F, K, and B together + +variable C11 equal v_F11+v_B[1]+v_K11 +variable C22 equal v_F22+v_B[2]+v_K22 +variable C33 equal v_F33+v_B[3]+v_K33 +variable C44 equal v_F44+v_B[4]+v_K44 +variable C55 equal v_F55+v_B[5]+v_K55 +variable C66 equal v_F66+v_B[6]+v_K66 + +variable C12 equal v_F12+v_B[7] +variable C13 equal v_F13+v_B[8] +variable C14 equal v_F14+v_B[9] +variable C15 equal v_F15+v_B[10] +variable C16 equal v_F16+v_B[11] + +variable C23 equal v_F23+v_B[12] +variable C24 equal v_F24+v_B[13] +variable C25 equal v_F25+v_B[14] +variable C26 equal v_F26+v_B[15] + +variable C34 equal v_F34+v_B[16] +variable C35 equal v_F35+v_B[17] +variable C36 equal v_F36+v_B[18] + +variable C45 equal v_F45+v_B[19] +variable C46 equal v_F46+v_B[20] + +variable C56 equal v_F56+v_B[21] diff --git a/examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in b/examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in new file mode 100644 index 0000000000..da5e9c452f --- /dev/null +++ b/examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in @@ -0,0 +1,24 @@ +# NOTE: This script can be modified for different pair styles +# See in.elastic for more info. + +# we must undefine any fix ave/time fix before using reset_timestep +if "$(is_defined(fix,avsigma))" then "unfix avsigma" +if "$(is_defined(fix,avsigmasq))" then "unfix avsigmasq" +reset_timestep 0 + +# Choose potential +pair_style sw +pair_coeff * * Si.sw Si + +# Setup neighbor style +neighbor 1.0 nsq +neigh_modify once no every 1 delay 0 check yes + +# Setup MD + +timestep ${timestep} +fix 4 all nve +if "${thermostat} == 1" then & + "fix 5 all langevin ${temp} ${temp} ${tdamp} ${seed}" + +