# -*-tcl-script-*-
# MD++ input file for BCC Ta
#
# compile MD++ by, e.g.
#  make fs build=R SYS=intel
#
# run MD++ by, e.g.
#  bin/fs_intel scripts/work/cnm_chap/ta_perfect_crystal.tcl 0
#
# the last command line argument sets the $status variable
#  which specifies the task for MD++ to do
#
# status = 0 : create perfect crystal and determine equilibrium lattice constant
# status = 1 : compute C11, C12, C44
# status = 2 : NVE MD simulation
# status = 3 : NVT MD simulation to compute thermal expansion coefficient
#
# For more information on MD++, go to
#  http://micro.stanford.edu/wiki/MD++_Manuals
#

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

#*******************************************
# Definition of procedures
#*******************************************
proc initmd { } {
#MD++ setnolog
MD++ setoverwrite
MD++ dirname = runs/ta-perf
MD++ NNM = 200
}

proc readpot { } { MD++ {
#--------------------------------------------
#Read in potential file
#
potfile = ../../potentials/ta_pot readpot
} }

proc create_crystal { n a } {
   
    # create perfect crystal of specified size n and lattice constant a
    # equilibrium lattice constant of Ta is 3.3058 Angstrom
    MD++ crystalstructure = body-centered-cubic
    MD++ latticesize = \[ 1 0 0 $n  0 1 0 $n  0 0 1 $n \]
    MD++ latticeconst = $a

    MD++ makecrystal 
    MD++ finalcnfile = perf.cn writecn #eval
}


proc setup_window { } { MD++ {
#--------------------------------------------
# Plot setting
#
atomradius = [1.0 0.78] bondradius = 0.3 bondlength = 0 #2.8285 #for Si
win_width=400 win_height=400
#atomradius = 0.9 bondradius = 0.3 bondlength = 0 
atomcolor = cyan highlightcolor = purple  bondcolor = red
fixatomcolor = yellow backgroundcolor = gray70
plot_atom_info = 1 # display reduced coordinates of atoms when clicked
#plot_atom_info = 2 # display real coordinates of atoms
#plot_atom_info = 3 # displayenergy of atoms
plotfreq = 10
#
rotateangles = [ 0 0 0 1.2 ]
} }

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

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

#--------------------------------------------
proc setup_md { } { MD++ {     
equilsteps = 0  timestep = 0.001 # (ps)
atommass = 180.94788 # (g/mol)
DOUBLE_T = 1
saveprop = 1 savepropfreq = 10 openpropfile #run
savecn = 1 savecnfreq = 1000 openintercnfile
savecfg = 1 savecfgfreq = 1000
plotfreq = 100 printfreq = 10
vt2 = 1e28  #1e28 2e28 5e28
wallmass = 2e4     # atommass * NP 
boxdamp = 1e-3     #
saveH # Use current H as reference (H0), needed for specifying stress
fixboxvec  = [ 0 1 1 
               1 0 1
               1 1 0 ]
output_fmt = "curstep EPOT KATOM Tinst HELM HELMP TSTRESSinMPa_xx TSTRESSinMPa_yy TSTRESSinMPa_zz H_11 H_22 H_33" 
} }
#end of proc setup_md


#--------------------------------------------
#proc datafile_process { filename index frac fracend operation } {
#   puts "total line $NUM \n"
#   set Nini [ expr round($NUM*$frac) ] 
#   set Nend [ expr round($NUM*$fracend) ]

proc datafile_process { filename index Nlast operation } {
   set fp [ open $filename r ]
   set data [ split [read $fp] \n]
   close $fp
   set NUM [ expr [ llength $data ] - 1 ]

   set Nend $NUM
   set Nini [ expr $Nend - $Nlast - 1 ]

   set Sum 0
   set k   0
   set Var 0
   set Std 0
   for { set i $Nini } { $i < $Nend } { incr i 1 } {
       set k [expr $k+1]
       set data1 [ lindex $data $i ]
       split $data1 
       set datum [ lindex $data1 $index ]
       if { $i == $Nini } { 
	set MAX [lindex $data1 $index]
 	set MIN [lindex $data1 $index]
       } elseif { $i > $Nini } {
        set MAX [ expr ($MAX>$datum)?$MAX:$datum ]
        set MIN [ expr ($MIN<$datum)?$MIN:$datum ]
       }
       set Sum [expr $Sum+$datum]
   }
   set Ave [ expr $Sum/$k]

   if { [ string match "*STD*" $operation ] || [ string match "*VAR*" $operation ] } {
      for { set i $Nini } { $i < $Nend } { incr i 1 } {
          set data1 [ lindex $data $i ]
          split $data1
          set datum [ lindex $data1 $index ]  
	  set Var [expr $Var+($datum-$Ave)*($datum-$Ave) ]
      }
      set Var [ expr $Var/$k ]
      set Std [ expr sqrt($Var) ]
   }
   split $operation
   set Nvar [ llength $operation ]
   for { set i 0 } { $i < $Nvar } { incr i 1 } {
      set var [ lindex $operation $i ]
      if { $var=="SUM" } { 
          lappend LIST $Sum
      } elseif { $var=="AVE" } {
          lappend LIST $Ave
      } elseif { $var=="MAX" } {
          lappend LIST $MAX
      } elseif { $var=="MIN" } {
          lappend LIST $MIN
      } elseif { $var=="VAR" } {
          lappend LIST $Var
      } elseif { $var=="STD" } {
          lappend LIST $Std
      }
   }
   return $LIST
}

