4.1.2. Split charge equilibrationΒΆ
The split charge equilibration (SQE) method postulates that partial charges \(q_i\) on atoms can only arise from transfer of charge \(p_{ij}\) between bonded atoms \(i\) and \(j\) only, according to:
\(q_i = \displaystyle{\frac{q_{tot}}{N}} + \sum\limits_{i,j}^{bonds} p_{ij}\)
where \(q_{tot}\) is the net charge on the compound and \(p_{ij} = -p_{ji}\). Although it is trivial to determine the partial charges \(q_i\) from the charge transfer \(p_{ij}\), the reverse is not necessarily true. To exemplify this we write the equations for methanol:
\(\begin{bmatrix}0 & p_{12} & p_{13} & p_{14} & p_{15} & 0 \\ -p_{12} & 0 & 0 & 0 & 0 & 0\\ -p_{13} & 0 & 0 & 0 & 0 & 0\\ -p_{14} & 0 & 0 & 0 & 0 & 0\\ -p_{15} & 0 & 0 & 0 & 0 & p_{56}\\ 0 & 0 & 0 & 0 & 0 -p_{56} & 0\\ \end{bmatrix} \begin{bmatrix}1\\ 1\\ 1 \\ 1 \\ 1 \\ 1 \end{bmatrix} = \begin{bmatrix}q_1\\ q_2\\ q_3 \\ q_4 \\ q_5 \\ q_6 \end{bmatrix}\)
With these equations we can express the \(p_{ij}\) in terms of the charges after manual Gaussian elimination as follows:
\(\begin{bmatrix}p_{12}\\ p_{13}\\ p_{14} \\ p_{15} \\ p_{56} \end{bmatrix} = \begin{bmatrix}-q_2\\ -q_3\\ -q_4\\ -q_5-q_6\\ -q_6\end{bmatrix}\)
with \(q_1\) by definition \(q_{tot}-q_2-q_3-q_4-q_5-q_6\) (with \(q_{tot}\) = 0). However, expressing the \(p_{ij}\) in \(q_i\) does not guarantee that the charge transfer through space is zero. Therefore it is necessary to do the reverse. As outlined by Chen et al. the problem can be solved by expressing the energy in terms of the charge transfer variables.
Verstraelen and co-workers used the following variant of the molecular energy:
\(E_{SQE} = E_{EEM} + \sum\limits_{i,j}^{M} \left(\frac{1}{2}\zeta_{ij}p_{ij}^2+\Delta\chi_{ij}(q_i-q_j)\right)\)
In this formulation, the energy is not just a function of the charges, but also of the charge transfer over \(M\) bonds. \(\zeta_{ij}\) is the bond hardness and \(\Delta\chi_{ij}\) is the electronegativity correction, both of these are bond-specific. We therefore replace the charge variable by charge transfer variables and substitute \(J_{ii} = \eta_i\), and for clarity we replace the summation over bonds by a summation over atom pairs, which is the same thing noting that only \(p_{ij}\) for chemical bonds are non-zero (note that we have to introduce a factor 1/2 to avoid double counting). If we introduce \(M_x\) as the number of bonds for species \(x\) we obtain
\(E_{SQE} = \sum\limits_{n=1}^{N} \left[\left(\frac{q_{tot}}{N} + \sum\limits_{m=1}^{M_n} p_{nm}\right)\left(\chi_n + \frac{1}{2} \eta_n \left(\frac{q_{tot}}{N} + \sum\limits_{m=1}^{M_n} p_{nm}\right) + \frac{1}{2}\sum\limits_{l\ne n}^{N}\left[\frac{q_{tot}}{N} + \sum\limits_{m=1}^{M_l} p_{lm}\right]J_{nl}\right)\right]+\sum\limits_{i,j}^{M}\left[\frac{1}{2}\zeta_{ij}p_{ij}^2+\Delta\chi_{ij}\left(\left[\frac{q_{tot}}{N} + \sum\limits_{m=1}^{M_i} p_{im}\right]-\left[\frac{q_{tot}}{N} + \sum\limits_{m=1}^{M_j} p_{jm}\right]\right)\right] \)
Our task now is to find the \(p_{ij}\) that minimize \(E_{SQE}\). Please note that the all summations are over atoms \(i,j,k,l\). Therefore, we take the derivative and equate it to zero:
\(0 = \displaystyle{\frac{\partial E_{SQE}}{\partial p_{ij}}} = \left(\chi_i - \chi_j+\frac{q_{tot}}{N}\left(\eta_i-\eta_j\right)+\sum\limits_{k=1}^{M_i}\Delta\chi_{ik}-\sum\limits_{k=1}^{M_j}\Delta\chi_{jk}\right) + \frac{1}{2}\left(\sum\limits_{l=1}^{N}J_{il}\left(\frac{q_{tot}}{N}+\sum\limits_{m=1}^{M_l} p_{lm}\right) - \sum\limits_{i=1}^{N}J_{li}\left(\frac{q_{tot}}{N}+\sum\limits_{k=1}^{M_i} p_{ik}\right)\right) + p_{ij}\zeta_{ij}\)
using the identity of the Coulomb-matrix elements (\(J_{ij} = J_{ji}\)) and where the terms involving the atomic hardness \(\eta_x\) are included on the diagonals of the \(J_{xy}\) matrix, except for the contribution of the total charge \(q_{tot}\). Note that the \(q_{tot}\) terms in the two sums cancel. The first term is the difference in electronegativity between atoms \(i\) and \(j\) sharing the bond plus the correction, the second term represents the difference in electrostatic potentials at the atoms, the third term is the interaction between \(i\) and \(j\) times the charge transfer. This leads to a coupled set of equations, written in matrix form as:
\(\mathbf{M}P = R\)
where \(\mathbf{M}\) is a square matrix of dimension number of bonds, \(P\) is a vector containing the charge transfers for all bonds and \(R\) the right hand side of the equations. If the summations are all written out we obtain the result published by Chen et al. for QTPIE and similar algorithms:
\(M_{ij,kl} = J_{ik}-J_{il}-J_{jk}+J_{jl}\)
for off-diagonal elements of the matrix. For the diagonal elements, the bond hardness \(\zeta_{ij}\) needs to be added. The right hand side is given by the electronegativity terms according to:
\(R_{ij} = \chi_j - \chi_i + \sum\limits_{k=1}^{M_j}\Delta\chi_{jk} - \sum\limits_{k=1}^{M_i}\Delta\chi_{ik}+\frac{q_{tot}}{N}\left(\sum\limits_{l=1}^N J_{jl}-\sum\limits_{l=1}^N J_{il}\right).\)
The equations are implemented in Python below.
import numpy as np
from scipy.linalg import lstsq
import math
from charge_utils import *
from test_systems import *
def calcRHS(names, bonds, JEEM, qtotal, model):
Natoms = len(names)
Nbond = len(bonds)
RHS = np.zeros(Nbond, dtype=float)
##### calc corrected chi ####
chi_corr = []
chi = get_chi(model)
for i in range(Natoms):
chi_corr.append(chi[names[i]])
for k in range(Nbond):
ai = names[bonds[k][0]]
aj = names[bonds[k][1]]
if bonds[k][0] == i:
chi_corr[i] += getDeltaChi(ai, aj, model)
if bonds[k][1] == i:
chi_corr[i] -= getDeltaChi(ai, aj, model)
for l in range(Natoms):
chi_corr[i] += JEEM[i][l]*(qtotal/Natoms)
#### calc EN diff for the bonds ####
for i in range(Nbond):
RHS[i] = chi_corr[bonds[i][1]] - chi_corr[bonds[i][0]]
return RHS
def solveSQE(names, coords, bonds, verbose, qtotal, model):
JEEM, JSQE = calcJSQE(names, coords, bonds, verbose, model)
rhs = calcRHS(names, bonds, JEEM, qtotal, model)
if verbose:
print("RHS: {}".format(rhs))
pij = np.linalg.solve(JSQE, rhs)
if verbose:
print("pij: {}".format(pij))
q = np.zeros(coords.shape[0], dtype=float)
for n in range(len(bonds)):
bi = bonds[n][0]
bj = bonds[n][1]
q[bi] += pij[n]
q[bj] -= pij[n]
Natoms = names.shape[0]
for i in range(Natoms):
q[i] += qtotal/Natoms
print("qSQE: {}".format(q))
return q
def run_compound(molname, verbose, shells):
mol = get_system(molname)
print("%s" % molname)
for bonds in mol["bonds"]:
q = solveSQE(mol["names"], mol["coords"], bonds, verbose, mol["qtotal"], "ACS-g")
printDipole(q, mol["coords"], mol["atomnr"], False)
verbose = False
shells = False
for compound in [ "methanol", "water", "acetate", "carbon-monoxide" ]:
run_compound(compound, verbose, shells)
methanol
qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]
qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]
qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]
qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]
qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]
Dipole = 2.07033 Debye
water
qSQE: [ 0.11909703 -0.23819405 0.11909703]
qSQE: [ 0.11909703 -0.23819405 0.11909703]
qSQE: [ 0.11909703 -0.23819405 0.11909703]
Dipole = 0.669411 Debye
acetate
qSQE: [-0.48654202 -0.0231567 -0.24282282 -0.24282009 -0.00136991 -0.00196462
-0.00132384]
qSQE: [-0.48654202 -0.0231567 -0.24282282 -0.24282009 -0.00136991 -0.00196462
-0.00132384]
Dipole = 1.41698 Debye
carbon-monoxide
No hardness information for bond c2-o
No deltaChi information for bond c2-o
No deltaChi information for bond c2-o
qSQE: [ 0.77205802 -0.77205802]
No hardness information for bond o-c2
No deltaChi information for bond o-c2
No deltaChi information for bond o-c2
qSQE: [ 0.77205802 -0.77205802]
Dipole = 0.445003 Debye
