kim init EAM_Dynamo_MendelevBorovikov_2020_FeNiCr__MO_922363340570_000 metal
include init.mod

# Compute initial state

variable thermostat equal 1
include potential.mod
run ${nequil}

if "${adiabatic} == 1" &
then "variable thermostat equal 0" &
else "variable thermostat equal 1"

print ${thermostat}

include potential.mod
run ${nrun}

variable pxx0 equal f_avp[1]
variable pyy0 equal f_avp[2]
variable pzz0 equal f_avp[3]

variable tmp equal lx
variable lx0 equal ${tmp}
variable tmp equal ly
variable ly0 equal ${tmp}
variable tmp equal lz
variable lz0 equal ${tmp}

# These formulas define the derivatives w.r.t. strain components
variable d1 equal -(v_pxx1-${pxx0})/(v_delta/v_len0)*${cfac}
variable d2 equal -(v_pyy1-${pyy0})/(v_delta/v_len0)*${cfac}
variable d3 equal -(v_pzz1-${pzz0})/(v_delta/v_len0)*${cfac}

# Write restart
write_restart restart.equil

# uxx Perturbation

variable dir equal 1
include displace.mod

# uyy Perturbation

variable dir equal 2
include displace.mod

# uzz Perturbation

variable dir equal 3
include displace.mod

# Output final values

variable C11all equal ${C11}
variable C12all equal 0.5*(${C12}+${C21})
variable C44all equal ${C44}


# Average moduli for cubic crystals

variable C11cubic equal (${C11all}+${C22all}+${C33all})/3.0
variable C12cubic equal (${C12all}+${C13all}+${C23all})/3.0
variable C44cubic equal (${C44all}+${C55all}+${C66all})/3.0

variable bulkmodulus equal (${C11cubic}+2*${C12cubic})/3.0
variable shearmodulus1 equal (${C11cubic}-${C12cubic+3.0* ${C44cubic})/2.0

print "========================================="
print "Components of the Elastic Constant Tensor"
print "========================================="

print "Elastic Constant C11all = ${C11all} ${cunits}"
print "Elastic Constant C12all = ${C12all} ${cunits}"
print "Elastic Constant C44all = ${C44all} ${cunits}"

print "========================================="
print "Average properties for a cubic crystal"
print "========================================="

print "Bulk Modulus = ${bulkmodulus} ${cunits}"
print "Shear Modulus 1 = ${shearmodulus} ${cunits}"