timing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2310a17c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cda1dd0",
   "metadata": {},
   "source": [
    "# Timing for the different functions of the code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cb53ce1",
   "metadata": {},
   "source": [
    "## Hermite polynomials\n",
    "\n",
    "Comparison of the recursive and scipy Hermite polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "14b61966",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import special\n",
    "\n",
    "def hermite_pol_recursive(x, n):\n",
    "    if n == 0:\n",
    "        return 1\n",
    "    elif n == 1:\n",
    "        return 2*x\n",
    "    else:\n",
    "        return 2*x*hermite_pol_recursive(x, n-1) - 2*(n - 1)*hermite_pol_recursive(x, n-2)\n",
    "\n",
    "def hermite_pol_scipy_1(x, n):\n",
    "    H = special.hermite(n)\n",
    "    return H(x)\n",
    "\n",
    "def hermite_pol_scipy_2(x, n):\n",
    "    return special.eval_hermite(n, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "51ca5c14",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Recursive Hermite polynomials\n",
      "3.68 µs ± 110 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      "128 µs ± 1.65 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "1.63 ms ± 38 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "Scipy Hermite polynomials (option 1)\n",
      "283 µs ± 5.96 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "361 µs ± 8.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "428 µs ± 7.15 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "Scipy Hermite polynomials (option 2)\n",
      "182 µs ± 4.8 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "214 µs ± 2.83 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "269 µs ± 11.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "N = 10000\n",
    "x = np.linspace(-10, 10, N)\n",
    "\n",
    "print(\"Recursive Hermite polynomials\")\n",
    "%timeit hermite_pol_recursive(x, 1)\n",
    "%timeit hermite_pol_recursive(x, 5)\n",
    "%timeit hermite_pol_recursive(x, 10)\n",
    "\n",
    "print(\"Scipy Hermite polynomials (option 1)\")\n",
    "%timeit hermite_pol_scipy_1(x, 1)\n",
    "%timeit hermite_pol_scipy_1(x, 5)\n",
    "%timeit hermite_pol_scipy_1(x, 10)\n",
    "\n",
    "print(\"Scipy Hermite polynomials (option 2)\")\n",
    "%timeit hermite_pol_scipy_2(x, 1)\n",
    "%timeit hermite_pol_scipy_2(x, 5)\n",
    "%timeit hermite_pol_scipy_2(x, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0fcf7bd",
   "metadata": {},
   "source": [
    "## Timing for the different functions of the code\n",
    "\n",
    "Write necessary functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "29b09bb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import lib\n",
    "import numpy as np\n",
    "from scipy import special\n",
    "import os\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "31702187",
   "metadata": {},
   "outputs": [],
   "source": [
    "lz = 0.5\n",
    "lx = 0.75\n",
    "\n",
    "OMEGA_X = 1/lx**2\n",
    "OMEGA_Z = 1/lz**2\n",
    "\n",
    "def Hermite_pol(x, n):\n",
    "    \"\"\"\n",
    "    Evaluates the Hermite polynomial of order n at x.\n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    x : float or np.ndarray\n",
    "        Values in which to evaluate the Hermite polynomial\n",
    "    n : int\n",
    "        Oder of the Hermite polynomial\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    float or np.ndarray\n",
    "    \"\"\"\n",
    "    return special.eval_hermite(n, x)\n",
    "\n",
    "def HO_wf(x, n, omega=OMEGA_X):\n",
    "    \"\"\"\n",
    "    Evaluates the n^th eigenfunction of the Harmonic Oscillator at x.\n",
    "    The units of x are sqrt(hbar/m*omega). \n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    x : float or np.ndarray(N)\n",
    "        Values in which to evaluate eigenfunction\n",
    "    n : int\n",
    "        Number of the eigenfunction\n",
    "    omega : float\n",
    "        Angular frequency of the Harmonic Oscillator\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    float or np.ndarray(N)\n",
    "    \"\"\"\n",
    "    #return omega**0.25 * Hermite_pol(np.sqrt(omega)*x, n) * np.exp(-omega*x**2 / 2) / np.sqrt(2**n * np.math.factorial(n) * np.sqrt(np.pi))\n",
    "    return omega**0.25 * Hermite_pol(np.sqrt(omega)*x, n) / np.sqrt(2**n * np.math.factorial(n) * np.sqrt(np.pi))\n",
    "\n",
    "def HO_wf_3D(x, y, z, nx, ny, nz, omega_x=OMEGA_X, omega_y=OMEGA_X, omega_z=OMEGA_Z):\n",
    "    \"\"\"\n",
    "    Evaluates the n^th eigenfunction of the (anisotropic) Harmonic Oscillator at x,y,z.\n",
    "    The units of x,y,z are sqrt(hbar/m*omega). \n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    x, y, z : float or np.ndarray(N)\n",
    "        Position in which to evaluate the wave function\n",
    "    nx, ny, nz : int\n",
    "        Number of the eigenfunction for each cartesian coordinate\n",
    "    omega_x, omega_y, omega_z : float\n",
    "        Harmonic osciallation constant for each cartesian coordinate\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    float or np.ndarray(N)\n",
    "    \"\"\"\n",
    "    return HO_wf(x, nx, omega=omega_x)*HO_wf(y, ny, omega=omega_y)*HO_wf(z, nz, omega=omega_z)\n",
    "\n",
    "def index_to_q_numbers(k):\n",
    "    \"\"\"\n",
    "    Returns the quantum numbers nx, ny, nz associated with the basis index k\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    k: int\n",
    "        Index of the basis from 0 to 13\n",
    "\n",
    "    Returns\n",
    "    ----------\n",
    "    nz ,ny, nz : int \n",
    "        Quantum numbers\n",
    "    \"\"\"\n",
    "\n",
    "    q_numbers = np.array([(0,0,0),(0,1,0),(1,0,0),(1,1,0),(2,0,0),(0,2,0),(2,1,0),(1,2,0),(3,0,0),(0,3,0),(2,2,0),(2,2,0),(3,1,0),(1,3,0)])\n",
    "    #q_numbers = np.array([(0,0,0),(0,0,1),(0,1,0),(1,0,0),(0,1,1),(1,0,1),(1,1,0),(1,1,1),(1,1,2),(1,2,1),(2,1,1),(1,2,2),(2,1,2),(2,2,1)])\n",
    "\n",
    "    return q_numbers[k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3229a164",
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrand_2(R, p, r, q, s):\n",
    "    r1 = R[:,0:3]\n",
    "    r2 = R[:,3:6]\n",
    "    r12 = np.sqrt(np.sum((r1 - r2)**2, axis=-1))\n",
    "    \n",
    "    x1 = R[:,0]\n",
    "    y1 = R[:,1]\n",
    "    z1 = R[:,2]\n",
    "    x2 = R[:,3]\n",
    "    y2 = R[:,4]\n",
    "    z2 = R[:,5]\n",
    "    \n",
    "    n_p = index_to_q_numbers(p-1)\n",
    "    n_r = index_to_q_numbers(r-1)\n",
    "    n_q = index_to_q_numbers(q-1)\n",
    "    n_s = index_to_q_numbers(s-1)\n",
    "    \n",
    "    result = 1/r12\n",
    "    \n",
    "    result *= HO_wf_3D(x1, y1, z1, n_p[0], n_p[1], n_p[2])\n",
    "    result *= HO_wf_3D(x2, y2, z2, n_r[0], n_r[1], n_r[2])\n",
    "    result *= HO_wf_3D(x1, y1, z1, n_q[0], n_q[1], n_q[2])\n",
    "    result *= HO_wf_3D(x2, y2, z2, n_s[0], n_s[1], n_s[2])\n",
    "  \n",
    "    \n",
    "    return result\n",
    "\n",
    "def norm_product(p, r, q, s):\n",
    "    return (np.pi/OMEGA_X)**2 * np.pi/OMEGA_Z\n",
    "\n",
    "def two_body_integrand(p, r, q, s):\n",
    "    f = lambda R: integrand_2(R, p, r, q, s)*norm_product(p, r, q, s)\n",
    "    \n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "48b0209d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytical_1(p,q):\n",
    "    n_p = index_to_q_numbers(p-1)\n",
    "    I = (p==q)*(OMEGA_X*(n_p[0] + n_p[1] + 1) + OMEGA_Z*(n_p[2] + 0.5))\n",
    "    \n",
    "    return I\n",
    "\n",
    "def f_cov(p, r, q, s):\n",
    "    \n",
    "    cov = 0.5*np.diag([1/OMEGA_X, 1/OMEGA_X, 1/OMEGA_Z, 1/OMEGA_X, 1/OMEGA_X, 1/OMEGA_Z])\n",
    "    \n",
    "    return cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7acd7ed0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_integrals(N_basis, N_points):\n",
    "    integrals_file = \"integrals_QD_timing.npy\"\n",
    "\n",
    "    normalized_wf = True\n",
    "    \n",
    "    # One- and Two-body integrals\n",
    "    MC_args = {\"f_cov\":f_cov, \"f_integrand\":two_body_integrand, \"N_points\":N_points}\n",
    "    integrals = lib.integral_master(N_basis)\n",
    "    integrals.calculate(integrals_file, analytical_1, MC_args = MC_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6175eba2",
   "metadata": {},
   "source": [
    "### Timing integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "79da8ec3",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Complete integral dictionary for N_points = 100000 and N_basis = 14\n",
      "Calculating two electron integrals: 1/15\n",
      "Calculating two electron integrals: 2/15\n",
      "Calculating two electron integrals: 3/15\n",
      "Calculating two electron integrals: 4/15\n",
      "Calculating two electron integrals: 5/15\n",
      "Calculating two electron integrals: 6/15\n",
      "Calculating two electron integrals: 7/15\n",
      "Calculating two electron integrals: 8/15\n",
      "Calculating two electron integrals: 9/15\n",
      "Calculating two electron integrals: 10/15\n",
      "Calculating two electron integrals: 11/15\n",
      "Calculating two electron integrals: 12/15\n",
      "Calculating two electron integrals: 13/15\n",
      "Calculating two electron integrals: 14/15\n",
      "Total time: 534.507654428482\n"
     ]
    }
   ],
   "source": [
    "# Time integrals\n",
    "N_points = 100000\n",
    "N_basis = 14\n",
    "\n",
    "print(\"Complete integral dictionary for N_points =\", N_points, \"and N_basis =\", N_basis)\n",
    "t0 = time.time()\n",
    "time_integrals(N_basis, N_points)\n",
    "t1 = time.time()\n",
    "print(\"Total time:\", t1-t0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00de1a9e",
   "metadata": {},
   "source": [
    "Total time for N_points = 100000 and N_basis = 14 is \n",
    "\n",
    "T = 557.6338496208191s = 9 min 18 s\n",
    "T = 534.507654428482 = 8 min 54 s"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9182d44",
   "metadata": {},
   "source": [
    "### Timing Self-Consistent Field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "27eb43f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_SCF(N_electrons, integrals, S, max_iter_SCF, eps_SCF):\n",
    "    lib.SCF(N_electrons, integrals, S, max_iter_SCF, eps_SCF)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f35c641e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral file already exists. Not computing the integrals. \n",
      "Self-Consistent Field for  N_electrons = 2\n",
      "E =  3.3149424095424243  | N(SCF) = 1\n",
      "E =  3.7574302193729383  | N(SCF) = 2\n",
      "E =  5.685716036011191  | N(SCF) = 3\n",
      "E =  6.1159252461496205  | N(SCF) = 4\n",
      "E =  6.2385004520760745  | N(SCF) = 5\n",
      "E =  6.252297204851878  | N(SCF) = 6\n",
      "SCF CONVERGED! E =  6.252297204851878\n",
      "Total time: 0.3181483745574951\n",
      "Self-Consistent Field for  N_electrons = 14\n",
      "E =  29.25181809359603  | N(SCF) = 1\n",
      "E =  28.02798604424418  | N(SCF) = 2\n",
      "E =  30.28372819345885  | N(SCF) = 3\n",
      "E =  30.125698598948564  | N(SCF) = 4\n",
      "E =  30.22125700085927  | N(SCF) = 5\n",
      "E =  30.050294331585178  | N(SCF) = 6\n",
      "E =  30.18343674001486  | N(SCF) = 7\n",
      "E =  30.233779114877905  | N(SCF) = 8\n",
      "E =  30.052845226384413  | N(SCF) = 9\n",
      "E =  30.260594841654726  | N(SCF) = 10\n",
      "E =  30.256984232807092  | N(SCF) = 11\n",
      "E =  30.23307767898799  | N(SCF) = 12\n",
      "E =  30.235347359011726  | N(SCF) = 13\n",
      "E =  29.92438185316049  | N(SCF) = 14\n",
      "E =  30.212080269159053  | N(SCF) = 15\n",
      "E =  30.111542884336878  | N(SCF) = 16\n",
      "E =  30.19751707610449  | N(SCF) = 17\n",
      "E =  30.275406620338742  | N(SCF) = 18\n",
      "E =  30.21436099116649  | N(SCF) = 19\n",
      "E =  30.214747382314393  | N(SCF) = 20\n",
      "E =  30.22506626218307  | N(SCF) = 21\n",
      "E =  30.292141819676665  | N(SCF) = 22\n",
      "E =  30.217029988988333  | N(SCF) = 23\n",
      "E =  30.17090025116955  | N(SCF) = 24\n",
      "E =  30.27516429952455  | N(SCF) = 25\n",
      "E =  30.17852111970671  | N(SCF) = 26\n",
      "E =  30.24475780559821  | N(SCF) = 27\n",
      "E =  29.978550205124186  | N(SCF) = 28\n",
      "E =  30.210299605047197  | N(SCF) = 29\n",
      "E =  30.160374488432875  | N(SCF) = 30\n",
      "E =  30.21530057361993  | N(SCF) = 31\n",
      "E =  29.99204100127889  | N(SCF) = 32\n",
      "E =  30.252725869411496  | N(SCF) = 33\n",
      "E =  30.306883674289438  | N(SCF) = 34\n",
      "E =  30.337221365875063  | N(SCF) = 35\n",
      "E =  30.017749990897695  | N(SCF) = 36\n",
      "E =  30.134435388218872  | N(SCF) = 37\n",
      "E =  30.219094859642524  | N(SCF) = 38\n",
      "E =  30.143420931254603  | N(SCF) = 39\n",
      "E =  30.199489379475402  | N(SCF) = 40\n",
      "E =  30.184214703439675  | N(SCF) = 41\n",
      "E =  30.201636863576944  | N(SCF) = 42\n",
      "E =  30.27298150967714  | N(SCF) = 43\n",
      "E =  30.281405260277403  | N(SCF) = 44\n",
      "E =  30.329515145629106  | N(SCF) = 45\n",
      "E =  30.112770180847797  | N(SCF) = 46\n",
      "E =  30.21941448229603  | N(SCF) = 47\n",
      "E =  30.318876074545972  | N(SCF) = 48\n",
      "E =  30.254043718636527  | N(SCF) = 49\n",
      "E =  30.23206284664751  | N(SCF) = 50\n",
      "E =  30.243896304062137  | N(SCF) = 51\n",
      "E =  29.956907269866637  | N(SCF) = 52\n",
      "E =  30.221222274232726  | N(SCF) = 53\n",
      "E =  30.255079167539662  | N(SCF) = 54\n",
      "E =  30.197866236437832  | N(SCF) = 55\n",
      "E =  30.152636138186235  | N(SCF) = 56\n",
      "E =  30.127982040184694  | N(SCF) = 57\n",
      "E =  30.056985263204023  | N(SCF) = 58\n",
      "E =  30.06663709792383  | N(SCF) = 59\n",
      "E =  30.18632918986566  | N(SCF) = 60\n",
      "E =  30.197486077867936  | N(SCF) = 61\n",
      "E =  30.13015037674533  | N(SCF) = 62\n",
      "E =  30.178512272591014  | N(SCF) = 63\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-15-5bd9990f023a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     21\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Self-Consistent Field for  N_electrons =\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mN_electrons\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     22\u001b[0m \u001b[0mt0\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 23\u001b[1;33m \u001b[0mtime_SCF\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mN_electrons\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mintegrals\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mS\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmax_iter_SCF\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0meps_SCF\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     24\u001b[0m \u001b[0mt1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     25\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Total time:\"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mt1\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mt0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-10-6cda67f0c5f0>\u001b[0m in \u001b[0;36mtime_SCF\u001b[1;34m(N_electrons, integrals, S, max_iter_SCF, eps_SCF)\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mtime_SCF\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mN_electrons\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mintegrals\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mS\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmax_iter_SCF\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0meps_SCF\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m     \u001b[0mlib\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSCF\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mN_electrons\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mintegrals\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mS\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmax_iter_SCF\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0meps_SCF\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m~\\Desktop\\Delft\\Computational Physics\\Project 3\\Project3_abermejillo_dbedialaunetar_mserraperalta\\lib.py\u001b[0m in \u001b[0;36mSCF\u001b[1;34m(N_electrons, integrals, S, max_iter_SCF, eps_SCF, max_delta_rho, C, print_E)\u001b[0m\n\u001b[0;32m     44\u001b[0m                         \u001b[0mn_iterations\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     45\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 46\u001b[1;33m                         \u001b[0mF\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcreate_F_matrix\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrho\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mintegrals\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     47\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     48\u001b[0m                         \u001b[1;32mif\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mS\u001b[0m \u001b[1;33m==\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0meye\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mN_basis\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mall\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;31m# orthonormal single basis\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\Desktop\\Delft\\Computational Physics\\Project 3\\Project3_abermejillo_dbedialaunetar_mserraperalta\\lib.py\u001b[0m in \u001b[0;36mcreate_F_matrix\u001b[1;34m(rho, integrals)\u001b[0m\n\u001b[0;32m    100\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    101\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mp\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mN_basis\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 102\u001b[1;33m                 \u001b[1;32mfor\u001b[0m \u001b[0mq\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mN_basis\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    103\u001b[0m                         \u001b[0mF\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mq\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mF\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mq\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mintegrals\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_1\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mq\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# add h matrix\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    104\u001b[0m                         \u001b[1;32mfor\u001b[0m \u001b[0mr\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mN_basis\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Time SCF\n",
    "N_points = 100000\n",
    "N_basis = 14\n",
    "\n",
    "MC_args = {\"f_cov\":f_cov, \"f_integrand\":two_body_integrand, \"N_points\":N_points}\n",
    "integrals_file = \"integrals_QD_timing.npy\"\n",
    "integrals = lib.integral_master(N_basis)\n",
    "integrals.calculate(integrals_file, analytical_1, MC_args = MC_args)\n",
    "S = np.eye(N_basis)\n",
    "max_iter_SCF = 500\n",
    "eps_SCF = 1e-4\n",
    "\n",
    "N_electrons = 2\n",
    "print(\"Self-Consistent Field for  N_electrons =\", N_electrons)\n",
    "t0 = time.time()\n",
    "time_SCF(N_electrons, integrals, S, max_iter_SCF, eps_SCF)\n",
    "t1 = time.time()\n",
    "print(\"Total time:\", t1-t0)\n",
    "\n",
    "N_electrons = 14\n",
    "print(\"Self-Consistent Field for  N_electrons =\", N_electrons)\n",
    "t0 = time.time()\n",
    "time_SCF(N_electrons, integrals, S, max_iter_SCF, eps_SCF)\n",
    "t1 = time.time()\n",
    "print(\"Total time:\", t1-t0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00c29e86",
   "metadata": {},
   "source": [
    "Total time for N_electrons = 2 is T = 5 s and for N_electrons = 14 is T = 4 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "19ee1db8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_SCF_average(N_max, N_av, integrals, S, max_iter_SCF, eps_SCF):\n",
    "    N_pairs = int(N_max/2)\n",
    "\n",
    "    E_i = np.zeros(N_av)\n",
    "    E_N = np.zeros(N_max + 1)\n",
    "    E_Nvar = np.zeros(N_max + 1)\n",
    "\n",
    "    for n in range(1, N_pairs + 1):\n",
    "        N_electrons = n*2\n",
    "        print ('Computing energy for N_electrons = ', N_electrons,'\\r')\n",
    "        for av in range(N_av):\n",
    "            E_i[av] = lib.SCF(N_electrons, integrals, S, max_iter_SCF, eps_SCF, print_E = False)\n",
    "\n",
    "        E_N[N_electrons] = np.average(E_i)\n",
    "        E_Nvar[N_electrons] = np.std(E_i) / np.sqrt(N_av)\n",
    "        print('E({})='.format(N_electrons), E_N[N_electrons])\n",
    "\n",
    "    print(E_N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e3e18c4f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral file already exists. Not computing the integrals. \n",
      "Self-Consistent Field averaging running\n",
      "Computing energy for N_electrons =  2 \n",
      "E(2)= 6.141950602156464\n",
      "Computing energy for N_electrons =  4 \n",
      "E(4)= 9.808849861384923\n",
      "Computing energy for N_electrons =  6 \n",
      "E(6)= 12.761459970048762\n",
      "Computing energy for N_electrons =  8 \n",
      "E(8)= 17.239099725630957\n",
      "Computing energy for N_electrons =  10 \n",
      "E(10)= 20.47351884546053\n",
      "Computing energy for N_electrons =  12 \n",
      "E(12)= 24.61423268061165\n",
      "Computing energy for N_electrons =  14 \n",
      "E(14)= 30.160809702750406\n",
      "[ 0.          0.          6.1419506   0.          9.80884986  0.\n",
      " 12.76145997  0.         17.23909973  0.         20.47351885  0.\n",
      " 24.61423268  0.         30.1608097 ]\n",
      "E =  26.35570230317822  | N(SCF) = 1\n",
      "E =  28.624753615026673  | N(SCF) = 2\n",
      "E =  30.417343763597454  | N(SCF) = 3\n",
      "E =  30.025647365725412  | N(SCF) = 4\n",
      "E =  30.184000278787234  | N(SCF) = 5\n",
      "E =  30.199150528017377  | N(SCF) = 6\n",
      "E =  30.207685310576586  | N(SCF) = 7\n",
      "SCF CONVERGED! E =  30.207685310576586\n",
      "Total time: 530.4971022605896\n"
     ]
    }
   ],
   "source": [
    "# Time SCF averaging\n",
    "N_basis = 14\n",
    "N_points = 100000\n",
    "\n",
    "MC_args = {\"f_cov\":f_cov, \"f_integrand\":two_body_integrand, \"N_points\":N_points}\n",
    "integrals_file = \"integrals_QD_timing.npy\"\n",
    "integrals = lib.integral_master(N_basis)\n",
    "integrals.calculate(integrals_file, analytical_1, MC_args = MC_args)\n",
    "S = np.eye(N_basis)\n",
    "max_iter_SCF = 500\n",
    "eps_SCF = 1e-4\n",
    "\n",
    "N_max = 14\n",
    "N_av = 20 #times we do the SCF to do statistics\n",
    "\n",
    "print(\"Self-Consistent Field averaging running\")\n",
    "t0 = time.time()\n",
    "time_SCF_average(N_max, N_av, integrals, S, max_iter_SCF, eps_SCF)\n",
    "t1 = time.time()\n",
    "print(\"Total time:\", t1-t0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55caa266",
   "metadata": {},
   "source": [
    "Total time for the averaging of the SCF is \n",
    "\n",
    "T = 530.497 s = 8 min 50 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb8c71d0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}