#-*-shell-script-*-
############################################################################### 
# MD++ simulation of graphene sheet using REBO potential
# 
#           Saad Bhamla and Yanming Wang 
#                     June, 2012
##############################################################################
# To compile: 
# make rebo build=R SYS=mac (or other specified system)
#
# To run: (examples)
# bin/rebo_mac scripts/ME346/graphene_md.tcl $dir $stat $temp
# 
# zero-temperature simulation of perfect graphene sheet and calculating 
# Hessian matrix 
# bin/rebo_mac scripts/ME346/graphene_md.tcl 0 0 
# 
# zero-temperature simulation of graphene sheet with vacancy
# bin/rebo_mac scripts/ME346/graphene_md.tcl 0 1
#
# zero-temperature simulation with various lattice constants
# bin/rebo_mac scripts/ME346/graphene_md.tcl 0 2
#
# finite temperature (NVT) simulation
# bin/rebo_mac scripts/ME346/graphene_md.tcl 0 3 300
#
# finite temperature (NVT) simulation adjusting lattice constants
# bin/rebo_mac scripts/ME346/graphene_md.tcl 0 3 300
#
#
#
#--------------------------------------------
# Definition of procedures
#--------------------------------------------
#initialize the simulation
#
proc initmd { n } {
MD++ setnolog
MD++ setoverwrite
MD++ dirname = ./runs/graphene-$n # specify run directory
}

#-------------------------------------------------------------
# read potential files
#
proc read_rebo { } {
MD++ {
rebofile = "~/Codes/MD++/potentials/AIREBO/CH.airebo"
rcut = 3.0 NNM=120
readREBO
}}

#-------------------------------------------------------------
#Create Perfect Lattice Configuration
#
proc create_graphene { a } {
MD++ crystalstructure = graphene 
MD++ latticeconst = \[ $a 0 6.708 \]  #(A) for graphene

MD++ {latticesize = [  1   0  0  12  # c1 = 12*[1 0 0]
                 0   1  0  12  # c2 = 12*[0 1 0]
                 0   0  1  1 ] # c3 =1*[0 0 1]
makecrystal 
NNM = 120 # neighbor list
input = [0 0 0.5] pbcshiftatom # move the sheet to the middle 
input = [3 3 0.5] shiftbox # enlarge the box along z-direction
finalcnfile = perf.cn writecn
finalcnfile = perf.cfg  writeatomeyecfg
}}

#-------------------------------------------------------------
# Create one vacancy
#
proc create_vacancy {} {
MD++ {
input = [ 1      # how many atoms will be removed  
          35 ]   # the number of the removed atoms
fixatoms_by_ID
removefixedatoms
finalcnfile = vacancy.cn writecn
finalcnfile = vacancy.cfg writeatomeyecfg
}}

#-------------------------------------------------------------
#Plot Configuration
#
proc openwindow {} {
MD++ {
atomradius = 1.0 bondradius = 0.3 bondlength = 0
atomcolor = cyan highlightcolor = purple backgroundcolor = gray
atomcolor = blue highlightcolor = purple backgroundcolor = white 
bondcolor = red fixatomcolor = yellow
plot_color_bar = [ 1 -3.37 -3.32 ]  highlightcolor = red
plot_atom_info = 3 plotfreq = 10
rotateangles = [ 0 0 0 1.0  ]
#
win_width = 600 win_height = 600
openwin alloccolors rotate saverot refreshnnlist plot
}}

#-------------------------------------------------------------
#Conjugate-Gradient relaxation
#
proc relaxation {} {
MD++ {
conj_ftol = 5e-5 conj_itmax = 1000 conj_fevalmax = 1000
conj_fixbox = 0    # the simulation box is allowed to change
relax finalcnfile = relaxed.cn writecn
finalcnfile = relaxed.cfg writeatomeyecfg
}}

