Add helium atom

Author: bjorgve

docs/notebooks/helium_atom.ipynb

@@ -0,0 +1,139 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Helium atom\n",
+    "\n",
+    "Once we have a solver for the Hydrogen Atom (1 atom), the next level off difficulty is the Helium atom (2 atoms). Here we'll show how\n",
+    "to make and SCF for the Helium atom within the framework of Hartree–Fock."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from vampyr import vampyr3d as vp\n",
+    "import numpy as np\n",
+    "\n",
+    "def laplace_operator(D, f_tree):\n",
+    "    return D(D(f_tree, 0), 0) + D(D(f_tree, 1), 1) + D(D(f_tree, 2), 2)\n",
+    "\n",
+    "def calculate_energy(phi_tree, V_tree):\n",
+    "    mra = phi_tree.MRA()\n",
+    "    d_oper = vp.ABGVDerivative(mra, 0.5, 0.5)\n",
+    "    return -0.5*vp.dot(laplace_operator(d_oper, phi_tree), phi_tree) + vp.dot(phi_tree, V_tree*phi_tree)\n",
+    "\n",
+    "def couloumb_potential(prec, phi_tree):\n",
+    "    mra = phi_tree.MRA()\n",
+    "    P = vp.PoissonOperator(mra, prec)\n",
+    "    return P(4.0*np.pi*phi_tree*phi_tree)\n",
+    "\n",
+    "# Analytic nuclear potential: f_nuc(r) = Z/|r|\n",
+    "def f_nuc(r):\n",
+    "    Z = 2.0                 # Nuclear charge\n",
+    "    R = np.sqrt(r[0]**2 + r[1]**2 + r[2]**2)\n",
+    "    return -Z / R\n",
+    "\n",
+    "# Analytic guess for wavefunction: f_phi(r) = exp(-r^2)\n",
+    "def f_phi(r):\n",
+    "    R2 = r[0]**2 + r[1]**2 + r[2]**2\n",
+    "    return np.exp(-R2)\n",
+    "\n",
+    "# Set target precision\n",
+    "precision = 1.0e-5\n",
+    "\n",
+    "# Define MRA basis and computational domain\n",
+    "k = 5                       # Polynomial order of MRA basis\n",
+    "world = [-20, 20]           # Computational domain [-L,L]^3 (a.u.)\n",
+    "MRA = vp.MultiResolutionAnalysis(order=k, box=world)\n",
+    "\n",
+    "# Define projector onto the MRA basis\n",
+    "P_mra = vp.ScalingProjector(mra=MRA, prec=precision)\n",
+    "\n",
+    "# Initialize the calculation\n",
+    "phi_n = P_mra(f_phi)           # Project analytic guess for wavefunction onto MRA\n",
+    "phi_n.normalize()              # Normalize the wavefunction guess\n",
+    "V_n = P_mra(f_nuc)               # Project analytic nuclear potential onto MRA\n",
+    "J = couloumb_potential(precision, phi_n)\n",
+    "E_n = calculate_energy(phi_n, V_n + J)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Loop parameters\n",
+    "iteration = 0                  # Iteration counter\n",
+    "max_iter = 30                  # Maximum iterations \n",
+    "thrs = precision # -1           # Convergence requirement. Set to -1 if you wish to limit using max_iter\n",
+    "update = 1.0                   # Initialize error measure (norm of wavefunction update)\n",
+    "# Minimization loop\n",
+    "while (update > thrs):\n",
+    "    if iteration > max_iter-1:\n",
+    "        break\n",
+    "    # Build Helmholtz operator from current energy\n",
+    "    mu = np.sqrt(-2*E_n)\n",
+    "    G = vp.HelmholtzOperator(mra=MRA, exp=mu, prec=precision)\n",
+    "    \n",
+    "    # Apply Helmholtz operator\n",
+    "    phi_np1 = -2*G((V_n + J)*phi_n)\n",
+    "    \n",
+    "    # Compute wavefunction and energy update\n",
+    "    d_phi_n = phi_np1 - phi_n\n",
+    "    update = d_phi_n.norm()\n",
+    "\n",
+    "    # Prepare for next iteration\n",
+    "    phi_n = phi_np1\n",
+    "    phi_n.normalize()\n",
+    "    J = couloumb_potential(precision, phi_n)\n",
+    "    E_n = calculate_energy(phi_n, V_n + J)\n",
+    "\n",
+    "    # Collect output\n",
+    "    print(iteration, \" |  E:\", E_n, \" |  d_phi:\", update)\n",
+    "    iteration += 1\n",
+    "E_tot = 2.0*E_n - vp.dot(phi_n*phi_n, couloumb_potential(precision, phi_n))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "documentation",
+   "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.10.8"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "a875265825dca88dc0cdca05b35ae15709f7b78b2bf941fac513dcdaf1078523"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}