# -*-shell-script-*-
# 3D Phase Field simulation
#-------------------------------------------

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

#*******************************************
# Definition of procedures
#*******************************************
proc init_single_phasefield { n size type } {
MD++ setnolog
MD++ setoverwrite
MD++ dirname = "runs/phasefield-$n"
MD++ srand48bytime 

#--------------------------------------------
# Initialization
#
MD++ num_fields = 1 gridsize = 0.4 model_type = 0
MD++ NX = $size  NY = $size   NZ = $size  input = $type  initphasefield
#incnfile = ../phasefield-0/phi$size.cn readcn
#--------------------------------------------
# Physical parameters
MD++ U_00 = 1  EPS_00 = 1   # interaction strength
MD++ K_00 = 1               # kinetic coefficients

# uncomment the following line to have anisotropic surface energy
MD++ EPS1_00 = 0.6993    EPS2_00 = -3.6595    EPS3_00 = -1.7384
} 

#--------------------------------------------
proc init_multi_phasefield { n size type } {
MD++ setnolog
MD++ setoverwrite
MD++ dirname = "runs/phasefield-multi-$n"
MD++ srand48bytime 

#--------------------------------------------
# Initializtion
#
MD++ num_fields = 3 gridsize = 0.4 model_type = 1
MD++ NX = $size  NY = $size   NZ = $size  input = $type  initphasefield
#incnfile = ../phasefield-multi-0/phi$size.cn readcn
#--------------------------------------------
# Physical parameters
MD++ K_12 = 1  K_21 = 1   # 1: solid,  2: vapor
MD++ K_20 = 1  K_02 = 1   # 2: vapor,  0: liquid
MD++ K_01 = 1  K_10 = 1   # 0: liquid, 1: solid

MD++ U_12 = 2  U_21 = 2   # 1: solid,  2: vapor
MD++ U_20 = 2  U_02 = 2   # 2: vapor,  0: liquid
MD++ U_01 = 2  U_10 = 2   # 0: liquid, 1: solid

MD++ EPS_12 = 1.7  EPS_21 = 1.7   # 1: solid,  2: vapor
MD++ EPS_20 = 2.0  EPS_02 = 2.0   # 2: vapor,  0: liquid
MD++ EPS_01 = 1.8  EPS_10 = 1.8   # 0: liquid, 1: solid

MD++ EPS1_12 = 0.6993    EPS2_12 = -3.6595    EPS3_12 = -1.7384 # 1: solid,  2: vapor
MD++ EPS1_21 = 0.6993    EPS2_21 = -3.6595    EPS3_21 = -1.7384 # 1: solid,  2: vapor

} 

#--------------------------------------------
proc openwindow { } { MD++ {
# Plot Configuration
#
atomradius = 1.0 win_width = 600 win_height = 600 rotateangles = [ 0 90 0 1.2]
color00 = pink color01 = blue backgroundcolor = gray80
#plot_threshold = [ 0.45 0.55  0.45 0.55 0.45 0.55 ]
plot_threshold = [ 0.43 0.58  0.43 0.58 0.499 0.501 ]
#plot_threshold = [ 0.40 0.60  0.40 0.60 0.40 0.60 ]
openwin alloccolors rotate saverot plot
} }
#end of proc openwindow





#*******************************************
# Main program starts here
#*******************************************
# status 0:prepare liquid by MD (NPT) simulation
#        1:MC equilibration
#        2:compute free energy by MC switching simulation
#
# 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 } {
  # Monte Carlo simulation (brute force)
  set size 50

  init_single_phasefield $n $size 4
  #init_single_phasefield $n $size 3
  #init_single_phasefield $n $size 2
  #init_single_phasefield $n $size -1

  MD++ eval

  openwindow
  #MD++ sleep quit

  MD++ timestep = 1e-3 dynamics_type = 2 # 0/2: non-conserved or 1: conserved
  #MD++ timestep = 1e-3 dynamics_type = 0 # 0/2: non-conserved or 1: conserved
  #MD++ timestep = 1e-3 dynamics_type = 1 # 0/2: non-conserved or 1: conserved
  #MD++ timestep = 1e-4 dynamics_type = 1 # 0/2: non-conserved or 1: conserved

  MD++ totalsteps = 1000000 finalcnfile = s$size.cn writecn 
  MD++ plotfreq = 500 
  MD++ saveprop = 1 savepropfreq = 100 printfreq = 500
  MD++ savecn = 1   savecnfreq   = 10000 openintercnfile

  # shorten simulation
  MD++ totalsteps = 20000 savecnfreq = 1000

  MD++ srand48bytime openpropfile run
  MD++ finalcnfile = phi$size-a.cn writecn
  
  MD++ sleep quit
  exitmd

} elseif { $status == 1 } {
  # Monte Carlo simulation (brute force)
  set size 75

  init_multi_phasefield $n $size 5

  openwindow
  #MD++ sleep quit

  MD++ timestep = 1e-3 dynamics_type = 2 # 0/2: non-conserved or 1: conserved
  #MD++ timestep = 1e-3 dynamics_type = 0 # 0/2: non-conserved or 1: conserved

  MD++ totalsteps = 1000000 finalcnfile = s$size.cn writecn 
  MD++ plotfreq = 100 #500 
  MD++ saveprop = 1 savepropfreq = 100 printfreq = 50 #500
  MD++ savecn = 1   savecnfreq   = 10000 openintercnfile

  # shorten simulation
  MD++ totalsteps = 20000 savecnfreq = 1000

  MD++ srand48bytime openpropfile run
  MD++ finalcnfile = phi$size-a.cn writecn
  
  MD++ sleep quit
  exitmd

} else {

  puts "unknown status = $status"
  exitmd

}