# -*-shell-script-*-
# test the MEAM potential of Au (TCL)
# run by meam-baskes
#------------------------------------------------------------

source "~/Codes/MD++/scripts/Examples/Tcl/startup.tcl"

#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } { MD++ {
setnolog
setoverwrite
dirname =~/Codes/MD++/runs/si-au
zipfiles = 0 writeall = 1
} }
#end of proc initmd

#------------------------------------------------------------
proc readmeam-baskes { } { MD++ {
#Read in MEAM potential (Baskes format)
meamfile = "~/Codes/MD++/Fortran/MEAM-Baskes/meamf"
meafile  = "~/Codes/MD++/Fortran/MEAM-Baskes/meafile_Au_orig"
ntypes = 1  ename0 = "Au1  " #leave spaces to be compatible with fortran
rcut = 6.0  kode0 = "library"
readMEAM
NNM = 300
} }

#------------------------------------------------------------
proc readmeam-lammps { } { MD++ {
#Read in MEAM potential (Baskes format)
#meamfile = "~/Codes/MD++/Fortran/MEAM-Lammps/meamf"
#nspecies = 1  element0 = "Au"  rcut = 6.0  readMEAM
#NNM = 300
meamfile = "~/Codes/MD++/potentials/MEAMDATA/meamf"
meafile  = "~/Codes/MD++/potentials/MEAMDATA/AuSi2nn.meam"
#nspecies = 2 element0="Si4" element1="AuBt" rcut = 6.0 readMEAM
nspecies = 2 element0="AuBt" element1="Si4" rcut = 6.0 readMEAM
NNM=400
} }

#-------------------------------------------------------------
proc readmeam { } { MD++ {
#Read in MEAM potential (MD++ format)
potfile = "~/Codes/MD++/potentials/MEAMDATA/meamdata.Au2" readMEAM
#potfile = "~/Codes/MD++/potentials/MEAMDATA/meamdata.Siz" readMEAM
#potfile = "~/Codes/MD++/potentials/MEAMDATA/meamdata.Si" readMEAM
#xzbl = -10 xzbl0 = -9 #turn off correction
} }

#-------------------------------------------------------------
proc readeam { } { MD++ {
#Read in EAM potential (MD++ format)
potfile = "~/Codes/MD++/potentials/EAMDATA/eamdata.AuFoiles" eamgrid = 500 readeam
NNM = 300
} }

#--------------------------------------------
proc create_FCC_Au { } { MD++ {
# Create Perfect Lattice Configuration
crystalstructure = face-centered-cubic  latticeconst = 4.07 #(A) FCC  structure
element0 = "Au"   element1 = "Si"  
atommass = [ 196.96655 28.0855 ] # Au, Si (g/mol)
latticesize = [ 1 0 0  4
                0 1 0  4
                0 0 1  4 ]
makecrystal finalcnfile = fcc-au-4x4x4.cn writecn
#makecrystal finalcnfile = fcc-au-3x3x3.cn writecn
#incnfile = fcc-au-3x3x3.cn readcn
} }        

#--------------------------------------------
proc create_BCC_Au { } { MD++ {
# Create Perfect Lattice Configuration
crystalstructure = body-centered-cubic  latticeconst = 3.0 #(A) BCC  structure
element0 = "Au"   element1 = "Si"  
atommass = [ 196.96655 28.0855 ] # Au, Si (g/mol)
latticesize = [ 1 0 0  5
                0 1 0  5
                0 0 1  5 ]
makecrystal finalcnfile = bcc-au-4x4x4.cn writecn
#incnfile = bcc-au-4x4x4.cn readcn
} }        

#--------------------------------------------
proc create_HCP_Au { } { MD++ {
# Create Perfect Lattice Configuration
crystalstructure = hexagonal-ortho  latticeconst = [ 2.87792 2.87792 4.6996 ] #(A) HCP  structure
element0 = "Au"   element1 = "Si"  
atommass = [ 196.96655 28.0855 ] # Au, Si (g/mol)
latticesize = [ 1 0 0  5
                0 1 0  5
                0 0 1  5 ]
makecrystal finalcnfile = hcp-au-5x5x5.cn writecn
#incnfile = hcp-au-5x5x5.cn readcn
} }        

