Problems with Global Charge Redistribution

Dear all,

we are currently trying to simulate thermal oxidation of silicon nanowires using ReaxFF and LAMMPS 10 Aug 2015.
During our simulation we discovered that oxygen molecules far away from the nanowire become negatively charged.
This results in long-ranging repulsive behaviour with oxygen atoms attached to the wire, which leads to an unphysical slow-down of the oxidation process.

We were able to replicate this unusual charge redistribution using the following input file. Si.reaxFF is a file containing parameters by either Buehler et al. (PRL 96, 2006) or Kulkarni et al. (J. Phys. Chem. C 117, 2013). Both yield the abnormal charges even though the exact numbers differ from each other.

---------------------- Input File ----------------------

units real
atom_style full

boundary p p p

lattice diamond 5.431 origin 0.1 0.1 0.1 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block 0 50 0 50 0 50 units lattice

create_box 2 box

region silicon block 0 0.1 0 0.1 0 0.1 units lattice

create_atoms 2 region silicon

region oxygen block 24 25 24 25 24 25 units lattice

create_atoms 1 region oxygen

mass 1 28.06
mass 2 15.999
pair_style reax/c NULL
pair_coeff * * Si.reaxFF Si O

neighbor 2 bin
neigh_modify every 10 delay 0 check no

compute type all property/atom type

timestep 0.05

fix qeq all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
fix nve all nve

dump custom all custom 1 output.custom id c_type q x y z vx vy vz

run 0

---------------------- Output File ----------------------

ITEM: TIMESTEP

ITEM: NUMBER OF ATOMS
9
ITEM: BOX BOUNDS pp pp pp
0 271.55
ITEM: ATOMS id c_type q x y z vx vy vz
1 2 -0.328226 0.5431 0.5431 0.5431 0 0 0
2 1 0.0719583 130.887 130.887 130.887 0 0 0
6 1 -0.0070383 132.245 132.245 132.245 0 0 0
5 1 0.0370514 133.603 133.603 130.887 0 0 0
9 1 0.0507171 134.96 134.96 132.245 0 0 0
4 1 0.0370514 133.603 130.887 133.603 0 0 0
8 1 0.0507171 134.96 132.245 134.96 0 0 0
3 1 0.0370514 130.887 133.603 133.603 0 0 0
7 1 0.0507171 132.245 134.96 134.96 0 0 0

As you can see the charges are redistributed instantaneously. There seems to be no spatial cut-off during the charge equilibrium (QEq).

We already checked the available documentation and mailing list and found a similar problem described in http://lammps.sandia.gov/threads/msg29916.html. Unfortunately, no solution was available back in 2012. Is there a way to fix this problem in 2017?

Kind regards,
Georg Heinze & Florian Fuchs

The current version of LAMMPS and its ReaxFF implementations do not include
the ACKS2 algorithm required for proper charge equilibration. This may be
the source of your problem.

Jim Kress

Hi Georg,

This is not unusual and unexpected. On the contrary, this is the exact behavior of the Electronegativity Equalization (QEq) approached developed by Goddard and others. The QEq approach works by first obtaining a scaler value of the global electronegativity that includes contributions from all of the atoms in the system. Then partial charge on each of the atoms are then optimized so that per-atom electronegativity is equal to the global electronegativity. Inevitably with this approach oxygen atoms beyond the pair cutoff exchange charges with the Si atoms.

To handle charge redistributions properly, several methods have been developed such as SQE, PTQIE, and ACKS2. Unfortunately none of these methods are implemented in LAMMPS. The fix qeq/variants would be the ideal parent class for new implementations.

Hope this helps,
Ray