Generating Shock by Shrinking Periodic Boundaries


I have been trying to conduct shock simulations by shrinking the periodic boundary along the longest dimension of a cuboid shaped domain (as used in this paper: The paper states:

A planar shock wave was initiated at time zero by shrinking the periodic computational cell length Lx according to Lx(t)=Lx(0)-2xUpxt ,while keeping the lengths Ly and Lz fixed.

This approach is expected to create two shock fronts that approach the center of the domain from the two periodic boundaries. My input script is as follows, where I am using ‘fix deform’ to shrink the periodic boundary at a certain velocity. However, no shock waves are created and the domain undergoes homogeneous shrinkage along z-axis at a constant rate. I understand that I am missing something fundamental here, since some kind of inhomogeneous deformation is needed to generate a shock front. Any help would be appreciated.

units real
boundary p p p

atom_style charge

read_restart restart.b4c.108000
pair_style reax/c lmp_control
pair_coeff * * ffield.reax C B B B B B B B B B B B C B C

neigh_modify every 5 delay 0 check yes

timestep 0.1

thermo 50

restart 4000 restart.b4c

compute p all stress/atom NULL
dump data all custom 200 id type mass x y z vx vy vz c_p[]*

fix 1 all qeq/reax 1 0.0 10.0 1e-6 reax/c
fix 2 all nve
fix 3 all deform 1 z vel -0.07
run 10000


Amith Adoor

The text you quote describes an instantaneous change of the box dimensions (at time t=0), thus fix deform is not the correct comand because it does a gradual change over the course of the trajectory.

Instead you should be using the change_box command after equilibration but before you do your production run.


Thank you for your input! I implemented your suggestion and found that it did initiate a shock wave emanating from the periodic boundary conditions. However, I also realized that ‘fix deform’ can be used to generate a shockwave. The command " fix 3 all deform 1 z vel -0.07" should be changed to " fix 3 all deform 1 z vel -0.07 remap v". The default setting for remap option is x. As per the documentation,

setting remap to x forces the atoms to deform via an affine transformation that exactly matches the box deformation. Note that though the atoms are effectively “moving” with the box over time, it is not due to their having a velocity that tracks the box change, but only due to the remapping

remap v is typically appropriate for fluids, where you want the atoms to respond to the change in box size/shape on their own and acquire a velocity that matches the box change, so that their motion will naturally track the box without explicit remapping of their coordinates.

1 Like

Please note that the change_box command also has a remap option. It is off by default, which makes its default behavior different from fix deform. Without remap all coordinates remain unchanged when the box dimensions are altered thus effectively slamming atoms on the opposite sides of the box into each other due to PBC. But for triggering a shock wave, this default behavior is what you want and I am confident that this is what the authors of the paper that you were quoting did.

1 Like