# -*-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 dislocation and relax structure and box # status = 1 : compute Peierls stress # status = 2 : NPT MD simulation to find equilibrium volume at finite T # status = 3 : NVT MD simulation to compute dislocation mobility # # 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 { n } { #MD++ setnolog MD++ setoverwrite MD++ dirname = runs/ta-edge-$n MD++ NNM = 200 } proc readpot { } { MD++ { #-------------------------------------------- #Read in potential file # potfile = ~/Codes/MD++/potentials/ta_pot readpot } } #-------------------------------------------- proc create_edge_dislocation { Nx Ny Nz a } { # create perfect crystal MD++ crystalstructure = body-centered-cubic # equilibrium lattice constant of Ta is 3.3058 Angstrom MD++ latticeconst = $a MD++ latticesize = \[ 1 1 1 $Nx -1 1 0 $Ny -1 -1 2 $Nz \] MD++ makecrystal #finalcnfile = perf.cn writecn #eval # Remove 1/4 of top and bottom layers of atoms set d_sy [expr 0.5/$Ny] set sy_bottom [expr -.5 + (0.125*2.0*$Ny - 1.0)/(2.0*$Ny) + .5*$d_sy] set sy_top [expr -.5 + (0.875*2.0*$Ny - 1.0)/(2.0*$Ny) + .5*$d_sy] puts "sy_bottom = $sy_bottom, sy_top = $sy_top" MD++ input = \[ 1 -1 1 -1 $sy_bottom -1 1 \] fixatoms_by_position MD++ input = \[ 1 -1 1 $sy_top 1 -1 1 \] fixatoms_by_position MD++ removefixedatoms MD++ vacuumratio = [expr 1.0 - ($sy_top - $sy_bottom)] puts "VACUUM Ratio = [expr 1.0 - ($sy_top - $sy_bottom)]" # dislocation position set x0 [expr 0.111 / [MD++_Get latticesize(3)] ] set y0 [expr -0.5+0.125/[MD++_Get latticesize(7)] ] set y1 [expr $y0 + 0.5] set sbx [expr 0.5 / [MD++_Get latticesize(3)] ] MD++ input = \[ 3 2 $sbx 0 0 $x0 $y0 $y1 0.305 -10 10 -10 10 1 0 0 0 0 \] MD++ makedipole MD++ finalcnfile = dipole.cn writecn } #-------------------------------------------- proc mark_groups { } { MD++ NCS = 8 eval calcentralsymmetry freeallatoms #identify top surface atom MD++ input = \[ 1 -10 10 0.25 10 -10 10 5 40 \] fixatoms_by_pos_topol set ntop [MD++_Get NPfixed] MD++ input = 1 setfixedatomsgroup freeallatoms #identify bottom surface atom MD++ input = \[ 1 -10 10 -10 -0.25 -10 10 5 40 \] fixatoms_by_pos_topol set nbot [MD++_Get NPfixed] MD++ input = 2 setfixedatomsgroup freeallatoms set natoms [list $ntop $nbot] return $natoms } #-------------------------------------------- proc relax_fixbox { } { MD++ { # Conjugate-Gradient relaxation conj_ftol = 1e-4 conj_itmax = 1000 conj_fevalmax = 1000 conj_fixbox = 1 relax } } #end of proc relax_fixbox #-------------------------------------------- proc relax_freebox { } { MD++ { # Conjugate-Gradient relaxation conj_ftol = 1e-4 conj_itmax = 1000 conj_fevalmax = 1000 conj_fixbox = 0 conj_fixboxvec = [ 0 1 1 1 1 1 1 1 0 ] relax } } #end of proc relax_fixbox proc setup_window { } { MD++ { #------------------------------------------------------------ #colors for Central symmetry view color00 = "red" color01 = "blue" color02 = "green" color03 = "magenta" color04 = "cyan" color05 = "purple" color06 = "gray80" color07 = "white" color08 = "orange" #-------------------------------------------- # Plot Configuration # 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 #2.8285 #for Si atomcolor = cyan highlightcolor = purple bondcolor = red fixatomcolor = yellow backgroundcolor = gray70 plot_color_axis = 2 plot_color_windows = [ 4 0.6 6 1 6 10 5 10 20 6 20 50 8 1 ] #plot_atom_info = 1 # reduced coordinates of atoms #plot_atom_info = 2 # real coordinates of atoms plot_atom_info = 3 # energy of atoms #plot_highlight = [ 0 0 1 2 3 4 5 6 7 8 9 ] plotfreq = 10 # rotateangles = [ 0 0 0 1.2 ] #openwin alloccolors rotate saverot plot } } 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 = 100 plotfreq = 100 printfreq = 100 vt2 = 1e28 #1e28 2e28 5e28 wallmass = 2e7 # atommass * NP = 14380 boxdamp = 1e-2 # optimal damping for 216 atoms and wallmass 1e-3 saveH # Use current H as reference (H0), needed for specifying stress fixboxvec = [ 0 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" writeall = 1 } } #end of proc setup_md #******************************************* # 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" if { $argc <= 1 } { set n 0 } elseif { $argc > 1 } { set n [lindex $argv 1] } puts "n = $n" if { $status == 0 } { # prepare dislocation and relax configuration to energy minimum initmd $n readpot create_edge_dislocation 30 40 2 3.3058 openwindow relax_fixbox relax_freebox MD++ eval MD++ finalcnfile = edge_relaxed.cn writecn MD++ finalcnfile = "edge_relaxed.cfg" writeatomeyecfg MD++ sleep exitmd } elseif { $status == 1 } { # read and relax dislocation structure initmd $n readpot MD++ incnfile = "edge_relaxed.cn" readcn # For less atoms stored in cfg files setup_window openwindow MD++ vacuumratio = 0.25 # obtain cross section area set Lx [MD++_Get H_11] set Ly [MD++_Get H_22] set Lz [MD++_Get H_33] set A [expr $Lx*$Lz ] puts "Lx = $Lx Lz = $Lz A = $A" # obtain number of top and bottom surface layer atoms set natoms [ mark_groups ] set ntop [lindex $natoms 0] set nbot [lindex $natoms 1] puts "ntop = $ntop nbot = $nbot" # normalize Fext to correspond to 1 MPa set F0 [expr $A / 160.2e3] set Fx [expr 1.0 * $F0] set minsigma 1 set dsigma 1 set maxsigma 200 set E0 [MD++_Get EPOT] set datafile "EPOT.dat" set fileID [open $datafile w ] puts $fileID "0 $E0 0 0" flush $fileID MD++ conj_ftol = 0.0 conj_dfpred = 1e-3 conj_fixbox = 1 conj_fevalmax = 100000 MD++ applyFext = 1 for { set stress $minsigma } { $stress < $maxsigma } { incr stress $dsigma } { set fx_top [expr $Fx*$stress/$ntop] set fz_top 0 set fx_bot [expr -$Fx*$stress/$nbot] set fz_bot 0 puts "fx_top = $fx_top fz_top = $fz_top fx_bot = $fx_bot fz_bot = $fz_bot" MD++ clearFext MD++ input = \[ 1 $fx_top 0 $fz_top \] addFext_to_group MD++ input = \[ 2 $fx_bot 0 $fz_bot \] addFext_to_group for { set repeat 0 } { $repeat < 3 } { incr repeat 1 } { MD++ relax eval set step [MD++_Get conj_step] # Estimate of dislocation motion distance (dr) set E [MD++_Get EPOT] set dE [expr $E-$E0 ] set dr [expr -0.5*$dE/($F0*$stress)] if { $step >= [MD++_Get conj_fevalmax] } { if { [expr abs($dr)] > 5 } { # If CGR reaches max eval and the disl moves long distance, quit break } } } MD++ finalcnfile = "ta-stress$stress.cn" writeall = 1 writecn MD++ finalcnfile = "ta-stress$stress.cfg" writeatomeyecfg puts $fileID "$stress $E $dE $dr" flush $fileID if { [expr abs($dr)] > 5 } { break } } close $fileID exitmd MD++ eval MD++ sleep exitmd } elseif { $status == 2 } { # NVT MD simulation to equilibrate at finite T initmd $n readpot # replicate the simulation cell along z MD++ allocmultiple = 5 incnfile = "edge_relaxed.cn" readcn MD++ input = \[ 3 5 0 \] extendbox MD++ finalcnfile = "edge_x5_relaxed.cn" writecn # allow thermal expansion set eps_T 0.00197 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 20 for { set iter 1 } { $iter <= $maxiter } { incr iter } { set sig_avg_xx 0 set sig_avg_yy 0 set sig_avg_zz 0 for { set subiter 1 } { $subiter <= 20 } { incr subiter} { MD++ continue_curstep = 1 MD++ totalsteps = 100 run set sig_xx [MD++_Get TSTRESSinMPa_xx] set sig_yy [MD++_Get TSTRESSinMPa_yy] set sig_zz [MD++_Get TSTRESSinMPa_zz] if { $subiter > 10 } { set sig_avg_xx [expr $sig_avg_xx + $sig_xx] set sig_avg_yy [expr $sig_avg_yy + $sig_yy] set sig_avg_zz [expr $sig_avg_zz + $sig_zz] } } set sig_avg_xx [expr $sig_avg_xx / 10] set sig_avg_yy [expr $sig_avg_yy / 10] set sig_avg_zz [expr $sig_avg_zz / 10] puts "sig_avg_xx = $sig_avg_xx sig_avg_yy = $sig_avg_yy sig_avg_zz = $sig_avg_zz" if { $iter < $maxiter } { set nu 0.3 set eps_xx [expr ($sig_avg_xx - $nu*$sig_avg_yy - $nu*$sig_avg_zz)*4e-6 ] set eps_yy [expr ($sig_avg_yy - $nu*$sig_avg_zz - $nu*$sig_avg_xx)*4e-6 ] set eps_zz [expr ($sig_avg_zz - $nu*$sig_avg_xx - $nu*$sig_avg_yy)*4e-6 ] puts "eps_xx = $eps_xx eps_yy = $eps_yy eps_zz = $eps_zz" MD++ input = \[ 1 1 $eps_xx \] changeH_keepS MD++ input = \[ 2 2 $eps_yy \] changeH_keepS MD++ input = \[ 3 3 $eps_zz \] changeH_keepS } } MD++ finalcnfile = "edge_x5_equil.cn" writeall = 1 writecn #MD++ sleep exitmd } elseif { $status == 3 } { # NVT MD simulation to test whether system has equilibrated at zero pressure initmd $n readpot MD++ incnfile = "equil/inter0030.cn" readcn 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 MD++ totalsteps = 10000 run MD++ eval MD++ finalcnfile = "edge_x5_equil.cn" writeall = 1 writecn #MD++ sleep exitmd } elseif { $status == 4 } { # NVT MD simulation to compute dislocation mobility initmd $n readpot set stress $n MD++ incnfile = "../ta-edge-0/edge_x5_equil.cn" readcn 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++ ensemble_type = "NVT" integrator_type = "VVerlet" NHMass = 2e-3 MD++ T_OBJ = 300 set natoms [ mark_groups ] set ntop [lindex $natoms 0] set nbot [lindex $natoms 1] puts "ntop = $ntop nbot = $nbot" #obtain cross section area set Lx [MD++_Get H_11] set Ly [MD++_Get H_22] set Lz [MD++_Get H_33] set A [expr $Lx*$Lz ] puts "Lx = $Lx Lz = $Lz A = $A" set F0 [expr $A / 160.2e3] set Fx [expr 1.0 * $F0] MD++ applyFext = 1 set fx_top [expr $Fx*$stress/$ntop] set fz_top 0 set fx_bot [expr -$Fx*$stress/$nbot] set fz_bot 0 puts "fx_top = $fx_top fz_top = $fz_top fx_bot = $fx_bot fz_bot = $fz_bot" MD++ clearFext MD++ input = \[ 1 $fx_top 0 $fz_top \] addFext_to_group MD++ input = \[ 2 $fx_bot 0 $fz_bot \] addFext_to_group MD++ savecnfreq = 10000 savecfgfreq = 1000 MD++ totalsteps = 1000000 run MD++ finalcnfile = "edge_x5_stress$stress.cn" writeall = 1 writecn #MD++ sleep exitmd } elseif { $status == 10 } { # visualization MD++ setnolog initmd "view" readpot set stress $n MD++ incnfile = "../ta-edge-$n/ta-edge-stress$stress.cn" readcn openwindow MD++ eval MD++ sleep exitmd } else { puts "unknown status = $status" exitmd }