# -*-shell-script-*- # Au nano wire is constructed # make fs build=R SYS=intel #source "scripts/Examples/Tcl/startup.tcl" #******************************************* # Definition of procedures #******************************************* proc initmd { n } { #MD++ setnolog MD++ setoverwrite MD++ dirname = runs/au-nw-torsion-$n MD++ NNM = 200 } proc initmd2 { n dphi } { #MD++ setnolog MD++ setoverwrite MD++ dirname = runs/au-nw-torsion-$n-dphi$dphi MD++ NNM = 200 } #------------------------------------------------------------ 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_orig" ntypes = 1 ename0 = "Au1 " #leave spaces to be compatible with fortran rcut = 6.0 kode0 = "library" readMEAM NNM = 300 } } #------------------------------------------------------------ proc readmeam-lammps { } { MD++ { #Read in MEAM potential (Baskes format) meamfile = "~/Codes/MD++/Fortran/MEAM-Lammps/meamf.lammps" meafile = "~/Codes/MD++/Fortran/MEAM-Lammps/Au.meam" nspecies = 1 element0 = "Au1" rcut = 6.0 readMEAM NNM = 300 } } #------------------------------------------------------------- proc readmeam { } { MD++ { #Read in MEAM potential (MD++ format) potfile = "~/Codes/MD++/potentials/MEAMDATA/meamdata.Au2" readMEAM #potfile = "~/Codes/MD++/potentials/MEAMDATA/meamdata.Siz" readMEAM #potfile = "~/Codes/MD++/potentials/MEAMDATA/meamdata.Si" readMEAM #xzbl = -10 xzbl0 = -9 #turn off correction } } #------------------------------------------------------------- proc readeam { } { MD++ { #Read in EAM potential (MD++ format) potfile = "~/Codes/MD++/potentials/EAMDATA/eamdata.AuFoiles" eamgrid = 500 readeam } } #------------------------------------------------------------- #-------------------------------------------- proc readmeam-according-to-myname { myname } { if { [ string match "*baskes*" $myname ] } { puts "readmeam-baskes" readmeam-baskes } elseif { [ string match "*lammps*" $myname ] } { puts "readmeam-lammps" readmeam-lammps } 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" } } #end of proc readmeam-according-to-myname proc makeperfcrystal { ND NL a } { #-------------------------------------------- # Create Perfect Lattice Configuration # MD++ crystalstructure = face-centered-cubic MD++ latticeconst = $a MD++ latticesize = \[ -1 -1 2 $ND 1 -1 0 $ND 1 1 1 $NL \] MD++ makecrystal #finalcnfile = perf.cn writecn #eval } proc make_nanowire { radius } { #-------------------------------------------- # make a cylinder # #input = [3 2 # cylinder axis, zind yind /* 1x/2y/3z */ # 0 0 # center of circle in reduced coordinate # $radius # radius of circle (A) # 0 ] # 0 remove all atoms outside a cylinder, 1 fix them MD++ input = \[3 2 0 0 $radius 0 \] MD++ makecylinder #finalcnfile = 5Dcylinder.cn writecn #eval } #-------------------------------------------- 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 = [ 1 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=600 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 #atomcolor = lightgrey highlightcolor = purple bondcolor = darkgrey #plot_color_windows = [ 5 # -5 -4 0 # -5.5 -5 8 #color00 = red # -6 -5.5 1 # -6.55 -6 2 # -6.72 -6.55 3 # -6.5 -6 6 #color06 = gray80 # -7.5 -6.5 4 # ] plot_color_axis = 2 plot_color_windows = [ 3 3 5 1 5 10 4 10 50 6 ] #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 = 1 rotateangles = [ 0 0 0 1.9 ] } } proc openwindow { } { setup_window #MD++ openwin alloccolors rotate saverot plot MD++ openwin alloccolors rotate saverot eval plot #MD++ plot_color_axis = 0 input = [ -8 -3 10] GnuPlotHistogram #MD++ plot_color_axis = 2 input = \[ 0.6 50 50 \] GnuPlotHistogram } #-------------------------------------------- proc exitmd { } { MD++ quit } #end of proc exitmd #-------------------------------------------- #-------------------------------------------- proc setup_md { } { MD++ { equilsteps = 0 timestep = 0.001 # (ps) atommass = 196.97 # (g/mol) DOUBLE_T = 1 saveprop = 1 savepropfreq = 100 openpropfile #run savecn = 0 savecnfreq = 5000 openintercnfile savecfg = 0 savecfgfreq = 500 plotfreq = 10 printfreq = 10 vt2 = 3e26 #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_zz P_com(0) P_com(1) P_com(2) F_com(0) F_com(1) F_com(2) Ptheta_com torquespec(2) Torque H_33" ensemble_type = "NVT" integrator_type = "VVerlet" implementation_type = 1 } } #end of proc setup_md #-------------------------------------------- proc setup_md_npt { } { MD++ { equilsteps = 0 timestep = 0.001 # (ps) atommass = 196.97 # (g/mol) DOUBLE_T = 0 saveprop = 1 savepropfreq = 100 openpropfile savecn = 0 savecnfreq = 10000 openintercnfile plotfreq = 100 vt2 = 3e26 #1e28 2e28 5e28 wallmass = 1e3 # boxdamp = 1e-2 # 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 = "NPT" integrator_type = "Gear6" } } #end of proc setup_md_npt #******************************************* # Main program starts here #******************************************* # status 0: # 1: # 2: # # 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" set myname [ MD++_Get "myname"] puts "myname = $myname" readmeam-according-to-myname $myname if { [ string match "*baskes*" $myname ] } { set myname "meam-baskes" set lattconst 4.07 } elseif { [ string match "*lammps*" $myname ] } { set myname "meam-lammps" set lattconst 4.07 } elseif { [ string match "*meam*" $myname ] } { set myname "meam" set lattconst 4.07 } elseif { [ string match "*eam*" $myname ] } { set myname "eam" set lattconst 4.08 } else { set myname "uh-oh" puts "I don't know the lattice constant" puts "I set it to 4.08" set lattconst 4.08 } if { $status == 0 } { # create nanowire, apply handle, relax, anneal and relax it again initmd $n # radius (in Angstroms) and aspect ratio of nanowire set R 25 set aspectratio 3.0 # Create the nanowire set ND [expr round(ceil(2*$R/$lattconst)) ] set NL [expr round(ceil($aspectratio*$R/$lattconst)) ] set D [format "%.0f" [expr $R*2/10]] MD++ allocmultiple = 2 makeperfcrystal $ND $NL $lattconst set NP_cryst [ MD++_Get NP ] make_nanowire $R set NP_nw [ MD++_Get NP ] MD++ input = \[ 1 1 1.0 \] changeH_keepR MD++ input = \[ 2 2 1.0 \] changeH_keepR set vacuumratio [ expr 1 - 1.0*$NP_nw/$NP_cryst/4 ] MD++ vacuumratio = $vacuumratio MD++ finalcnfile = "nw-perf.cn" writecn openwindow #setup_window # Create handles for applying t-PBC set dphi 0.02 set dz 0.1 set phi 0 MD++ torsionsim = 1 torquespec = \[ $dphi $dz $phi \] maketorquehandle MD++ finalcnfile = "nw-perf.image" writeimagefile set vacuumratio [ expr 1 - 1.0*$NP_nw/$NP_cryst/4/(1+4*$dz) ] MD++ vacuumratio = $vacuumratio MD++ eval # Conjugate gradient relaxation relax_fixbox MD++ finalcnfile = "handle-relaxed.cn" writeall = 1 writecn MD++ eval # Molecular Dynamics simulation setup_md MD++ T_OBJ = 1000 totalsteps = 1000 initvelocity run MD++ finalcnfile = "handle-annealed.cn" writeall = 1 writecn MD++ eval # Conjugate gradient relaxation (again) relax_fixbox MD++ finalcnfile = "handle-annealed-relaxed.cn" writeall = 1 writecn MD++ eval exitmd } elseif { $status == 1 } { # create nanowire, relax, anneal and relax it again # (using PBc, to compare with status 0) initmd $n # radius (in Angstroms) and aspect ratio of nanowire set R 25 set aspectratio 3.0 # Create the nanowire set ND [expr round(ceil(2*$R/$lattconst)) ] set NL [expr round(ceil($aspectratio*$R/$lattconst)) ] set D [format "%.0f" [expr $R*2/10]] #MD++ allocmultiple = 2 #(not necessary in PBC) makeperfcrystal $ND $NL $lattconst set NP_cryst [ MD++_Get NP ] make_nanowire $R set NP_nw [ MD++_Get NP ] MD++ input = \[ 1 1 1.0 \] changeH_keepR MD++ input = \[ 2 2 1.0 \] changeH_keepR set vacuumratio [ expr 1 - 1.0*$NP_nw/$NP_cryst/4 ] MD++ vacuumratio = $vacuumratio MD++ finalcnfile = "nw-perf.cn" writecn openwindow #setup_window MD++ eval # Conjugate gradient relaxation relax_fixbox MD++ finalcnfile = "nw-relaxed.cn" writeall = 1 writecn MD++ eval # Molecular Dynamics simulation setup_md MD++ T_OBJ = 1000 totalsteps = 1000 initvelocity run MD++ finalcnfile = "nw-annealed.cn" writeall = 1 writecn MD++ eval # Conjugate gradient relaxation (again) relax_fixbox MD++ finalcnfile = "nw-annealed-relaxed.cn" writeall = 1 writecn MD++ eval exitmd } elseif { $status == 2 } { initmd $n # radius (in Angstroms) and aspect ratio of nanowire set R 25 set aspectratio 3.0 # Create the nanowire set ND [expr round(ceil(2*$R/$lattconst)) ] set NL [expr round(ceil($aspectratio*$R/$lattconst)) ] set D [format "%.0f" [expr $R*2/10]] MD++ allocmultiple = 2 makeperfcrystal $ND $NL $lattconst set NP_cryst [ MD++_Get NP ] make_nanowire $R set NP_nw [ MD++_Get NP ] MD++ input = \[ 1 1 1.0 \] changeH_keepR MD++ input = \[ 2 2 1.0 \] changeH_keepR set vacuumratio [ expr 1 - 1.0*$NP_nw/$NP_cryst/4 ] MD++ vacuumratio = $vacuumratio MD++ finalcnfile = "nw-perf.cn" writecn openwindow #setup_window # Create handles for applying t-PBC set dphi 0.02 set dz 0.1 set phi 0 MD++ torsionsim = 1 torquespec = \[ $dphi $dz $phi \] maketorquehandle set vacuumratio [ expr 1 - 1.0*$NP_nw/$NP_cryst/4/(1+4*$dz) ] MD++ vacuumratio = $vacuumratio # Conjugate gradient relaxation relax_freebox MD++ finalcnfile = "handle-relaxed-freebox.cn" writeall = 1 writecn MD++ eval # Molecular Dynamics simulation setup_md_npt MD++ T_OBJ = 300 NHMass = 1e-2 totalsteps = 30000 initvelocity run MD++ finalcnfile = "handle-annealed-NPT.cn" writeall = 1 writecn MD++ eval # Conjugate gradient relaxation (again) relax_freebox MD++ finalcnfile = "handle-annealed-relaxed-freebox.cn" writeall = 1 writecn MD++ eval exitmd } elseif { $status == 3 } { # create nanowire, relax, anneal and relax it again # (using PBc, to compare with status 0) initmd $n # radius (in Angstroms) and aspect ratio of nanowire set R 25 set aspectratio 3.0 # Create the nanowire set ND [expr round(ceil(2*$R/$lattconst)) ] set NL [expr round(ceil($aspectratio*$R/$lattconst)) ] set D [format "%.0f" [expr $R*2/10]] MD++ allocmultiple = 2 makeperfcrystal $ND $NL $lattconst set NP_cryst [ MD++_Get NP ] make_nanowire $R set NP_nw [ MD++_Get NP ] MD++ input = \[ 1 1 1.0 \] changeH_keepR MD++ input = \[ 2 2 1.0 \] changeH_keepR set vacuumratio [ expr 1 - 1.0*$NP_nw/$NP_cryst/4 ] MD++ vacuumratio = $vacuumratio MD++ finalcnfile = "nw-perf.cn" writecn #MD++ incnfile = "nw-relaxed.cn" readcn openwindow #setup_window MD++ eval # Conjugate gradient relaxation relax_freebox MD++ finalcnfile = "nw-relaxed-freebox.cn" writeall = 1 writecn MD++ eval # Molecular Dynamics simulation setup_md_npt MD++ T_OBJ = 300 NHMass = 1e-2 totalsteps = 30000 initvelocity run MD++ finalcnfile = "nw-annealed-NPT.cn" writeall = 1 writecn MD++ eval # Conjugate gradient relaxation (again) relax_freebox MD++ finalcnfile = "nw-annealed-relaxed-freebox.cn" writeall = 1 writecn MD++ eval exitmd } elseif { $status == 4 } { # starting from NPT equilibrated nanowire, set length, run NVT MD initmd $n # Create handles for applying t-PBC set dphi 0.02 set dz 0.1 set phi 0 MD++ torsionsim = 1 torquespec = \[ $dphi $dz $phi \] MD++ incnfile = "../au-nw-torsion-2/handle-annealed-NPT.cn" readcn set NP0 15637 MD++ NP0 = $NP0 NIMAGES = [expr [MD++_Get NP] - $NP0] # mass and timestep need to be set to interprete velocity and compute stress correctly MD++ timestep = 0.001 atommass = 196.97 MD++ vacuumratio = 0.963766201747 MD++ H_33 = [ expr 133.20*(1+4*$dz) ] SHtoR eval openwindow # Molecular Dynamics simulation setup_md MD++ { output_fmt = "curstep EPOT KATOM Tinst HELM HELMP TSTRESSinMPa_xx TSTRESSinMPa_yy TSTRESSinMPa_zz H_11 H_22 H_33" } MD++ T_OBJ = 300 totalsteps = 10000 run MD++ finalcnfile = "handle-annealed.cn" writeall = 1 writecn MD++ eval # Conjugate gradient relaxation (again) relax_fixbox MD++ finalcnfile = "handle-annealed-relaxed.cn" writeall = 1 writecn MD++ eval exitmd } elseif { $status == 5 } { # starting from NPT equilibrated nanowire, set length, run NVT MD # (using PBc, to compare with status 0) initmd $n MD++ incnfile = "../au-nw-torsion-3/nw-annealed-NPT.cn" readcn # mass and timestep need to be set to interprete velocity and compute stress correctly MD++ timestep = 0.001 atommass = 196.97 MD++ vacuumratio = 0.949272682446 MD++ H_33 = 133.20 SHtoR eval openwindow # Molecular Dynamics simulation setup_md MD++ { output_fmt = "curstep EPOT KATOM Tinst HELM HELMP TSTRESSinMPa_xx TSTRESSinMPa_yy TSTRESSinMPa_zz H_11 H_22 H_33" } MD++ T_OBJ = 300 totalsteps = 10000 run MD++ finalcnfile = "nw-annealed.cn" writeall = 1 writecn MD++ eval # Conjugate gradient relaxation (again) relax_fixbox MD++ finalcnfile = "nw-annealed-relaxed.cn" writeall = 1 writecn MD++ eval exitmd } elseif { $status == 6 } { # load configs from previous PBC run, eval, apply t-PBC and eval again # (results should be independent of whether handle is applied before or after annealing) initmd $n MD++ allocmultiple = 2 #MD++ incnfile = "../au-nw-torsion-3/nw-annealed-NPT.cn" readcn MD++ incnfile = "../au-nw-torsion-5/nw-annealed.cn" readcn # mass and timestep need to be set to interprete velocity and compute stress correctly MD++ timestep = 0.001 atommass = 196.97 MD++ vacuumratio = 0.949272682446 # potential call in PBC MD++ eval # Create handles for applying t-PBC set dphi 0.02 set dz 0.1 set phi 0 MD++ torsionsim = 1 torquespec = \[ $dphi $dz $phi \] MD++ incnfile = "../au-nw-torsion-0/nw-perf.image" readimagefile fix_imageatoms MD++ input = \[ 3 3 [expr 4*$dz] \] changeH_keepR copytorqueatoms MD++ vacuumratio = 0.963766201747 # potential call in t-PBC MD++ eval # Molecular Dynamics simulation setup_md MD++ { output_fmt = "curstep EPOT KATOM Tinst HELM HELMP TSTRESSinMPa_xx TSTRESSinMPa_yy TSTRESSinMPa_zz H_11 H_22 H_33" } MD++ T_OBJ = 300 totalsteps = 10000 run MD++ finalcnfile = "handle-annealed.cn" writeall = 1 writecn MD++ eval exitmd } elseif { $status == 7 } { # continue a previous NVT-MD run under t-PBC initmd $n MD++ incnfile = "../au-nw-torsion-4/handle-annealed.cn" readcn # mass and timestep need to be set to interprete velocity and compute stress correctly MD++ timestep = 0.001 atommass = 196.97 # set handle information set dphi 0.02 set dz 0.1 set phi 0 MD++ torsionsim = 1 torquespec = \[ $dphi $dz $phi \] MD++ incnfile = "../au-nw-torsion-0/nw-perf.image" readimagefile MD++ vacuumratio = 0.963766201747 # potential call in t-PBC MD++ eval # Molecular Dynamics simulation setup_md MD++ { output_fmt = "curstep EPOT KATOM Tinst HELM HELMP TSTRESSinMPa_xx TSTRESSinMPa_yy TSTRESSinMPa_zz H_11 H_22 H_33" } MD++ T_OBJ = 300 totalsteps = 10000 run MD++ finalcnfile = "handle-annealed.cn" writeall = 1 writecn MD++ eval exitmd } elseif { $status == 8 } { # apply torsion on configuration already equilibrated under NVT/t-PBC # set handle information set dphi 0.005 set dz 0.1 set phi 0 set total_twist_steps 400 initmd2 $n $dphi MD++ incnfile = "../au-nw-torsion-7/handle-annealed.cn" readcn # mass and timestep need to be set to interprete velocity and compute stress correctly MD++ timestep = 0.001 atommass = 196.97 MD++ torsionsim = 1 torquespec = \[ $dphi $dz $phi \] MD++ incnfile = "../au-nw-torsion-0/nw-perf.image" readimagefile MD++ vacuumratio = 0.963766201747 #openwindow setup_window # potential call in t-PBC MD++ eval # Molecular Dynamics simulation setup_md MD++ T_OBJ = 300 totalsteps = 10000 for { set x 1 } { $x <= $total_twist_steps } { incr x 1 } { MD++ addtorque run copytorqueatoms MD++ finalcnfile = "run$x.cfg" writeatomeyecfg MD++ finalcnfile = "run$x.cn" writeall = 1 writecn MD++ eval } exitmd } elseif { $status == 9 } { # apply torsion on configuration already equilibrated under NVT/PBC # need to apply handle to it before apply twist angle # set handle information set dphi 0.005 set dz 0.1 set phi 0 set total_twist_steps 400 initmd2 $n $dphi MD++ allocmultiple = 2 MD++ incnfile = "../au-nw-torsion-5/nw-annealed.cn" readcn # mass and timestep need to be set to interprete velocity and compute stress correctly MD++ timestep = 0.001 atommass = 196.97 MD++ vacuumratio = 0.949272682446 # potential call in PBC MD++ eval # Create handles for applying t-PBC MD++ torsionsim = 1 torquespec = \[ $dphi $dz $phi \] MD++ incnfile = "../au-nw-torsion-0/nw-perf.image" readimagefile fix_imageatoms MD++ input = \[ 3 3 [expr 4*$dz] \] changeH_keepR copytorqueatoms MD++ vacuumratio = 0.963766201747 #openwindow setup_window # potential call in t-PBC MD++ eval # Molecular Dynamics simulation setup_md MD++ T_OBJ = 300 totalsteps = 10000 for { set x 1 } { $x <= $total_twist_steps } { incr x 1 } { MD++ addtorque run copytorqueatoms MD++ finalcnfile = "run$x.cfg" writeatomeyecfg MD++ finalcnfile = "run$x.cn" writeall = 1 writecn MD++ eval } exitmd } else { puts "unknown status = $status" exitmd }