{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Electronegativity equalization method\n",
    "The electronegativity equalization method (EEM) is a second order expansion of the molecular energy $E_{EEM}$ in terms of the partial charges $q_i$:\n",
    "\n",
    "$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]$\n",
    "\n",
    "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:\n",
    "\n",
    "$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}$\n",
    "\n",
    "\n",
    "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:\n",
    "\n",
    "$\\begin{bmatrix}\\eta_1 & J_{12} & J_{13} & ... & J_{1N} & -1 \\\\ J_{21} & \\eta_2 & J_{23} & ... & J_{2N} &  -1 \\\\\n",
    "... & ... & ...  &... & ... & ... \\\\\n",
    "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}$\n",
    "\n",
    "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](https://doi.org/10.1021/j100161a070). \n",
    "\n",
    "In the code below we test the method on a methanol molecule, using parameters from [Verstraelen et al.](http://dx.doi.org/10.1063/1.3187034) as well as experimental data from [Cordero et al.](https://doi.org/10.1039/B801115J).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "carbon-monoxide\n",
      "q[c2] = 0.100795\n",
      "q[o] = -0.100795\n",
      "Dipole = 0.0580969 Debye\n",
      "\n",
      "water\n",
      "q[hw] = 0.130157\n",
      "q[ow] = -0.260313\n",
      "q[hw] = 0.130157\n",
      "Dipole = 0.731574 Debye\n",
      "\n",
      "methanol\n",
      "q[oh] = -0.252355\n",
      "q[hp] = 0.13162\n",
      "q[c3] = -0.151608\n",
      "q[h1] = 0.0878756\n",
      "q[h1] = 0.087886\n",
      "q[h1] = 0.0965812\n",
      "Dipole = 1.22594 Debye\n",
      "\n",
      "acetate\n",
      "q[c3] = -0.254055\n",
      "q[cm] = -0.108471\n",
      "q[om] = -0.322013\n",
      "q[om] = -0.322018\n",
      "q[hc] = 0.00241518\n",
      "q[hc] = 0.00166854\n",
      "q[hc] = 0.00247338\n",
      "Dipole = 0.841977 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 solveEEM(names, coords, qtotal, model, verbose=False):\n",
    "    # Compute Coulomb\n",
    "    J = calcJEEM(names, coords, model)\n",
    "    N = coords.shape[0]\n",
    "\n",
    "    # Right hand side of the equation\n",
    "    rhs = np.zeros(N+1,dtype=float)\n",
    "    chi = get_chi(model)\n",
    "    for i in range(N):\n",
    "        if names[i] in chi:\n",
    "            rhs[i] = -chi[names[i]]\n",
    "        else:\n",
    "            print(\"No chi for %s\" % names[i])\n",
    "            exit(1)\n",
    "    rhs[N] = qtotal\n",
    "    q = np.linalg.solve(J, rhs)\n",
    "    if verbose:\n",
    "        print(\"J = \\n{}\".format(J))\n",
    "        print(\"rhs = {}\".format(rhs))\n",
    "        print(\"q = {}\".format(q))\n",
    "        y = np.dot(J,q)\n",
    "        print(\"y = {}\".format(y))\n",
    "    return q\n",
    "\n",
    "def run_compound(molname, verbose, alexandria=False):\n",
    "    mol = get_system(molname)\n",
    "    if not mol:\n",
    "        print(\"No test system %s defined\" % molname)\n",
    "    else:\n",
    "        print(\"\\n%s\" % molname)\n",
    "        q = solveEEM(mol[\"names\"], mol[\"coords\"], mol[\"qtotal\"], verbose, alexandria)\n",
    "        for i in range(mol[\"coords\"].shape[0]):\n",
    "            print(\"q[%s] = %g\" % (mol[\"names\"][i], q[i]))\n",
    "        printDipole(q, mol[\"coords\"], mol[\"atomnr\"], False)\n",
    "\n",
    "verbose = False\n",
    "\n",
    "for compound in [ \"carbon-monoxide\", \"water\", \"methanol\", \"acetate\" ]:\n",
    "    run_compound(compound, \"ACM-g\", verbose)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the positive charge on O and negative charge on one of the H with the EEM chi and eta values from Verstraelen.\n",
    "These charges are too small since  the experimental dipole for methanol is about 1.69 Debye."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}