Dear shocksimulation community,
I am quite new in atomistic level (MD) shock compression simulation, but I have experience in LAMMPS and as well as MD. However I have set the methodology successfully in LAMMPS. I have read through all the previous discussions regarding shock wave simulation by Ray Shan, Oscar G and Steve along wit the literature.
What I am doing:
I am setting the so called “momentum mirror method” by utilizing “wall/reflect” velocity “set” with keyword “sum” “yes” and doing NVE to study metallic and polymer system. Basically i’m setting a velocity of each atoms towards the mirror. I think my methodology is correct as I am seeing that the shock wave propagates through the bulk, by computing cna/atom and doing centrosymetry analysis we confirmed.
I am computing the density and velocity (after ) per region basis (I have defined several special regions by group and region command) by using “reduce/region” ,“ave/spatial” etc.
What my confusions are ?

As I intend to calculate the shock velocity(Us), I am plugging the two densities (before and after shock passed) and particle velocity in the Rankin Hugoniot jump condition ( momentum conservation eqn.). But here which particle velocity (Up) value I should plug the velocity I am setting by velocity “set” or the velocity of the particles i have computed after subtracting off the COM bias ?
Is it a correct way to compute shock velocity(Us) or any other method I should use, please share your experience.
(i have gone through “Budzien, Thompson and Zybin, JPC B 113, 13142 2009” suggested by Ray in some previous discussion regarding shock)

The second thing is that I compute temperature from kinetic energy( I have computed after subtracting off the COM velocity) from the eqn E=3/2 KT and it’s a good agreement wit the existing experimental values and as well as simulation literature, but the way in which I should compute pressure of the shock region is a bit confusion to me.
Would you please discuss the way I should follow to get the pressure from the velocity or kinetic energy I have computed for the shock region of the sample ?
Thank you for your attention.
Thanks
Anirban Dhar
Tata Institute of Fundamental Research(TIFR)
Mumbai, India
 As I intend to calculate the shock velocity(Us), I am plugging the two densities (before and after shock passed) and particle velocity in the Rankin Hugoniot jump condition ( momentum conservation eqn.). But here which particle velocity (Up) value I should plug the velocity I am setting by velocity “set” or the velocity of the particles i have computed after subtracting off the COM bias ?
Is it a correct way to compute shock velocity(Us) or any other method I should use, please share your experience. (i have gone through “Budzien, Thompson and Zybin, JPC B 113, 13142 2009” suggested by Ray in some previous discussion regarding shock)
COMMENT:
Up is the piston velocity (impact speed) and it can be set to any arbitrary value. Given UP you find a US , you plot UP vs US to calculate the hugoniot curve and from the rankle conditions you can find pressure, density and energy , but not temperature .
There are many ways to calculate (US) , For example the reference from Lan Ha et al “Molecular dynamics simulations of shock wave in oriented nitromethane single crystals” find the shock velocity by calculating the center of mass position for each layer along the shock direction and the layer with maximum kinetic energy (please read the reference). http://link.aip.org/link/doi/10.1063/1.3561397
Other methods for tracking of a shockwave are by Mottsmith profile or by means of calculating the density shifts . Please read :
http://www.ifd.mavt.ethz.ch/research/group_tr/projects/proj_schlamp/PDFs/torstenhofmann.pdf page 13 ) .
There is a rich and variety of reference available for the people interested in Shock wave phenomena. Please take the time to browser around the web…
2.The second thing is that I compute temperature from kinetic energy( I have computed after subtracting off the COM velocity) from the eqn E=3/2 KT and it’s a good agreement wit the existing experimental values and as well as simulation literature, but the way in which I should compute pressure of the shock region is a bit confusion to me. Would you please discuss the way I should follow to get the pressure from the velocity or kinetic energy I have computed for the shock region of the sample ?
COMMENT : I dont undertand quite well what your saying, but of you have the velocity (US) then you can find the pressure by using the rankie condition
P = densityUSUP
NOTE: this kind of questions are usually answered by looking at references sources by yourself. However i dont mind to use a little of my time to help others .
Since many People Are asking , i am tempted to make a brief tutorial to explain how to model Shock/Uniaxial and Quasiisentropic compression using lammps (please stay tuned for more info)
A salute
Oscar Guerrero
Hi Anirban,
I’ll just add the following to Carlos’ nice comments.

In this kind of simulations, particle velocity is the piston velocity or the impact velocity, which is the velocity you “set”.

If you wish to estimate stress tensor, instead of pressure, you can use the following:
Ray,
It wasn’t me the one who wrote the comments. It was Oscar who has a Spanish name as well but as far as I know is not my twin or clone…
Carlos
Pressed ‘enter’ too fast …
 If you wish to estimate stress tensor along shock direction, you can use the following:
Stress_xx = [ Sum(Force_xx of i) / (LyLz) ] + [ Sum(mass_i * Vx* Vx) / (LxLyLz)],
where Vx = (V_x of i minus Vcom_x), and Lx, Ly, Lz are spatial decomposed bins of the system (which can be thin slices), Vcom_x is the COM of the bin, and i runs through all atoms in a bin.
Cheers,
Ray
Right! I was ahead of myself. Sorry, Carlos and Oscar; my bad.
Cheers,
Ray
thanks a lot Ray and Oscar for your valuable attention and useful comments, i’m trying to implement what you are suggesting.
Oscar it will be a great help for all the persons specially for the learners of “shock simulation community”, if you share your experience likely with some documentation, a advance thanks.
Dear Ray and oscar,
I have read the reference you suggested and have learn how to compute shock velocity by different method. Although this is not related to lammps problem directly yet i thought that you would be so pleased to share your knowledge.
the things, I have confusion is that from the several literature I’ve seen that they are predicting the time evolution of pressure and temperature by averaging over the “fixed volume movable” bin and generally they take the position of the bin in the middle half of the entire sample to exclude the surface or near surface atoms which are not in equilibrium during shock run, but as i intend to track the shock front and want to study the phenomenon in the shock front then it is not sure that the bin is always just behind the front, when the front leads a appreciable distance from the user defined bin then is it conceptually correct to predict time evolution of P and T due to shock ?
Ray,
I was trying to implement the formula you suggested, Stress_xx = [ Sum(Force_xx of i) / (LyLz) ] + [ Sum(mass_i * Vx* Vx) / (LxLyLz)],
I have computed stress_xx, stress_yy and stress_zz and doing average ( i.e P = 1/3 (stress_xx +stress_yy + stress_zz )) to see if it same as I computed from the Rankin Hugonoit jump condt. P= density* Us * Up but it differs significantly, how can you explain that ?
I have calculated the “Force_xx” and “mass_i” of the region I required, by the fix " ave/spatial " and devided that two by area and volume of the bin respectively. Is there anything, I have missed. Please make your comment and thanking you for your valuable time.