#---------------------------------------------------------------------------
# This is a Tcl script for computing melting points of pure elements in MD++
#
# This file can be run with various empirical potentials developed in MD++ such as
#
#  meam-lammps, meam-baskes, fs, eam, sw, swge, sworig, tersoff
#
# It has been tested for various elements such as Si, Au, etc.
# For details, see Seunghwa Ryu and Wei Cai, 
#  "Comparison of Thermal Properties Predicted by Interatomic Potential Models", 
# Modelling and Simulation in Materials Science and Engineering, 16, 085006, (2008).
#
#
# How to use:
#
#   Step 1.  Compile MD++ with the potential you want to use, for example
#
#             make sworig build=R
#                                      or
#             make meam-lammps build=R
#  
#            see http://micro.stanford.edu/wiki/MD++_Manuals
#            for more instructions on compiling MD++
#
#   Step 2.  Make a new subdirectory 
#
#             mkdir runs/Single_Elem
#
#   Step 3.  Run MD++.
#
#            Because there are many steps required to compute melting point,
#             we need to run the same tcl script over and over
#             with different input arguments.
#            The following jobs need to be run (using sworig as an example):
#  
#  (Determine equilibrium volume of liquid and solid and compute hessian matrix of solid)
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  1 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  7 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  2 1 SW_Si
#
#  (Determine free energy of solid)
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  3 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  4 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  5 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  6 1 SW_Si
#
#  (Determine free energy of liquid)
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  8 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  9 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl 10 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl 11 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl 12 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl 13 1 SW_Si
#
#  (Generate matlab files from simulation data to determine melting point
#      14: Hessian_Siz.m, 15: Free_Energy_of_Siz.m, 16: tar_Siz )
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl 14 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl 15 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl 16 1 SW_Si
#
#            If you have a computer cluster, you can write a automatic script
#             that submit these jobs simultaneously. But you need to pay attention
#             to the inter-dependence of some of the jobs.
#
#            Describe the meaning of the three input arguments.
#
#               First argument is switching step, i = 1,2,...,16. 
#
#               Second argument is job index, x
#                   for i=3,4,5,6 and i=8,9,10,11,12,13, 
#                   you can run multiple jobs with different x
#                   simultaneously to collect more statistics, for example,
#                   the following will run 4 jobs simultaneously for step 3:
#
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  3 1 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  3 2 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  3 3 SW_Si
#             bin/sworig_gpp scripts/Examples/Tcl/melting.tcl  3 4 SW_Si
#
#                Third one is material name, which has to follow a convention
#             to be described below.
#
#            Describe the inter-dependence among steps: 
#            
#            1 -> 2, 7    
#            2 -> 3, 4
#            3 -> 5, 6 
#            7 -> 8, 9
#            8 -> 12, 13
#            9 -> 10, 11
#            
#            Before starting right-side jobs,
#            some results of left-side jobs are needed.
#            Until left-side jobs give signal to 'flag' file,
#            right-side jobs will wait for the results.
#            e.g) 0K solid property from step 1 will be used for
#                 initial setup for step 7.
#                 T1 K solid property from step 1 will be used for
#                 quasi-harmonic crystal hessian matrix 
#                 computation in step 2.        
#           
#            Describe which step [1-16] :
#	        step 1 : determine thermal expansion of solid at T1 K 
#		step 2 : compute hessian matrix of quasi-harmonic crystal(QHC)
#		step 3 : adiabatic switching from real solid to QHC
#		step 4 : adiabatic switching from QHC to real solid
#		step 5 : reversible scaling from T1 to T1/lambda1
#               step 6 : reversible scaling from T1/lambda to T1
#            
#               step 7 : determine thermal expansion of liquid at T2 K
#		step 8 : adiabatic switching from real liquid to gaussian fluid
#		step 9 : adiabatic switching from gaussian fluid to real liquid
#		step 10: switching from gaussian fluid to ideal gas
#		step 11: switching from ideal gas to gaussian fluid
#		step 12: reversible scaling from T2 to T2/lambda2
#		step 13: reversible scaling from T2/lambda2 to T2
#
#               e.g) bin/sworig_gpp scripts/Examples/Tcl/melting.tcl 3 1 SW_Si
#   First argument 3 means "Run step 3 calculation"
#   Second argument 1 means "Save the results in ~~~_1.dat or ~~_1.out" 
#   Third arugement SW_Si means "Use initial condition assigned for SW_Si"
#
#
#    Finally, to get the melting point, you need to run the following in Matlab:
#
#    First, you may need to transfer important data
#       to another computer that has Matlab, 
#
#    You can do so by running bash script generated by step 16.
#
#      tar_Siz 
#
#    which will generate inter_Siz.tar
#
#    Transfer this file to another computer 
#    and untar it in your favorite directory.
#
#    Launch Matlab and run 
#
#    Hessian_Siz.m 
#    Free_Energy_of_Si.m
#  
#     They will generate a plot of solid and liquid free energy
#     as functions of temperature
#     Melting point will be printed to the screen.
#     More data will be saved into file SW_Si_Publication_Info.dat, 
#     including entropy, latent heat, error analysis, and etc.
# 
#---------------------------------------------------------------------------------------
#
# Advanced topics:
#
#    When running meam-lammps the last argument has to be the same
#             element type in the meamf file, such as
#
#             bin/meam-lammps single_element_MEAM_free_energy_wcr.tcl 3 1 Siz
#
#    To modify initial conditions such as T1, T2, reversible scaling factor,
#    change initial conditions specified at the begginning of function
#    definitions .
#
#    If you want to run with potentials other than sw/sworig, you will need to 
#     to uncomment some lines from line 783 to line 798 for details
#
#    Naming convention for the third argument is the following 
#
#      Ta, Mo, W, Siz, Cu, Au, Ag, Ge   => meam materials
#      ta, mo, w 			=> fs materials
#      auf, cuf			        => eam materials
#      SW_Si, SW_Ge			=> sw materials
#
#
#-----------------------------------------------------------------------



#**********************************************
# Function Definitions
#**********************************************

# default values
set material "Siz"
# Input On or Off for Graphic & Log  
set graphic "Off"

set material [ lindex $argv 2 ]

