# -*-shell-script-*-
# KMC YSZ 3D simulation

#*******************************************
# Definition of procedures
#*******************************************

proc initmd { n a b } {
#MD++ setnolog
MD++ setoverwrite
#MD++ dirname = "runs/YSZ_1V/Cube$a-$b-$n"
MD++ dirname = "runs/YSZ_Rcut_test/Cube$a-$b-$n"
}
#end of proc initmd

proc createconfig { } { MD++ {
# Number of Neighbors and Atoms in a Cell
NNM = 500
# Create Perfect Lattice Configuration
crystalstructure = CaF2
latticeconst = 5.14 #in Angstrom
latticesize = [ 1  0  0  3     #(x)
                0  1  0  3     #(y)
                0  0  1  3  ]  #(z)
makecrystal
writeall = 1 finalcnfile = zro2.cn  writecn
} }
#end of proc createconfig

proc setelements { } { MD++ {
nspecies = 4
element0 = "Zr" element1 = "O" element3 = "Y"
atommass = [  91.22400   15.99940  0.0   88.90585]
atomcharge = [ 4.0  -2.0  0.0  3.0 ]
#element0 = "Ce" element1 = "O" element3 = "Gd" 
#atommass = [ 140.11600   15.99940  0.0  157.25000] 
#atomcharge = [ 4.0  -2.0  0.0  3.0 ]
} }
#end of proc setelements

proc dopingY { } { MD++ {
input = [ 2 158 169 ]
fixatoms_by_ID
input = [3]
setfixedatomsspecies
freeallatoms
} }
#end of proc dopingY

proc createV { } { MD++ {
#input = [ 1 162 ]
#input = [ 1 165 ]
#input = [ 1 124 ]
#input = [ 1 129 ]
#input = [ 1 237 ]
input = [ 1 244 ]
fixatoms_by_ID
removefixedatoms
freeallatoms
} }
#end of proc createV

proc setR_cutoff { } { MD++ {
SKIN  = 0.75
RLIST = 9.75 #in Angstrom
} }
#end of proc setR_cutoff

proc runEwald { Rcutinewald } {
MD++ Ewald_Rcut="$Rcutinewald"
MD++ {
Ewald_CE_or_PME = 0
Ewald_option_Alpha = 1
Ewald_volume_change = 1
Ewald_precision = 4.05
Ewald_Alpha = 0.45
Ewald_time_ratio = 2.5
PME_K1 =  75
PME_K2 =  75
PME_K3 =  75
PME_bsp_n = 10
Ewald_init
} }
#end of proc runEwald

proc CG_relax { } { MD++ {
conj_ftol = 1.0e-8
conj_itmax = 30000
conj_fevalmax = 30000
conj_fixbox = 1
conj_fixboxvec = [ 1 1 1
                   1 1 1
                   1 1 1 ]
relax eval
writeall = 1 finalcnfile = relaxed.cn writecn
} }
#end of proc CG_relax

proc runMDnow { } { MD++ {
equilsteps = 0
timestep = 0.00005
DOUBLE_T = 0
saveprop = 1
savepropfreq = 10
openpropfile 
savecn = 1
savecnfreq = 20
openintercnfile
#savecfg = 1 savecfgfreq = 50 intercfgfile = ysztrial.cfg
vt2 = 4e28  #1e28 2e28 5e28
wallmass = 2.50e2     # atommass * NP
boxdamp = 7.50e-2     # damping 
saveH # Use current H as reference (H0), needed for specifying stress
fixboxvec  = [ 1 1 1 
               1 1 1
               1 1 1 ]
stress = [    0     0     0
              0     0     0
              0     0     0 ]
output_fmt = "curstep EPOT KATOM Tinst HELM HELMP TSTRESS_xx TSTRESS_yy TSTRESS_zz H_11 H_22 H_33" 
T_OBJ = 1800
initvelocity
ensemble_type = "NPT"
integrator_type = "Gear6"
totalsteps = 400
run eval
writeall = 1 finalcnfile = md_final.cn writecn
} }
#end of proc runMDnow


proc openwindow { } { MD++ {
# Plot Configuration
atomradius = [ 0.5 0.3 0 0.6] bondradius = 0.1 bondlength = 2.5
highlightcolor = red  bondcolor = red  backgroundcolor = white
atomcolor0 = orange atomcolor1 = red
atomcolor2 = cyan   atomcolor3 = blue
fixatomcolor = yellow
#-------------------------------------------------------------
plot_map_pbc = 1
plot_atom_info = 1
plotfreq = 10
rotateangles = [ 0 0 0 1.2 ] 
#openwin  alloccolors rotate saverot eval plot
#openwin alloccolors rotate saverot plot
eval
} }
#end of proc openwindow

proc getfilecounter { path precurs suffix fmt } {
set file_id 0
set yesno 1
while {$yesno != 0} {
    set file_id [expr {$file_id + 1}]
    set filename "$path/$precurs[format "$fmt" \ $file_id].$suffix"
    set yesno [file exists $filename]
#    puts "filenameis $filename and yesno is $yesno"
#    puts $file_id
}
return $file_id
}
#end of proc getfilecounter

proc setrandseed { } { MD++ {
srand48bytime
} }
#end of proc setrandseed

proc exitmd { } { MD++ quit }
#end of proc exitmd

proc sleepmd { } { MD++ sleep }
#end of proc sleepmd


#*******************************************
# Main program starts here
#*******************************************
# Assign directory
if { $argc == 0 } {
 set status 0
} elseif { $argc > 0 } {
 set status [lindex $argv 0]
 set Rcutinewald [lindex $argv 1]
}

set dbid [expr $status]
#set dbid [expr $status-20] 
#set dbid 1
set sc 3
set dop 0.93

initmd $status $sc $dop

# Set randseed by time
setrandseed    

# Create atomic configuration
createconfig
setelements
puts "createconfig done"

# Make YSZ - doping Y ions and create vacancy
dopingY
puts "dopingY done"   
createV
puts "createV done"

# Cut-off Radius   
setR_cutoff
puts "setR_cutoff done" 

# Set Ewald parameters and run it
runEwald $Rcutinewald
puts "runEwald done"

# Activate the graphics - based on Atomeye
openwindow
puts "openwindow started"

######### Execute dynamics ############

# Conjugate-Gradient Relaxation
MD++ writeall = 1 finalcnfile = initial.cn  writecn
#CG_relax
#puts "CG_relax done"
# Run MD   
#runMDnow       

# Close out MD++
#sleepmd
exitmd