[lammps-users] Instability of reaxFF for ZrO2

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS-users (and Aidan in particular!),

I'm trying to do some simulations of cubic / tetragonal ZrO2 using the
reaxFF potential, but LAMMPS crashes after 10 fs with the message
Too many bonds on atom. Increase MBONDDEF.
I've tried various tricks to stablize it, and in the best case have
managed to reach 10 ps, but inevitably LAMMPS crashes with the same
error message.

The mailing list suggests the following possible errors, but none of
them solve the problem:
- Unphysical configuration: I've simulated this configuration (cubic or
tetragonal unit cell of ZrO2) with a Buckingham potential and get
sensible answers.
- Minimize the unit cell: I've tried "minimize" for 200 iterations
before starting NVT, but sooner or later I get the same error message
- Reduce the number of atoms per processor: The simulation cell is 3x3x3
unit cells, each of which contains 12 atoms, and I've tried both 1 and 8
processors.
- Remove the drag term from NVT and NPT: The system crashes in either
case. Incidentally, I'm using LAMMPS downloaded on 23 March 2010, before
NVT chains were added, but surely this shouldn't be so critical?

As a desparate attempt I increased the lattice parameters by 50% (which
is unphysical) and performed NVT for 875 fs followed by NPT for 875 fs.
Things looked quite stable, but then the simulation crashed with the
same error message.

Below is my input script. At the beginning of the simulation the
temperature tends to rise dramatically (to 9000 K), so I suppress the
fluctuations using temp/rescale. (I know it's not advisable to do this
together with NVT, but it's certainly not the cause of the above error
message.) The force field parameters are taken from Adri C.T. van Duin
et al., J. Phys. Chem. A 112, 139 (2008), by simply copying the
Supplementary Information into a text file. The timestep of 0.035 fs was
recommended by these authors in a subsequent paper, from which I also
took the initial charges of O and Zr.

I'd be very grateful for any suggestions,

Thanks,
Philip

boundary p p p
units real
atom_style charge
lattice custom 6 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.01 &
              origin 0.0 0.0 0.0 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
basis 0.5 0.0 0.5 basis 0.0 0.5 0.5 &
              basis 0.25 0.25 0.25 basis 0.75 0.25 0.25 basis 0.25 0.75
0.25 basis 0.75 0.75 0.25 basis 0.25 0.25 0.75 basis 0.75 0.25 0.75
basis 0.25 0.75 0.75 basis 0.75 0.75 0.75
region mybox block 0 3 0 3 0 3
create_box 2 mybox
create_atoms 2 box basis 1 2 basis 2 2 basis 3 2 basis 4 2 &
              basis 5 1 basis 6 1 basis 7 1 basis 8 1 &
              basis 9 1 basis 10 1 basis 11 1 basis 12 1

group O type 1
group Zr type 2
set group O charge -1
set group Zr charge 2
mass 1 16
mass 2 91.22

neigh_modify one 20000 page 200000
timestep 0.035
pair_style reax 15 1.0e-6
pair_coeff * * ffield.reax 3 5
thermo_style custom step temp press pxx pyy pzz pe etotal evdwl ecoul
lx ly lz
thermo_modify flush yes
thermo 10

velocity all create 2000 891286 mom yes rot yes dist Gaussian
fix RESCALE all temp/rescale 1 2000 200070 1
fix NVT all nvt 2000 2000 0.35
run 25000
unfix NVT
unfix RESCALE

fix NPT all npt 2000 2000 0.35 aniso 0.0 0.0 0.0 0.0 0.0 0.0
10
run 25000
unfix NPT

2

I'll let Aidan respond, but there should be no reason to make MBONDDEF
unnaturally large. If you exceed that limit, I think that would
indicate you have a poorly
parameterized potential and the geometry of the atoms is doing something
odd.

Steve

3

You did not say what value of MBONDDEF you were using. For some metal
systems, the number of "bonds" can be much larger than for covalently bonded
materials. Since you have a very small system, make NATDEF small and
MBONDDEF big. You should also try comparing your minimized structures and
energies with those in the paper, to see if you are really running the same
ReaxFF potential.

