Velocity zeros out for atoms in MD with ReaxFF

LAMMPS
1

I’m trying to run an MD simulation of a gaseous (with various densities) system with some metal and non-metal atoms included using ReaxFF potential. I start with individual atoms randomly distributed within periodic cell. Overtime atoms are expected to aggregate into clusters.

The problem I encountered was that after some clusters were formed, which happened pretty quickly, the remaining individual atoms stopped moving (speed was reduced dramatically, or it was zeroed out completely). I have tried different initial conditions (coordinates, velocities, etc.), various thermostats, different types of atoms, and even different parameter sets with no change in behavior.
This behavior is not apparent if one works with nanoparticles instead of gas, as atoms within clusters continue vibrating.

I also tried several versions of LAMMPS (latest 2024, LAMMPS-64bit-2Aug2023,LAMMPS-64bit-15Sep2022), all on Windows from precompiled executables downloaded from the official LAMMPS repository.
In the older version, I also tried the reax/c package instead.

I’m not sure where the problem can originate.

Below is the LAMMPS script I used (one of the iterations, as I mentioned, many things were already tried):

units real
atom_style charge
dimension 3
boundary p p p
atom_modify sort 0 0
atom_modify map array
neigh_modify one 10000
neigh_modify page 1000000

read_data data.ClusterForm

pair_style reaxff NULL
pair_coeff * * …/ffield.reax.Ni2010 Ni C
fix _qeq all qeq/reax 1 0.0 10.0 1.0e-6 reaxff
set atom * charge 0

thermo 500
thermo_modify flush yes
dump _trj all custom 500 ClusterForm.lammpstrj id element x y z
dump_modify _trj element Ni C
dump_modify _trj sort id

timestep 0.05

velocity all create 700 59267 dist gaussian units box

fix _dyn all nvt temp 700 700 1.0
run 500000
write_data data.ClusterForm_fin nocoeff
unfix _dyn

#-----------------------------------------
clear

2

I don’t know what your problem is, but a time constant of 1 fs for the Nose-Hoover thermostat seems like a crazy value. Any reason for that choice?

3

That was part of the thermostat setting tests. I wasn’t sure if the system’s dynamics were not resolved properly, causing this weird behavior.
Using general guidelines from the LAMMPS documentation:

A good choice for many models is a Tdamp of around 100 timesteps.

I tested ts (timestep) = (0.05, 0.1, 0.25, 0.5) and Tdamp = (20ts, 50ts, 100ts, 150ts, 200ts) for each ts value.
This particular script was for ts=0.05 and Tdamp=20ts=1.

4

And when you say:

what about the temperature of the system? Is it not 700K ?

5

What happens if you don’t have a thermostat, i.e. only fix nve? If the system blows up, e.g. lost atoms then you know something is very wrong. A thermostat can mask issues by putting on the brakes so hard that everything else slows down.

1 Like
6

what about the temperature of the system? Is it not 700K ?

The temperature oscillates around 700.
Individual atoms stop/slow down, but that is compensated by atomic vibrations within clusters.

7

What happens if you don’t have a thermostat, i.e. only fix nve ?

You are onto something. No atoms were lost, but…
I have tested this for two densities: (1) 120 atoms within a 30x30x30 periodic cell and (2) 120 atoms in a 60x60x60 cell. In the first case, the temperature went crazy immediately (as soon as particles started interacting), and I terminated the computation at 10000K. For the second case, it was delayed due to lower density and less frequent atomic collisions, but eventually went in the same way (max temp of 3500K).
Does it tell us something?

8

It means your system is unstable. You could have overlapping atom positions that need to be relaxed, e.g. by minimization minimize command — LAMMPS documentation, or bad force field parameters. The thermostat works by taking out energy from the system by reducing the velocities. If some atoms are exploding apart, then the thermostat slows down everything, not just the fast atoms.

9

See this related post: Why non-reacting atoms (part ethylene) move much slowly in the whole?

10

It means your system is unstable. You could have overlapping atom positions that need to be relaxed…

There are no overlapping atoms. The original configuration was used for VASP optimization and MD with no problems. I also tested it with many different starting coordinates. I even ran an equilibration run with repulsive Buckingham potential to make sure that velocities were decently distributed. The same problems occur every time.

…or bad force field parameters.

I also tested it with several ReaxFF parameter sets, including:
(1) https://doi.org/10.1021/jp9035056 , using Ni+C atoms
(2) Redirecting , using Zn+N, Zn+H
(3) https://doi.org/10.1021/acscatal.5b01766 , Fe+S, Cr+C, Cr+S
The same behavior every time.

11

You don’t seem to be conserving energy for some reason. Perhaps your timestep is too large? When all else fails, visualizing the simulation can be very helpful and is what I would recommend.

12

When all else fails, visualizing the simulation can be very helpful and is what I would recommend.

I agree, and this is always my first sanity check.
All atoms start moving normally, but as soon as small clusters start forming, all individual atoms stop, while atom vibrations within clusters continue. From your questions, I concluded that the kinetic energy of vibrations overwhelms the thermostat. But there is no clue what is causing it.
I will try an even smaller timestep, but 0.05fs seems on the small side already.

13

The thermostat time constant corresponds to the coupling frequency.

I fail to see the justification of a thermostat when simulating a gaseous(?) system and particularly a Nose-Hoover thermostat which was designed for bulk liquids.
What bulk do your atoms/molecules couple to? Also the mechanism of exchanging kinetic energy is different between gas and liquid. Gas is far more compressible and molecules/atoms are expected to move freely and undisturbed between collisions (which makes a langevin thermostat and even worse choice than Nose-Hoover for a gas after the system is properly equilibrated).

14

I think for a ReaxFF based system where atoms start clustering an increase of the potential energy is to be expected as you have on average more atoms within the cutoff radius for the clusters than for the distributed system and you don’t have accounting of energy for atoms outside the cutoff (unlike with long-range coulomb).

15

I fail to see the justification of a thermostat when simulating a gaseous(?) system and particularly a Nose-Hoover thermostat which was designed for bulk liquids.

Thank you for your clarification.
What would be your recommendation for the system where we expect and look for the cluster formation, so gas → liquid/solid transformation? Essentially, by the end of the simulation, the entirety of the system is a single nanoparticle.

16

This is tricky business since you are not doing a simulation in equilibrium.
When your atoms clump together, their kinetic energy has to go somewhere. Also you may free some potential energy due to the assembly of the particles. That may translate into even larger kinetic energy of the particles so they accelerate, or so much internal kinetic energy that the particles may break apart again.

In the gas phase there is no way to get rid of that energy unless you have some collisions with a material that can absorb the energy (and then dissipate it into its bulk).

If you have just one large particle, then you could end up with a so-called “flying-icecube” scenario.

You have to look at what happens in equivalent experiments. How do the (large) nano particles behave after they formed? How long do they stay around before they make contact?

I am certain that there is some published literature on similar scenarious.

1 Like