#--------------------------------------------
proc create_DC_Si { } { MD++ {
# Create Perfect Lattice Configuration
crystalstructure = diamond-cubic  latticeconst = 5.4309529817532409 #(A) for Si
element0 = "Au"   element1 = "Si"  
atommass = [ 196.96655 28.0855 ] # Au, Si (g/mol)
latticesize = [ 1 0 0  3
                0 1 0  3
                0 0 1  3 ]
makecrystal finalcnfile = dc-si-3x3x3.cn 
fixallatoms input = 1 setfixedatomsspecies freeallatoms
writeall = 1 writecn
#incnfile = dc-si-2x2x2.cn readcn
} }        


#--------------------------------------------
proc create_B1_AuSi { } { MD++ {
# Create Perfect Lattice Configuration
crystalstructure = NaCl  latticeconst = 5.400  #(A) for AuSi (B1)
element0 = "Au"   element1 = "Si"
atommass = [ 196.96655 28.0855 ] # Au, Si (g/mol)
latticesize = [ 1 0 0  4
                0 1 0  4
                0 0 1  4 ]
makecrystal finalcnfile = b1-ausi-4x4x4.cn
#fixallatoms input = 1 setfixedatomsspecies freeallatoms
writeall = 1 writecn
} }


#-------------------------------------------------------------
proc openwindow { } { MD++ {
# Plot Configuration
#
atomradius = [ 0.67 0.34 ] bondradius = 0.3 bondlength = 1.0
atomcolor = orange bondcolor = red backgroundcolor = gray70 fixatomcolor = yellow
plotfreq = 100  rotateangles = [ 0 0 0 1.25 ]
color00 = "orange"  color01 = "red"    color02 = "green"
color03 = "magenta" color04 = "cyan"   color05 = "purple"
color06 = "gray80"  color07 = "white"  color08 = "blue"
#
#plot_color_axis = 2  #2: use CSD (default 0: use local energy)
#
plot_color_windows = [ 0
                      1.5  3  8 #color08 = blue
                       3  20  6 #color06 = gray80
		       20 40  7
                       40 100 1
                      ]
#plot_limits = [ 1 -10 10 -0.05 10 -10 10 ]
plot_atom_info = 1  plotfreq = 10 
#
openwin  alloccolors rotate saverot eval plot
} }
#end of proc openwindow

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

#--------------------------------------------
proc relax_freebox { } { MD++ {
# Conjugate-Gradient relaxation
conj_ftol = 2e-6 conj_itmax = 1000 conj_fevalmax = 10000
conj_fixbox = 0  conj_fixboxvec = [ 0 1 1
                                    1 0 1
                                    1 1 0 ]
relax
} }
#end of proc relax_freebox
#--------------------------------------------
proc relax_fixbox { } { MD++ {
# Conjugate-Gradient relaxation
conj_ftol = 2e-6 conj_itmax = 1000 conj_fevalmax = 10000
conj_fixbox = 1  conj_fixboxvec = [ 0 1 1
                                    1 0 1
                                    1 1 0 ]
relax
} }
#end of proc relax_freebox

#--------------------------------------------
proc setup_md { } { MD++ {     
equilsteps = 0  timestep = 0.001 # (ps)
atommass = [ 196.96655 28.0855 ] # Au, Si (g/mol)
DOUBLE_T = 0
saveprop = 1 savepropfreq = 100 openpropfile #run
savecn = 1 savecnfreq = 1000 openintercnfile
plotfreq = 100
vt2 = 2e28  #1e28 2e28 5e28
wallmass = 2e3     # atommass * NP = 14380
boxdamp = 1e-3     # optimal damping for 216 atoms and wallmass 1e-3
saveH # Use current H as reference (H0), needed for specifying stress
fixboxvec  = [ 0 1 1 
               1 0 1
               1 1 0 ]
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" 
} }
#end of proc setup_md



