Nan Values on First Time Step of Shear Strain

Dear LAMMPS users,

I’m having trouble with a shear strain (in xy direction) simulation of an anisotropic crystalline unit cell (~70 atoms) in ReaxFF. This system has spent a total of 5.15 ns in ReaxFF (Liu parameter set) with 4 ns spent as an orthogonal box and 1.15 ns spent equilibrating in NPT at 300 K and 1 atm after changing the box over to triclinic. I noticed that nan values would occur for all thermodynamics except potential energy after the zero time step when I attempted to strain the system. The fixes in use are fix 1 all npt temp 300 300 100 xz 1 1 1000 yz 1 1 1000 and fix 2 all deform 1 xy erate 0.0000002. The time step is 0.1 fs. The version of LAMMPS I’m using is 02/29/2016.

I changed the thermo output frequency to output at every time step and the nan values begin on the first time step. The potential energy on the zero time step is -9366 and on the subsequent steps it’s -82. I read on the mailing list and in the manual that nan values are usually caused by a bad model (overlapping atoms, etc.). The potential energy value staying the same seems to indicate that something is wrong with the system. I don’t think it’s a bad model before applying the shear strain because it equilibrated for 1.15 ns with a triclinic box with no issues (thermo output was normal). Maybe the first little application of shear (fix deform xy) is pushing atoms on top of each other? Is there a way I can ensure that atoms do not overlap during the shearing? Is there another shearing approach I could try?

Thanks,
Will

Dear LAMMPS users,

I'm having trouble with a shear strain (in xy direction) simulation of an
anisotropic crystalline unit cell (~70 atoms) in ReaxFF. This system has
spent a total of 5.15 ns in ReaxFF (Liu parameter set) with 4 ns spent as an
orthogonal box and 1.15 ns spent equilibrating in NPT at 300 K and 1 atm
after changing the box over to triclinic. I noticed that nan values would
occur for all thermodynamics except potential energy after the zero time
step when I attempted to strain the system. The fixes in use are fix 1 all
npt temp 300 300 100 xz 1 1 1000 yz 1 1 1000 and fix 2 all deform 1 xy erate
0.0000002. The time step is 0.1 fs. The version of LAMMPS I'm using is
02/29/2016.

I changed the thermo output frequency to output at every time step and the
nan values begin on the first time step. The potential energy on the zero
time step is -9366 and on the subsequent steps it's -82. I read on the
mailing list and in the manual that nan values are usually caused by a bad
model (overlapping atoms, etc.). The potential energy value staying the same
seems to indicate that something is wrong with the system. I don't think
it's a bad model before applying the shear strain because it equilibrated
for 1.15 ns with a triclinic box with no issues (thermo output was normal).
Maybe the first little application of shear (fix deform xy) is pushing atoms
on top of each other? Is there a way I can ensure that atoms do not overlap
during the shearing? Is there another shearing approach I could try?

will,

have you tried the same simulation protocol with a different pair
style (on a similar system)? that way you can rule out whether there
is an issue with your protocol. also, 70 atoms sounds a bit smallish.
have you tried a larger system? finally, please check if the issue
persists when use the very latest patchlevel of LAMMPS? which of the
two reaxff implementations in LAMMPS are you using?

thanks,
    axel.

Will,

You might try switching from NPT to NVT to see if it’s related to your deform command or your NPT command. You might also try increasing your pressure coupling coefficients a few orders of magnitude to see if the box size is exploding/imploding.

-Ben

Will,

You might try switching from NPT to NVT to see if it's related to your
deform command or your NPT command. You might also try increasing your
pressure coupling coefficients a few orders of magnitude to see if the box
size is exploding/imploding.

ben,

this is a *very* important point. well spotted. there must not be a
fix npt and fix deform modifying the same box dimensions at the same
time. that will *definitely* lead to NaN pressure and everything else.
this should have been spotted with the first check that i suggested,
as this would happen with *any* pair style.

axel.

Ben and Axel,

Thinking on what you both said made me realize that my NPT command only set the xz and yz pressures when I should have set xz, yz, x, y, and z (when deforming in the xy direction). My simulation will take several days to run to completion, but it’s looking like the nan value issue has been fixed. Thank you both very much!

Thanks,
Will

I did some quick tests using npt and deform in the examples/melt example. It turns out that even without deform, any attempt to use fix npt with a shear component and without at least one active axial component triggers the recently added error message:

ERROR on proc 0: Non-numeric atom coords - simulation unstable (…/domain.cpp:510)

This is not something that we had tested before, but is a valid use case. There are a few parameters that are defined as averages over the active axial components, and these should default to zero when no axial components are active. After fixing these, the expected behavior is recovered. Patch is appended below and will be in the next release.

Aidan

[[email protected]…5834… src]$ svn diff
Index: fix_nh.cpp