update notebook

Author: PaulWAyers

book/problems/EEM.ipynb

@@ -259,13 +259,106 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 18,
    "id": "ca743a8b",
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Atomic charges: [-0.28375011 -0.28374982 -0.01416517  0.18020443  0.1505785   0.15057809\n",
+      "  0.06419838  0.0641997  -0.00754719 -0.01027312 -0.0102737 ]\n",
+      "Equilibrated chemical potential: 0.24684627271641876\n",
+      "Sum of charges: -5.204170427930421e-17\n"
+     ]
+    }
+   ],
    "source": [
-    "# Set up the equations to solve for the charges and print out the charges.  "
+    "# Set up the equations to solve for the charges and print out the charges.\n",
+    "\n",
+    "def eem_charges(atnums, atcoords, mu, eta, mol_charge, alpha=0.5):\n",
+    "    \"\"\"Compute EEM atomic charges.\n",
+    "\n",
+    "    Builds and solves the (P+1) x (P+1) EEM linear system A·δ = m, where\n",
+    "    the diagonal contains the atomic hardnesses, the off-diagonal atom-atom\n",
+    "    blocks contain the erfgau-screened Coulomb interaction J_AB, the border\n",
+    "    row and column enforce the charge constraint, and the right-hand side\n",
+    "    contains the negative chemical potentials and the total molecular charge.\n",
+    "\n",
+    "    Parameters\n",
+    "    ----------\n",
+    "    atnums : array-like of int, shape (P,)\n",
+    "        Atomic numbers (used only to determine P; properties come from mu/eta).\n",
+    "    atcoords : numpy.ndarray, shape (P, 3)\n",
+    "        Atomic coordinates in bohr.\n",
+    "    mu : numpy.ndarray, shape (P,)\n",
+    "        Atomic chemical potentials (negative electronegativities).\n",
+    "    eta : numpy.ndarray, shape (P,)\n",
+    "        Atomic chemical hardnesses.\n",
+    "    mol_charge : float\n",
+    "        Total molecular charge Q.\n",
+    "    alpha : float\n",
+    "        Range-separation parameter for the erfgau Coulomb screening.\n",
+    "\n",
+    "    Returns\n",
+    "    -------\n",
+    "    charges : numpy.ndarray, shape (P,)\n",
+    "        EEM atomic charges q_A = -Delta N_A.\n",
+    "    mu_total : float\n",
+    "        Equilibrated chemical potential (Lagrange multiplier).\n",
+    "    \"\"\"\n",
+    "    atcoords = np.asarray(atcoords)\n",
+    "    mu = np.asarray(mu)\n",
+    "    eta = np.asarray(eta)\n",
+    "    P = len(mu)\n",
+    "\n",
+    "    # Build the (P+1) x (P+1) EEM matrix A\n",
+    "    A = np.zeros((P + 1, P + 1))\n",
+    "\n",
+    "    # Diagonal: atomic hardnesses\n",
+    "    A[:P, :P] = np.diag(eta)\n",
+    "\n",
+    "    # Off-diagonal atom-atom blocks: erfgau-screened Coulomb interaction\n",
+    "    for i in range(P):\n",
+    "        for j in range(i + 1, P):\n",
+    "            R_ij = np.linalg.norm(atcoords[i] - atcoords[j])\n",
+    "            J_ij = erfgau(R_ij, alpha)\n",
+    "            A[i, j] = J_ij\n",
+    "            A[j, i] = J_ij\n",
+    "\n",
+    "    # Border row and column: charge constraint\n",
+    "    A[:P, P] = 1.0\n",
+    "    A[P, :P] = 1.0\n",
+    "    # A[P, P] = 0 (already zero)\n",
+    "\n",
+    "    # Build the right-hand side m\n",
+    "    m = np.empty(P + 1)\n",
+    "    m[:P] = -mu\n",
+    "    m[P] = -mol_charge\n",
+    "\n",
+    "    # Solve A·delta = m  where delta = [Delta_N_1, ..., Delta_N_P, mu_total]\n",
+    "    delta = np.linalg.solve(A, m)\n",
+    "\n",
+    "    charges = -delta[:P]    # q_A = -Delta N_A\n",
+    "    mu_total = delta[P]\n",
+    "\n",
+    "    return charges, mu_total\n",
+    "\n",
+    "\n",
+    "charges, mu_total = eem_charges(mol0.atnums, mol0.atcoords, mu, eta, mol_charge=0)\n",
+    "print(\"Atomic charges:\", charges)\n",
+    "print(\"Equilibrated chemical potential:\", mu_total)\n",
+    "print(\"Sum of charges:\", charges.sum())\n"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3cd69598",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {