Problems with charged systems on pair style reax

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I’ve been trying to perform some calculations on charged systems using the pair style reax, but the potential seems to be ignoring the charges I input via the data file. Despite the value of the initial charge on my input the charges are calculated as if there was no net charge on the system. Can anyone help me on this? My data file (in charge format) reads like this:
162 atoms
0 bonds
0 angles
0 dihedrals
0 impropers
2 atom types
0 bond types
0 angle types
0 dihedral types
0 improper types
-50.419500 49.580500 xlo xhi
-50.839500 49.160500 ylo yhi
-0.903500 3.803500 zlo zhi

# Pair Coeffs

2

ReaxFF performs charge equilibration, so it resets the
charge on each atom every timestep. Thus your initial
settings are effectively ignored.

Steve

3

Steve,

Thanks for your reply. Is it then impossible for ReaxFF to perform calculations in systems with non zero net charge? Isn’t there anyway to equilibrate the extra charge?

Eric

4

Aidan might know the answer to this - I think it is an ill-defined problem.
I believe the charge equil step in ReaxFF will reset the net charge
to 0.0 every time it is invoked.

Steve

5

Yes, it is possible to do charge equilibration with net charge, but the
pair_style reax implementation in LAMMPS does not support it. I don't
think it would be too hard to turn it on, but it is not something that
people generally want. You could also try using pair_style reax/c with
qeq/reax, which might handle charged systems.

6

Thanks for the help. I tried using reax/c with qeq but it also overran my charge input.

7

pair reax and reax/c formulate the energy minimization of the
electrostatic energy using two different linear sovles that are
mathematically equivalent. In both cases, a charge conservation equation
is included, which constrains the total charge to match some target value,
within the tolerance specified. This target value is currently zero for
both reax and reax/c, but I think it would not be hard to change the
target to equal the initial value. In fact, in the case of pair_reax.cpp,
some unused code to do this is already in place:

  double qtot;
  MPI_Allreduce(&qsum,&qtot,1,MPI_DOUBLE,MPI_SUM,world);
  elcvec[nlocal+nghost] = 0.0;
  ch[nlocal+nghost] = chpot;

It may be sufficient to change this to:

  double qtot;
MPI_Allreduce(&qsum,&qtot,1,MPI_DOUBLE,MPI_SUM,world);
elcvec[nlocal+nghost] = qtot;
ch[nlocal+nghost] = chpot;

I found that the latter version sometimes led to numerical problems for
charge-neutral systems.

Aidan