Simulating shock in LAMMPS - about changing BCs and fix qeq/reax

Hello all,

I’m using LAMMPS-17Nov16 to simulate shock propagation in an all-atom polyurea model, with preliminary minimization and equilibration using an NPT ensemble, and have a few questions:

  • I apply periodic boundary conditions (BCs) in all directions during equilibration
    (‘boundary p p p’),
    but subsequently apply non-periodic BCs along the direction of shock propagation
    (‘change_box all boundary p p sf remap units box’).
    The original all-atom model was generated by someone else using Materials Studio, and I’m confident there are no bonds crossing the cell boundaries at the start (before minimization). Is there something I should check before changing the BCs?

  • The log file for the simulation contains the message “WARNING: Fix qeq/reax CG convergence failed after 200 iterations at 50047 step
    (…/fix_qeq_reax.cpp:722)” for multiple timesteps of shock propagation.
    I will try changing the velocity of the wall in the ‘fix wall/piston’ command (it’s 1.5 A/fs), and check the force-field using a different simulation (constant engineering strain), but is there a way to change the number of iterations mentioned?

  • This question is more general. In the “real” system of units in LAMMPS, the unit of pressure would be (1 kcal/mol-A / 1 A^2) or approx. 69,500 bar. Since the unit of pressure is called atm, do I understand I should divide a reported pressure of, say, 1800 atm by 69,500 to get the actual value?

I’d appreciate any help on these questions; I’m new to LAMMPS.

Regards,
Ramaswami.

  1. If the system has been equilibrated with NPT ensemble with no other reflective walls to prevent bonds from crossing PBC, it is hard to believe that “there are no bonds crossing the cell boundaries”. You should check the system carefully and possibly re-do the NPT equilibration with proper walls.

  2. If qeq/reax fails to converge with 200 steps, increasing the number of steps wouldn’t help at all. I believe the convergence problem is a result of dangling bonds at the boundaries.

  3. Pressure is a property of the entire system/box not per-atom. It doesn’t make sense to normalize the pressure. You can use fix ave/chuck to get localized pressures, and that can tell you the pressure near the piston and those ahead/behind the shock front.

Ray

Hello Ray,

Thank you for replying so fast.

  1. You’re right; the system expands during equilibration, comparing the box bounds before and after this step, so dangling bonds must be cut off when I change the boundary conditions (BCs). Before starting equilibration, I will try using ‘fix wall/lj93’ with the walls enclosing a space slightly larger than that indicated by the final box bounds.

  2. “Pressure is a property of the entire system/box not per-atom. It doesn’t make sense to normalize the pressure.”

Maybe I was unclear earlier; I’m writing thermodynamic data to the log file, including the pressure in each direction, and did understand pressure to be a property of the entire system/box. But I will see how to use ‘fix ave/chunk’ as you suggest.

Regards,

Ramaswami.