###################################### Perfect #####################################

clear
log perfect.lammps
units		metal
dimension	3
boundary	p p p
atom_style	atomic
variable lp equal 3.52395413816153

variable factor equal 1


lattice fcc ${lp} orient x -1 1 0 orient y -1 -1 2 orient z 1 1 1

variable xdim equal ${lp}*sqrt(2)
variable ydim equal ${lp}*sqrt(6)
variable zdim equal (${lp}*sqrt(3)*${factor})


variable atoms_end equal ${lp}*sqrt(3)*36

variable length equal (${zdim}*4)  


region full_box block 0 ${xdim} 0 ${ydim} 0 ${zdim} units box
region atoms block 0 ${xdim} 0 ${ydim} 0 ${atoms_end} units box



create_box 1 full_box
create_atoms 1 region full_box


#replicate 6 3 36 bbox

variable x_length equal "lx"
variable y_length equal "ly"
variable z_length equal "lz"



#region top block 0 ${x_length} 0 ${y_length} ${length} ${z_length} units box
#region bottom block 0 ${x_length} 0 ${y_length} 0 ${length} units box

#group top_atoms region top
#group bottom_atoms region bottom

#region vacuum block 0 ${x_length} 0 ${y_length} ${atoms_end} ${z_length} units box
#delete_atoms region vacuum
# --------------------- EAM POTENTIAL ---------------------

pair_style	eam/alloy
 pair_coeff * * "/home/labuser/snap/snapd-desktop-integration/178/potential_files/FeNiCr_Foster_Sills.setfl" Ni 
neighbor 2.0 bin 
# --------------------- SETTINGS ---------------------

compute peratom all pe/atom
compute total_energy all reduce sum c_peratom


####################### thermal properties output on screen ###################### 

thermo 10		        
thermo_style custom step pe lx ly lz press pxx pyy pzz c_total_energy   


#dump dump1 all cfg 1 dump_perfect/perfect_*.cfg mass type xs ys zs c_peratom 

# Dump for top atoms
#dump dump2 top_atoms cfg 1 dump_perfect/perfect_top_*.cfg mass type xs ys zs c_peratom 

# Dump for bottom atoms
#dump dump3 bottom_atoms cfg 1 dump_perfect/perfect_bottom_*.cfg mass type xs ys zs c_peratom


dump dump4 all custom 1 perfect.dat id type xs ys zs c_peratom 
min_style cg
minimize 1e-10 1e-10 5000 10000  
variable tengperfect equal "c_total_energy"
variable natomsperfect equal "count(all)"
variable ecohperfect equal "v_tengperfect/v_natomsperfect"

variable ecohperfect_final equal "v_ecohperfect"

variable areaperfect equal "lx*ly"

print "cohesive energy of perfect ${ecohperfect}"
print "total energy of perfect ${tengperfect}"
print "perfect area equals ${areaperfect}"
print "No. of atoms in perfect ${natomsperfect}"

write_data pure_data.txt nocoeff nofix


###################################### Fault #####################################

clear
log fault.lammps
units		metal
dimension	3
boundary	p p p
atom_style	atomic
variable lp equal 3.52395413816153

variable factor equal 1


read_data pure_data.txt

variable x_length equal "lx"
variable y_length equal "ly"
variable z_length equal "lz"


#Shear 1
variable length_1 equal (${zdim}*24)
region top_1 block 0 ${x_length} 0 ${y_length} ${length_1} ${z_length} units box
region bottom_1 block 0 ${x_length} 0 ${y_length} 0 ${length_1} units box


group top_atoms_1 region top_1
group bottom_atoms_1 region bottom_1

variable y_displace equal -1*(${lp}*sqrt(6))/6
displace_atoms top_atoms_1 move 0.0 ${y_displace} 0.0 units box
fix 1 all setforce 0 0 NULL

#Shear 2
variable length_2 equal (${zdim}*12)
region top_2 block 0 ${x_length} 0 ${y_length} ${length_2} ${z_length} units box
region bottom_2 block 0 ${x_length} 0 ${y_length} 0 ${length_2} units box

group top_atoms_2 region top_2
group bottom_atoms_2 region bottom_2

variable y_displace equal -1*(${lp}*sqrt(6))/6
displace_atoms top_atoms_2 move 0.0 ${y_displace} 0.0 units box
fix 2 all setforce 0 0 NULL

#region vacuum block 0 ${x_length} 0 ${y_length} ${atoms_end} ${z_length} units box
#delete_atoms region vacuum




# --------------------- EAM POTENTIAL ---------------------

pair_style	eam/alloy
pair_coeff * * "/home/labuser/snap/snapd-desktop-integration/178/potential_files/FeNiCr_Foster_Sills.setfl" Ni 
neighbor 2.0 bin 
# --------------------- SETTINGS ---------------------

compute peratom all pe/atom
compute total_energy all reduce sum c_peratom


####################### thermal properties output on screen ###################### 

thermo 10		        #'thermo' specifies the output during minimization at every 10 steps
thermo_style custom step pe lx ly lz press pxx pyy pzz c_total_energy   #'thermo_style' specifies what type of output is shown to screen eg: pe : potential energy;lx,ly,lz : simulation cell dimensions; press : pressure etc..


dump dump1 all cfg 10 dump_fault/fault_*.cfg mass type xs ys zs c_peratom 

dump dump4 all custom 1 fault.dat id type xs ys zs c_peratom 


# Dump for top atoms
dump dump2 top_atoms_1 cfg 10 dump_fault/fault_top_*.cfg mass type xs ys zs c_peratom 

# Dump for bottom atoms
dump dump3 bottom_atoms_1 cfg 10 dump_fault/fault_bottom_*.cfg mass type xs ys zs c_peratom
min_style cg
minimize 1e-10 1e-10 5000 10000  
variable tengfault equal "c_total_energy"
variable natomsfault equal "count(all)"
variable ecohfault equal "v_tengfault/v_natomsfault"

variable ecohfault_final equal "v_ecohfault"

variable areafault equal "lx*ly"


print "cohesive energy of fault ${ecohfault}"
print "total energy of fault ${tengfault}"
print "Fault area equals ${areafault}"
print "No. of atoms in fault ${natomsfault}"