# for given job index x, we will have 1+iter1 iteration for each step
# iter2 is iteration number for ideal gas - Gaussian fluid step.
# iter1 is iteration number for other switching steps.
set iter1 1 
set iter2 3
# Initial Parameter Set
if {$material=="SW_Si"} {
# lattice constant a : this will be recomputed in step 1
   set a 5.431
# size of cubic cell L : size of L x L x L simulation cell will be used. 
   set L 4
# initial temperature for solid free energy computation.
   set T_solid_0 1600
# elastic strain at T_solid_0 : This will be computed in step 1.
# leave this for 
   set eT_solid_0 0.002
# by reversible scaling, free energy of solid will be obtained in the range of
# T=T_solid_0 ~ T=T_solid_0/factor_solid_RS
   set factor_solid_RS 0.8
# initial temperature for liquid free energy computation.
   set T_liquid_0 2000
# strain at T_liquid_0 : This will be computed in step 7.
   set eT_liquid_0 0.02
# by reversible scaling, free energy of liquid will be obtained in the range of
# T=T_liquid_0 ~ T=T_liquid_0/factor_liquid_RS
   set factor_liquid_RS 1.25
# in the intermediate step for liquid free energy computation, 
# we have ideal gas - Gaussian potential switching step.
# gaussian potential will be changed from U to factor_gauss*U in 1 switching.
# total six switching will lead (0.1)^6*U which is practically zero potential 
   set factor_gauss 0.1
# atomic mass of the material
   set atom_mass 28.086
# crystal structure of the material
   set CS "DC"
}
if {$material=="SW_Ge"} {
   set a 5.6575
   set L 4
   set T_solid_0 2500
   set eT_solid_0 0.002
   set factor_solid_RS 0.8
   set T_liquid_0 3125
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.25
   set factor_gauss 0.1
   set atom_mass 72.64
   set CS "DC"
}
if {$material=="auf"} {
   set a 4.078
   set L 6
   set T_solid_0 730
   set eT_solid_0 0.002
   set factor_solid_RS 0.2
   set T_liquid_0 1170
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.3
   set factor_gauss 0.1
   set atom_mass 196.9665
   set CS "FCC"
}
if {$material=="cuf"} {
   set a 3.615
   set L 6
   set T_solid_0 1080
   set eT_solid_0 0.00
   set factor_solid_RS 0.8
   set T_liquid_0 1350
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.25
   set factor_gauss 0.1
   set atom_mass 63.54
   set CS "FCC"
}
if {$material=="w"} {
   set a 3.165
   set L 6
   set T_solid_0 2400
   set eT_solid_0 0.002
   set factor_solid_RS 0.5
   set T_liquid_0 4500
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.25
   set factor_gauss 0.1
   set atom_mass 183.85
   set CS "BCC"
}
if {$material=="ta"} {
   set a 3.165
   set L 6
   set T_solid_0 3600
   set eT_solid_0 0.002
   set factor_solid_RS 0.8
   set T_liquid_0 4500
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.25
   set factor_gauss 0.1
   set atom_mass 180.948
   set CS "BCC"
}
if {$material=="mo"} {
   set a 3.165
   set L 6
   set T_solid_0 1600
   set eT_solid_0 0.002
   set factor_solid_RS 0.5
   set T_liquid_0 3400
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.3
   set factor_gauss 0.1
   set atom_mass 95.94
   set CS "BCC"
}
if {$material=="Siz"} {
   set a 5.431
   set L 4
   set T_solid_0 1600
   set eT_solid_0 0.002
   set factor_solid_RS 0.8
   set T_liquid_0 2000
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.25
   set factor_gauss 0.1
   set atom_mass 28.086
   set CS "DC"
}
if {$material=="TR_Si"} {
   set a 5.431
   set L 4
   set T_solid_0 2210
   set eT_solid_0 0.002
   set factor_solid_RS 0.85
   set T_liquid_0 2760
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.15
   set factor_gauss 0.1
   set atom_mass 28.086
   set CS "DC"
}
if {$material=="Ta"} {
   set a 3.303
   set L 6
   set T_solid_0 800
   set eT_solid_0 0.002
   set factor_solid_RS 0.2
   set T_liquid_0 4600
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.3
   set factor_gauss 0.1
   set atom_mass 180.948
   set CS "BCC"
}
if {$material=="Mo"} {
   set a 3.15
   set L 6
   set T_solid_0 800
   set eT_solid_0 0.002
   set factor_solid_RS 0.2
   set T_liquid_0 3400
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.4
   set factor_gauss 0.1
   set atom_mass 95.94
   set CS "BCC"
}
if { $material=="Au" } {
   set a 4.078
   set L 5
   set T_solid_0 840
   set eT_solid_0 0.002636
   set factor_solid_RS 0.7
   set T_liquid_0 1170
   set eT_liquid_0 -0.02486
   set factor_liquid_RS 1.3
   set factor_gauss 0.1
   set atom_mass 196.96655
   set CS "FCC"
}
if {$material=="Ge"} {
   set a 5.6575
   set L 4
   set T_solid_0 1200
   set eT_solid_0 0.0168
   set factor_solid_RS 0.8
   set T_liquid_0 1500
   set eT_liquid_0 -0.00388
   set factor_liquid_RS 1.25
   set factor_gauss 0.1
   set atom_mass 72.64
   set CS "DC"
}
if { $material == "Ag" } {   
   set a 4.08
   set L 5
   set T_solid_0 840
   set eT_solid_0 0.002
   set factor_solid_RS 0.7
   set T_liquid_0 1200
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.4
   set factor_gauss 0.1
   set atom_mass 107.87
   set CS "FCC"
}
if {$material=="Cu"} {
   set a 3.6
   set L 5
   set T_solid_0 1000
   set eT_solid_0 0.002
   set factor_solid_RS 0.7
   set T_liquid_0 1400
   set eT_liquid_0 0.02
   set factor_liquid_RS 1.4
   set factor_gauss 0.1
   set atom_mass 63.54
   set CS "FCC"
}

# initiate md
#------------------------------------------------------------
proc initmd { material n } { 
MD++ setoverwrite

if { $n < 7 } {
MD++ dirname = runs/Single_Elem_Tm/${material}_Solid 
} elseif { $n <14 } {
MD++ dirname = runs/Single_Elem_Tm/${material}_Liquid
} elseif { $n >=14 } { 
MD++ dirname = runs/Single_Elem_Tm
} 
MD++  zipfiles = 0 writeall = 1

} 

proc get_struc { CS } {
	if { $CS=="L1_2" || $CS=="NaCl" || $CS=="CsCl" } {
	return $CS
	} elseif { $CS=="ZB" } {
	return zinc-blende
	} elseif { $CS=="FCC" } {
	return face-centered-cubic
	} elseif { $CS=="BCC" } {
	return body-centered-cubic
	} elseif { $CS=="DC" } {
	return diamond-cubic
	} elseif { $CS=="HCP" } {
	return hexagonal-ortho
	}
}

# potential options
#------------------------------------------------------------
proc readmeam-baskes { } { MD++ {
#Read in MEAM potential (Baskes format)
#meamfile = "~/Codes/MD++/Fortran/MEAM-Baskes/meamf"
#meafile  = "~/Codes/MD++/Fortran/MEAM-Baskes/meafile_Au_Si"
meamfile = "/home/shryu/Codes/MD++/Fortran/MEAM-Baskes/meamf"
meafile  = "/home/shryu/Codes/MD++/Fortran/MEAM-Baskes/meafile_Au_Si_Stanford"
ntypes = 2  ename0 = "Au1  " ename1 = "Si4  " #leave spaces to be compatible with fortran
rcut = 6.0  kode0 = "library"
readMEAM    NNM = 300
} }

#------------------------------------------------------------
proc readmeam-lammps { material } { MD++ {
#Read in MEAM potential (Baskes format)
meamfile = "/home/shryu/Codes/MD++/potentials/meamf"
#meafile = "/home/shryu/Codes/MD++/potentials/AuSi.meam" 
nspecies = 1 }
MD++  element0 = $material
MD++ {
rcut = 6.0  readMEAM
NNM = 300
} }

#-------------------------------------------------------------
proc readeam { material } { 
#Read in MEAM potential (MD++ format)
MD++ potfile = "~/Codes/MD++/potentials/EAMDATA/Eamdata_${material}oiles.txt"
MD++ eamgrid = 500 readeam NNM=300
} 

# set species 
#----------------------------------------------------------------------
proc setspecies {material} { MD++ {
	nspecies=2
	element0="Au" element1="Si"
}
if { $material=="Au" } {
MD++	input=\[ 1 -10 10 -10 10 -10 10 \] fixatoms_by_position
MD++	input=0 setfixedatomsspecies freeallatoms
}
if { $material=="Si" } {
MD++ 	input=\[ 1 -10 10 -10 10 -10 10 \] fixatoms_by_position
MD++	input=1 setfixedatomsspecies freeallatoms
}
MD++	atommass=\[196.96655 28.0855\]
}

# set plot options
#----------------------------------------------------------------------
proc openplot {} { MD++ {
	atomradius = [ 0.8 0.5 ]
	win_width = 500 win_height = 500
	atomcolor1 = [ orange gray ]
        highlightcolor = purple
	fixedatomcolor = yellow backgroundcolor = gray70
	plot_atom_info = 2
	plotfreq=100
	rotateangles = [ 0 0 0 1.2 ]
	openwin alloccolors rotate saverot plot
}}

# set create_crystal
#----------------------------------------------------------------------
proc create_crystal { material CS a L } {

MD++     latticestructure=[get_struc $CS ]
MD++     latticeconst=$a
MD++     latticesize=\[ 1 0 0 $L 0 1 0 $L 0 0 1 $L \]
MD++     makecrystal

#	setspecies $material
#        openplot
}

# set create_liquid
#--------------------------------------------------------------------
proc create_liquid { material CS a L } { 

create_crystal $material $CS $a $L
MD++ srand48bytime randomposition 
}

# relaxation method
#--------------------------------------------------------------------
proc CGR_fixbox {} { MD++ {
        conj_ftol=2e-6 conj_itmax=2000 conj_fevalmax=20000
        conj_fixbox=1 conj_fixboxvec= [ 0 1 1
                                        1 0 1
                                        1 1 0 ]
        relax
}}

proc CGR_freebox {} { MD++ {
        conj_ftol=2e-6 conj_itmax=2000 conj_fevalmax=20000
        conj_fixbox=0 conj_fixboxvec= [ 0 1 1
                                        1 0 1
                                        1 1 0 ]
        relax
}}