Aidan

4

Thanks for the very fast advice. I'm using MBONDDEF = 20 and will try increasing it (unfortunately I don't have compile rights on the cluster I'm using, so I can't do this immediately).

However, I suspect there's something more fundamentally wrong, since in some circumsatances the the total energy of my starting configuration depends on the number of processors I use. For the simulation with 3x3x3 unit cells in the script I posted there's no difference, but when I simulate 8x8x8 unit cells the total energy differs by about 3% depending on whether I use 1 or 8 processors. I've checked explicitly that the atom positions are the same in each case, and there's no dependence on the neighbor list settings, so surely the initial energy should be the same?

Could you suggest a diagnosis for this problem? Thanks for your help,
Philip

5

LAMMPS, including ReaxFF, is generally very consistent w.r.t. number of
processors. Given an initial geometry, the processor count should only
affect the forces and energies in the last couple of digits. This is due to
reordering of arithmetic operations. Subsequent timesteps will show larger
difference, as the trajectories will slowly diverge.

The charge equilibration is an iterative operation, so it can introduce
additional variation, if the small arithmetic errors cause a difference in
iteration counts.

Feel free to send me an example of this problem.

Aidan

6

Aidan,

Thanks for your offer, which I'll gladly accept. I'm familiar with LAMMPS' usual near-independence of the number of processors, which is why I find this result so disturbing. Attached you find my input script, the ffield.reax (taken from the Supplementary Information to Adri C.T. van Duin et al., J. Phys. Chem. A 112, 139 (2008)) and the results of running it on 1 or 8 processors. LAMMPS was compiled using the default settings for reaxFF.

The script generates cubic ZrO2 at T=0 and changes the volume (without relaxing the atomic positions), which is related to Fig.4 of the above paper (although I don't normalize the scale on the axes). This should be continuous, but if you plot the energy vs. the lattice constant you'll see that it's not.
  (for gnuplot use
  set xrange [5.18:5.32]
  plot "log.lammps.n1" u 7:4, "log.lammps.n8" u 7:4
This might indicate that the reax cut-off is too small, but when I increase it, say from 15 to 18, LAMMPS generates the error message
  Setting up run ...
  ^@^@-^A^@Fatal: Unknown bond in molecule

This is the first time I've used reaxFF so I've got no experience of its convergence behavior. One question: reaxFF adapts the charges etc. to the local environment, but does this happen self-consistently within each timestep or do the values at one timestep reflect the configuration of the previous timestep?

Looking forward to hearing your comments,
Philip

log.lammps.n1 (16.9 KB)

log.lammps.n8 (16.9 KB)

ffield.reax (17.5 KB)

minimize_ZrO2.run (2.44 KB)

7

Howard,

Thank you for the nice example. The calculated energy per atom (use
thermo_modify norm yes to see this easily) for the initial geometry is
always the same, independent of system size and processor count, as it
should be....

...except for the case you mention, 8x8x8, 1 processor, and a very large
value of hbcut=15.0. The number of Verlet pairs (related to number of LAMMPS
neighbors, but bigger) was very large, and exceeded the statically allocated
array. Some memory checking that would detect this was disabled when we
ported ReaxFF to LAMMPS. I have put those checks back in, so now, when you
run this problem, you will see:

ERROR on proc 0: reax_defs.h::NNEIGHMAXDEF too small

I also put in a check on NATDEF that was missing. This will be released as a
patch shortly. In the mean time, stay away from those huge neighbor lists
:-).

Aidan

8

Aidan,

Thanks for sorting out the problem. I chose these parameters based
on experience with other potentials, but clearly they were
pathological. When I follow your advice it works fine.

According to Adri van Duin's documentation of reaxFF, hbcut determines
the hydrogen-bond cut-off. I'd assumed that LAMMPS uses it as some
sort of global cut-off for pair interactions, but this appears not
to be the case: I'm studying ZrO2 (i.e. no hydrogen), and my results
are completely independent of hbcut. Apparently I can set hbcut=1
with no ill-effects, which presumably reduces the size of arrays
which are allocated but never actually used because I have no
hydrogen. Is this correct? If so, does this mean that I can only
influence the precision and hence the computational cost of the force
calculations via the PPPM accuracy?

Philip

9

According to Adri van Duin's documentation of reaxFF, hbcut determines
the hydrogen-bond cut-off. I'd assumed that LAMMPS uses it as some
sort of global cut-off for pair interactions, but this appears not
to be the case: I'm studying ZrO2 (i.e. no hydrogen), and my results
are completely independent of hbcut. Apparently I can set hbcut=1
with no ill-effects, which presumably reduces the size of arrays
which are allocated but never actually used because I have no
hydrogen. Is this correct?

Yes.

If so, does this mean that I can only
influence the precision and hence the computational cost of the force
calculations via the PPPM accuracy?

The second parameter is the charge equilibration precision. It defaults to
1.0e-6. It controls how many iterations are performed by the CG solver every
timestep. You can experiment with larger or smaller values, the goal being
to find the largest value that still produces accurate results.

Aidan

Philip

From: Thompson, Aidan [mailto:[email protected]]
Sent: Wednesday, June 30, 2010 6:09 PM
To: Howell, Philip Clissold; Steve Plimpton
Cc: [email protected]
Subject: Re: [lammps-users] Instability of reaxFF for ZrO2

Howard,

Thank you for the nice example. The calculated energy per atom (use
thermo_modify norm yes to see this easily) for the initial geometry is
always the same, independent of system size and processor count, as it
should be....

...except for the case you mention, 8x8x8, 1 processor, and a very large
value of hbcut=15.0. The number of Verlet pairs (related to number of LAMMPS
neighbors, but bigger) was very large, and exceeded the statically allocated
array. Some memory checking that would detect this was disabled when we
ported ReaxFF to LAMMPS. I have put those checks back in, so now, when you
run this problem, you will see:

ERROR on proc 0: reax_defs.h::NNEIGHMAXDEF too small

I also put in a check on NATDEF that was missing. This will be released as a
patch shortly. In the mean time, stay away from those huge neighbor lists
:-).

