# -*-shell-script-*-
# Analyze the LAMMPS simulation dump files for Au-Si
# and compute the number of crystalline Si atoms

source "scripts/Examples/Tcl/startup.tcl"

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

#-------------------------------------------------------------
proc initmd { status } { 
  #MD++ setnolog
  MD++ setoverwrite
  MD++	dirname = runs/si_au_test_status-$status
  MD++ zipfiles = 0 writeall = 1
}
#end of proc initmd

#-------------------------------------------------------------
proc readeam {} { MD++ {
potfile=~/Codes/MD++/potentials/EAMDATA/eamdata.AuFoiles
eamgrid = 500 readeam NNM = 3000
} }
#end of proc readeam

#------------------------------------------------------------
proc readmeam-lammps {} {
    MD++ {
        meamfile = "~/WORKSPACE/yanmingw/Codes/MD++/potentials/MEAMDATA/meamf"
        meafile  = "~/WORKSPACE/yanmingw/Codes/MD++/potentials/MEAMDATA/AuSi2nn.meam"
        nspecies=2 element0="AuBt" element1="Si4" rcut = 4.5 readMEAM
        NNM=3000
    }
}
#end of proc readmeam-lammps

#-------------------------------------------------------------
proc setup_window { } { MD++ {
# Plot Configuration
#
	atomradius=[ 0.8 0.5 ] bondradius = 0.3
	win_width= 500 win_height = 500
	bondlength=0 atomcolor = orange
	atomcolor1 = gray90 highlightcolor = purple
	bondcolor = green backgroundcolor = gray70
	plot_atom_info = 2
	plotfreq= 100   printfreq = 10

        fixatomcolor = yellow
        color00 = "orange"  color01 = "purple" color02 = "green"
        color03 = "magenta" color04 = "cyan"   color05 = "red"
        color06 = "gray80"  color07 = "white"
        #plot_color_axis = 8  #(probably should use this option)
        #plot_color_axis = 7
        plot_color_windows = [ 0
                              -10   -1.9999999  6
                              -1.999999 10   0
                              ]
	rotateangles = [ 0 0 0 1.15 ]
} }

#------------------------------------------------------------
proc openwindow { } {
setup_window
MD++ openwin alloccolors rotate saverot
MD++ refreshnnlist eval
MD++ plot
}

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

#--------------------------------------------

#---------------------------------------------------
#       Conjugate-Gradient relaxation
#---------------------------------------------------
proc relax_freebox { } { MD++ {
conj_ftol = 5e3 conj_itmax = 500 conj_fevalmax = 1000
conj_fixbox = 0
conj_fixboxvec = [ 0 1 1   #box remain rectangular
                   1 0 1
                   1 1 0 ]
relax
}}

proc relax_freeboxy { } { MD++ {
conj_ftol = 5e3 conj_itmax = 500 conj_fevalmax = 1000
conj_fixbox = 0
conj_fixboxvec = [ 1 1 1   #box remain rectangular
                   1 0 1
                   1 1 1 ]
relax
}}

proc relax_fixbox { } { MD++ {
conj_ftol = 5e3 conj_itmax = 500 conj_fevalmax = 1000
conj_fixbox = 1
conj_fixboxvec = [ 0 1 1   #box remain rectangular
                   1 0 1
                   1 1 0 ]
relax
}}

#--------------------------------------------
proc setup_QLM { } {
    MD++ {
	UMB_order_param = 103
        #setting for 
        #103: q3, for detecting diamond cubic structure,
        #111:HCPonly atoms,
        #112:FCConly atoms,
        #113:HCP+FCC atoms,
	Rc_for_QLM = 3.48
	QLM_cutoff = -0.7
	N_lgst_index = -2
    }
}
#end of proc setup_QLM