#---------------------------------------------
# Compute elastic constants for cubic crystal
proc compute_elast_const { datafile structname argv0 } {    
relax_freebox
MD++ finalcnfile = relaxed.cn writecn

set nx [MD++_Get latticesize  3]
set ny [MD++_Get latticesize  7]
set nz [MD++_Get latticesize 11]

set Lx [MD++_Get H 0]
set Ly [MD++_Get H 4]
set Lz [MD++_Get H 8]

set NP [MD++_Get NP]

set a0 [expr $Lx/$nx]

set nx0 [expr $Lx/$a0]
set ny0 [expr $Ly/$a0]
set nz0 [expr $Lz/$a0]

#-----------------------
# Compute bulk modulus
#-----------------------
MD++ saveH conj_fixbox = 1

set fileID [open $datafile w ]

puts $fileID "%This file is generated by MD++ ($argv0)\n\n"
puts $fileID "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
puts $fileID "% $structname "
puts $fileID "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
puts $fileID "$structname = struct();"
puts $fileID "$structname.N = $NP;"
puts $fileID "$structname.a0 = $a0;"
puts $fileID "$structname.C1 = \[ $nx0 0 0 \];"
puts $fileID "$structname.C2 = \[ 0 $ny0 0 \];"
puts $fileID "$structname.C3 = \[ 0 0 $nz0 \];"

puts $fileID "$structname.data = \[ %lattice_const.(A)  Epot(eV)"

for {set x 0.997} {$x<=1.003} {set x [expr $x+0.001]} {

    MD++ restoreH input = $x scaleH relax eval

    puts $fileID "[format %.3f%3d%24.13e $x 1 [MD++_Get EPOT]]"
}

puts $fileID "];"
puts $fileID "$structname.data(:,1) = $structname.data(:,1)*$structname.a0;\n%"
puts $fileID "$structname = fit_aEB($structname);\n%"
puts $fileID "%disp(sprintf('$structname: a0 = %.6f A, Ecoh = %.6f eV, B = %.6f GPa', ..."
puts $fileID "%             $structname.a0, $structname.Ecoh, $structname.B));\n"

#-----------------------
# Compute C11 C12
#-----------------------
puts $fileID "$structname.data = \[ %lattice_const.(A)  Epot(eV)"

for {set x -0.003} {$x<=0.003} {set x [expr $x+0.001]} {

    MD++ restoreH input = \[ 1 1 $x \] shiftbox relax eval

    puts $fileID "[format %.3f%3d%24.13e [expr $x+1] 1 [MD++_Get EPOT]]"
}

puts $fileID "];"
puts $fileID "$structname = fit_C11($structname);\n%"
puts $fileID "%disp(sprintf('$structname: B = %.6f GPa, C11 = %.6f GPa, C12 = %.6f GPa', ..."
puts $fileID "%             $structname.B, $structname.C11, $structname.C12));\n"



#-----------------------
# Compute C44
#-----------------------
# rotate crystal by 45 degrees around z axis
MD++  {
    restoreH relax
    finalcnfile = perf.cn writeall = 1 writecn
    incnfile = perf.cn input = [ 1 0 ] splicecn
    input = [ 2 1 0.5 ] redefinepbc
    input = [ 1 2 -1  ] redefinepbc  reorientH 
    relax eval saveH
    finalcnfile = perf2.cn writecn
    incnfile = perf2.cn readcn
}

puts $fileID "$structname.N = [expr $NP*2];"
puts $fileID "$structname.C1 = \[  $nx0 -$ny0 0 \];"
puts $fileID "$structname.C2 = \[  $nx0 $ny0 0 \];"
puts $fileID "$structname.C3 = \[  0 0 $nz0 \];"

puts $fileID "$structname.data = \[ %lattice_const.(A)  Epot(eV)"

MD++ fixallatoms  #Note: only for SC and BCC crystals!

for {set x -0.003} {$x<=0.003} {set x [expr $x+0.001]} {

    MD++ restoreH input = \[ 1 1 $x \] shiftbox relax eval

    puts $fileID "[format %.3f%3d%24.13e [expr $x+1] 1 [MD++_Get EPOT]]"
}

puts $fileID "];"
puts $fileID "$structname = fit_C44($structname);\n%"
puts $fileID "disp(sprintf('$structname: a0 = %.6f A, Ecoh = %.6f eV, B = %.6f GPa', ..."
puts $fileID "             $structname.a0, $structname.Ecoh, $structname.B));"
puts $fileID "disp(sprintf('        C11 = %.6f GPa, C12 = %.6f GPa, C44 = %.6f GPa', ..."
puts $fileID "             $structname.C11, $structname.C12, $structname.C44));\n"


close $fileID
}     



