# ---------- Initialize Simulation --------------------- clear units metal dimension 3 boundary p p p atom_style atomic atom_modify map array variable a loop 9 # ---------- define variables --------------------- variable Ts equal 300 variable alat equal 3.304 variable myseed equal 12345 variable xm equal 10 variable ym equal 10 variable zm equal 420 #LJ variable dt equal 0.001 variable t_eq equal 10000 variable t_sh equal 50000 variable vpis index 0.1 0.3 0.5 1.0 1.3 1.8 2.2 2.5 3 variable Nevery equal 10 variable Nrepeat equal 5 variable Nfreq equal 500 variable dz equal 10 variable atomrate equal 5000 variable tdamp equal "v_dt*100" variable pdamp equal "v_dt*1000" variable Up equal "10*v_vpis" timestep ${dt} # ---------- Create Atoms --------------------- lattice bcc ${alat} origin 0.0 0.0 0.0 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1 #tmd 27g/NA 6e23 *4 atoms/cell/(4.05e-8cm)**3 = 2.7 g/cm**3 # define size of the simulation box region sim_box block 0 ${xm} 0 ${ym} 0 ${zm} units lattice create_box 2 sim_box # define atoms in sim_box region atom_box block 0 ${xm} 0 ${ym} 0 ${zm} create_atoms 1 region atom_box # define a group for the atom_box region group atom_box region atom_box region piston block INF INF INF INF INF 4 region bulk block INF INF INF INF 4 INF group piston region piston group bulk region bulk set group piston type 1 set group bulk type 2 # ---------- Define Interatomic Potential --------------------- # for lj, epsilon (eV) = heat of sublimation/atom # for lj, sigma (ang)= 2**(-1/6)*nearest neighbour distance # epsilon is 3.34 eV (cohesive energy per atom from Avinc) # sigma is 2**(-1/6) * (4.05**2+4.05**2)**.5 /2 angstrom == 2.55 Angstrom * 2.5 == 6.38 mass 1 180.95 mass 2 180.95 #pair_style lj/cut 6.38 #pair_coeff * * 3.34 2.55 pair_style eam/alloy pair_coeff * * Ta.eam.alloy Ta Ta # al1.eam.fs is file available from http://www.ctcms.nist.gov/potentials # pair_style eam/fs # pair_coeff * * al1.eam.fs Al Al # al99.eam.alloy is file available from http://www.ctcms.nist.gov/potentials # pair_style eam/alloy # pair_coeff * * Al99.eam.alloy Al Al #mass 1 180.95 #mass 2 180.95 compute myCN bulk cna/atom 3.0 compute myKE bulk ke/atom compute myPE bulk pe/atom compute myCOM bulk com compute peratom bulk stress/atom NULL compute vz bulk property/atom vz compute vorvol bulk voronoi/atom # ------------ Equilibrate --------------------------------------- reset_timestep 0 # Now, assign the initial velocities using Maxwell-Boltzmann distribution velocity atom_box create ${Ts} ${myseed} rot yes dist gaussian fix equilibration bulk npt temp ${Ts} ${Ts} ${tdamp} iso 0 0 ${pdamp} drag 1 variable eq1 equal "step" variable eq2 equal "pxx/10000" variable eq3 equal "pyy/10000" variable eq4 equal "pzz/10000" variable eq5 equal "lx" variable eq6 equal "ly" variable eq7 equal "lz" variable eq8 equal "temp" variable eq9 equal "etotal" variable eq11 equal "density" fix output1 bulk print 10 "${eq1} ${eq2} ${eq3} ${eq4} ${eq5} ${eq6} & ${eq7} ${eq8} ${eq9} ${eq11}" file run.out.${vpis} screen no thermo 10 thermo_style custom step pxx pyy pzz lx ly lz temp etotal density run ${t_eq} unfix equilibration unfix output1 # -------------- Shock ------------------------------------------- change_box all boundary p p s reset_timestep 0 # WE CREATE A PISTON USING A FEW LAYERS OF ATOMS AND THEN WE GIVE IT # A CONSTANT POSTIVE SPEED. YOU COULD ALSO USE LAMMPS' FIX WALL/PISTON COMMAND fix 1 all nve velocity piston set 0 0 v_Up sum no units box fix 2 piston setforce 0.0 0.0 0.0 # WE CREATE BINS IN ORDER TO TRACK THE PASSING OF THE SHOCKWAVE compute 3 bulk chunk/atom bin/1d z lower ${dz} units box fix density_profile bulk ave/chunk ${Nevery} ${Nrepeat} ${Nfreq} 3 density/mass file denz.${vpis} variable temp atom c_myKE/(1.5*8.61e-5) fix temp_profile bulk ave/chunk ${Nevery} ${Nrepeat} ${Nfreq} 3 v_temp file temp.${vpis} # meanpress was originally pressure *volume per atom (cubic distance units) # manual suggests use compute voronoi/atom to estimate atomic volume variable meanpress atom -(c_peratom[1]+c_peratom[2]+c_peratom[3])/3/(c_vorvol[1]+1e-99) fix pressure_profile bulk ave/chunk ${Nevery} ${Nrepeat} ${Nfreq} 3 v_meanpress file pressure.${vpis} fix velZ_profile bulk ave/chunk ${Nevery} ${Nrepeat} ${Nfreq} 3 c_vz file velocityZcomp.${vpis} variable eq1 equal "step" variable eq2 equal "pxx/10000" variable eq3 equal "pyy/10000" variable eq4 equal "pzz/10000" variable eq5 equal "lx" variable eq6 equal "ly" variable eq7 equal "lz" variable eq8 equal "temp" variable eq9 equal "etotal" variable eq10 equal "c_myCOM[3]" variable eq11 equal "density" fix shock bulk print 10 "${eq1} ${eq2} ${eq3} ${eq4} ${eq5} ${eq6} ${eq7} ${eq8} ${eq9} ${eq10} ${eq11}" file run.${Ts}K.${vpis} screen no thermo 10 thermo_style custom step pxx pyy pzz lx ly lz temp etotal c_myCOM[3] density #Use cfg(?) for AtomEye dump 1 all cfg ${atomrate} dump._*.${vpis}.cfg mass type xs ys zs id c_myPE c_myKE c_myCN dump_modify 1 element Ta Ta #dump 2 all image 1000 image.*.jpg type type axes yes 0.8 0.02 view 120 -90 run ${t_sh} clear next vpis next a jump Tashock.in