# Shock Modeling in Lammps

hi
I am studying the shock response of material by FEM and as well as SPH , now i intend to study this type of mechanics by MD,so I want to do some shock simulation by Lammps,as per i go through the manual(shock package) and discussion about this topics in this mail list by steve, oscar G and Ray ,I come to know that I can generate shock in different ways ,for this i tried to use fix_ msst although i go through the paper by Reed et al(2003) i have one question-
to use msst ppp bc is necessary so when i send some shock to +Z direction, from video it is showing as if it’s propagating from each direction(+Z & -Z dirn.) simultaneously,is it reflecting from opposite dirn ?
for this difficulty i am using wall/piston option according to steve suggestion , here do i have to group some region of my geometry ? e.g named as PISTON and then

``````"fix walls PISTON wall/piston zlo pos 0.0 0.0 0.0 vel 0.0 0.0 0.0 units box" ?

best rgds
oxem

``````

Aidan may be able to help.

Steve

Hi Selesta,

Assumption of MSST is that the shock wave is homogeneous within the
small representative sample. Therefore the sample is compressed
uni-axially from both ends along the shock orientation when shock wave
arrives. You may get gradient of some sort with direct NEMD
simulations when looking at a slice of same size, but MSST has no such

I think with wall/piston, you should bound the whole simulation box by
applying to group-all.

Ray

Read the MSST papers that we cite on the doc page.

Hi…

Thanks for your interest in my earlier post =>

I will copy and paste some of the reference in regard shock wave :

In NEMD simulations, there are three principal ways to generate a shock wave: (1) As in laboratory planar-
shockwave experiments, one can computationally “hurl” a flyer plate toward a stationary target at a velocity of
2up, where up is the ”piston” or particle velocity. This is equivalent to slamming the two plates together at &up;
in this symmetric-impact case, a pair of shock waves move out from the interface at fu8 (il8 is the shockwave
velocity). (2) The symmetric impact can be generated by inhomogeneously shrinking the longitudinal periodic
length. This is useful, particularly for fluid shockwave simulations, in order to eliminate entirely the effect of
any free surfaces ]. (3) Material can be pushed by an infinitely massive piston moving at velocity up; all parti-
cles coming in contact with the piston face are specularly reflected (“momentum mirror,).

Usually for my MD simulations in Copper and Tantalum I use the “Momentum Mirror” although you coul also set a piston using the command fixforce and then the command move…

correction : There are also other ways to generare shockwaves (Mimics)
Oscar G.

Hi oscar & Ray

thanks a lot to you for the earlier so valuable feedback.Now i am using classical method i.e" momentum mirror " method to produce shock.I am setting “-up” to push my system towards a massive piston by “wall/piston”.But after getting momentum towards -Z,my system shrinks and after few steps it explodes towards +Z direction and as well as temp is increasing very rapidly although i am using “temp” in “wall/piston” for langevin kicks.Furthermore my simulation box also shifted towards +Z dirn.need some discussion please. My input file is here–

units metal
dimension 3
boundary p p s
atom_style atomic
lattice fcc 4.05
region box block 0 10 0 10 0 30

create_box 1 box
create_atoms 1 box
mass 1 27
pair_style eam/alloy
pair_coeff * * Al99.eam.alloy Al
velocity all create 300.0 9999
fix 2 all nvt temp 300.0 300.0 100
timestep .001
thermo 100
run 1000
unfix 2
fix walls all wall/piston zlo pos 0.0 0.0 0.0 vel 0.0 0.0 0.0 temp 10 .001 10 40 units box
velocity all set 0 0 -10.00
fix 1 all nve
thermo 500
timestep .001

compute cna all cna/atom 4.05
compute csym all centro/atom fcc
compute peratom all pe/atom
compute cluster all cluster/atom 4.05

thermo_style custom step lx ly lz press pe temp

dump 1 all custom 500 dump.comp_*.cfg id type xs ys zs c_cna c_csym c_peratom c_cluster fx fy fz
dump_modify 1 element Al
run 50000

output is as follows-----

Step Temp E_pair E_mol TotEng Press Volume
0 300 -40798.541 0 -40325.488 1549.7896 199330.23
100 171.04149 -40595.109 0 -40325.404 2720.0803 199330.23

900 153.37401 -40567.235 0 -40325.388 452.0537 199330.23
1000 155.78787 -40571.04 0 -40325.387 172.40572 199330.23