Aidan

--
      Aidan P. Thompson
      01435 Multiscale Dynamic Materials Modeling
      Sandia National Laboratories
      PO Box 5800, MS 1322 Phone: 505-844-9702
      Albuquerque, NM 87185 Fax : 505-845-7442
      E-mail:[email protected] Cell : 505-550-2614

From: "Howell, Philip Clissold" <[email protected]...>
Date: Mon, 28 Jun 2010 10:05:41 -0600
To: Aidan Thompson <[email protected]>, Steve Plimpton <[email protected]>
Cc: <[email protected]>
Subject: RE: [lammps-users] Instability of reaxFF for ZrO2

Aidan,

Thanks for your offer, which I'll gladly accept. I'm familiar with LAMMPS'
usual near-independence of the number of processors, which is why I find this
result so disturbing. Attached you find my input script, the ffield.reax
(taken from the Supplementary Information to Adri C.T. van Duin et al., J.
Phys. Chem. A 112, 139 (2008)) and the results of running it on 1 or 8
processors. LAMMPS was compiled using the default settings for reaxFF.

The script generates cubic ZrO2 at T=0 and changes the volume (without
relaxing the atomic positions), which is related to Fig.4 of the above paper
(although I don't normalize the scale on the axes). This should be
continuous,
but if you plot the energy vs. the lattice constant you'll see that it's not.
  (for gnuplot use
  set xrange [5.18:5.32]
  plot "log.lammps.n1" u 7:4, "log.lammps.n8" u 7:4
This might indicate that the reax cut-off is too small, but when I increase
it, say from 15 to 18, LAMMPS generates the error message
  Setting up run ...
  ^@^@-^A^@Fatal: Unknown bond in molecule

This is the first time I've used reaxFF so I've got no experience of its
convergence behavior. One question: reaxFF adapts the charges etc. to the
local environment, but does this happen self-consistently within each
timestep
or do the values at one timestep reflect the configuration of the previous
timestep?

Looking forward to hearing your comments,
Philip

From: Thompson, Aidan [mailto:[email protected]]
Sent: Monday, June 28, 2010 5:04 PM
To: Howell, Philip Clissold; Steve Plimpton
Cc: [email protected]
Subject: Re: [lammps-users] Instability of reaxFF for ZrO2

LAMMPS, including ReaxFF, is generally very consistent w.r.t. number of
processors. Given an initial geometry, the processor count should only
affect the forces and energies in the last couple of digits. This is due to
reordering of arithmetic operations. Subsequent timesteps will show larger
difference, as the trajectories will slowly diverge.

The charge equilibration is an iterative operation, so it can introduce
additional variation, if the small arithmetic errors cause a difference in
iteration counts.

Feel free to send me an example of this problem.

Aidan

--
      Aidan P. Thompson
      01435 Multiscale Dynamic Materials Modeling
      Sandia National Laboratories
      PO Box 5800, MS 1322 Phone: 505-844-9702
      Albuquerque, NM 87185 Fax : 505-845-7442
      E-mail:[email protected] Cell : 505-550-2614

From: "Howell, Philip Clissold" <[email protected]...>
Date: Mon, 28 Jun 2010 08:55:27 -0600
To: Aidan Thompson <[email protected]>, Steve Plimpton <[email protected]>
Cc: <[email protected]>
Subject: RE: [lammps-users] Instability of reaxFF for ZrO2

Thanks for the very fast advice. I'm using MBONDDEF = 20 and will try
increasing it (unfortunately I don't have compile rights on the cluster I'm
using, so I can't do this immediately).

However, I suspect there's something more fundamentally wrong, since in some
circumsatances the the total energy of my starting configuration depends on
the number of processors I use. For the simulation with 3x3x3 unit cells in
the script I posted there's no difference, but when I simulate 8x8x8 unit
cells the total energy differs by about 3% depending on whether I use 1 or 8
processors. I've checked explicitly that the atom positions are the same in
each case, and there's no dependence on the neighbor list settings, so
surely
the initial energy should be the same?

Could you suggest a diagnosis for this problem? Thanks for your help,
Philip

From: Thompson, Aidan [mailto:[email protected]]
Sent: Friday, June 25, 2010 5:12 PM
To: Steve Plimpton; Howell, Philip Clissold
Cc: [email protected]
Subject: Re: [lammps-users] Instability of reaxFF for ZrO2

You did not say what value of MBONDDEF you were using. For some metal
systems, the number of "bonds" can be much larger than for covalently bonded
materials. Since you have a very small system, make NATDEF small and
MBONDDEF big. You should also try comparing your minimized structures and
energies with those in the paper, to see if you are really running the same
ReaxFF potential.

Aidan

--
      Aidan P. Thompson
      01435 Multiscale Dynamic Materials Modeling
      Sandia National Laboratories
      PO Box 5800, MS 1322 Phone: 505-844-9702
      Albuquerque, NM 87185 Fax : 505-845-7442
      E-mail:[email protected] Cell : 505-550-2614

From: Steve Plimpton <[email protected]>
Date: Fri, 25 Jun 2010 08:45:23 -0600
To: "Howell, Philip Clissold" <[email protected]...>
Cc: <[email protected]>, Aidan Thompson
<[email protected]>
Subject: Re: [lammps-users] Instability of reaxFF for ZrO2

I'll let Aidan respond, but there should be no reason to make MBONDDEF
unnaturally large. If you exceed that limit, I think that would
indicate you have a poorly
parameterized potential and the geometry of the atoms is doing something
odd.

Steve

Dear LAMMPS-users (and Aidan in particular!),

I'm trying to do some simulations of cubic / tetragonal ZrO2 using the
reaxFF potential, but LAMMPS crashes after 10 fs with the message
Too many bonds on atom. Increase MBONDDEF.
I've tried various tricks to stablize it, and in the best case have
managed to reach 10 ps, but inevitably LAMMPS crashes with the same
error message.

The mailing list suggests the following possible errors, but none of
them solve the problem:
- Unphysical configuration: I've simulated this configuration (cubic or
tetragonal unit cell of ZrO2) with a Buckingham potential and get
sensible answers.
- Minimize the unit cell: I've tried "minimize" for 200 iterations
before starting NVT, but sooner or later I get the same error message
- Reduce the number of atoms per processor: The simulation cell is 3x3x3
unit cells, each of which contains 12 atoms, and I've tried both 1 and 8
processors.
- Remove the drag term from NVT and NPT: The system crashes in either
case. Incidentally, I'm using LAMMPS downloaded on 23 March 2010, before
NVT chains were added, but surely this shouldn't be so critical?

As a desparate attempt I increased the lattice parameters by 50% (which
is unphysical) and performed NVT for 875 fs followed by NPT for 875 fs.
Things looked quite stable, but then the simulation crashed with the
same error message.

Below is my input script. At the beginning of the simulation the
temperature tends to rise dramatically (to 9000 K), so I suppress the
fluctuations using temp/rescale. (I know it's not advisable to do this
together with NVT, but it's certainly not the cause of the above error
message.) The force field parameters are taken from Adri C.T. van Duin
et al., J. Phys. Chem. A 112, 139 (2008), by simply copying the
Supplementary Information into a text file. The timestep of 0.035 fs was
recommended by these authors in a subsequent paper, from which I also
took the initial charges of O and Zr.

I'd be very grateful for any suggestions,

Thanks,
Philip

boundary p p p
units real
atom_style charge
lattice custom 6 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.01 &
             origin 0.0 0.0 0.0 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
basis 0.5 0.0 0.5 basis 0.0 0.5 0.5 &
             basis 0.25 0.25 0.25 basis 0.75 0.25 0.25 basis 0.25 0.75
0.25 basis 0.75 0.75 0.25 basis 0.25 0.25 0.75 basis 0.75 0.25 0.75
basis 0.25 0.75 0.75 basis 0.75 0.75 0.75
region mybox block 0 3 0 3 0 3
create_box 2 mybox
create_atoms 2 box basis 1 2 basis 2 2 basis 3 2 basis 4 2 &
             basis 5 1 basis 6 1 basis 7 1 basis 8 1 &
             basis 9 1 basis 10 1 basis 11 1 basis 12 1

group O type 1
group Zr type 2
set group O charge -1
set group Zr charge 2
mass 1 16
mass 2 91.22

neigh_modify one 20000 page 200000
timestep 0.035
pair_style reax 15 1.0e-6
pair_coeff * * ffield.reax 3 5
thermo_style custom step temp press pxx pyy pzz pe etotal evdwl ecoul
lx ly lz
thermo_modify flush yes
thermo 10

velocity all create 2000 891286 mom yes rot yes dist Gaussian
fix RESCALE all temp/rescale 1 2000 200070 1
fix NVT all nvt 2000 2000 0.35
run 25000
unfix NVT
unfix RESCALE

fix NPT all npt 2000 2000 0.35 aniso 0.0 0.0 0.0 0.0 0.0 0.0
10
run 25000
unfix NPT

__________________________

Dr. Philip Howell

Siemens AG
Corporate Technology
CT T DE HW 2
Otto-Hahn-Ring 6
81739 Munich, Germany
Tel.: +49 (89) 636-80225
Fax: +49 (89) 636-48131
mailto:philip.howell@…528…

Siemens Aktiengesellschaft: Chairman of the Supervisory Board: Gerhard
Cromme; Managing Board: Peter Loescher, Chairman, President and Chief
Executive Officer; Wolfgang Dehen, Heinrich Hiesinger, Joe Kaeser,
Barbara Kux, Hermann Requardt, Siegfried Russwurm, Peter Y. Solmssen;
Registered offices: Berlin and Munich, Germany; Commercial registries:
Berlin Charlottenburg, HRB 12300, Munich, HRB 6684; WEEE-Reg.-No. DE
23691322

Important notice: This e-mail and any attachment thereof contain
corporate proprietary information. If you have received it by mistake,
please notify us immediately by reply e-mail and delete this e-mail and
its attachments from your system. Thank you.

----------------------------------------------------------------------------->>