# with relaxation in this code

##################################### 
# definition of functions
##################################### 
proc write_POSCAR { posfile a0 ax ay az bx by bz cx cy cz } {

  #set posfile "POSCAR"
  set fileID [open $posfile w ]

  puts $fileID "POSCAR for FCC Au (created by tcl)"
  puts $fileID "$a0"
  puts $fileID [format "%15.12f %15.12f %15.12f" $ax $ay $az]
  puts $fileID [format "%15.12f %15.12f %15.12f" $bx $by $bz]
  puts $fileID [format "%15.12f %15.12f %15.12f" $cx $cy $cz]
  puts $fileID "1"
  puts $fileID "Cartesian  (real coordinates r)"
  puts $fileID "0    0    0"

  close $fileID
} 

proc read_POSCAR { posfile } {

  #set posfile "POSCAR"
  set fileID [open $posfile r ]

  # first line is comment
  gets $fileID line

  # second line is a0
  gets $fileID line
  scan $line "%f" a0

  # now comes three lines specifying the cell
  gets $fileID line
  scan $line "%f %f %f" ax ay az

  gets $fileID line
  scan $line "%f %f %f" bx by bz

  gets $fileID line
  scan $line "%f %f %f" cx cy cz

  # skip the rest of the file
  close $fileID

  # return cell geometry as a list
  return [list $a0 $ax $ay $az $bx $by $bz $cx $cy $cz]
} 

proc get_stress { } {

  set fileID [open "get_stress" w ]
  puts $fileID "grep \"in kB\" OUTCAR | sed \'s/in kB //\' "
  close $fileID
  file attributes get_stress -permissions u+x

  exec ./get_stress

}

proc run_vasp { ncpu } {

  #if { $ncpu == 1 } { 
  #    exec ./vasp
  #} else {
      set fileID [open "run_vasp" w ]
      #puts $fileID "mpiexec --comm=ib -np $ncpu ../../../../bin/vasp.2008-11-15"
      puts $fileID "/opt/mpiexec/bin/mpiexec --comm=pmi -np $ncpu ../../../../bin/vasp.caiwei.mva2"
      close $fileID
      file attributes run_vasp -permissions u+x
      exec ./run_vasp
  #}

}