# running condition setup
#-------------------------------------------------------------------
proc run_setup { n nn atom_mass } { 
MD++        atommass=$atom_mass
MD++ {
        equilsteps=0 timestep=0.0001
        DOUBLE_T=0
        savecfg=1 savecfgfreq=999999
        savecn=1 savecnfreq=999999
        saveprop=1 savepropfreq=100 }
    if { $n==1 || $n==7 } {
    MD++ savecnfreq=99999
    MD++ savecfgfreq=99999
    }
    MD++ outpropfile="MEAM_F${n}_${nn}.out"
    MD++ intercnfile="MEAM_F${n}_.cn"
    MD++ intercfgfile="MEAM_F${n}_.cfg"
    MD++ {
        openpropfile
        openintercnfile
        openintercfgfile
#       savecn=1 savecnfreq=10000 openintercnfile
        wallmass=2e3
        vt2=2e28
        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 HELM Tinst TSTRESS_xx TSTRESS_yy TSTRESS_zz H_11 H_22 H_33 OMEGA"
        writeall=1
}}

# scale supercell
#-----------------------------------------------------------------
proc scale_H { eT } {
        MD++ input=\[ 1 1 $eT \] shiftbox
        MD++ input=\[ 2 2 $eT \] shiftbox
        MD++ input=\[ 3 3 $eT \] shiftbox
}

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

# call potential
#--------------------------------------------
proc readmeam-according-to-myname { myname material } {
 if { [ string match "*baskes*" $myname ] } {
  puts "readmeam-baskes"
  readmeam-baskes 
 } elseif { [ string match "*lammps*" $myname ] } {
  puts "readmeam-lammps"
  readmeam-lammps $material
 } elseif { [ string match "*meam*" $myname ] } {
  puts "readmeam"
  readmeam
 } elseif { [ string match "*eam*" $myname ] } {
  puts "readeam"
  readeam
 } else {
  puts "not an eam potential, not reading any files"
 }
}

proc make_flag { n } {
   set fp [ open "flag" a+ ]
   puts $fp "$n"
   close $fp 
}

proc wait_flag { n } {
   set key -1
   if { $n==1 || $n==7 } {
   set key 1 
   set fp [ open "info.dat" w ]
   puts $fp "a0 L eT Ecoh V T NP"
   close $fp 
   }
   while { 1 } {
      set fp [ open "flag" a+ ]
      close $fp
      set fp [ open "flag" r ]
      set tags [ split [ read $fp ] \n ]
      set NUM [ llength $tags ]
      for { set i 0 } { $i < $NUM } {incr i } {
         if { $n == [ lindex $tags $i ] } {
	    set key 1
  	    break
         }
      }
      close $fp
      if { $key > 0 } { break }
      exec sleep 60
      
   }
} 

proc check_127 { n } {
      set key -1
      set fp [ open "flag" a+ ]
      close $fp
      set fp [ open "flag" r ]
      set tags [ split [ read $fp ] \n ]
      set NUM [ llength $tags ]

      for { set i 0 } { $i < $NUM } { incr i } {
         if { $n == [ lindex $tags $i ] } {
            set key 1
            break
         }
      }

      if { $key > 0 } { 
      close $fp
      MD++ quit
      } else {
      close $fp
      make_flag $n 
      }

}

proc get_data { type } {
   set fp [ open "info2.dat" r ]
   set data [ split [ read $fp ] ]
   close $fp
   if { $type=="a0" } { 
   return [ lindex $data 0 ]
   } elseif { $type=="L" } {
   return [ lindex $data 1 ]
   } elseif { $type=="eT" } {
   return [ lindex $data 2 ]
   } elseif { $type=="Ecoh" } {
   return [ lindex $data 3 ]
   } elseif { $type=="V" } {
   return [ lindex $data 4 ]
   } elseif { $type=="T" } {
   return [ lindex $data 5 ]
   } elseif { $type=="NP" } {
   return [ lindex $data 6 ]
   }
}   

proc fileclose { n } {
   exec mv A.log A$n.log
   MD++ quit
}

proc datafile_process { filename index frac fracend operation } {
   set fp [ open $filename r ]
   set data [ split [read $fp] \n]
   close $fp
   set NUM [ llength $data ]
   set NUM [ expr $NUM-1 ]
   set Nini [ expr round($NUM*$frac) ] 
   set Nend [ expr round($NUM*$fracend) ]
#   puts "total line $NUM \n"
   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 ]
       }
#       puts "$datum"
       set Sum [expr $Sum+$datum]
   }
   set Ave [ expr $Sum/$k]
#   puts "$Sum $Ave"
#   puts "$MAX $MIN"

    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
     }
   }
#   puts "$Std $Var" 
#   puts "$LIST"
   return $LIST
}
#}      
   

#****************************************
# MAIN PROGRAM START
#****************************************

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

if { $graphic=="On"} {
MD++ setnolog
openplot
}

initmd $material $n

# wait until any computation need for step n is done #
if { $n < 14 } { 
wait_flag $n
}
if { $n==1 || $n==7 } {
check_127 $n
}

#### Use following lines for different potentials ####

### 1 meam-lammps ###
set myname [ MD++_Get "myname" ]
readmeam-according-to-myname $myname $material

### 2 fs potential ###
#MD++ potfile = ~/Codes/MD++/potentials/${material}_pot readpot NNM=200

### 3 eam potentail ###
#readeam $material

### 4 SW potential ###
# comment out all the lines above 

### 
#openplot

### set up initial temperature T1 and T2 for solid and liquid ### 
if { $n<=6 } {
    set T $T_solid_0
    set factor $factor_solid_RS
} else {
    set T $T_liquid_0
    set factor $factor_liquid_RS
}


