electrostatics in the KIM API: plan and open questions

OpenKIM
21

Hi,

Here are my comments on the "QEq" part:

* First things first. It bothers me (and many others that work on density functional theory) that the term QEq is so often used for models with variable charges. The pioneering work was in fact done by Mortier et al (in 1985-1986), long before the QEq paper came out (1991). QEq is mostly a reinvention of Mortier's work, with some minor improvements and without all the nice insights from Density Functional Theory. For that reason, I prefer to use Mortier's terminology: the Electronegativity Equalization Method (EEM) instead of Charge Equilibration (QEq). Any model with variable charges is essentially a variant of EEM. It would be nice to reflect this in the KIM API, especially because this community cares so much about provenance. (I'll do this consistently below.)

* It is not so useful to combine plain EEM with long-range electrostatics. One major deficiency/feature of plain EEM is that the entire system exhibits a metallic polarizability scaling. This means that the polarizability scales cubically with the system size. As a consequence, the charges always evolve such that all long-range electrostatics are completely screened away. The current implementation of ReaxFF is therefore fine without long-range electrostatics. This special behavior disappears when restraints or constraints are imposed on the variable charges (except for the usual constraint on the total charge.) As I explain it here, it may seem that metallic polarizability scaling is a good thing but for most simulations, e.g. on dielectric systems, it results in completely erratic electrostatic interactions. Hence, there are several attempts to fix this problem. Alas, if a model fixes it, costly long-range electrostatics become relevant again. For example, at the SES meeting, it was mentioned that COMB uses some restraints in the charge-solver (not in the model). New EEM variants like Nistor's SQE and our ACKS2 fix this polarizability scaling issue in the model.

* The API for EEM-ish models can be simplified. The Model only needs to provide the charge-gradient of the energy as function of the carges, for the part of the energy that is computed in the model. In the case of short-range electrostatics, that gradient can entirely be computed by the model. In case of long-range electrostatics, the long-range part of the gradient must be computed by the Simulator. I think any existing atomistic simulator converges the charges with some sort of self-consistent, conjugate-gradient or minres algorithm. All these algorithms can be implemented in the Simulator, given an API to get (part of) the gradient from the model. The details of the sort-range contributions to that gradient (electronegativity, hardness, off-diagonal terms, ... ) are irrelevant for the API. They completely live in the model space. The constraint on the total charge is better taken care of in the Simulator because (i) it requires some communication in parallel codes and (ii) it is usually not included in the charge-gradient.

* For split-charge equilibration, the situation is more complicated. In that model, the variables are not the charges, but split charges (i.e. charge-transfers) between some (not all!) pairs of atoms. To find the optimal split charges, the model must be able to compute the derivative of the energy toward the split charges. The simulator can then iteratively update the split charges to zero this gradient. This is all very similar to the EEM API, except for two points:

1) The model and the simulator need to agree on a list of atom pairs between which the split-charges are present. This list may not be constant throughout an MD simulation, e.g. when bonds break. This complexity is comparable to the topology in a typical covalent force field.

2) It is computationally easy to transform the charge-gradient to a split-charge gradient, but not the other way around. Hence, in the case of long-range electrostatics, the Simulator needs to carry out a transformation in that way to obtain the long-range contribution to the split-charge gradient.

* One of the reasons we derived the ACKS2 model, is to avoid all the bookkeeping of the split charges in SQE. In ACKS2, every atom has two variables (in the minimal case): an atomic charge and an effective potential. (The effective potential is different from the electrostatic potential. It acts as a genuine variable.) The set of possible ACKS2 models is mathematically speaking a superset of the possible SQE models, so ACKS2 can be used as an "alternative implementation" of SQE with only atomic fluctuating degrees of freedom. This may be more bearable for the KIM API. There are other advantages of ACKS2 as well, but I won't bother you with these. You can read all about it in our recent JCP paper (http://dx.doi.org/10.1063/1.4901513, just published, so free download) and in our first paper on the ACKS2 model (http://dx.doi.org/10.1063/1.4791569, harder to digest than the newest paper, contains proof of equivalence with SQE).

Best Regards,

Toon

22

Toon raises some very good points that need to be considered in the API. I
have two specific comments on those below. I addition, I will add the
general comment that variable charge methods are much less standardized
than fixed charged methods. Hence, it is much more difficult to define an
API that includes most of the methods out there. I suggest starting with a
simple prototype API that handles just one or two methods. The lessons
learned from getting that to work can be used to decide how other methods
should be handled.

Aidan¹

23

Hi,

I'd like to put my above comment a little into perspective. I started wondering how much precision loss one could get in a typical test case. In the following simple test, it turns out not to be a big issue, as opposed to what I believed before.

I took a snapshot from a 32-water molecule simulation with the flexible TIP3P CHARMM force field. This is a simple case where I was expecting significant precision loss. The energy and RMSD of the forces are as follows:

                            Energy RMSD Force

Sum off all terms: -0.5012681286 0.013524530470
     Bond stretch: 0.0452453653 0.015609513945
     Angle bending: 0.0297902428 0.006708814297
     Lennard-Jones: 0.0801189458 0.004467270550
     Electrostatics: -0.6564226826 0.013226552963

