# Find minimum energy fcc configuration
# Mark Tschopp, 2010
# This requires the variable latconst to be input via the command line
# e.g., lmp_win_no-mpi -var latconst 4 < calc_fcc_ver1.in

# ---------- Initialize Simulation --------------------- 
clear 
units metal 
dimension 3 
boundary p p p 
atom_style atomic 
atom_modify map array
timestep 0.0005




variable latconst equal 3.2506

# ---------- Create Atoms ---------------------
region simbox block -30 30 -20 85 -50 80
#Creating tool


lattice 	custom ${latconst}  a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.5



#region tool union 3 red blue yellow
region 1 prism -2 2 21 26 -15 -7 0 0 -3
region 2 prism -2 2 24 80 -9 -6 0 0 -10
region 3 intersect 2 1 2
region tool union 2 1 2



create_box 2 simbox

create_atoms 2 region tool basis 1 2 basis 2 1


lattice 	custom ${latconst}  a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.5

#region workpiece block -20 0 10 20 0 20

variable theta equal 4*dt*step


region boundary cylinder y 0 0 6 -5 -3 rotate v_theta 0 0 0 0 1 0
region thermw cylinder y 0 0 6 -3 -1 rotate v_theta 0 0 0 0 1 0
region newton cylinder y 0 0 6 -1 20 rotate v_theta 0 0 0 0 1 0

region workpiece union 3 boundary thermw newton



create_atoms 2 region workpiece basis 1 2 basis 2 1
delete_atoms region 3
delete_atoms region 2

group tool region tool

group boundary region boundary
group thermw region thermw
group newton region newton

group workpiece union boundary thermw newton

group moving union newton thermw







# ---------- Define Interatomic Potential --------------------- 
pair_style meam 
pair_coeff * * library.meam Zr Cu ZrCu.meam Zr Cu

#pair_style eam/fs 
#pair_coeff * * Cu-Zr_2.eam.fs Zr Cu

#minimize 1.0-e4  1.0e-4 1000 100000



neighbor 2.0 bin 
neigh_modify delay 10 check yes 
 
# ---------- Define Settings --------------------- 
#compute eng all pe/atom 
#compute eatoms all reduce sum c_eng 

dump		1 all atom 100 dump.structure


# ---------- Run Minimization --------------------- 
#reset_timestep 0 
#fix 1 all box/relax iso 0.0 vmax 0.001
#thermo 10 
#thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms 
#min_style cg 

# initialize temperature, except for fixed atoms
velocity moving create 293 58142275 units box







fix thermw thermw langevin 293.0 293.0 0.05 832743
fix tempw2 newton  langevin 293.0 293.0 0.05 928349

run 15000


variable natoms equal "count(all)" 
variable teng equal "c_eatoms"
variable length equal "lx"
variable ecoh equal "v_teng/v_natoms"



unfix thermw
unfix tempw2

variable A equal 2
variable T equal 0.1
variable omega equal 2.0*PI/$T
variable x equal swiggle(100*dt*step,$A,$T)
variable v equal v_omega*($A-cwiggle(0.0,$A,$T))

variable vmax equal 500
variable sin equal abs(sin(0.5*v_omega*dt*step))
variable velocity equal (v_vmax*v_sin)

fix md moving nve
fix 1 tool move variable NULL NULL NULL NULL NULL v_velocity

#fix 3 boundary move rotate 0 0 0 0 1 0 ${T}







fix 4 tool setforce 0 0 0



compute ke all ke/atom

variable fx equal fcm(tool,x)
variable fy equal fcm(tool,y)
variable fz equal fcm(tool,z)


thermo 1000
thermo_style custom step temp pe ke v_fx v_fy v_fz



dump  2  all custom 10 v_ht.*.lammpstrj id type xu yu zu vx vy vz c_ke 
run 5000
#print "Total energy (eV) = ${teng};"
#print "Number of atoms = ${natoms};"
print "Lattice constant (Angstoms) = ${length};"
#print "Cohesive energy (eV) = ${ecoh};"



print "All done!"