NaN after long simulation times

Just wondering if I could get some ideas/feedback on what might be going on here. I'm in the process of setting up my initial simulations using lammps version 15May15 to look at the MSD of poly(dimethylsiloxane) using PCFF and employing fix shake to constrain hydrogen with a timestep of 1 fs. I initially perform an NVT equilibration cycle which runs for 100,000 timesteps and the energies, temperature, etc. appear to be behaving appropriately. Then I begin a production run to get the MSD, again using an NVT ensemble, but now for 2,000,000 timesteps. Once again energies, temperature, etc. initially appear to be behaving normally without any issues, however after a significant simulation time I begin to see NaN in the log output (see below). This seems to occur at random, i.e. not necessarily at the same timestep if I rerun the simulation, but usually occur after ~750,000 timesteps. Below are portions of the input deck and log file. Any help is much appreciated, hopefully it’s not some trivial issue that went over my head, if so I apologize in advance.

Thanks for the feedback,
Stefan

^^^^^^^^^^^^^^^^^
Input deck:

units real
atom_style full
boundary p p p

neighbor 1.0 bin
neigh_modify delay 10 check yes

# Using PCFF
pair_style lj/class2/coul/long 12.0
pair_modify mix sixthpower tail yes
bond_style class2
angle_style class2
dihedral_style class2
improper_style class2
special_bonds lj 0 0 1 coul 0 0 1
kspace_style ewald 1.0e-6

read_data data.pdms

#Equilibrate
timestep 1.00
thermo 10000
thermo_style custom step temp press vol density etotal pe emol evdwl ecoul elong
velocity all create 298 152558 sum no rot yes mom yes dist gaussian loop local
fix f0 all shake 1.0e-6 100 10000 m 1.00 a 11
fix f1 all nvt temp 298 298 100.0 tchain 5 tloop 5
run 100000
unfix f1

#Production run
reset_timestep 0
thermo 10000
compute c1 all chunk/atom molecule
compute msd1 all msd/chunk c1
fix fmsd1 all ave/time 1000 10 10000 c_msd1 file pdms-msd.dat mode vector
dump d1 all custom 10000 dump.pdms id mol type x y z xu yu zu
fix f3 all nvt temp 298 298 100.0
run 2000000

^^^^^^^^^^^^^^^^
Log file output:

#Equilibration
.....................
.....................
Neighbor list info ...
  1 neighbor list requests
  update every 1 steps, delay 10 steps, check yes
  master list distance cutoff = 13
SHAKE stats (type/ave/delta) on step 0
  3 1.101 1.3507e-06
  6 1.101 5.40316e-07
  11 107.66 6.33878e-05
Memory usage per processor = 20.4906 Mbytes
Step Temp Press Volume Density TotEng PotEng E_mol E_vdwl E_coul E_long
       0 377.18698 -1232.7704 78314.941 0.89216771 -41585.029 -46696.2 -10322.608 -2895.2825 -16248.132 -17230.177
SHAKE stats (type/ave/delta) on step 10000
  3 1.101 1.19264e-07
  6 1.101 5.74244e-08
  11 107.66 4.13521e-06
   10000 295.38762 -531.14887 78314.941 0.89216771 -42700.538 -46703.266 -10275.729 -2912.8268 -16282.927 -17231.782
SHAKE stats (type/ave/delta) on step 20000
  3 1.101 1.02566e-06
  6 1.101 1.68047e-07
  11 107.66 1.24696e-05
   20000 300.70857 91.982684 78314.941 0.89216771 -42698.822 -46773.652 -10347.48 -2901.048 -16293.048 -17232.077
...................
...................
...................
SHAKE stats (type/ave/delta) on step 90000
  3 1.101 3.52165e-07
  6 1.101 1.90757e-07
  11 107.66 2.37887e-05
   90000 296.80783 -107.37605 78314.941 0.89216771 -42675.612 -46697.584 -10266.078 -2884.8488 -16315.522 -17231.135
SHAKE stats (type/ave/delta) on step 100000
  3 1.101 4.04882e-07
  6 1.101 2.53903e-07
  11 107.66 2.48655e-05
  100000 302.49349 28.625901 78314.941 0.89216771 -42486.528 -46585.545 -10179.984 -2868.1142 -16306.536 -17230.911
Loop time of 1512.78 on 24 procs (24 MPI x 1 OpenMP) for 100000 steps with 5755 atoms

Have you tried running plain NVE to check for energy conservation over an equivalent simulation time? Let’s hope one response will help you to gather some more from other participants with more potential insight on your problem.
Carlos

Just wondering if I could get some ideas/feedback on what might be going
on here. I'm in the process of setting up my initial simulations using
lammps version 15May15 to look at the MSD of poly(dimethylsiloxane) using
PCFF and employing fix shake to constrain hydrogen with a timestep of 1 fs.
I initially perform an NVT equilibration cycle which runs for 100,000
timesteps and the energies, temperature, etc. appear to be behaving
appropriately. Then I begin a production run to get the MSD, again using an
NVT ensemble, but now for 2,000,000 timesteps. Once again energies,
temperature, etc. initially appear to be behaving normally without any
issues, however after a significant simulation time I begin to see NaN in
the log output (see below). This seems to occur at random, i.e. not
necessarily at the same timestep if I rerun the simulation, but usually
occur after ~750,000 timesteps. Below are portions of the input deck and
log file. Any help is much appreciated, hopefully it’s not some trivial
issue that went over my head, if so I apologize in advance.

​the fist thing that i would worry about, is whether you get "dangerous
builds", because your neighbor list settings with a rather small neighbor
list skin​ of 1.0 (for typical size atoms in real units), and your pretty
large neighbor list rebuild delay of 10 time steps.

what you describe is something that can happen when atoms get too close.
that can be due to too large a time step, too soft a potential or atoms
"suddenly" appearing inside the cutoff range where they already have
significant interactions. a typo in the potential parameters or incorrectly
converted units is what has happened to me most frequently in the past.

axel

I’ll add that you could try re-running starting at step 1320000 (good)

for 10K steps to 1330000 when things are bad, and print

out thermo info every step. You may well see things blow

up within a small number of steps at some intermediate

time.

Steve