Step Lx Ly Lz Press PotEng Temp
1000 40.5 40.5 122.59395 147491.33 -40571.04 17744.502
WARNING: Too many common neighbors in CNA 6 times (compute_cna_atom.cpp:353)
1500 40.5 40.5 101.58173 398357.92 -30907.723 11676.172
SHOCK: step 2000 t 1 zpos 0 vz 0 az 4.94066e-324 zlo -0.0119286
WARNING: Too many common neighbors in CNA 21 times (compute_cna_atom.cpp:353)
2000 40.5 40.5 81.208517 740164.67 -20951.113 5365.5948
WARNING: Too many common neighbors in CNA 23 times (compute_cna_atom.cpp:353)
2500 40.5 40.5 84.344177 716846.23 -19406.908 4388.7737
SHOCK: step 3000 t 2 zpos 0 vz 0 az 4.94066e-324 zlo -0.0121474
WARNING: Too many common neighbors in CNA 44 times (compute_cna_atom.cpp:353)
3000 40.5 40.5 105.2804 408054.93 -27669.71 9634.3831
WARNING: Too many common neighbors in CNA 186 times (compute_cna_atom.cpp:353)
3500 40.5 40.5 126.51777 220966.02 -34787.445 14149.444
SHOCK: step 4000 t 3 zpos 0 vz 0 az 4.94066e-324 zlo -0.00225957
WARNING: Too many common neighbors in CNA 40 times (compute_cna_atom.cpp:353)
4000 40.5 40.5 147.51142 97152.835 -36148.263 15012.589
WARNING: Too many common neighbors in CNA 8 times (compute_cna_atom.cpp:353)
4500 40.5 40.5 169.31071 51149.705 -35194.725 14407.951
SHOCK: step 5000 t 4 zpos 0 vz 0 az 4.94066e-324 zlo 9.06546
WARNING: Too many common neighbors in CNA 6 times (compute_cna_atom.cpp:353)
5000 40.5 40.5 184.40538 41607.615 -34607.247 14035.392
WARNING: Too many common neighbors in CNA 5 times (compute_cna_atom.cpp:353)
5500 40.5 40.5 194.71686 43615.175 -34430.182 13923.095
SHOCK: step 6000 t 5 zpos 0 vz 0 az 4.94066e-324 zlo 27.8177
WARNING: Too many common neighbors in CNA 8 times (compute_cna_atom.cpp:353)
6000 40.5 40.5 206.20151 43720.753 -34355.677 13875.829
WARNING: Too many common neighbors in CNA 13 times (compute_cna_atom.cpp:353)
6500 40.5 40.5 213.50777 47332.831 -34490.981 13961.618
SHOCK: step 7000 t 6 zpos 0 vz 0 az 4.94066e-324 zlo 49.7798
WARNING: Too many common neighbors in CNA 21 times (compute_cna_atom.cpp:353)
7000 40.5 40.5 219.05479 52555.997 -34627.552 14048.204
WARNING: Too many common neighbors in CNA 19 times (compute_cna_atom.cpp:353)
7500 40.5 40.5 227.14221 53848.569 -34679.184 14080.941
.
.

as per Oscar comment to Swain,I have checked with different “up”,but the scenario are nearly same. This seems like my system behaves like a fluid, i.e its satisfy Oscar’s 2nd point i.e shock_fluid (see oscar’s previous comment).But i want to track shock to study system behavior,bt how can it be solved? please comment.

best
Oxem

Hi Selesta,

velocity all set 0 0 -10.00
as per Oscar comment to Swain,I have checked with different "up",but the
scenario are nearly same.

I was wondering what values of Up have you also tried? -10 A/ps is
-1000 km/s, which is roughly 2 orders of magnitude larger than usual.
This explains why your temperature was flying high, and the shock wave
travels through the sample so quickly. In these kinds of simulations,
the simulation ends when the shock wave hits the +Z free surface,
which in your case ended in ~2000th step. Try with smaller velocities
and smaller timesteps as well.

Really 10 A/ps is 1000 km/s ? … Sighs… I just made a terrible algebra converson error …=(

Thanks
Oscar Guerrero

No wait, I was wrong. 10 A/ps is 1000 "m"/s = 1 km/s. Maybe all
Selesta needs is a smaller timestep.

Also try without the temp keyword.

Ray

Maybe nigel may wish to comment in regard why you obtain this abruptly change of temperature in your simulation .

Thanks
Oscar G.

Also am… You need to read more references… And temperature should increase (since you apply a shockwave)… sighs…
Oscar G.

Selesta ,

I saw you script and i saw a terrible error in your simulation . Your using :
“velocity all set 0 0 -10.00” and that is wrong .

You must use : “velocity atom_box set 0.0 0.0 -10 sum no units box” , where
-10 represents -10 A/ps

The units option is used by set and ramp. If units = box, the velocities and coordinates specified in the velocity command are in the standard units described by the units command (e.g. Angstroms/fmsec for real units). If units = lattice, velocities are in units of lattice spacings per time (e.g. spacings/fmsec) and coordinates are in lattice spacings

The option defaults is units = lattice

So you must use “unit box” to have the correct units…

Oscar G.

Hi Oscar & Ray

really it was a bad mistake, thanks goes to you for your qualified suggestion.
BTW now temperature does not increasing so rapidly and as well as now I can visualize the shock propagation, i am attaching one .avi file consisting some timesteps, please see.To see what happen i introduced one void into it and it now it is collapsing when the shock propagates.
But another problem occurs that after some steps, may be due to high reflected momentum (Up = 1.2 Km/S) my total simulation box shifted gradually along shock direction.How you explain this ?

secondly, here it is shown that before passing the shock, partially the void becomes distorted (i think it should not be happen as i give same velocity to all ),please discuss about it.

best
Selesta

jpg.avi (337 KB)

Selesta,

But another problem occurs that after some steps, may be due to high
reflected momentum (Up = 1.2 Km/S) my total simulation box shifted gradually
along shock direction.How you explain this ?

If you throw a ball against a wall, it bounces back to you. Same
thing happened here.

secondly, here it is shown that before passing the shock, partially the void
becomes distorted (i think it should not be happen as i give same velocity