proc create_droplet { SiFrac orien T } {

MD++ {
allocmultiple = 3 #anticipating enlargement of system size
#--------------------------
# Create Perfect Crystal Au
crystalstructure = face-centered-cubic latticeconst = 4.08
latticesize = [  -1   0   1   14
                  0   1   0   13
                  1   0   1   14
              ]
#latticesize = [ -1  0  1  9
#                 0  1  0  13
#                 1  0  1  4
#              ]
#latticesize = [ -1  0  1  7
#                 0  1  0  6
#                 1  0  1  7
#              ]

makecrystal
input = [ 1 -10 10 -10 10 -10 10 ] fixatoms_by_position
# set all atoms to be Au
input = 0 setfixedatomsspecies freeallatoms
}

set SiFracequil [ GetEquilSiFrac $T ]

assignSiFrac $SiFracequil

MD++ nspecies = 2

MD++ srand48bytime randomposition

MD++ { input = [ 3 3 0.00165441176 ] shiftbox }

if { $orien == 111 } {
MD++ { input = [ 1 1 0.02470239 ] shiftbox }
} elseif { $orien == 110 } {
MD++ { input = [ 1 1 0.0458325 ] shiftbox }
} elseif { $orien == 112 } {
MD++ { input = [ 1 1 0.08686101 ] shiftbox }
} else {
MD++ { input = [ 1 1 0.03532135 ] shiftbox }
}

relax_freeboxy
relax_fixbox
relax_freeboxy
relax_fixbox

}
#--------------------------------------------


proc create_SiNW { SiFrac orien } {
#--------------------------
# Create Perfect Crystal Si
#
MD++ crystalstructure = diamond-cubic latticeconst = 5.431

if { $orien == 111 } {
MD++ {
latticesize = [
                 1  -2   1   6
                 1   1   1   2
                -1   0   1   10
              ]
makecrystal
}
} elseif { $orien == 110 } {
MD++ {

#latticesize = [  0  1  0 15
#                 1  0  1 5
#                -1  0  1 10 ]
#latticesize = [ 0  1  0  10
#                1  0  1  2
#               -1  0  1  3 ]
latticesize = [ 0  1  0  15
                1  0  1  2
               -1  0  1  10 ]



makecrystal
}

} elseif { $orien == 112 } {
MD++ {
latticesize = [  1  1  1 6
                 1 -2  1 1
                -1  0  1 3 ]
makecrystal
}

} else {
MD++ {
latticesize = [ 1 0 1 7
                0 1 0 3
               -1 0 1 3 ]
makecrystal
#MD++ { input = [ 2 2 0.1 ] redefinepbc }
 }
}

MD++ { input=[ 1 1 0.01 ] pbcshiftatom }


MD++ {

input = [ 1 -10 10 -10 10 -10 10 ] fixatoms_by_position
input = 1 setfixedatomsspecies freeallatoms
eval
}

MD++ Fold_into_Unitcell

}


#-------------------------------------------------------

proc GetEquilSiFrac { T } {
if { $T == 600 } {
return 22
} elseif { $T == 630 } {
return 23
} elseif { $T == 650 } {
return 24
} elseif { $T == 700 } {
return 26
} elseif { $T == 800 } {
return 28
} elseif { $T == 900 } {
return 32
} elseif { $T == 1000 } {
return 36
} elseif { $T == 1100 } {
return 40
} elseif { $T == 1200 } {
return 48
} elseif { $T == 1300 } {
return 56
} elseif { $T == 1400 } {
return 65
} elseif { $T == 1500 } {
return 75
} elseif { $T == 1600 } {
return 85
}

}


#-------------------------------------------------------


proc assignSiFrac { SiFrac } {

set NP [MD++_Get NP]
set N_Si [ expr $NP*$SiFrac/100.0]
for { set x 0 } { $x < $N_Si } { incr x 1 } {
MD++ input=\[ 1 $x \]
MD++ fixatoms_by_ID
}

#-------------------------------------------------------

proc AdjustSiFrac { dN_Si NP } {
puts "NP = ${NP}"
set NP_ [ expr $NP - 1 ]
set count 0
 if { $dN_Si == 0 } {
   return 0
 } elseif { $dN_Si > 0 } {
    for { set x 0 } { $x < $NP } { incr x 1 } {
        set species [MD++_Get species $x ]
	if { $species == 0 } {
	    set count [ expr $count + 1 ]
	    MD++_Set species($x) 1
	    if { $dN_Si == $count } {
	        break;
	    }
	}
    }
 } else {
    for { set x $NP_ } { $x >= 0 } { set x [ expr $x-1 ] } {
	set species [ MD++_Get species $x ]
	if { $species == 1 } {
	    set count [ expr $count -1 ]
            MD++_Set species($x) 0
            if { $dN_Si == $count } {
                break;
            }
	}
    }
 }
 puts "count = $count"
 MD++ initvelocity
}
MD++ input=1 setfixedatomsspecies freeallatoms
}



