LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
[lammps-users] Shock propagation , script not running
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[lammps-users] Shock propagation , script not running


From: anuradha singla <asanusingla@...24...>
Date: Tue, 19 Dec 2017 11:20:46 +0530

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_dt*100"
variable pdamp equal "v_dt*1000"

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/cm**3

# 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.05**2+4.05**2)**.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]