units           metal
dimension       3
boundary        p       p       p
atom_style      atomic
variable latparam equal 3.1652

# ----------------------- ATOM DEFINITION ----------------------------
lattice         bcc ${latparam}
region          whole block 0 31.5 0 31.5  0 31.5
create_box      1 whole
lattice         bcc ${latparam} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms    1 region whole








pair_style       eam/alloy
pair_coeff      * * MoTaWTiZr.set  W

min_style       cg
min_modify      dmax 0.1
thermo          10000
thermo_style    multi
minimize        1e-15 1e-15 1000000 1000000

compute csym all centro/atom bcc
compute peratom all pe/atom
reset_timestep  0

timestep 0.001
velocity all create 4000 123456 dist gaussian








fix 1 all npt temp 4000 4000 1 iso 0 0 1 drag 1

thermo 1000

thermo_style custom step temp press ke pe etotal vol density
run 90000

compute myRDFall_1 all rdf 50
fix 116 all ave/time 100 1 100 c_myRDFall_1[*] file tmp1_all.rdf mode vector






compute         vacf1_all all vacf



fix             5 all vector 1 c_vacf1_all[4]
variable        vacf_all equal dt*trap(f_5)




thermo_style    custom step temp  v_vacf_all  c_msd1_all[4]

thermo          5

run 5000

unfix 1




fix 10 all npt temp 4000 3300 1 iso 0 0 1 drag 1



thermo 10000

thermo_style custom step temp press ke pe etotal vol density
run 2000000
unfix 10







fix 11 all npt temp 3300 3300 1 iso 0 0 1 drag 1


thermo 1000

thermo_style custom step temp press ke pe etotal vol density
run 50000





compute myRDFall_2 all rdf 50
fix 216 all ave/time 100 1 100 c_myRDFall_2[*] file tmp2_all.rdf mode vector








compute         msd2_all all msd com no


compute         vacf2_all all vacf





fix             9 all vector 1 c_vacf2_all[4]
variable        vacf2_all equal dt*trap(f_9)


thermo_style    custom step temp v_vacf2_all  c_msd2_all[4]

thermo          5

run 1000


unfix 11








fix 18 all npt temp 3300 2883.8 1 iso 0 0 1 drag 1



thermo 10000

thermo_style custom step temp press ke pe etotal vol density
run 2000000
unfix 18






fix 19 all npt temp 2883.8 2883.8 1 iso 0 0 1 drag 1


thermo 1000

thermo_style custom step temp press ke pe etotal vol density
run 50000




compute myRDFall_3 all rdf 50
fix 316 all ave/time 100 1 100 c_myRDFall_3[*] file tmp3_all.rdf mode vector








compute         msd3_all all msd com no



compute         vacf3_all all vacf









fix             9 all vector 1 c_vacf3_all[4]
variable        vacf_all equal dt*trap(f_9)


thermo_style    custom step temp v_vacf_all  c_msd3_all[4]

thermo          5

run 1000


unfix 19




fix 20 all npt temp 2883.8 2500 1 iso 0 0 1 drag 1


thermo 10000

thermo_style custom step temp press ke pe etotal vol density
run 2000000
unfix 20




fix 21 all npt temp 2500 2500 1 iso 0 0 1 drag 1


thermo 1000

thermo_style custom step temp press ke pe etotal vol density
run 50000




compute myRDFall_4 all rdf 50
fix 416 all ave/time 100 1 100 c_myRDFall_4[*] file tmp4_all.rdf mode vector








compute         msd4_all all msd com no


compute         vacf4_all all vacf




fix             9 all vector 1 c_vacf4_all[4]
variable        vacf_all equal dt*trap(f_9)


thermo_style    custom step temp v_vacf_all  c_msd4_all[4]

thermo          5

run 1000


unfix 21






fix 22 all npt temp 2500 1100 1 iso 0 0 1 drag 1


thermo 10000

thermo_style custom step temp press ke pe etotal vol density
run 2000000
unfix 22



fix 23 all npt temp 1100 1100 1 iso 0 0 1 drag 1


thermo 1000

thermo_style custom step temp press ke pe etotal vol density
run 50000



compute myRDFall_5 all rdf 50
fix 516 all ave/time 100 1 100 c_myRDFall_5[*] file tmp5_all.rdf mode vector







compute         msd5_all all msd com no


compute         vacf5_all all vacf




fix             9 all vector 1 c_vacf5_all[4]
variable        vacf_all equal dt*trap(f_9)


thermo_style    custom step temp v_vacf_all  c_msd5_all[4]

thermo          5

run 1000


unfix 23






fix 24 all nve
thermo 1000
thermo_style custom step temp press ke pe etotal vol density
run 10000
unfix 24

print "equi done"

#store final cell length for strain calculations
variable tmp equal "lx"
variable L0 equal ${tmp}
print "Initial Length, L0: ${L0}"

variable a loop 5
label loop

fix 25 all npt temp 1100 1100 1 iso 0 0 1 drag 1
thermo 1000
thermo_style custom step temp press ke pe etotal vol density
run 90000
unfix 25

fix 26 all nve
thermo 1000
thermo_style custom step temp press ke pe etotal vol density
run 50000
unfix 26


fix             27 all npt temp 1100 1100 1 y 0 0 1 z 0 0 1 drag 1
variable srate equal 1.0e10
variable srate1 equal "v_srate/1.0e12"
fix             28 all deform 1 x erate ${srate1} units box remap x

# Output strain and stress info to file
# for units metal, pressure is in [bars] = 100 [kPa] = 1/10000 [GPa]
# p2, p3, p4 are in MPa
variable strain equal "(lx - v_L0)/v_L0"
variable p1 equal "v_strain"
variable p2 equal "-pxx/10000"
variable p3 equal "-pyy/10000"
variable p4 equal "-pzz/10000"
fix def1 all print 1000 "${p1} ${p2} ${p3} ${p4}" file load.$a screen yes



run            20000

unfix 28
unfix 27

next a
jump input_W.txt loop
######################################
# SIMULATION DONE
print "All done"