#--------------------------------------------
proc setup_md { } { MD++ {
   atommass= [ 196.96655 28.0855 ] equilsteps=0 timestep=0.0005
   DOUBLE_T=0
   savecn=1 savecnfreq=10000
   savecfg=1 savecfgfreq=10000
   saveprop=1 savepropfreq=100 

   openpropfile openintercnfile openintercfgfile
   ensemble_type = "NVTC" integrator_type = "Gear6" NHChainLen = 4 NHMass=[1e-1 1e-2 1e-2 1e-2 ]
   #ensemble_type = "NVT" integrator_type = "VVerlet" NHMass = 1e-1 

   boxdamp=1e-3 saveH
   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 TSTRESSinMPa_xx TSTRESSinMPa_yy TSTRESSinMPa_zz N_lgst_cluster"
   writeall=1
}}


#*******************************************
# 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"

if { $argc <= 1 } {
 set n 0
} elseif { $argc > 1 } {
 set n [lindex $argv 1]
}    
puts "n = $n"

if { $status == 0 } {
  # convert LAMMPS dump files into .cn and .cfg files 
  # and compute number of crystalline Si atoms
  MD++ setnolog
  initmd $status

  set lammps_folder $n

  set dump_freq 1000000
  set num_files 80
  #set num_files 63

  readmeam-lammps
  set fileID [open "analysis.out" w]

  setup_window
  #openwindow

  setup_QLM

  for {set x 0} {$x <= $num_files} {incr x} {
     set y [expr ($dump_freq*$x)]
     MD++ incnfile= "$lammps_folder/dump_$y.md"
     MD++ input=\[ $y	0 \]
     MD++ readLAMMPS

     # We need to manually input the box size here because
     #  LAMMPS does not print enough digits for the box in the dump file
     MD++ H_11 = 8.3543607e+01
     MD++ H_22 = 1.17028139e+02
     MD++ H_33 = 8.1664514e+01
     MD++ RHtoS


     #MD++ incnfile=config_$y.cn readcn

     #MD++ finalcnfile = config_$y.lammps writeLAMMPS
     #MD++ finalcnfile = force_$y.log writeFORCE

     MD++ allocUMBorder eval
     MD++ calcrystalorder
     puts "N_lgst_cluster =  [ MD++_Get N_lgst_cluster ]"

     MD++ writeall = 1 
     MD++ finalcnfile=config_$y.cn writecn
     MD++ finalcnfile = config_$y.cfg writeatomeyecfg
     MD++ curstep = $y
     puts $fileID "[format %12.8e [MD++_Get curstep]] [format %12.8e [MD++_Get EPOT]] [format %12.8e [MD++_Get Tinst]] [format %12.8e [MD++_Get TSTRESSinMPa_xx]] [format %12.8e [MD++_Get TSTRESSinMPa_yy]] [format %12.8e [MD++_Get TSTRESSinMPa_zz]] [format %12.8e [MD++_Get N_lgst_cluster]]"
     flush $fileID
  }
  close $fileID
  exitmd

} elseif { $status == 1 } {
  # read MD++ inter .cn files and compute number of crystalline Si atoms
  MD++ setnolog
  initmd $status

  set mdpp_folder $n

  #set num_files 63
  set num_files 326
  #set num_files 3

  readmeam-lammps
  set fileID [open "analysis.out" w]

  setup_md
  setup_QLM
  setup_window

  if { 0 } {
     MD++ incnfile = "$mdpp_folder/init.cn" readcn nspecies = 2
     MD++ allocUMBorder eval calcrystalorder
     puts "N_lgst_cluster =  [ MD++_Get N_lgst_cluster ]"
     exitmd
  }

  for {set x 1} {$x <= $num_files} {incr x} {
     set y [format "%04d" $x]
     MD++ incnfile= "$mdpp_folder/inter$y.cn"  readcn  nspecies = 2

     #MD++ finalcnfile = config_$y.cfg writeatomeyecfg
     #MD++ finalcnfile = config_$y.lammps writeLAMMPS
     #MD++ finalcnfile = force_$y.log writeFORCE

     MD++ finalcnfile = newinter$y.cn writeall = 1 writecn

     if { $x == 1 }  {
       MD++ allocUMBorder 
     }

     MD++ eval calcrystalorder
     puts "N_lgst_cluster =  [ MD++_Get N_lgst_cluster ]"

     MD++ curstep = $y
     puts $fileID "[format %12.8e [MD++_Get curstep]] [format %12.8e [MD++_Get EPOT]] [format %12.8e [MD++_Get Tinst]] [format %12.8e [MD++_Get N_lgst_cluster]]"
     flush $fileID

     #if { $x == 1 } {
     #   openwindow
     #} else {
     #   MD++ eval plot
     #}

     #MD++ totalsteps = 1000 T_OBJ = 800 
     #MD++ run

     #MD++ sleepseconds = 2 sleep
  }
  close $fileID

  MD++ finalcnfile = last_intercn.lammps writeLAMMPS
  MD++ sleepseconds = 100 sleep
  exitmd

} elseif { $status == 2 } {
  # read LAMMPS dump file and continue with MD simulation
  MD++ setnolog
  initmd $status-$n

  #set lammps_folder $n
  #set lammps_folder .
  set lammps_folder ../lammps-convert-$status

  set dump_freq 100000
  set file_id   0

  readmeam-lammps

  setup_QLM
  setup_md
  setup_window

  set y [expr ($dump_freq*$file_id)]
  MD++ incnfile= "$lammps_folder/dump_$y.md" input=\[ $y 0 \] readLAMMPS 

  #MD++ input = \[ 1 -10 10 -10 -0.20 -10 10 \] fixatoms_by_position 
  MD++ input = \[ 1 -10 10 -10 -0.28 -10 10 \] fixatoms_by_position 
  MD++ write_all = 1 finalcnfile=init.cn writecn

  MD++ allocUMBorder eval calcrystalorder
  puts "N_lgst_cluster =  [ MD++_Get N_lgst_cluster ]"

  MD++ totalsteps = 2000000 T_OBJ = 800 

  #openwindow

  # Randomize initial velocities
  MD++ srand48bytime
  MD++ initvelocity_type = Gaussian  initvelocity
  MD++ zerototmom zerorotation

  MD++ run

  exitmd

} elseif { $status ==3 } {
  #create new configurations 
  if { $argc <= 3 } {
  MD++ quit
} elseif { $argc > 3 } {
 set p [lindex $argv 2]
 set temp [lindex $argv 3]
}    


 #MD++ setnolog
  initmd $status-$n
  
  readmeam-lammps  

  set TEMP [expr $temp]
  set SiFrac [expr $p]
  set orien 111

  create_droplet $SiFrac $orien $TEMP
  

  relax_fixbox
  MD++ Fold_into_Unitcell
  
  set NP_Droplet [MD++_Get NP]
  set EquilFrac [ GetEquilSiFrac $TEMP ]
  set N_Siequil [expr int($NP_Droplet*$EquilFrac/100.0) ]
  puts "N_Siequil = ${N_Siequil}"

  MD++ {
  writeall = 1 
  finalcnfile = droplet.cn writecn
  finalcnfile = droplet.cfg writeatomeyecfg
 }
 #MD++ sleep 

create_SiNW $SiFrac $orien
if { $orien == 111 } {
MD++ { input = [ 0.0 -0.05 0.0 ] pbcshiftatom 
input = [ 1 -10 10 -10 -0.20 -10 10 ] fixatoms_by_position }
MD++ input = \[ 2 2 0.12 \] redefinepbc
} elseif { $orien == 100 } {
MD++ { input = [ 0.0 -0.05 0.0 ] pbcshiftatom 
input = [ 1 -10 10 -10 -0.20 -10 10 ] fixatoms_by_position }
MD++ input = \[ 2 2 0.07 \] redefinepbc
} elseif { $orien == 110 } {
MD++ { input = [ 0.0 -0.05 0.0 ] pbcshiftatom 
input = [ 1 -10 10 -10 -0.20 -10 10 ] fixatoms_by_position }
MD++ input = \[ 2 2 0.12 \] redefinepbc
} elseif { $orien == 112 } {
MD++ { input = [ 0.0 -0.1 0.0 ] pbcshiftatom 
input = [ 1 -10 10 -10 -0.20 -10 10 ] fixatoms_by_position }
MD++ input = \[ 2 2 0.15 \] redefinepbc
}

MD++ finalcnfile = SiNW.cn writecn

MD++ {
incnfile = droplet.cn input = [ 2 0 ] splicecn
input = [ 2 2 0.5 ] redefinepbc
}

set therm_exp [ expr 1.3E-5*$TEMP ]
### take care of thermal expansion here
MD++ input = \[ 1 1 ${therm_exp} \] shiftbox
MD++ input = \[ 2 2 ${therm_exp} \] shiftbox
MD++ input = \[ 3 3 ${therm_exp} \] shiftbox

relax_fixbox
setup_md
setup_window
#openwindow

MD++ conj_step = 0


MD++ initvelocity_type = Gaussian
MD++ T_OBJ = $TEMP initvelocity
MD++ zerototmom zerorotation

MD++ totalsteps = 20000 run
setup_QLM
MD++ allocUMBorder eval
set Nmax [ MD++_Get N_lgst_cluster ]
puts "$Nmax"
set NP_TOT [MD++_Get NP]
set NP [ expr $NP_TOT - $Nmax ]  #number of atoms in liquid
set N_Si [expr int($NP*$SiFrac/100.0) ] #objective Si atoms in liquid
set dN_Si [ expr ${N_Si} - ${N_Siequil} -(${NP} - ${NP_Droplet})] #the difference
puts "N_Si = ${N_Si}" 
puts "N_Siequil = ${N_Siequil}"
puts "dN_Si = ${dN_Si}"
AdjustSiFrac ${dN_Si} ${NP}
MD++ nspecies = 2
MD++ finalcnfile = initstruc_$n.lammps writeLAMMPS
MD++ finalcnfile = initstruc_$n.cfg writeatomeyecfg
MD++ continue_curstep = 1
#MD++ totalsteps = 5000000 run
MD++ sleepseconds = 100 sleep
exitmd

} elseif { $status == 4 } {
  #run MD simulation of Si crystal
  MD++ setnolog
  initmd $status-$n
   
  readmeam-lammps
  
  set SiFrac 100
  set TEMP 700
  set orien 110
 
  create_SiNW $SiFrac $orien

 
 
MD++  finalcnfile = SiNW.cn writecn
if { $n == 0 } {
  set therm_exp [ expr 1.004565476E-9*$TEMP*$TEMP + 1.2060036310E-6*$TEMP + 4.2183422619E-4 ]
} elseif { $n ==1 } {
  set therm_exp [ expr 1.36E-5*$TEMP ]
}

MD++ input = \[ 1 1 ${therm_exp} \] shiftbox
#MD++ input = \[ 2 2 ${therm_exp} \] redefinepbc
MD++ input = \[ 2 2 0.5 \] redefinepbc
MD++ input = \[ 3 3 ${therm_exp} \] shiftbox

relax_fixbox
setup_md
setup_window
#openwindow

MD++ initvelocity_type =Gaussian
MD++ T_OBJ = $TEMP initvelocity
MD++ zerototmom zerorotation
MD++ totalsteps = 100000 run
MD++ finalcnfile = SiNW.cfg writeatomeyecfg
exitmd

} elseif { $status == 5 } {
  # read cn file and continue with MD simulation
  MD++ setnolog
  initmd $status-$n

  readmeam-lammps

  setup_QLM
  setup_md
  setup_window

  MD++ incnfile= "~/Codes/MD++.svn/inter0001.cn" readcn

  #MD++ input = \[ 1 -10 10 -10 -0.20 -10 10 \] fixatoms_by_position 
  MD++ input = \[ 1 -10 10 -10 -0.28 -10 10 \] fixatoms_by_position 
  MD++ write_all = 1 finalcnfile=init.cn writecn

  MD++ allocUMBorder eval calcrystalorder
  puts "N_lgst_cluster =  [ MD++_Get N_lgst_cluster ]"

  MD++ totalsteps = 20000 T_OBJ = 630 

  #openwindow

  # Randomize initial velocities
  MD++ srand48bytime
  MD++ initvelocity_type = Gaussian  initvelocity
  MD++ zerototmom zerorotation

  MD++ run

  exitmd

} elseif { $status == 10 } {
  # Visualization
  MD++ setnolog
  initmd $status

  set lammps_folder $n

  set dump_freq 100000
  set file_start 0 
  set file_end   63

  readmeam-lammps
  set fileID [open "analysis.out" w]

  setup_window
  #openwindow

  setup_QLM

  for {set x $file_start} {$x <= $file_end} {incr x} {
     set y [expr ($dump_freq*$x)]
     MD++ incnfile= "$lammps_folder/dump_$y.md"
     MD++ input=\[ $y	0 \]
     MD++ readLAMMPS

     #MD++ writeall = 1 finalcnfile=config_$y.cn writecn
     #MD++ finalcnfile = config_$y.cfg writeatomeyecfg

     MD++ allocUMBorder eval

     if { $x == $file_start} {
        openwindow
     } else {
        MD++ eval plot
     }
     MD++ sleepseconds = 2 sleep
  }

  MD++ sleepseconds = 100 sleep
  exitmd

} else {

  puts "unknown status = $status"
  exitmd

}