# -*-shell-script-*- # testing MEAM potential for Si source "scripts/Examples/Tcl/startup.tcl" #******************************************* # Definition of procedures #******************************************* proc initmd { } { MD++ setnolog MD++ setoverwrite MD++ dirname = runs/meam-si-c } #end of proc initmd proc create_dc_crystal { nx ny nz a } { #-------------------------------------------- # Create Perfect Lattice Configuration MD++ crystalstructure = diamond-cubic MD++ latticeconst = $a MD++ latticesize = \[ 1 0 0 $nx 0 1 0 $ny 0 0 1 $nz \] MD++ makecrystal MD++ writeall = 1 zipfiles = 1 finalcnfile = "dc-perf.cn" writecn MD++ incnfile = "dc-perf.cn" readcn } #end of proc create_dc_crystal proc create_zb_crystal { nx ny nz a } { #-------------------------------------------- # Create Perfect Lattice Configuration MD++ crystalstructure = zinc-blende MD++ latticeconst = $a MD++ latticesize = \[ 1 0 0 $nx 0 1 0 $ny 0 0 1 $nz \] MD++ makecrystal MD++ writeall = 1 zipfiles = 1 finalcnfile = "zb-perf.cn" writecn MD++ incnfile = "zb-perf.cn" readcn } #end of proc create_zb_crystal #------------------------------------------------------------ proc readmeam-lammps-Si { } { MD++ { #Read in MEAM potential (Baskes format) meamfile = "../../scripts/work/meam/meamf" nspecies = 1 element0 = "Siz" rcut = 6.0 readMEAM NNM = 300 } } #------------------------------------------------------------ proc readmeam-Si { } { MD++ { #Read in MEAM potential (MD++ format) potfile = "../../scripts/work/meam/meamdata.Siz" readMEAM rcut = 6.0 NNM = 300 } } #------------------------------------------------------------ proc readmeam-lammps-Si-C { } { MD++ { #Read in MEAM potential (Baskes format) meamfile = "../../scripts/work/meam/meamf" meafile = "../../scripts/work/meam/SiC.meam" nspecies = 2 element0 = "Si" element1 = "C" rcut = 4.5 readMEAM NNM = 300 } } #------------------------------------------------------------- proc openwindow { } { MD++ { # Plot Configuration # atomradius = [ 0.67 0.50 ] bondradius = 0.3 bondlength = 2 #2.8285 #for Si atomcolor0 = orange atomcolor1 = magenta bondcolor = red backgroundcolor = gray70 fixatomcolor = yellow plotfreq = 100 rotateangles = [ 0 0 0 1.25 ] openwin alloccolors rotate saverot eval plot } } #end of proc openwindow #-------------------------------------------- proc exitmd { } { MD++ quit } #end of proc exitmd #-------------------------------------------- proc relax_freebox { } { MD++ { # Conjugate-Gradient relaxation conj_ftol = 2e-6 conj_itmax = 1000 conj_fevalmax = 3000 conj_fixbox = 0 conj_fixboxvec = [ 1 1 1 1 1 1 1 1 0 ] relax } } #end of proc relax_freebox #-------------------------------------------- proc relax_fixbox { } { MD++ { # Conjugate-Gradient relaxation conj_ftol = 2e-6 conj_itmax = 1000 conj_fevalmax = 3000 conj_fixbox = 1 conj_fixboxvec = [ 0 1 1 1 0 1 1 1 0 ] relax } } #end of proc relax_freebox #-------------------------------------------- proc setup_md { } { MD++ { equilsteps = 0 timestep = 0.0001 # (ps) atommass = [ 28.0855 12.011 ] # (g/mol) DOUBLE_T = 0 saveprop = 1 savepropfreq = 100 openpropfile savecn = 1 savecnfreq = 1000 openintercnfile plotfreq = 100 vt2 = 2e28 #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 stress = [ 0 0 0 0 0 0 0 0 0 ] output_fmt = "curstep EPOT KATOM Tinst HELM HELMP H_12 H_22 H_32 TSTRESSinMPa_xy TSTRESSinMPa_yy TSTRESSinMPa_zy " ensemble_type = "NVT" integrator_type = "VVerlet" implementation_type = 1 } } #end of proc setup_md #******************************************* # Main program starts here #******************************************* # status 0:prepare liquid, 1: prepare liquid-solid interface, 2,3: measure Tm # 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 { $status == 0 } { # MD simulation of perfect Si crystal (NVT) ensemble initmd readmeam-lammps-Si # for comparison purpose (run with meam_$(SYS) ) #readmeam-Si create_dc_crystal 5 5 5 5.431 openwindow MD++ plot setup_md MD++ ensemble_type = "NVT" integrator_type = "VVerlet" implementation_type = 1 MD++ T_OBJ = 1000 totalsteps = 10000 initvelocity run MD++ finalcnfile = "si-equil.cn" writeall = 1 writecn relax_fixbox MD++ finalcnfile = "si-relaxed.cn" writeall = 1 writecn MD++ sleep exitmd } elseif { $status == 1 } { # MD simulation of perfect Si crystal (NPT) ensemble initmd readmeam-lammps-Si # for comparison purpose (run with meam_$(SYS) ) #readmeam-Si create_dc_crystal 5 5 5 5.431 openwindow MD++ plot setup_md MD++ ensemble_type = "NPT" integrator_type = "Gear6" implementation_type = 0 # fix three box components to prevent arbitrary rotation MD++ { fixboxvec = [ 0 0 0 1 0 1 1 0 0 ] } # specify applied stress components MD++ { stress = [ 0 1 0 1 0 0 0 0 0 ] } MD++ stressmul = 1000 # in MPa MD++ T_OBJ = 1000 totalsteps = 10000 initvelocity run MD++ finalcnfile = "si-equil.cn" writeall = 1 writecn relax_fixbox MD++ finalcnfile = "si-relaxed.cn" writeall = 1 writecn MD++ sleep exitmd } elseif { $status == 2 } { # MD simulation of perfect SiC crystal (NPT) ensemble initmd readmeam-lammps-Si-C create_zb_crystal 5 5 5 4.3596 openwindow MD++ plot setup_md MD++ ensemble_type = "NPT" integrator_type = "Gear6" implementation_type = 0 # fix three box components to prevent arbitrary rotation MD++ { fixboxvec = [ 0 0 0 1 0 1 1 0 0 ] } # specify applied stress components MD++ { stress = [ 0 1 0 1 0 0 0 0 0 ] } MD++ stressmul = 1000 # in MPa MD++ T_OBJ = 1000 totalsteps = 10000 initvelocity run MD++ finalcnfile = "sic-equil.cn" writeall = 1 writecn relax_fixbox MD++ finalcnfile = "sic-relaxed.cn" writeall = 1 writecn MD++ sleep exitmd } else { puts "unknown status = $status" exitmd }