### equilibrate solid to measure thermal strain eT at T1 ###  
if { $n==1 } { 

# get 0K properties for solid and save into info1.dat #
    create_crystal $material $CS $a $L
    CGR_freebox
    MD++ eval
    set NP [MD++_Get NP] 
    set EPOT0 [MD++_Get EPOT]
    set Ecoh0 [expr $EPOT0/$NP]
    set V0 [MD++_Get OMEGA]
    set a0 [expr pow($V0,1.0/3)/$L]     
    set fp [ open "info1.dat" w ]
    puts $fp "$a0 $L 0 $Ecoh0 $V0 0 $NP"
    close $fp

# run MD simulation of solid at T1 #
    run_setup $n $nn $atom_mass
    MD++ NHMass=\[ 1e-2 1e-2 1e-2 1e-2 \] T_OBJ=$T
    MD++ ensemble_type="NPTC" integrator_type="Gear6"
    MD++ NHChainLen = 4 initvelocity totalsteps=1000000
    MD++ run

# compute average volume at T1 #
    set V [datafile_process "MEAM_F1_${nn}.out" 11 0.7 1 "AVE" ]
    set a [expr pow($V,1.0/3)/$L]
    set eT [expr $a/$a0-1]

# compute thermal strain eT and Ecoh at T1 and save into info2.dat #      
    create_crystal $material $CS $a $L
    CGR_fixbox
    MD++ eval
    set EPOT [ MD++_Get EPOT ]
    set Ecoh [ expr $EPOT/$NP ]
    set fp [ open "info2.dat" w ]
    puts $fp "$a0 $L $eT $Ecoh $V $T $NP"  
    close $fp
# activate step 2 #
    make_flag 2
    MD++ quit
} elseif { $n==2 } {
# step 2 is for 'computing hessian matrix for quasi-harmonic crystal #
# read info2.dat and make crystal with strain eT #
    create_crystal $material $CS [get_data a0] [get_data L]
    scale_H [get_data eT]
    CGR_fixbox

# compute hessian matrix to get quasi-harmonic crystal(QHC) #
    MD++ input=0 calHessian
    exec mv hessian.out hessian$T.out
    MD++ setoverwrite
# activate step 3 and 4 #
    make_flag 3
    make_flag 4
    MD++ quit
} elseif { $n==3 } {
# step 3 is for 'switching from real solid potential to QHC' #
# create crystal with strain eT #
    create_crystal $material $CS [get_data a0] [get_data L]
    scale_H [get_data eT]

    MD++ finalcnfile="${T}K_perf_crystal.cn" writecn setconfig1
    MD++ incnfile="hessian$T.out"
    MD++ readHessian

# equilibrate solid before switching run #
    run_setup $n $nn $atom_mass
    MD++ srand48bytime NHMass=\[ 1e-2 1e-2 1e-2 1e-2\] totalsteps=100000
    MD++ ensemble_type="NVTC" integrator_type="Gear6" NHChainLen=4
    MD++ T_OBJ=$T initvelocity saveprop=0 
    MD++ run
    MD++ finalcnfile="${T}K_equil_crystal.cn" writeall=1 writecn
    MD++ incnfile="${T}K_equil_crystal.cn"

# activate step 5 and step 6
    make_flag 5
    make_flag 6
# switching setup
    MD++ Ecoh=[get_data Ecoh]
    MD++ switchfreq=10 saveprop=1 savepropfreq=100 openpropfile printfreq=5000
    MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33"


for { set x 0 } {$x<=$iter1} {incr x 1} {
	set y [format %04d $x]
	MD++ readcn totalsteps=50000 saveprop=0

# equilibration run
	if { $iter1>0 } {
	MD++ initvelocity run
        }
# switching run  
	MD++ lambda0=0 lambda1=1 refpotential=1
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH
	set datafile "MEAM_F$n.dat"
	set fileID [ open $datafile a+]
	puts $fileID "$x [MD++_Get Wavg] [MD++_Get Wtot]"
	close $fileID        
	MD++ finalcnfile=f$n$y.cn writeall=1 writecn
}
MD++ quit
}  elseif { $n==4 } {
# step 4 is for 'switching from QHC to real solid' # 
create_crystal $material $CS [get_data a0] [get_data L]

scale_H [get_data eT]
MD++ finalcnfile="${T}K_perf_crystal.cn" writeall=1 writecn setconfig1
MD++ incnfile="hessian$T.out"
MD++ readHessian

# equilibrate QHC at T1 # 
run_setup $n $nn $atom_mass
MD++ srand48bytime NHMass=\[ 1e-2 1e-2 1e-2 1e-2 \] totalsteps=100000
MD++ ensemble_type="NVTC" integrator_type="Gear6" NHChainLen=4
MD++ T_OBJ=$T initvelocity saveprop=0 
MD++ Ecoh=0 switchfreq=10 ecspring=0
MD++ lambda0=1 lambda1=1 refpotential=1
MD++ runMDSWITCH
MD++ finalcnfile="${T}K_equil_ha_crystal.cn" writeall=1 writecn
MD++ incnfile="${T}K_equil_ha_crystal.cn"

# switching setup
MD++ Ecoh=[ get_data Ecoh ]
MD++ switchfreq=10 saveprop=1 savepropfreq=100 openpropfile printfreq=5000
MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33"


for { set x 0 } {$x<=$iter1} {incr x 1} {
	set y [format %04d $x]
	MD++ readcn totalsteps=50000 saveprop=0

# equilibration run
	MD++ lambda0=1 lambda1=1 refpotential=1
	if { $iter1>0 } {
	MD++ initvelocity runMDSWITCH
	}
# switching run
	MD++ lambda0=1 lambda1=0 refpotential=1
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH
	set datafile "MEAM_F$n.dat"
	set fileID [open $datafile a+]
	puts $fileID "$x [MD++_Get Wavg] [MD++_Get Wtot]"
	close $fileID
	MD++ finalcnfile=f$n$y.cn writeall=1 writecn
}
MD++ quit
} elseif { $n==5 } {
# step 5 is reversible scaling from T1 to T1/solid_factor_RS #
MD++ incnfile="${T}K_equil_crystal.cn" readcn setconfig1

run_setup $n $nn $atom_mass
MD++ srand48bytime NHMass=\[ 1e-2 1e-2 1e-2 1e-2 \] 
MD++ ensemble_type="NPTC" integrator_type="Gear6" NHChainLen=4
MD++ T_OBJ=$T initvelocity

# switching setup

MD++ switchfreq=10 saveprop=1 savepropfreq=100 openpropfile printfreq=5000
MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33"


for { set x 0 } {$x<=$iter1} {incr x 1} {
	set y [format %04d $x]

# equilibration run
	if { $iter1 > 0 } {
	MD++ readcn totalsteps=50000 saveprop=0 initvelocity run
	}

# switching run
	MD++ lambda0=1 lambda1=$factor switchfunc=1 refpotential=4
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH
	set datafile "MEAM_F$n.dat"
	set fileID [open $datafile a+]
	puts $fileID "$x [MD++_Get Wavg] [MD++_Get Wtot]"
	close $fileID
	MD++ finalcnfile=f$n$y.cn writeall=1 writecn
}
# save every independent switching work into MEAM_F5.out file #
exec cat MEAM_F${n}_${nn}.out >> MEAM_F${n}.out
MD++ quit
} elseif { $n==6 } {
# step 6 is reversible scaling from T1/solid_factor_RS to T1 #
set T1 [expr $T/$factor]

MD++ incnfile="${T}K_equil_crystal.cn" readcn setconfig1

# equilibrate with potential scaled by solid_factor_RS #
run_setup $n $nn $atom_mass
MD++ srand48bytime NHMass=\[ 1e-2 1e-2 1e-2 1e-2 \] totalsteps=100000
MD++ ensemble_type="NPTC" integrator_type="Gear6" NHChainLen=4
MD++ T_OBJ=$T initvelocity saveprop=0
MD++ lambda0=$factor lambda1=$factor refpotential=4 switchfunc=1
MD++ runMDSWITCH finalcnfile="${T1}K_equil_crystal.cn" writeall=1 writecn
MD++ incnfile="${T1}K_equil_crystal.cn"

# switching setup

MD++ switchfreq=10 saveprop=1 savepropfreq=100 openpropfile printfreq=5000
MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33"

for { set x 0 } {$x<=$iter1} {incr x 1} {
	set y [format %04d $x]

# equilibration run
	MD++ readcn totalsteps=50000 saveprop=0
	MD++ lambda0=$factor lambda1=$factor refpotential=4 switchfunc=1
	if { $iter1> 0 } {
	MD++ initvelocity runMDSWITCH
	}

# switching run
	MD++ lambda0=$factor lambda1=1 switchfunc=1 refpotential=4
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH

	set datafile "MEAM_F$n.dat"
	set fileID [open $datafile a+]	
	puts $fileID "$x [MD++_Get Wavg] [MD++_Get Wtot]"
	close $fileID
	MD++ finalcnfile=f$n$y.cn writeall=1 writecn
}
# save every indepedent switching work into MEAM_F6.out #
exec cat MEAM_F${n}_${nn}.out >> MEAM_F${n}.out
MD++ quit

} elseif { $n==7 } {
# equilibrate liquid at T2 #
# read 0K solid properties from info1.dat #
exec sleep 10
set fp [ open "../${material}_Solid/info1.dat" r ]
set data [ split [ read $fp ] ]
set a0 [ lindex $data 0 ]
close $fp

# using lattice constant from info1.dat, #
# construct a supercell with randomized atoms #
create_liquid $material $CS $a0 $L
CGR_fixbox

MD++ finalcnfile="0K_liquid.cn" writeall=1 writecn

# equilibrate the liquid at T2 #
run_setup $n $nn $atom_mass
MD++ T_OBJ=$T
MD++ initvelocity ensemble_type="NPT" integrator_type="Gear6"
MD++ totalsteps=1000000 run
MD++ finalcnfile="${T}K_equil_liquid_eT.cn" writecn
MD++ eval
set Ecoh [MD++_Get Ecoh]
set NP [MD++_Get NP]
set V [ datafile_process "MEAM_F7_${nn}.out" 11 0.5 1 "AVE" ]
set a [ expr pow($V,1.0/3)/$L]
set eT [ expr $a/$a0-1 ]

# save equilibrated liquid properties into info2.dat #
set datafile "info2.dat"
set fp [ open $datafile w ]
puts $fp "$a0 $L $eT $Ecoh $V $T $NP"
close $fp

# activate step 8 and step 9 # 
make_flag 8
make_flag 9
MD++ quit
  
} elseif { $n==8 } {
# step 8 is for 'switching from real liquid to gaussian fluid #
create_liquid $material $CS [get_data a0] [get_data L]
scale_H [get_data eT]
CGR_fixbox

# equilbrating run of real liquid
run_setup $n $nn $atom_mass
MD++ srand48bytime totalsteps=200000
MD++ ensemble_type="NVT" integrator_type="VVerlet" 
MD++ T_OBJ=$T initvelocity saveprop=0
MD++ run finalcnfile="${T}K_equil_liquid.cn" writeall=1 writecn setconfig1
MD++ incnfile="${T}K_equil_liquid.cn"

# activate step 12 and 13 for reversible scaling steps #
make_flag 12
make_flag 13
# switch setup
MD++ switchfreq=10 saveprop=1 openpropfile printfreq=5000
MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33 TSTRESS_xx TSTRESS_yy TSTRESS_zz"
MD++ Gauss_rc=3.771 Gauss_sigma=0.7 Gauss_epsilon=50


for { set x 0 } {$x<=$iter1} {incr x 1} {
	set y [format %04d $x]
	MD++ readcn totalsteps=100000 saveprop=0
	# equilibration run
	if { $iter1 > 0 } {
	MD++ initvelocity run
	}
	# switching run
	MD++ lambda0=0 lambda1=1 refpotential=5
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH

	set datafile "MEAM_F$n.dat"
	set fileID [ open $datafile a+ ]
	puts $fileID "$x [MD++_Get Wavg] [MD++_Get Wtot]"
	close $fileID
	
	MD++ finalcnfile=f$n$y.cn writecn
}


MD++ quit
} elseif { $n==9 } {
# step 9 is for 'switching from gaussian fluid to real liquid #
create_liquid $material $CS [get_data a0] [get_data L]
scale_H [get_data eT]
CGR_fixbox

#MD++ Ecoh=-4.3335
# equilbrating run of gaussian fluid
run_setup $n $nn $atom_mass
MD++ srand48bytime totalsteps=200000 setconfig1
MD++ ensemble_type="NVT" integrator_type="VVerlet" 
MD++ T_OBJ=$T initvelocity saveprop=0
MD++ switchfreq=10
MD++ lambda0=1 lambda1=1 refpotential=5
# gaussian potential parameter setup #
MD++ Gauss_rc=3.771 Gauss_sigma=0.7 Gauss_epsilon=50
MD++ runMDSWITCH finalcnfile="${T}K_equil_ga_liquid.cn" writeall=1 writecn
MD++ incnfile="${T}K_equil_ga_liquid.cn"
# activate step 10 and 11 for Gaussian fluid -ideal gas run #
make_flag 10
make_flag 11

# switch setup
MD++ switchfreq=10 saveprop=1 openpropfile printfreq=5000
MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33 TSTRESS_xx TSTRESS_yy TSTRESS_zz"


for { set x 0 } {$x<=$iter1} {incr x 1} {
	set y [format %04d $x]
	MD++ readcn totalsteps=100000 saveprop=0
	MD++ lambda0=1 lambda1=1 refpotential=5
	# equilibration run
	if { $iter1 > 0 } {
	MD++ initvelocity runMDSWITCH
	}

	# switching run
	MD++ lambda0=1 lambda1=0 refpotential=5
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH

	set datafile "MEAM_F$n.dat"
	set fileID [ open $datafile a+ ]
	puts $fileID "$x [MD++_Get Wavg] [MD++_Get Wtot]"
	close $fileID
	
	MD++ finalcnfile=f$n$y.cn writecn
}

MD++ quit
} elseif { $n==10 } {
# step 10 is for 'switching from gaussian fluid to ideal gas #
set factor 0.1

MD++ incnfile="${T}K_equil_ga_liquid.cn" readcn setconfig1

# switch setup
run_setup $n $nn $atom_mass
MD++ srand48bytime totalsteps=200000
MD++ ensemble_type="NVT" integrator_type="VVerlet" 
MD++ T_OBJ=$T initvelocity saveprop=0
MD++ switchfreq=10 saveprop=1 openpropfile printfreq=5000
MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33 TSTRESS_xx TSTRESS_yy TSTRESS_zz"
MD++ Gauss_rc=3.771 Gauss_sigma=0.7 Gauss_epsilon=50


for { set x 0 } {$x<=$iter2} {incr x 1} {
	set y [format %04d $x]
	MD++ readcn totalsteps=100000 saveprop=0
	MD++ lambda0=1 lambda1=1 refpotential = 5 switchfunc=1

	# equilibration run of gaussian fluid
	if { $iter2 > 0 } {
	MD++ initvelocity runMDSWITCH
	}

	# switching run
	MD++ lambda0=1 lambda1=$factor refpotential=6
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH
	set W1 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,1)] lambda1=[expr pow($factor,2)]
	MD++ runMDSWITCH
	set W2 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,2)] lambda1=[expr pow($factor,3)]
	MD++ runMDSWITCH
	set W3 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,3)] lambda1=[expr pow($factor,4)]
	MD++ runMDSWITCH
	set W4 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,4)] lambda1=[expr pow($factor,5)]
	MD++ runMDSWITCH
	set W5 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,5)] lambda1=[expr pow($factor,6)]