#*******************************************
# Main program starts here
#*******************************************
# status 0:
#        1:
#        2:
#
# read status from command line argument

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


# Do different tasks depending on  
if { $status == 0 } {
  # prepare perfect crystal and compute energy as a function of a
  initmd 
  readpot

  set nx 5

  # open data file to store results
  set fid [open "EPOT_a_coarse.dat" "w"]
  for { set i -10 } { $i <= 10 } { incr i 1 } {

       set a [expr 3.306 + 0.001*$i]
       create_crystal $nx $a
       MD++ eval
       set EPOT [MD++_Get "EPOT"]
       set NP   [MD++_Get "NP"  ]
       set E0   [expr $EPOT/$NP ]
       puts $fid "$a $NP $E0"

  }
  close $fid

  set fid [open "EPOT_a_fine.dat" "w"]
  for { set i -10 } { $i <= 10 } { incr i 1 } {

       set a [expr 3.3058 + 0.00001*$i]
       create_crystal $nx $a
       MD++ eval
       set EPOT [MD++_Get "EPOT"]
       set NP   [MD++_Get "NP"  ]
       set E0   [expr $EPOT/$NP ]
       puts $fid "$a $NP $E0"

  }
  close $fid

  # visualization of perfect crystal
  openwindow
  MD++ sleep

  # to plot the data and to fit a0, Ecoh, B
  # run fit_a0EB('EPOT_a_coarse.dat',1,2) or
  #     fit_a0EB('EPOT_a_fine.dat',  1,2) in Matlab or Octave
  # you can get fit_a0EB.m by the following command
  #  wget http://micro.stanford.edu/mediawiki/images/2/26/Fit_a0EB.m.txt -O fit_a0EB.m

  exitmd

} elseif { $status == 1 } {
  # determine elastic constants C11, C12, C44
  initmd
  readpot

  set nx 5
  set a0 3.3058
  create_crystal $nx $a0
  MD++ saveH eval

  # store volume at equilibrium
  set vol0 [MD++_Get OMEGA]

  # open data file to store results
  set fid [open "C11_C12.dat" "w"]
  for { set i -10 } { $i <= 10 } { incr i 1 } {

       # apply strain to simulation cell
       set eps_xx [expr 1e-5*$i]
       MD++ restoreH input = \[ 1 1 $eps_xx \] changeH_keepS clearR0
       
       # relaxation is needed only if basis has more than one atom
       # relax_fixbox

       # record stress values
       MD++ eval
       set sig_xx [expr (-1)*[MD++_Get TSTRESSinMPa_xx]]
       set sig_yy [expr (-1)*[MD++_Get TSTRESSinMPa_yy]]
       set EPOT   [MD++_Get EPOT]
    
       puts $fid "[format %18.14e $eps_xx]\t \
                  [format %18.14e $sig_xx]\t \
                  [format %18.14e $sig_yy]\t \
                  [format %18.14e [expr $EPOT/$vol0]]"

  }
  close $fid

  set fid [open "C44.dat" "w"]
  for { set i -10 } { $i <= 10 } { incr i 1 } {

       # apply strain to simulation cell
       set eps_xy [expr 1e-5*$i]
       MD++ restoreH input = \[ 1 2 $eps_xy \] changeH_keepS
       
       # relaxation is needed only if basis has more than one atom
       # relax_fixbox

       # record stress values
       MD++ eval
       set sig_xy [expr (-1)*[MD++_Get TSTRESSinMPa_xy]]
       set EPOT   [MD++_Get EPOT]
    
       puts $fid "[format %18.14e $eps_xy]\t \
                  [format %18.14e $sig_xy]\t \
                  [format %18.14e [expr $EPOT/$vol0]]"

  }
  close $fid

  # to plot the data and to fit C11, C12, C44
  # run fit_C11C12C44('C11_C12.dat','C44.dat') 
  # you can get fit_C11C12C44.m by the following command
  #  wget http://micro.stanford.edu/mediawiki/images/2/26/Fit_C11C12C44.m.txt -O fit_C11C12C44.m

  exitmd
  
} elseif { $status == 2 } {
  # NVE MD simulations
  initmd
  readpot

  set nx 5
  set a0 3.3058
  create_crystal $nx $a0

  setup_md
  MD++ ensemble_type = "NVE" integrator_type = "VVerlet"

  setup_window
  #openwindow

  MD++ T_OBJ = 300 DOUBLE_T = 1 randseed = 12345 srand48 initvelocity
  MD++ totalsteps = 10000 run
  
  
  #MD++ sleep
  exitmd
  
} elseif { $status == 3 } {
  # NPT MD simulation to determine thermal expansion coefficient
  initmd 
  readpot
 
  set nx 5
  set a0 3.3058
  create_crystal $nx $a0

  setup_md

  MD++ ensemble_type = "NPTC" integrator_type = "Gear6" timestep = 0.0001 #ps
  MD++ NHChainLen = 4 NHMass = \[ 2e-3 2e-6 2e-6 2e-6 \]
  MD++ savepropfreq = 100 savecnfreq = 10000 savecfgfreq = 10000 printfreq = 100

  setup_window
  #openwindow

  MD++ T_OBJ = 300 DOUBLE_T = 0 randseed = 67890 srand48 initvelocity
  MD++ totalsteps = 1000000 run
  
  #MD++ sleep
  exitmd

} elseif { $status == 4 } {
  # NVT MD simulations that adjusts volume to reach zero stress
  initmd 
  readpot
 
  set nx 5
  set a0 3.3058
  create_crystal $nx $a0
  set Lx0 [MD++_Get H_11]

  setup_md

  # allow thermal expansion
  #set eps_T 0.00197
  #set eps_T 0.00200
  #set eps_T 0.00203
  #set eps_T 0.00205
  #MD++ input = \[ 1 1 $eps_T \] changeH_keepS
  #MD++ input = \[ 2 2 $eps_T \] changeH_keepS
  #MD++ input = \[ 3 3 $eps_T \] changeH_keepS
 
  setup_window
  #openwindow

  setup_md

  MD++ ensemble_type = "NVTC" integrator_type = "VVerlet"
  MD++ NHChainLen = 4 NHMass = \[ 2e-3 2e-6 2e-6 2e-6 \]
  MD++ T_OBJ = 300 DOUBLE_T = 0 randseed = 12345 srand48 initvelocity

  set maxiter  100
  set mdsteps  1000
  set propfreq 10

  set fid [open "P_L.dat" "w"]

  for { set iter 1 } { $iter <= $maxiter } { incr iter 1 } {
      MD++ continue_curstep = 1
      MD++ totalsteps = $mdsteps savepropfreq = $propfreq run
      set Nlast [expr $mdsteps/$propfreq - 10]
      set sig_xx [datafile_process "prop.out" 6 $Nlast "AVE" ]
      set sig_yy [datafile_process "prop.out" 7 $Nlast "AVE" ]
      set sig_zz [datafile_process "prop.out" 8 $Nlast "AVE" ]
      set pressure [ expr ($sig_xx + $sig_yy + $sig_zz) / 3.0 ]
      set Lx [MD++_Get H_11]
      puts "iter = $iter  p = $pressure  Lx = $Lx  strain = [expr ($Lx - $Lx0)/$Lx0]"
      puts $fid "$iter $pressure $Lx [expr ($Lx - $Lx0)/$Lx0]"
      flush $fid

      if { $iter < $maxiter } { 
         set eps_T [ expr $pressure * 5e-7 ]
         MD++ input = \[ 1 1 $eps_T \] changeH_keepS
         MD++ input = \[ 2 2 $eps_T \] changeH_keepS
         MD++ input = \[ 3 3 $eps_T \] changeH_keepS
      }
  }
  close $fid
  
  #MD++ sleep
  exitmd

} elseif { $status == 5 } {
  # NVT MD simulation to generate stress flucutation data

    MD++ timestep = 0.001 incnfile = equil1.cn readcn eval

    MD++ outpropfile = prop2.out
    setup_md
    
    MD++ T_OBJ = $T ensemble_type = "NVT" 
    MD++ integrator_type = "VVerlet" implementation_type = 2
    MD++ Dexx = 1e-4 
    
    MD++ {output_fmt = "curstep EPOT KATOM Tinst HELMP \
                        TSTRESS_xx TSTRESS_yy TSTRESS_zz \
                        TSTRESS_xy TSTRESS_xz TSTRESS_yz \
                        DVIRIAL(0) DVIRIAL(1) DVIRIAL(2) \
                        DVIRIAL(3) DVIRIAL(4) DVIRIAL(5) \
                        DVIRIAL(6) DVIRIAL(7) DVIRIAL(8) \
                        VIRIAL(0) VIRIAL(1) VIRIAL(2) \
                        VIRIAL(3) VIRIAL(4) VIRIAL(5) \
                        VIRIAL(6) VIRIAL(7) VIRIAL(8) \
                        PSTRESS(0) PSTRESS(1) PSTRESS(2) \
                        PSTRESS(3) PSTRESS(4) PSTRESS(5) \
                        PSTRESS(6) PSTRESS(7) PSTRESS(8) " }
    
    MD++ vt2 = 1e28 totalsteps = 500000 timestep = 0.001 
    MD++ savecn = 1 savecnfreq = 100000 
    MD++ continue_curstep = 1 curstep = 100000 
    MD++ run
    MD++ finalcnfile = equil-500000.cn writecn
    MD++ finalcnfile = equil-500000.cfg writeatomeyecfg
    exitmd

} else {
        
 puts "unknown status = $status"
 exitmd 

}