Created a Born matrix example for silicon
This commit is contained in:
Aidan Thompson
1
examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw
Symbolic link
1
examples/ELASTIC_T/BORN_MATRIX/Silicon/Si.sw
Symbolic link
@ -0,0 +1 @@
|
||||
../../../../potentials/Si.sw
|
||||
76
examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic
Normal file
76
examples/ELASTIC_T/BORN_MATRIX/Silicon/in.elastic
Normal file
@ -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}"
|
||||
38
examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in
Normal file
38
examples/ELASTIC_T/BORN_MATRIX/Silicon/init.in
Normal file
@ -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
|
||||
|
||||
|
||||
139
examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in
Normal file
139
examples/ELASTIC_T/BORN_MATRIX/Silicon/output.in
Normal file
@ -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]
|
||||
24
examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in
Normal file
24
examples/ELASTIC_T/BORN_MATRIX/Silicon/potential.in
Normal file
@ -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}"
|
||||
|
||||
|
||||
Reference in New Issue
Block a user