Hii,
I am trying to study shock wave propagation in Aluminium using LAMMPS. To start LJ potential is being used. I am using the following script. But it stops running giving an error
ERROR: Invalid compute ID in variable formula (…/variable.cpp:903).
I feel the problem lies in the value of Nevery Nrepeat and Nfreq. But I can’t figure out the problem. I am using LAMMPS-2013 version where ave/spatial is working. The script is
---------- Initialize Simulation ---------------------
clear
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array
---------- define variables ---------------------
variable Ts equal 300
variable alat equal 4.05
variable xm equal 4
variable ym equal 4
variable zm equal 42 #LJ
#variable zm equal 140 #EAM
variable dt equal 0.001
variable t_eq equal 10000
variable t_sh equal 50000
variable vpis equal 0.10
variable Nevery equal 1000
variable Nrepeat equal 5
variable Nfreq equal 10000
variable dz equal 20
variable atom equal 100
variable tdamp equal “v_dt100"
variable pdamp equal "v_dt1000”
variable Up equal “10*v_vpis”
timestep ${dt}
---------- Create Atoms ---------------------
lattice fcc ${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/cm3
define size of the simulation box
region sim_box block 0 {xm} 0 {ym} 0 ${zm}
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.052+4.052)**.5 /2 angstrom == 2.55 Angstrom * 2.5 == 6.38
pair_style lj/cut 6.38
pair_coeff * * 3.34 2.55
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 27.0
mass 2 27.0
#compute myCN bulk cna/atom 3.45708
compute myCN bulk cna/atom 6.38
compute myKE bulk ke/atom
compute myPE bulk pe/atom
compute myCOM bulk com
compute peratom bulk stress/atom virial
compute vz bulk property/atom vz
------------ Equilibrate ---------------------------------------
reset_timestep 0
Now, assign the initial velocities using Maxwell-Boltzmann distribution
velocity atom_box create {Ts} 12345 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”
fix output1 bulk print 10 “{eq1} {eq2} {eq3} {eq4} {eq5} {eq6} &
{eq7} {eq8} ${eq9}” file run.out screen no
thermo 10
thermo_style custom step pxx pyy pzz lx ly lz temp etotal
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
fix density_profile bulk ave/spatial {Nevery} {Nrepeat} {Nfreq} z lower {dz} density/mass file denz.profile units box
variable temp atom c_myKE/(1.5*8.61e-5)
fix temp_profile bulk ave/spatial {Nevery} {Nrepeat} {Nfreq} z lower {dz} v_temp file temp.profile units box
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_vol[1]
fix pressure_profile bulk ave/spatial {Nevery} {Nrepeat} {Nfreq} z lower {dz} v_meanpress units box file pressure.profile
fix velZ_profile bulk ave/spatial {Nevery} {Nrepeat} {Nfreq} z lower {dz} c_vz units box file velocityZcomp.profile
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]”
fix shock bulk print 10 “{eq1} {eq2} {eq3} {eq4} {eq5} {eq6} {eq7} {eq8} {eq9} {eq10}” file run.${Ts}K.out screen no
thermo 10
thermo_style custom step pxx pyy pzz lx ly lz temp etotal c_myCOM[3]