(All values in atomic units.) When I compute the electrostatic term without the exclusion lists, i.e. adding all interactions between point charges, including those within the water molecules, I get the following:

Full Electrostatics: -10.7677968129 0.048682429176

The full electrostatics is about two orders of magnitude larger than the one where intra-molecular terms are excluded. This implies than Ryan's suggestion (to compute the full point-charge electrostatics in the Simulator and to compensate short-range effects in the Model) leads to a precision loss of only a few digits in the energy. This is much less than I expected and I guess this is acceptable for most applications? For the forces, there is even less precision loss.

Best Regards,

Toon

24

* First things first. It bothers me (and many others that work on
density functional theory) that the term QEq is so often used for
models with variable charges. The pioneering work was in fact done
by Mortier et al (in 1985-1986), long before the QEq paper came out
(1991). QEq is mostly a reinvention of Mortier's work, with some
minor improvements and without all the nice insights from Density
Functional Theory. For that reason, I prefer to use Mortier's
terminology: the Electronegativity Equalization Method (EEM)
instead of Charge Equilibration (QEq). Any model with variable
charges is essentially a variant of EEM. It would be nice to
reflect this in the KIM API, especially because this community
cares so much about provenance. (I'll do this consistently below.)

I have a problem with specific terms like EEM and QEq. Some people
use them to refer to a general class of similar methods, while others
use them to make distinctions between these similar methods. This
makes it hard to know what meaning is intended. I suggest that KIM
avoid this difficulty by picking a generic descriptive term like
Variable Charge, as Toon did above.

Hi,

This is a good idea. "Variable charges" is the most clear and self-explaining name and the subtle differences between EEM and QEq should not matter for the API.

* The API for EEM-ish models can be simplified. The Model only
needs to provide the charge-gradient of the energy as function of
the carges, for the part of the energy that is computed in the
model. In the case of short-range electrostatics, that gradient can
entirely be computed by the model.

This is true, but it might not be the most convenient choice for
some Model/Simulator combinations.

Could you explain why? I'd like to understand this.

In case of long-range electrostatics, the long-range part of the
gradient must be computed by the Simulator. I think any existing
atomistic simulator converges the charges with some sort of
self-consistent, conjugate-gradient or minres algorithm.

Note that some potentials, such as COMB, use damped dynamics to
relax charges, but I think it is still true that only charge-gradient
is required.

I agree. When fluctuating charges are propagated dynamically with an extended Lagrangian (cfr. Car-Parrinello MD), the charge-gradient is sufficient. I assume the damped dynamics in COMB is almost the same as the extended Lagrangian formalism.

Best Regards,

Toon

25

questions

* The API for EEM-ish models can be simplified. The Model only
needs to provide the charge-gradient of the energy as function of
the carges, for the part of the energy that is computed in the
model. In the case of short-range electrostatics, that gradient can
entirely be computed by the model.

This is true, but it might not be the most convenient choice for
some Model/Simulator combinations.

Could you explain why? I'd like to understand this.

No, I can¹t really explain it, that¹s why I said ³might². I just know that
the variable charge method in ReaxFF explicitly requires the elements of a
sparse matrix A, where each non-zero off-diagonal element is roughly
1/r_ij. Those elements can not be obtained from the charge gradient. I
think the charge gradient would be sufficient to calculate the
matrix-vector product A.q, but that is not how some simulators are coded.

Aidan

26

Hi,

I can provide some more background here. (I've coded an ACKS2-ish model in the ReaxFF code of van Duin.)

The matrix-vector product you mention is the charge-gradient itself. The sparse matrix is constructed to allow for a rapid evaluation of the charge-gradient for different charge vectors. In ReaxFF, the energy is a quadratic function of the charge vector, so the charge-gradient is a linear function of the charge vector.

This sparse matrix contains all the short-range electrostatics (i.e. 1/r) on the off-diagonals and the atomic hardness parameters on the diagonal elements. It is therefore very similar to a neighborlist. For the KIM API, it is fine to construct this matrix inside a model. For the sake of efficiency, it should be possible for the simulator to let the model know that the charges changed without a change in the atomic positions, such that the sparse matrix doesn't have to be reconstructed inside the loop of the charge-solver.

Best Regards,

Toon

27

Hi All,

I just wanted to give a short note of thanks for the interesting discussion that is continuing here. I'm following the thread, but it will be some time before I'm able to get back to this and carefully consider all the issues. It is great to see all of this discussed and recorded on the mailing list. I think it is really useful for the community, in general.

Thanks again.

Ryan

28

Hello,

One more correction to my previous comment came to mind. The self-consistent update used in variable-charge models can be reformulated as a newton method to minimize the energy as function of the charges, where only the diagonal of the charge-Hessian is used and the off-diagonals are approximated to be zero. (The gradient is still computed properly.) If one wants to support the self-consistent update in a Simulator, the API must be extended such that the simulator can ask the diagonal of the charge-Hessian from the model. If a variable-charge model is not simply quadratic in the charges, e.g. as in QEq, this diagonal depends on the current value of the charges and it can, in that case, not be communicated once in the beginning of the simulation.

Best Regards,

Toon