#-------------------------------------------------------------
#Setup MD simulation
#
proc setup_md { } { MD++ {     
equilsteps = 0  timestep = 0.0005 # (ps)
atommass = 12.0107
DOUBLE_T = 1
saveprop = 1 savepropfreq = 100 openpropfile #run
savecn = 1 savecnfreq = 100 openintercnfile
savecfg = 1 savecfgfreq = 2500
plotfreq = 10 printfreq = 10
vt2 = 1e28  #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  = [ 1 1 1 
               1 1 1
               1 1 0 ]
stress = [ 0 0 0
           0 0 0
           0 0 0 ]
output_fmt = "curstep EPOT KATOM Tinst HELM HELMP TSTRESSinMPa_xx TSTRESSinMPa_yy TSTRESSinMPa_zz H_11 H_22 H_33"
ensemble_type = "NVT" integrator_type = "VVerlet" implementation_type = 1
zerorot = "all" # zero out angular momentum when initvelocity
writeall = 1 # write velocity and species information into file
}}

#-------------------------------------------------------------
# Process the output datafile to find average, maximum, minimum, 
# variation and std for the data
#
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
##############################################################

if { $argc == 0 } {
set n 1
set status 0
set T_OBJ 300
} elseif {$argc == 1 } {
set n [lindex $argv 0]
set status 0
set T_OBJ 300
} elseif {$argc == 2 } {
set n [lindex $argv 0]
set status [lindex $argv 1]
set T_OBJ 300
} elseif {$argc > 2 } {
set n [lindex $argv 0]
set status [lindex $argv 1]
set T_OBJ [lindex $argv 2]
}

#-------------------------------------------------------------
if { $status == 0 } {
initmd $n
read_rebo
set a 2.420 
create_graphene $a
openwindow
MD++ eval
MD++ input = 0 
MD++ calHessian
MD++ sleep quit
}

#-------------------------------------------------------------
elseif { $status == 1 } {
initmd $n
read_rebo
set a 2.420
create_graphene $a
openwindow
relaxation
MD++ eval
MD++ sleep quit
}

#-------------------------------------------------------------
elseif { $status == 2 } {
initmd $n
read_rebo
set fileID [open "graphene.dat" w]

for { set i 0} { $i <= 240 } { incr i 1 } {
set a [expr (2.3 + $i * (0.001))]
create_graphene $a
#openwindow
MD++ eval
 set EPOT     [MD++_Get EPOT]
  set N        [MD++_Get NP]
  set VOL      [MD++_Get OMEGA]
  set ELAT     [expr $EPOT / $N]
  set rho      [expr $N / double($VOL)]
  set AtomVol  [expr double($VOL) / $N]

  puts $fileID "[format %18.12e $a]\t \
                [format %18.12e $rho] \t \
                [format %18.12e $AtomVol]\t \
                [format %18.12e $EPOT]\t \
                [format %18.12e $ELAT]"
}
close $fileID
MD++ sleep quit
}

#-------------------------------------------------------------
elseif { $status == 3 } {
initmd $n
read_rebo
set a 2.420
create_graphene $a
openwindow
relaxation
setup_md
MD++ T_OBJ = $T_OBJ  initvelocity
MD++ totalsteps = 10000 run
MD++ finalcnfile = equil.cn writeall = 1 writecn
MD++ sleep quit
}

#-------------------------------------------------------------
elseif { $status == 4 } {
initmd $n
read_rebo
set a 2.420
create_graphene $a
setup_md
MD++ T_OBJ = $T_OBJ  initvelocity
set Lx0 [MD++_Get H_11]
 
  set maxiter  500
  set mdsteps  8000
  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) / 2.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 * 2e-8 ]
         MD++ input = \[ 1 1 $eps_T \] changeH_keepS
         MD++ input = \[ 2 2 $eps_T \] changeH_keepS
 
      }
  }
  close $fid
 
MD++ finalcnfile = equil.cn writeall = 1 writecn
MD++ sleep quit
}

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