#*******************************************
# Main program starts here
#*******************************************
# read in status from command line argument
if { $argc == 0 } {
 set status 0
} elseif { $argc > 0 } {
 set status [lindex $argv 0]
}
puts "status = $status"

set myname  [ MD++_Get "myname"]
puts "myname = $myname"

initmd

if { [ string match "*baskes*" $myname ] } {
  puts "readmeam-baskes"
  readmeam-baskes
} elseif { [ string match "*lammps*" $myname ] } {
  puts "readmeam-lammps"
  readmeam-lammps
} elseif { [ string match "*meam*" $myname ] } {
  puts "readmeam"
  readmeam
} elseif { [ string match "*eam*" $myname ] } {
  puts "readeam"
  readeam
}

if { $status == 0 } {
  #computing Ecoh of diamond-cubic Si
  create_DC_Si
  #openwindow
  MD++ eval

  set N  [MD++_Get "NP"]
  set E1 [MD++_Get "EPOT"]

  puts "N = $N  E1 = $E1"

  MD++ fixallatoms  
  relax_freebox  

  MD++ eval
  set E0 [MD++_Get "EPOT"]
  set Ecoh [expr $E0/$N]
  set Lx [MD++_Get H_11]
  set Nx [MD++_Get latticesize(3)]
  set a0 [expr $Lx/$Nx]

  puts "Lx = $Lx Nx = $Nx a0 = $a0"
  puts "E0 = $E0  Ecoh = $Ecoh"

  #MD++ sleep

} elseif { $status == 1 } {
  #computing Ecoh of FCC Au
  create_FCC_Au
  #openwindow
  MD++ eval

  set N  [MD++_Get "NP"]
  set E1 [MD++_Get "EPOT"]

  puts "N = $N  E1 = $E1"

  MD++ fixallatoms  
  relax_freebox  

  MD++ eval
  set E0 [MD++_Get "EPOT"]
  set Ecoh [expr $E0/$N]
  set Lx [MD++_Get H_11]
  set Nx [MD++_Get latticesize(3)]
  set a0 [expr $Lx/$Nx]

  puts "Lx = $Lx Nx = $Nx a0 = $a0"
  puts "E0 = $E0  Ecoh = $Ecoh"

  #MD++ sleep
  

} elseif { $status == 2 } {
  #computing Ecoh of B1 AuSi
  create_B1_AuSi
  #openwindow
  MD++ eval

  set N  [MD++_Get "NP"]
  set E1 [MD++_Get "EPOT"]

  puts "N = $N  E1 = $E1"

  MD++ fixallatoms
  relax_freebox

  MD++ eval
  set E0 [MD++_Get "EPOT"]
  set Ecoh [expr $E0/$N]
  set Lx [MD++_Get H_11]
  set Nx [MD++_Get latticesize(3)]
  set a0 [expr $Lx/$Nx]

  puts "Lx = $Lx Nx = $Nx a0 = $a0"
  puts "E0 = $E0  Ecoh = $Ecoh"

  #MD++ sleep

} elseif { $status == 3 } {
  #computing energy of Si impurity in Au
  create_FCC_Au
  #openwindow
  MD++ eval

  set N  [MD++_Get "NP"]
  set E1 [MD++_Get "EPOT"]

  puts "N = $N  E1 = $E1"

  MD++ fixallatoms  
  relax_freebox  

  MD++ eval
  set E0 [MD++_Get "EPOT"]
  set Ecoh [expr $E0/$N]
  set Lx [MD++_Get H_11]
  set Nx [MD++_Get latticesize(3)]
  set a0 [expr $Lx/$Nx]

  puts "Lx = $Lx Nx = $Nx a0 = $a0"
  puts "E0 = $E0  Ecoh = $Ecoh"

  # turn a Au atom into Si
  MD++ species(0) = 1

  relax_freebox
  MD++ eval

  set E1 [MD++_Get "EPOT"]
  set Ecoh_si -4.63
  set Eimp [expr $E1 - ($N-1.0)/$N*$E0 - $Ecoh_si]
  puts "Eimp = $Eimp  (eV) Si atom in Au fcc" 
  #MD++ sleep
  
} elseif { $status == 4 } {
  #computing energy of Au impurity in Si
  create_DC_Si
  #openwindow
  MD++ eval

  set N  [MD++_Get "NP"]
  set E1 [MD++_Get "EPOT"]

  puts "N = $N  E1 = $E1"

  MD++ fixallatoms
  relax_freebox

  MD++ eval
  set E0 [MD++_Get "EPOT"]
  set Ecoh [expr $E0/$N]
  set Lx [MD++_Get H_11]
  set Nx [MD++_Get latticesize(3)]
  set a0 [expr $Lx/$Nx]

  puts "Lx = $Lx Nx = $Nx a0 = $a0"
  puts "E0 = $E0  Ecoh = $Ecoh"

  # turn a Si atom into Au
  MD++ species(0) = 0

  relax_freebox
  MD++ eval

  set E1 [MD++_Get "EPOT"]
  set Ecoh_au -3.93
  set Eimp [expr $E1 - ($N-1.0)/$N*$E0 - $Ecoh_au]
  puts "Eimp = $Eimp  (eV) Au atom in Si DC"
  #MD++ sleep

#####

} elseif { $status == 10 } {
  #computing unrelaxed vacancy energy
  create_FCC_Au
  MD++ eval
  set E0 [MD++_Get "EPOT"]
  set N  [MD++_Get "NP"]
  set Ecoh [expr $E0/$N]

  MD++ input = \[ 1 0 \] fixatoms_by_ID removefixedatoms
  #openwindow
  MD++ eval
  set N1 [MD++_Get "NP"]
  set E1 [MD++_Get "EPOT"]
  set Ev [expr $E1 - $E0 + $Ecoh]

  puts "N = $N  E0 = $E0 E1 = $E1  Ecoh = $Ecoh  Ev = $Ev"

  relax_fixbox  
  set E2 [MD++_Get "EPOT"]
  set Ev_relaxed [expr $E2 - $E0 + $Ecoh]
  puts "Ev_relaxed = $Ev_relaxed"

  #MD++ sleep

} elseif { $status == 20 } {
  #BCC structure
  create_BCC_Au
  #openwindow
  MD++ eval

  set N  [MD++_Get "NP"]
  set E1 [MD++_Get "EPOT"]

  puts "N = $N  E1 = $E1"

  MD++ fixallatoms  
  relax_freebox  

  MD++ eval
  set E0 [MD++_Get "EPOT"]
  set Ecoh [expr $E0/$N]
  set Lx [MD++_Get H_11]
  set Nx [MD++_Get latticesize(3)]
  set a0 [expr $Lx/$Nx]

  puts "Lx = $Lx Nx = $Nx a0 = $a0"
  puts "E0 = $E0  Ecoh = $Ecoh"


  #MD++ sleep

} elseif { $status == 30 } {
  #HCP structure
  create_HCP_Au
  #openwindow
  #MD++ sleep
  MD++ eval

  set N  [MD++_Get "NP"]
  set E1 [MD++_Get "EPOT"]

  puts "N = $N  E1 = $E1"
  
  MD++ fixallatoms  
  relax_freebox  

  MD++ eval
  set E0 [MD++_Get "EPOT"]
  set Ecoh [expr $E0/$N]
  set Lx [MD++_Get H_11]
  set Nx [MD++_Get latticesize(3)]
  set a0 [expr $Lx/$Nx]
  set Lz [MD++_Get H_33]
  set Nz [MD++_Get latticesize(11)]
  set c0 [expr $Lz/$Nz]
  set covera [expr $c0/$a0]
  puts "Lx = $Lx Nx = $Nx a0 = $a0"
  puts "Lz = $Lz Nz = $Nz c0 = $c0"
  puts "E0 = $E0  Ecoh = $Ecoh  covera = $covera"


  #MD++ sleep

}  elseif { $status == view } {
  #load structure to view

 set inputdir ../au-test-s0-md-300K

 MD++ incnfile = "relaxed-2x2.cn" readcn
 openwindow
  
 MD++ plot_color_axis = 2 eval plot
 
 MD++ sleep
}

exitmd