4.1.1. Electronegativity equalization method

The electronegativity equalization method (EEM) is a second order expansion of the molecular energy \(E_{EEM}\) in terms of the partial charges \(q_i\):

\(E_{EEM}(q_1,q_2,...,q_N) = \sum\limits_{i=1}^{N} \left[ \chi_i q_i + \frac{1}{2}\eta_i q_i^2 + \frac{1}{2}\sum\limits_{l\ne i}^{N}q_iq_lJ_{il}\right]\)

where \(N\) is the number of atoms, \(\chi_i\) are the atomic electronegativities, \(\eta_i\) the atomic hardness and \(J_{il}\) the Coulomb interaction between atoms. The factor 1/2 before the Coulomb matrix is to avoid double counting. In equilibrium, the chemical potential is the same for each atom, \(\chi_{eq}\). The EEM can be implemented by taking the derivative of \(E_{EEM}\) with respect to \(q_i\) and equating it to zero:

\(0 = \displaystyle{\frac{\partial E_{EEM}}{\partial q_i}} = \chi_{eq} = \chi_i + q_i\eta_i + \frac{1}{2}\sum\limits_{l\ne i}^{N} q_l J_{il}\)

which leads to a set of \(N+1\) equations that is linear in the \(q_i\) and that can be solved using linear algebra tools. After taking the derivative with respect to \(q_i\) a matrix equation is obtained:

\(\begin{bmatrix}\eta_1 & J_{12} & J_{13} & ... & J_{1N} & -1 \\ J_{21} & \eta_2 & J_{23} & ... & J_{2N} & -1 \\ ... & ... & ... &... & ... & ... \\ J_{N1} & J_{N2} & J_{N3} & ... & \eta_N & -1 \\ 1 & 1 & 1 & ... & 1 & 0\end{bmatrix} \begin{bmatrix}q_1\\ q_2\\ ... \\ q_N \\ \chi_{eq} \end{bmatrix} = \begin{bmatrix}-\chi_1\\-\chi_2\\...\\ -\chi_N\\ q_{tot}\end{bmatrix}\)

Note the last column in the matrix is there to make sure that the electronegativity for all atoms is the same (\(\chi_{eq}\)), while the last row is there to make sure the total charge \(q_{tot}\) is maintained. A good reference for the method is the classic paper by Rappe and Goddard.

In the code below we test the method on a methanol molecule, using parameters from Verstraelen et al. as well as experimental data from Cordero et al..

import numpy as np
from scipy.linalg import lstsq
import math
from charge_utils import *
from test_systems import *

def solveEEM(names, coords, qtotal, model, verbose=False):
    # Compute Coulomb
    J = calcJEEM(names, coords, model)
    N = coords.shape[0]

    # Right hand side of the equation
    rhs = np.zeros(N+1,dtype=float)
    chi = get_chi(model)
    for i in range(N):
        if names[i] in chi:
            rhs[i] = -chi[names[i]]
        else:
            print("No chi for %s" % names[i])
            exit(1)
    rhs[N] = qtotal
    q = np.linalg.solve(J, rhs)
    if verbose:
        print("J = \n{}".format(J))
        print("rhs = {}".format(rhs))
        print("q = {}".format(q))
        y = np.dot(J,q)
        print("y = {}".format(y))
    return q

def run_compound(molname, verbose, alexandria=False):
    mol = get_system(molname)
    if not mol:
        print("No test system %s defined" % molname)
    else:
        print("\n%s" % molname)
        q = solveEEM(mol["names"], mol["coords"], mol["qtotal"], verbose, alexandria)
        for i in range(mol["coords"].shape[0]):
            print("q[%s] = %g" % (mol["names"][i], q[i]))
        printDipole(q, mol["coords"], mol["atomnr"], False)

verbose = False

for compound in [ "carbon-monoxide", "water", "methanol", "acetate" ]:
    run_compound(compound, "ACM-g", verbose)

carbon-monoxide
q[c2] = 0.100795
q[o] = -0.100795
Dipole = 0.0580969 Debye

water
q[hw] = 0.130157
q[ow] = -0.260313
q[hw] = 0.130157
Dipole = 0.731574 Debye

methanol
q[oh] = -0.252355
q[hp] = 0.13162
q[c3] = -0.151608
q[h1] = 0.0878756
q[h1] = 0.087886
q[h1] = 0.0965812
Dipole = 1.22594 Debye

acetate
q[c3] = -0.254055
q[cm] = -0.108471
q[om] = -0.322013
q[om] = -0.322018
q[hc] = 0.00241518
q[hc] = 0.00166854
q[hc] = 0.00247338
Dipole = 0.841977 Debye

Note the positive charge on O and negative charge on one of the H with the EEM chi and eta values from Verstraelen. These charges are too small since the experimental dipole for methanol is about 1.69 Debye.