##################################### 
# beginning of main program
##################################### 

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

  # modulus in kilobar (100 MPa)
  set C11 2020
  set C12 1690
  set C44  450

  if { 0 } {
     # compute nu, E using isotropic approximation
     set lambda $C12
     set mu [expr ($C11-$C12)/2.0 ] 
     set nu [expr $lambda / (2.0*($lambda + $mu)) ]
     # lambda = 2*mu*nu / (1-2*nu)
     set YoungsModulus [expr 2*$mu*(1+$nu) ]
     puts "mu = $nu nu = $nu  E = $YoungsModulus"
  } else {
     # directly set nu, E
     set YoungsModulus 780
     set nu 0.44
     puts "nu = $nu  E = $YoungsModulus"
  }

  set S11 [expr 1/$YoungsModulus ]
  set S12 [expr -$nu/$YoungsModulus ]

  puts "S11 = $S11 S12 = $S12"

  set factor 0.4
  set sig_tol 0.1

  set a0 4.0615
  #set a0 4.062

  set ax0 [expr sqrt(6)/4]
  set ay0 0
  set az0 [expr -sqrt(2)/4]

  set bx0 [expr sqrt(6)/6]  
  set by0 [expr sqrt(3)/3]  
  set bz0 0

  set cx0 [expr sqrt(6)/4]
  set cy0 0
  set cz0 [expr sqrt(2)/4]

  write_POSCAR "POSCAR.0" $a0 $ax0 $ay0 $az0 $bx0 $by0 $bz0 $cx0 $cy0 $cz0
  #run_vasp $ncpu

  puts "bx0 = $bx0"

  set datafile "stress_rlx.dat"
  set rlxfileID [open $datafile w ]

  set max_dx_iter  50
  set stop_dx_iter 50
  set start_dx_iter 1
  set reset_poscar  1
  set max_rlx_iter 20

  for { set dx_iter $start_dx_iter } { $dx_iter <= $stop_dx_iter } { incr dx_iter 1 } {

      # set cell dimension: bx
      set dx [expr sqrt(6)/6 * $dx_iter/$max_dx_iter]
      set bx [expr $bx0 + $dx]
      puts "dx = $dx bx = $bx"

      # set cell dimensions other than bx
     if { $dx_iter == $start_dx_iter } {
      if { $reset_poscar } {
          # initialize cell value
          set ax $ax0
          set cx $cx0
          set by $by0
          set az $az0
          set cz $cz0
          write_POSCAR "POSCAR"          $a0 $ax $ay0 $az $bx $by $bz0 $cx $cy0 $cz
          write_POSCAR "POSCAR.$dx_iter" $a0 $ax $ay0 $az $bx $by $bz0 $cx $cy0 $cz

      } else {
          # read cell value from POSCAR 
          set cell [ read_POSCAR "POSCAR" ]
          set a0 [lindex $cell 0 ]
          set ax [lindex $cell 1 ]
          set ay [lindex $cell 2 ]
          set az [lindex $cell 3 ]
          # set bx [lindex $cell 4 ]
          set by [lindex $cell 5 ]
          set bz [lindex $cell 6 ]
          set cx [lindex $cell 7 ]
          set cy [lindex $cell 8 ]
          set cz [lindex $cell 9 ]

          puts "read cell geometry:"
          puts " $a0"
          puts " $ax $ay $az"
          puts " $bx $by $bz"
          puts " $cx $cy $cz"
      }
     }

      # adjusting box to relax normal stress components to zero
      for { set rlx_iter 0 } { $rlx_iter <= $max_rlx_iter } { incr rlx_iter 1 } {

          run_vasp $ncpu

          set stress [ get_stress ]
          scan $stress "%f %f %f %f %f %f" sig_xx sig_yy sig_zz sig_xy sig_yz sig_zx

          puts $rlxfileID "$dx $rlx_iter $stress"
          flush $rlxfileID

          if { $rlx_iter < $max_rlx_iter } {

              if { 0 } {
                  # account for Poisson's effect
                  set e_xx [expr $sig_xx * $S11 + $sig_yy * $S12 + $sig_zz * $S12]
                  set e_yy [expr $sig_xx * $S12 + $sig_yy * $S11 + $sig_zz * $S12]
                  set e_zz [expr $sig_xx * $S12 + $sig_yy * $S12 + $sig_zz * $S11]

                  set bx [expr $e_yy * ($bx0+$dx) * $factor + $bx ]
              } else {
                  # do not account for Poisson's effect
                  set e_xx [expr $sig_xx / $C11 ]
                  set e_yy [expr $sig_yy / $C11 ]
                  set e_zz [expr $sig_zz / $C11 ]
              }

              puts "sig_xx = $sig_xx  sig_yy = $sig_yy  sig_zz = $sig_zz"
              puts "e_xx = $e_xx  e_yy = $e_yy  e_zz = $e_zz"

              set ax [expr $e_xx * $ax0 * $factor + $ax ]
              set cx [expr $e_xx * $cx0 * $factor + $cx ]
              set by [expr $e_yy * $by0 * $factor + $by ]
              set az [expr $e_zz * $az0 * $factor + $az ]
              set cz [expr $e_zz * $cz0 * $factor + $cz ]

              write_POSCAR "POSCAR" $a0 $ax $ay0 $az $bx $by $bz0 $cx $cy0 $cz
              write_POSCAR "POSCAR.$dx_iter.[expr $rlx_iter+1]" $a0 $ax $ay0 $az $bx $by $bz0 $cx $cy0 $cz
          }
      }
  }  
  close $rlxfileID

# End of program