#	MD++ lambda0=[expr pow($factor,6)] lambda1=0
	MD++ runMDSWITCH
	set W6 [ MD++_Get Wavg ]

	set datafile "MEAM_F$n.dat"
	set fileID [ open $datafile a+ ]
	puts $fileID "$x $W1 $W2 $W3 $W4 $W5 $W6"
	close $fileID
	
	MD++ finalcnfile=f$n$y.cn writecn
}

MD++ quit

} elseif { $n==11 } {
# step 11 for 'switching from idea gas to gaussian fluid #
set factor 0.1
set T1 [expr $T/pow($factor,6)]

MD++ incnfile="${T}K_equil_ga_liquid.cn" readcn setconfig1

# switch setup
run_setup $n $nn $atom_mass
MD++ srand48bytime totalsteps=200000
MD++ ensemble_type="NVT" integrator_type="VVerlet" 
MD++ T_OBJ=$T initvelocity saveprop=0
MD++ switchfreq=10 saveprop=0 openpropfile printfreq=5000
MD++ lambda0=[expr pow($factor,6)] lambda1=[expr pow($factor,6)]
MD++ refpotential=6 switchfunc=1
MD++ Gauss_rc=3.771 Gauss_sigma=0.7 Gauss_epsilon=50
MD++ runMDSWITCH finalcnfile="${T1}K_equil_ga_liquid.cn" writecn
MD++ incnfile="${T1}K_equil_ga_liquid.cn"

MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33 TSTRESS_xx TSTRESS_yy TSTRESS_zz"



for { set x 0 } {$x<=$iter2} {incr x 1} {
	set y [format %04d $x]
	MD++ readcn totalsteps=100000 saveprop=0
	MD++ lambda0=[expr pow($factor,6)] lambda1=[expr pow($factor,6)]
	MD++ refpotential = 6 switchfunc=1

	# equilibration run
	if { $iter2>0 } {
	MD++ initvelocity runMDSWITCH
	}

	# switching run
	MD++ lambda0=[expr pow($factor,6)] lambda1=[expr pow($factor,5)]
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH
	set W6 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,5)] lambda1=[expr pow($factor,4)]
	MD++ runMDSWITCH
	set W5 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,4)] lambda1=[expr pow($factor,3)]
	MD++ runMDSWITCH
	set W4 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,3)] lambda1=[expr pow($factor,2)]
	MD++ runMDSWITCH
	set W3 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,2)] lambda1=[expr pow($factor,1)]
	MD++ runMDSWITCH
	set W2 [ MD++_Get Wavg ]

        MD++ lambda0=[expr pow($factor,1)] lambda1=1
#	MD++ lambda0=[expr pow($factor,6)] lambda1=0
	MD++ runMDSWITCH
	set W1 [ MD++_Get Wavg ]

	set datafile "MEAM_F$n.dat"
	set fileID [ open $datafile a+ ]
	puts $fileID "$x $W6 $W5 $W4 $W3 $W2 $W1"
	close $fileID
	
	MD++ finalcnfile=f$n$y.cn writecn
}

