ReaxFF/C charge equilibration for charged systems

Hi,

I have found several messages in the mailing list archive saying that
the charge equilibration for ReaxFF/C currently only supports neutral
systems. As I am interested in running charged systems and it was also
mentioned that generalizing this would not be too big of a change, I
would like to ask if it would be possible to do this? I would be happy
to help with coding (though I would need some guidance) and definitely
with testing.

Thank you for considering this,
Ondrej Marsalek

Hi Ondrej,

It is not recommended to have a non-zero net charge with the QEq method. From the fix qeq/variant doc page: “IMPORTANT NOTE: To avoid the evaluation of the derivative of charge with respect to position, which is typically ill-defined, the system should have a zero net charge.”

If you would like to try this regardless, you can look at pair_reax.cpp in the REAX package, particularly the “qtot” variable. That is how the fortran version of ReaxFF includes the non-zero net charge. You just have to change the elcvec value. You can do the same with the C version of ReaxFF.

— pair_reax.cpp 2014-11-10 15:11:19.000000000 -0700

+++ REAX/pair_reax.cpp 2014-11-10 15:11:28.000000000 -0700

@@ -887,7 +887,7 @@

double qtot;

MPI_Allreduce(&qsum,&qtot,1,MPI_DOUBLE,MPI_SUM,world);

  • elcvec[nlocal+nghost] = qtot;
  • elcvec[nlocal+nghost] = 0.0;

ch[nlocal+nghost] = chpot;

Ray

Hi Ray,

thank you for the recommendation. I went by what was said in these emails:

http://lammps.sandia.gov/threads/msg23146.html
http://lammps.sandia.gov/threads/msg23164.html

so I thought it should be possible. The original ReaxFF implementation
apparently does it and also any electronic structure method will
"equilibrate charge", so in principle this should work. I am not sure
I understand what exactly the note refers to, but perhaps it is a
technical/implementation limitation, rather than a fundamental one?

Ondrej

As an additional thought on the possible numerical instability
mentioned in the 2011 email, would it perhaps be possible and helpful
to specify the target value of charge as an argument to qeq/reax,
rather than calculating it by summing the charges?

Ondrej

Ondrej,

Please refer to the attachement. If the system is non-neutral, then evaluating the derivative of charge wrt to position is inevitable; and this derivative is ill-defined. Some variations of ReaxFF code does it but it does not mean it has considered this and the electronic strucutre method is another story. It is a fundamental and a mathematical issue, and this is how “numerical instability” is originated. Unless “partial q over partial r” is carefully and rigorously evaluated, I don’t expect there to be an easy workaround.

Ray

The attachment.

Ray

min_w_qeq.pdf (48.6 KB)

Ray,

am I not correct that at charge equilibrium partial derivatives of E with respect to q_i should be 0 (by the concept of minimum)? They are equal to \mu_i at zero i-th charge (eq. 2 in Rappe and Goddard paper), not at equilibrium.

Oleg

6:17, 11 ноября 2014 г., Ray Shan <rayshan819@…24…>:

No, I’m not.

Oleg

Oleg Sergeev <seoman@…567…> 13 ноября 2014 г. 16:23:34 написал: