{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Split charge equilibration\n", "The split charge equilibration (SQE) method postulates that partial charges $q_i$ on atoms \n", "can only arise from transfer of charge $p_{ij}$ between bonded atoms $i$ and $j$ only, according to:\n", "\n", "$q_i = \\displaystyle{\\frac{q_{tot}}{N}} + \\sum\\limits_{i,j}^{bonds} p_{ij}$\n", "\n", "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:\n", "\n", "$\\begin{bmatrix}0 & p_{12} & p_{13} & p_{14} & p_{15} & 0 \\\\ \n", "-p_{12} & 0 & 0 & 0 & 0 & 0\\\\\n", "-p_{13} & 0 & 0 & 0 & 0 & 0\\\\\n", "-p_{14} & 0 & 0 & 0 & 0 & 0\\\\\n", "-p_{15} & 0 & 0 & 0 & 0 & p_{56}\\\\\n", "0 & 0 & 0 & 0 & 0 -p_{56} & 0\\\\\n", "\\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}$\n", "\n", "With these equations we can express the $p_{ij}$ in terms of the charges after manual [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) as follows:\n", "\n", "$\\begin{bmatrix}p_{12}\\\\ p_{13}\\\\ p_{14} \\\\ p_{15} \\\\ p_{56} \\end{bmatrix} =\n", "\\begin{bmatrix}-q_2\\\\ -q_3\\\\ -q_4\\\\ -q_5-q_6\\\\ -q_6\\end{bmatrix}$\n", "\n", "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.](https://doi.org/10.1063/1.3021400) the problem can be solved by expressing the energy in terms of the charge transfer variables.\n", "\n", "[Verstraelen and co-workers](http://dx.doi.org/10.1063/1.3187034) used the following variant of the molecular energy:\n", "\n", "$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)$\n", "\n", "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\n", "\n", "$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] $\n", "\n", "\n", "\n", "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:\n", "\n", "$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}$\n", "\n", "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.\n", "This leads to a coupled set of equations, written in matrix form as:\n", "\n", "$\\mathbf{M}P = R$\n", "\n", "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 \n", "result published by [Chen et al.](https://aip.scitation.org/doi/10.1063/1.3021400) for QTPIE and similar algorithms:\n", "\n", "$M_{ij,kl} = J_{ik}-J_{il}-J_{jk}+J_{jl}$\n", "\n", "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:\n", "\n", "$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).$\n", "\n", "The equations are implemented in Python below." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "methanol\n", "qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]\n", "qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]\n", "qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]\n", "qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]\n", "qSQE: [-0.50802387 0.33691938 -0.33742533 0.16436908 0.16438559 0.17977516]\n", "Dipole = 2.07033 Debye\n", "water\n", "qSQE: [ 0.11909703 -0.23819405 0.11909703]\n", "qSQE: [ 0.11909703 -0.23819405 0.11909703]\n", "qSQE: [ 0.11909703 -0.23819405 0.11909703]\n", "Dipole = 0.669411 Debye\n", "acetate\n", "qSQE: [-0.48654202 -0.0231567 -0.24282282 -0.24282009 -0.00136991 -0.00196462\n", " -0.00132384]\n", "qSQE: [-0.48654202 -0.0231567 -0.24282282 -0.24282009 -0.00136991 -0.00196462\n", " -0.00132384]\n", "Dipole = 1.41698 Debye\n", "carbon-monoxide\n", "No hardness information for bond c2-o\n", "No deltaChi information for bond c2-o\n", "No deltaChi information for bond c2-o\n", "qSQE: [ 0.77205802 -0.77205802]\n", "No hardness information for bond o-c2\n", "No deltaChi information for bond o-c2\n", "No deltaChi information for bond o-c2\n", "qSQE: [ 0.77205802 -0.77205802]\n", "Dipole = 0.445003 Debye\n" ] } ], "source": [ "import numpy as np\n", "from scipy.linalg import lstsq\n", "import math\n", "from charge_utils import *\n", "from test_systems import *\n", "\n", "def calcRHS(names, bonds, JEEM, qtotal, model):\n", " Natoms = len(names)\n", " Nbond = len(bonds)\n", " RHS = np.zeros(Nbond, dtype=float)\n", " ##### calc corrected chi ####\n", " chi_corr = []\n", " chi = get_chi(model)\n", " for i in range(Natoms):\n", " chi_corr.append(chi[names[i]])\n", " for k in range(Nbond):\n", " ai = names[bonds[k][0]]\n", " aj = names[bonds[k][1]]\n", " if bonds[k][0] == i:\n", " chi_corr[i] += getDeltaChi(ai, aj, model)\n", " if bonds[k][1] == i:\n", " chi_corr[i] -= getDeltaChi(ai, aj, model)\n", " for l in range(Natoms):\n", " chi_corr[i] += JEEM[i][l]*(qtotal/Natoms)\n", " #### calc EN diff for the bonds ####\n", " for i in range(Nbond):\n", " RHS[i] = chi_corr[bonds[i][1]] - chi_corr[bonds[i][0]]\n", " return RHS\n", "\n", "def solveSQE(names, coords, bonds, verbose, qtotal, model):\n", " JEEM, JSQE = calcJSQE(names, coords, bonds, verbose, model)\n", " rhs = calcRHS(names, bonds, JEEM, qtotal, model)\n", " if verbose:\n", " print(\"RHS: {}\".format(rhs))\n", " pij = np.linalg.solve(JSQE, rhs)\n", " if verbose:\n", " print(\"pij: {}\".format(pij))\n", " q = np.zeros(coords.shape[0], dtype=float)\n", " for n in range(len(bonds)):\n", " bi = bonds[n][0]\n", " bj = bonds[n][1]\n", " q[bi] += pij[n]\n", " q[bj] -= pij[n]\n", " Natoms = names.shape[0]\n", " for i in range(Natoms):\n", " q[i] += qtotal/Natoms\n", " print(\"qSQE: {}\".format(q))\n", " return q\n", "\n", "def run_compound(molname, verbose, shells):\n", " mol = get_system(molname)\n", " print(\"%s\" % molname)\n", " for bonds in mol[\"bonds\"]:\n", " q = solveSQE(mol[\"names\"], mol[\"coords\"], bonds, verbose, mol[\"qtotal\"], \"ACS-g\")\n", " printDipole(q, mol[\"coords\"], mol[\"atomnr\"], False)\n", "\n", "verbose = False\n", "shells = False\n", "\n", "for compound in [ \"methanol\", \"water\", \"acetate\", \"carbon-monoxide\" ]:\n", " run_compound(compound, verbose, shells)\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.1" } }, "nbformat": 4, "nbformat_minor": 4 }