MD++ quit
} elseif { $n==12 } { 
# step 12 is for reversible scaling from T to T/liquid_factor_RS #
MD++ incnfile="${T}K_equil_liquid.cn" readcn setconfig1

# switch setup
run_setup $n $nn $atom_mass
MD++ srand48bytime totalsteps=200000
MD++ ensemble_type="NPT" integrator_type="Gear6" 
MD++ T_OBJ=$T initvelocity saveprop=0
MD++ switchfreq=10 saveprop=1 openpropfile printfreq=5000 switchfunc=1
MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33 TSTRESS_xx TSTRESS_yy TSTRESS_zz"

for { set x 0 } {$x<=$iter1} {incr x 1} {
	set y [format %04d $x]
	MD++ readcn totalsteps=100000 saveprop=0
	# equilibration run
	if { $iter1>0 } {
	MD++ initvelocity run
	}

	# switching run
	MD++ lambda0=1 lambda1=$factor refpotential=4
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH
	set W [ MD++_Get Wavg ]

	set datafile "MEAM_F$n.dat"
	set fileID [ open $datafile a+ ]
	puts $fileID "$x $W"
	close $fileID	
	
	MD++ finalcnfile=f$n$y.cn writecn
}

exec cat MEAM_F${n}_${nn}.out >> MEAM_F${n}.out
MD++ quit
} elseif { $n==13 } {
# step 13 is for reversible scaling form T/liquid_factor_RS to T #

set T1 [expr $T/$factor]

MD++ incnfile="${T}K_equil_liquid.cn" readcn setconfig1

 # equilibration 
run_setup $n $nn $atom_mass
MD++ srand48bytime totalsteps=200000
MD++ ensemble_type="NPT" integrator_type="Gear6" 
MD++ T_OBJ=$T initvelocity saveprop=0
MD++ switchfreq=10 printfreq=5000 refpotential=4 switchfunc=1
MD++ lambda0=$factor lambda1=$factor runMDSWITCH
MD++ finalcnfile="${T1}K_equil_liquid.cn" writeall=1 writecn
MD++ incnfile="${T1}K_equil_liquid.cn"
MD++ output_fmt="curstep EPOT KATOM HELMP Tinst dEdlambda Wtot Wavg H_11 H_22 H_33 TSTRESS_xx TSTRESS_yy TSTRESS_zz"


for { set x 0 } {$x<=$iter1} {incr x 1} {
	set y [format %04d $x]
	MD++ readcn totalsteps=100000 saveprop=0
	# equilibration run
	MD++ lambda0=$factor lambda1=$factor
	if { $iter1>0 } {
	MD++ initvelocity runMDSWITCH
	}
	# switching run
	MD++ lambda0=$factor lambda1=1 refpotential=4
	MD++ totalsteps=1000000 saveprop=1 runMDSWITCH
	set W [ MD++_Get Wavg ]

	set datafile "MEAM_F$n.dat"
	set fileID [ open $datafile a+ ]
	puts $fileID "$x $W "
	close $fileID
	
	MD++ finalcnfile=f$n$y.cn writecn
}

exec cat MEAM_F${n}_${nn}.out >> MEAM_F${n}.out
MD++ quit

} elseif { $n==14 } {
# step 14 is for making Hessian_XX.m file computing QHC free energy #
        set fp [ open "Hessian_${material}.m" w ]
   	puts $fp "infile = '${material}_Solid/hessian${T_solid_0}.out';"
	puts $fp "T=${T_solid_0};"
        puts $fp "temp=load('${material}_Solid/info2.dat');"
        puts $fp "Vol=temp(5);"
        puts $fp "EPOT=temp(4);"
        puts $fp "NP=temp(7);"
        puts $fp "atomic_mass=${atom_mass}"
        puts $fp "kB=8.617343*10^(-5); % eV/K"
        puts $fp "hbar=6.58211899*10^(-16); %eV s"
        puts $fp "Na=6.02214179*10^23;"
        puts $fp "data = load(infile);"
        puts $fp "data=reshape(-data',NP*3,NP*3);"
        puts $fp "data=(data+data')/2;"
        puts $fp "D1=eig(data);"
        puts $fp "w=sqrt(D1*16.022*Na/(atomic_mass*0.001));"
        puts $fp "w=w(4:end);"
        puts $fp "F1=-kB*T*sum(log(kB*T/hbar./w));"
        puts $fp "gam=sqrt(6.626068^2*6.022/(1.3806503*atomic_mass*0.001*T*2*pi))*10^(-1);"
	puts $fp "%F_ideal=-kB*T*log(Vol/gam^3);"
 	puts $fp "F=F1;%+F_ideal;"
        puts $fp "F_per=F/NP;"
        puts $fp "disp(sprintf('F=%f (eV) dFhar= %f (eV) F_per=%f(eV)',F+EPOT*NP,F,F_per));"
        puts $fp "fid=fopen('${material}_Solid_FreeE','w');"
  	puts $fp "fprintf(fid,'%15.8f',F+EPOT*NP);"
        puts $fp "fclose(fid);"
	puts $fp "fid=fopen('${material}_Intermediate_Values','w');"
	puts $fp "fprintf(fid,'eT_solid at ${T_solid_0} = %15.8f\\n',temp(3));"
	puts $fp "fprintf(fid,'EO for %d = %15.8f \\nF = %15.8f\\n',NP,EPOT*NP,F);"
	puts $fp "fclose(fid);" 
	close $fp
	MD++ quit
} elseif { $n==15 } {
# step 15 is for Free_Energy_of_XX.m file computing melting point #
     	set fp [ open "Free_Energy_of_${material}.m" w]
	puts $fp " data3=load('${material}_Solid/MEAM_F3.dat');"
        puts $fp " data4=load('${material}_Solid/MEAM_F4.dat');"
        puts $fp " data5=load('${material}_Solid/MEAM_F5.out');"
        puts $fp " data6=load('${material}_Solid/MEAM_F6.out');"
        puts $fp " data8=load('${material}_Liquid/MEAM_F8.dat');"
        puts $fp " data9=load('${material}_Liquid/MEAM_F9.dat');"
        puts $fp "data10=load('${material}_Liquid/MEAM_F10.dat');"
        puts $fp "data11=load('${material}_Liquid/MEAM_F11.dat');"
        puts $fp "data12=load('${material}_Liquid/MEAM_F12.out');"
        puts $fp "data13=load('${material}_Liquid/MEAM_F13.out');"
        puts $fp "kB=8.616343*10^(-5); %eV/K"
        puts $fp "temp=load('${material}_Liquid/info2.dat');"
      	puts $fp "NP=temp(7);"
        puts $fp "Vliq=temp(5);"
        puts $fp "Tliq=${T_liquid_0};"
	puts $fp "atomic_mass=${atom_mass}"
        puts $fp "Fha_${T_solid_0}=load('${material}_Solid_FreeE');"
if { $iter1 >0 } {
	puts $fp "free3=mean(data3(:,2));"
	puts $fp "free3_e=std(data3(:,2));" 
	puts $fp "free4=mean(data4(:,2));"
	puts $fp "free4_e=std(data4(:,2));"
} else {
	puts $fp "free3=(data3(:,2));"
	puts $fp "free3_e=std(data3(:,2));" 
	puts $fp "free4=(data4(:,2));"
	puts $fp "free4_e=std(data4(:,2));"
}
	puts $fp "dFre_Fha=(free4-free3)/2;"
	puts $fp "dFre_Fha_e=(free3_e+free4_e)/2;"
 	puts $fp "Fre_${T_solid_0}=Fha_${T_solid_0}+dFre_Fha;"
        set d [ expr $factor_solid_RS-1.0 ]
	puts $fp "lambda=1:${d}/10000:${factor_solid_RS};"
	puts $fp "iter=floor(length(data5(:,8))/10001);"
	puts $fp "temp1=reshape(data5(1:10001*iter,8),10001,iter);"
	puts $fp "iter=floor(length(data6(:,8))/10001);"
	puts $fp "temp2=reshape(data6(1:10001*iter,8),10001,iter);"
if {$iter1>0} {
        puts $fp "w1=mean(temp1');"
	puts $fp "w2=mean(temp2');"
} else {
        puts $fp "w1=(temp1');"
	puts $fp "w2=(temp2');"
}
	puts $fp "for i=1:10001"
	puts $fp "w1_inv(i)=w2(10001)-w2(10002-i);"
	puts $fp "end"
	puts $fp "for i=1:10001"
	puts $fp "w(i)=(w1(i)-w1_inv(i))/2;"
	puts $fp "end"
	puts $fp "T=zeros(1,10001);"
	puts $fp "T=${T_solid_0}./lambda;"
	puts $fp "Fs=zeros(1,10001);"
	puts $fp "for i=1:10001"
	puts $fp "Fs(i)=((Fre_${T_solid_0}+w(i))/${T_solid_0}-3/2*(NP-1)*kB*log(T(i)/${T_solid_0}))*T(i);"
	puts $fp "end"
	puts $fp "% Liquid Free E"
        set fp1 [ open "logfactorial.m" w]
	puts $fp1 "function y = logfactorial (x)"
	puts $fp1 "y=(x+1/2)*log(x)-x+1/2*log(2*pi);"
	close $fp1
	puts $fp "logfac=logfactorial(NP-1);"
	puts $fp "gam=sqrt(6.626068^2*6.022/(1.3806503*atomic_mass*0.001*Tliq*2*pi))*10^(-1);"
	puts $fp "Fid=-kB*Tliq*((NP-1)*log(Vliq/gam^3)-logfac);"
if {$iter1>0} {
	puts $fp "free8=mean(data8(:,2));"
	puts $fp "free9=mean(data9(:,2));"
} else {
	puts $fp "free8=(data8(:,2));"
	puts $fp "free9=(data9(:,2));"
}
	puts $fp "dFsw_Fga=(free9-free8)/2;"
if {$iter2>0 } {
	puts $fp "free10=mean(sum(data10(:,2:7)'));"
	puts $fp "free11=mean(sum(data11(:,2:7)'));"
} else {
	puts $fp "free10=(sum(data10(:,2:7)'));"
	puts $fp "free11=(sum(data11(:,2:7)'));"
}
	puts $fp "dFga_Fid=(free11-free10)/2;"
	puts $fp "Fsw_${T_liquid_0}=Fid+dFga_Fid+dFsw_Fga;"
	set d [ expr ${factor_liquid_RS}-1.0]
	puts $fp "lambda1=1:${d}/10000:${factor_liquid_RS};"
	puts $fp "iter=floor(length(data12(:,8))/10001);"
	puts $fp "temp3=reshape(data12(1:10001*iter,8),10001,iter);"
	puts $fp "iter=floor(length(data13(:,8))/10001);"
	puts $fp "temp4=reshape(data13(1:10001*iter,8),10001,iter);"
if { $iter1>0 } {      
        puts $fp "w3=mean(temp3');"
	puts $fp "w4=mean(temp4');"
} else {
        puts $fp "w3=(temp3');"
	puts $fp "w4=(temp4');"
}
	puts $fp "for i=1:10001"
	puts $fp "w3_inv(i)=w4(10001)-w4(10002-i);"
	puts $fp "end"
	puts $fp "for i=1:10001"
	puts $fp "ww(i)=(w3(i)-w3_inv(i))/2;"
	puts $fp "end"
	puts $fp "T1=zeros(1,10001);"
	puts $fp "T1=${T_liquid_0}./lambda1;"
	puts $fp "Fl=zeros(1,10001);"
	puts $fp "for i=1:10001"
	puts $fp "Fl(i)=((Fsw_${T_liquid_0}+ww(i))/${T_liquid_0}-1.5*(NP-1)*kB*log(T1(i)/${T_liquid_0}))*T1(i);"
	puts $fp "end"
	puts $fp "figure(1)"
	puts $fp "plot(T,Fs,T1,Fl)"
	puts $fp "temp=load('${material}_Liquid/info2.dat');"
	puts $fp "fid=fopen('${material}_Intermediate_Values','a');"
	puts $fp "fprintf(fid,'Fso-Fha at ${T_solid_0} = %15.8f\\n',dFre_Fha);"
	puts $fp "fprintf(fid,'Fso(%15.8f) = %15.8f, Fso(%15.8f)=%15.8f\\n',T(1),Fs(1),T(end),Fs(end));"
	puts $fp "fprintf(fid,'eT_liquid at ${T_liquid_0} = %15.8f\\n',temp(3));"
	puts $fp "fprintf(fid,'Fideal at ${T_liquid_0} = %15.8f\\n',Fid);"
	puts $fp "fprintf(fid,'Fga-Fid at ${T_liquid_0} = %15.8f\\n',dFga_Fid);"
	puts $fp "fprintf(fid,'Fli-Fga at ${T_liquid_0} = %15.8f\\n',dFsw_Fga);"
	puts $fp "fprintf(fid,'Fli(%15.8f) = %15.8f, Fli(%15.8f) = %15.8f\\n',T1(1),Fl(1),T1(end),Fl(end));"
	puts $fp "fclose(fid);"

	puts $fp "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
	puts $fp "%%% For publication data process"
	puts $fp "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"

	puts $fp "N3=length(data3(:,2));"
	puts $fp "N4=length(data4(:,2));"
	puts $fp "N5=length(temp1(1,:));"
	puts $fp "N6=length(temp2(1,:));"
	puts $fp "N8=length(data8(:,2));"
	puts $fp "N9=length(data9(:,2));"
	puts $fp "N10=length(data10(:,2));"
	puts $fp "N11=length(data11(:,2));"
	puts $fp "N12=length(temp3(1,:));"
	puts $fp "N13=length(temp4(1,:));"


	puts $fp "temp_sol=load('${material}_Solid/info2.dat');"
	puts $fp "a_0=temp_sol(1);"
	puts $fp "e_s=temp_sol(3);"
	puts $fp "e_l=temp(3);"
	puts $fp "E_coh=temp_sol(4)*NP;"
	puts $fp "E_0=temp_sol(4)*NP;"
	puts $fp "F_ha=load('${material}_Solid_FreeE');"
	puts $fp "F_w=F_ha-E_0;"
	puts $fp "F2_fow=free3;"
	puts $fp "F2_fow_std=free3_e/sqrt(N3);"
	puts $fp "F2_rev=free4;"
	puts $fp "F2_rev_std=free4_e/sqrt(N4);"
	puts $fp "dF2=(F2_rev-F2_fow)/2;"

	puts $fp "temp5=load('${material}_Solid/MEAM_F5.dat');"
	puts $fp "temp6=load('${material}_Solid/MEAM_F6.dat');"

	puts $fp "F3_fow=mean(temp5(:,2));"
	puts $fp "F3_fow_std=std(temp5(:,2))/sqrt(N5);"
	puts $fp "F3_rev=mean(temp6(:,2));"
	puts $fp "F3_rev_std=std(temp6(:,2))/sqrt(N6);"

	puts $fp "F5_1_fow=free8;"
	puts $fp "F5_1_fow_std=std(data8(:,2))/sqrt(N8);"
	puts $fp "F5_1_rev=free9;"
	puts $fp "F5_1_rev_std=std(data9(:,2))/sqrt(N9);"

	puts $fp "F5_2_fow=free10;"
	puts $fp "F5_2_fow_std=std(sum(data10(:,2:7)'))/sqrt(N10);"
	puts $fp "F5_2_rev=free11;"
	puts $fp "F5_2_rev_std=std(sum(data11(:,2:7)'))/sqrt(N11);"

	puts $fp "temp12=load('${material}_Liquid/MEAM_F12.dat');"
	puts $fp "temp13=load('${material}_Liquid/MEAM_F13.dat');"

	puts $fp "F6_fow=mean(temp12(:,2));"
	puts $fp "F6_fow_std=std(temp12(:,2))/sqrt(N12);"
	puts $fp "F6_rev=mean(temp13(:,2));"
	puts $fp "F6_rev_std=std(temp13(:,2))/sqrt(N13);"

	puts $fp "fid=fopen('${material}_Publication_Info.dat','w');"
	puts $fp "fprintf(fid,'Solid Free Energy Info\\n');"
	puts $fp "fprintf(fid,'lattice const = %10.5f\\n',a_0);"
	puts $fp "fprintf(fid,'strain = %10.5f\\n',e_s);"
	puts $fp "fprintf(fid,'E_0 = %10.5f eV ( %10.5f eV/atom)\\n',E_0, E_0/NP);"
	puts $fp "fprintf(fid,'F_w = %10.5f eV ( %10.5f eV/atom)\\n',F_w,F_w/NP);"
	puts $fp "fprintf(fid,'F_ha = E_0 + F_w = %10.5f eV ( %10.5f eV/atom)\\n',F_ha,F_ha/NP);"
	puts $fp "fprintf(fid,'W_re_to_ha = %10.5f eV ( %10.5f eV/atom)\\n',F2_fow,F2_fow/NP);"
	puts $fp "fprintf(fid,'W_ha_to_re = %10.5f eV ( %10.5f eV/atom)\\n',F2_rev,F2_rev/NP);"
	puts $fp "fprintf(fid,'F_re - F_ha = (W_ha_to_re - W_re_to_ha)/2 = %10.5f eV ( %10.5f eV/atom)\\n',dFre_Fha,dFre_Fha/NP);"
	puts $fp "fprintf(fid,'F_re at %f = %10.5f eV ( %10.5f eV/atom)\\n',T(1),Fs(1),Fs(1)/NP);"
	puts $fp "fprintf(fid,'W_T1_to_T11 = %10.5f eV ( %10.5f eV/atom)\\n',F3_fow,F3_fow/NP);"
	puts $fp "fprintf(fid,'W_T11_to_T1 = %10.5f eV ( %10.5f eV/atom)\\n',F3_rev,F3_rev/NP);"
	puts $fp "fprintf(fid,'dW1 = (W_T1_to_T11 - W_T11_to_T1)/2 = %10.5f eV ( %10.5f eV/atom)\\n',(F3_fow-F3_rev)/2,(F3_fow-F3_rev)/2/NP);"
	puts $fp "fprintf(fid,'F_re at %f = %10.5f eV ( %10.5f eV/atom)\\n',T(end),Fs(end),Fs(end)/NP);"
	puts $fp "fprintf(fid,'\\n');"
	puts $fp "fprintf(fid,'Liquid Free Energy Info\\n');"
	puts $fp "fprintf(fid,'strain = %10.5f\\n',e_l);"
	puts $fp "fprintf(fid,'F_ideal = %10.5f ( %10.5f eV/atom)\\n',Fid, Fid/NP);"
	puts $fp "fprintf(fid,'W_re_to_ga = %10.5f eV ( %10.5f eV/atom)\\n',F5_1_fow,F5_1_fow/NP);"
	puts $fp "fprintf(fid,'W_ga_to_re = %10.5f eV ( %10.5f eV/atom)\\n',F5_1_rev,F5_1_rev/NP);"
	puts $fp "fprintf(fid,'F_re - F_ga = (W_ga_to_re - W_re_to_ga)/2 = %10.5f eV ( %10.5f eV/atom)\\n',dFsw_Fga,dFsw_Fga/NP);"
	puts $fp "fprintf(fid,'W_ga_to_id = %10.5f eV ( %10.5f eV/atom)\\n',F5_2_fow,F5_2_fow/NP);"
	puts $fp "fprintf(fid,'W_id_to_ga = %10.5f eV ( %10.5f eV/atom)\\n',F5_2_rev,F5_2_rev/NP);"
	puts $fp "fprintf(fid,'F_ga - F_id = (W_id_to_ga - W_ga_to_id)/2 = %10.5f eV ( %10.5f eV/atom)\\n',dFga_Fid,dFga_Fid/NP);"
	puts $fp "fprintf(fid,'F_re at %f = %10.5f eV ( %10.5f eV/atom)\\n',T1(1),Fl(1),Fl(1)/NP);"
	puts $fp "fprintf(fid,'W_T2_to_T22 = %10.5f eV ( %10.5f eV/atom)\\n',F6_fow,F6_fow/NP);"
	puts $fp "fprintf(fid,'W_T22_to_T2 = %10.5f eV ( %10.5f eV/atom)\\n',F6_rev,F6_rev/NP);"
	puts $fp "fprintf(fid,'dW2 = (W_T2_to_T22 - W_T22_to_T2)/2 = %10.5f eV ( %10.5f eV/atom)\\n',(F6_fow-F6_rev)/2,(F6_fow-F6_rev)/2/NP);"
	puts $fp "fprintf(fid,'F_re at %f = %10.5f eV ( %10.5f eV/atom)\\n',T1(end),Fl(end),Fl(end)/NP);"
	puts $fp "fprintf(fid,'\\n');"
	puts $fp "fprintf(fid,'Error Analysis\\n\\n');"
	puts $fp "fprintf(fid,'Solid Side\\n');"
	puts $fp "fprintf(fid,'Err_re_to_ha = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F2_fow_std,F2_fow_std/NP,F2_fow_std/NP/kB);"
	puts $fp "fprintf(fid,'Err_ha_to_re = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F2_rev_std,F2_rev_std/NP,F2_rev_std/NP/kB);"
	puts $fp "fprintf(fid,'Err_T1_to_T11 = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F3_fow_std,F3_fow_std/NP,F3_fow_std/NP/kB);"
	puts $fp "fprintf(fid,'Err_T11_to_T1 = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F3_rev_std,F3_rev_std/NP,F3_rev_std/NP/kB);"
	puts $fp "fprintf(fid,'\\n');"
	puts $fp "fprintf(fid,'Liquid Side\\n');"
	puts $fp "fprintf(fid,'Err_re_to_ga = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F5_1_fow_std,F5_1_fow_std/NP,F5_1_fow_std/NP/kB);"
	puts $fp "fprintf(fid,'Err_ga_to_re = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F5_1_rev_std,F5_1_rev_std/NP,F5_1_rev_std/NP/kB);"
	puts $fp "fprintf(fid,'Err_ga_to_id = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F5_2_fow_std,F5_2_fow_std/NP,F5_2_fow_std/NP/kB);"
	puts $fp "fprintf(fid,'Err_id_to_ga = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F5_2_rev_std,F5_2_rev_std/NP,F5_2_rev_std/NP/kB);"
	puts $fp "fprintf(fid,'Err_T2_to_T22 = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F6_fow_std,F6_fow_std/NP,F6_fow_std/NP/kB);"
	puts $fp "fprintf(fid,'Err_T22_to_T2 = %10.5f eV ( %10.8f eV/atom)(%6.3f K)\\n',F6_rev_std,F6_rev_std/NP,F6_rev_std/NP/kB);"
	puts $fp "dT=sqrt(F2_fow_std^2+F2_rev_std^2+F3_fow_std^2+F3_rev_std^2+F5_1_fow_std^2+F5_1_rev_std^2+F5_2_fow_std^2+F5_2_rev_std^2+F6_fow_std^2+F6_rev_std^2)/NP/kB;"
	puts $fp "fprintf(fid,'Max dT = %6.3f\\n',dT);"
	puts $fp "fclose(fid);"

	puts $fp "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
	puts $fp "%%%% End of Data Process"
	puts $fp "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"


	puts $fp "Tmin=max(T(1),T1(end));"
	puts $fp "Tmax=min(T(end),T1(1));"
	puts $fp "T2=Tmin:(Tmax-Tmin)/10000:Tmax;"
	puts $fp "Fsol=interp1(T,Fs,T2);"
	puts $fp "Fliq=interp1(T1,Fl,T2);"
	puts $fp "\[Fcoin,num\]=min(abs(Fsol-Fliq));"
	puts $fp "Tm=T2(num);"
	puts $fp "Fcoin=Fsol(num);"
	puts $fp "num2=min(num+100,length(T2));"
	puts $fp "num1=max(num-100,1);"
	puts $fp "Ss=-(Fsol(num2)-Fsol(num1))/(T2(num2)-T2(num1));"
	puts $fp "Sl=-(Fliq(num2)-Fliq(num1))/(T2(num2)-T2(num1));"
	puts $fp "L=Tm*(Sl-Ss);"
	puts $fp "ss=Ss/NP;sl=Sl/NP;"
	puts $fp "dt=dT*kB/(sl-ss);"
	puts $fp "disp(sprintf('Tm = %f, Fm = %f',Tm,Fcoin));"
	puts $fp "coef_sol=polyfit(T2(:),Fsol(:),2);"
	puts $fp "Cp_sol=-2*coef_sol(1)*1.6*10^(-19)*6.023*10^23/atomic_mass*Tm;"
	puts $fp "coef_liq=polyfit(T2(:),Fliq(:),2);"
	puts $fp "Cp_liq=-2*coef_liq(1)*1.6*10^(-19)*6.023*10^23/atomic_mass*Tm;"
	puts $fp "fid=fopen('${material}_Publication_Info.dat','a+');"
	puts $fp "fprintf(fid,'Tm= %10.5f (K), err=%6.3f Fm = %10.5f (eV/atom)\\n',Tm,dt,Fcoin/NP);"
	puts $fp "fprintf(fid,'L=%20.10f (eV) Ss=%f Sl=%f\\n',L,Ss,Sl);"
	puts $fp "fprintf(fid,'L=%20.10f (eV/atom) ss=%f sl=%f\\n',L/NP,ss,sl);"
	puts $fp "fprintf(fid,'L=%20.10f (J/mol) ss=%f sl=%f\\n',L/NP*1.6*10^(-19)*6.02*10^23,ss*1.6*10^(-19)*6.02*10^23,sl*1.6*10^(-19)*6.02*10^23);"
	puts $fp "fprintf(fid,'L=%20.10f (J/kg) ss=%f sl=%f\\n',L/NP*1.6*10^(-19)*6.02*10^23*1000/atomic_mass,ss*1.6*10^(-19)*6.02*10^23*1000/atomic_mass,sl*1.6*10^(-19)*6.02*10^23*1000/atomic_mass);"
	puts $fp "fprintf(fid,'Cp_sol=%20.10f (J/gK) Cp_liq=%20.10f (J/gK)\\n',Cp_sol,Cp_liq);"
	puts $fp "fclose(fid);"
	puts $fp "fid=fopen('Free_Energy_${material}.dat','a+');"
	puts $fp "for i=1:length(T2)"
	puts $fp "fprintf(fid,'%10.5f %10.5f %10.5f\\n',T2(i),Fsol(i)/NP,Fliq(i)/NP);"
	puts $fp "end"
	puts $fp "fclose(fid);"
	puts $fp "figure(2)"
	puts $fp "plot(T2,Fsol,T2,Fliq)"
	close $fp
	MD++ quit
} elseif { $n==16 } {
# step 16 generates executable file 'tar_XX' zipping data into inter_XX.tar #
	set fp [ open "tar_${material}" w ]
	puts $fp "tar -cvf inter_${material}.tar \*.m ${material}_Solid/\*.dat ${material}_Solid/hessian${T_solid_0}.out ${material}_Liquid/\*.dat ${material}_Solid/MEAM_F5.out ${material}_Solid/MEAM_F6.out ${material}_Liquid/MEAM_F12.out ${material}_Liquid/MEAM_F13.out"
	close $fp
	exec chmod 755 tar_${material}
}