check_convergence.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2310a17c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cda1dd0",
   "metadata": {},
   "source": [
    "# Checks for the convergence of the two electron integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "14b61966",
   "metadata": {},
   "outputs": [],
   "source": [
    "# LOAD PARAMETERS FOR MONTE CARLO INTEGRATION\n",
    "\n",
    "import MC_integration as mc\n",
    "import basis_set as bs\n",
    "\n",
    "L_size = 5\n",
    "N_walkers = 400\n",
    "N_steps = 10000\n",
    "N_skip = 1000\n",
    "trial_move = 0.6\n",
    "\n",
    "integrand = bs.two_body_integrand\n",
    "dimension = 6\n",
    "sampling = bs.sampling_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "51ca5c14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GENERATE RANDOM STEPS\n",
    "\n",
    "indices = np.array([1,1,1,1])\n",
    "\n",
    "init_points = mc.rand_init_point(L_size, dimension, N_walkers)\n",
    "steps, _, acceptance_ratio = mc.random_walkers(sampling, N_steps, init_points, trial_move)\n",
    "steps = steps[N_skip:, :, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "2643e371",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of distances smaller than 1E-5:  0\n"
     ]
    }
   ],
   "source": [
    "# CHECK IF THERE ARE DISTANCES SMALLER THAN A CERTAIN VALUE (FOR CONVERGENCE OF 1/r12)\n",
    "\n",
    "r1, r2 = steps[:,:,0:3], steps[:,:,3:6]\n",
    "dist = np.sum((r1 - r2)**2, axis=2)\n",
    "dist = dist.flatten() # 1D vector\n",
    "\n",
    "print(\"Number of distances smaller than 1E-5: \", len(np.where(dist <= 1E-5)[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "20c96b3d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[inf inf inf inf inf inf]\n"
     ]
    }
   ],
   "source": [
    "# CHECK IF THERE ARE DIVERGENCES IN E_LOCAL()\n",
    "\n",
    "E_alpha_walkers = mc.MC_average_walkers(integrand, steps, indices)\n",
    "\n",
    "print(E_alpha_walkers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "972406aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of infinities:  66\n"
     ]
    }
   ],
   "source": [
    "# CHECK IF THERE ARE DIVERGENCES IN E_LOCAL()\n",
    "\n",
    "E_local = integrand(steps, indices)\n",
    "\n",
    "print(\"Number of infinities: \", len(np.where(E_local == np.inf)[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b57d773a",
   "metadata": {},
   "source": [
    "Therefore, the error is in the expression of the integrand (two_body_integrand)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "893989a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PREVIOUS INTEGRAND\n",
    "def two_body_integrand(R, indices):\n",
    "\t\"\"\"\n",
    "\tReturns the value of the <pr|g|qs> integrand evaluated at R\n",
    "\n",
    "\tParameters\n",
    "\t----------\n",
    "\n",
    "\tR: np.ndarray(6,N)\n",
    "\t\tPositions of the two electrons x1,y1,z1,x2,y2,z2 \n",
    "\t\tat N diferent coordinates (to evaluate random walkers simultaneously)\n",
    "\tindices: np.ndarray(4)\n",
    "\t\tIndices of the basis functions \n",
    "\n",
    "\tReturns\n",
    "\t----------\n",
    "\tI : np.ndarray(N)\n",
    "\t\t\n",
    "\t\"\"\"\n",
    "\tp,q,r,s = indices[0],indices[1],indices[2],indices[3]\n",
    "\n",
    "\talpha, beta, gamma, delta = basis_radius(p), basis_radius(q), basis_radius(r), basis_radius(s)\n",
    "\n",
    "\ta = alpha + beta\n",
    "\tb = gamma + delta\n",
    "\n",
    "\tR[0:3] = a**(-0.5)*R[0:3]\n",
    "\tR[3:6] = b**(-0.5)*R[3:6]\n",
    "\n",
    "\tr1 = np.sqrt(R[0]**2 + R[1]**2 + R[2]**2)\n",
    "\tr2 = np.sqrt(R[3]**2 + R[4]**2 + R[5]**2)\n",
    "\tr12 = abs(r1-r2)\n",
    "\tI = (a*b)**(-1.5)*HO_wf_3D_basis(R[0:3],p)*HO_wf_3D_basis(R[3:6],r)*(1/r12)*HO_wf_3D_basis(R[0:3],q)*HO_wf_3D_basis(R[3:6],s)/sampling_function(R[0:6])\n",
    "\n",
    "\treturn I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ed3e529b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.5066282746309887, 1.5957691216057306, 1.5957691216057306, 1.5957691216057306]\n"
     ]
    }
   ],
   "source": [
    "# CHECK THE BASIS RADIUS FUNCTION\n",
    "\n",
    "r = [basis_radius(i+1) for i in range(4)]\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "4c3b31f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# NEW INTEGRAND\n",
    "from basis_set import *\n",
    "\n",
    "def two_body_integrand_new(R, indices):\n",
    "\n",
    "    p,q,r,s = indices\n",
    "    alpha, beta, gamma, delta = basis_radius(p), basis_radius(q), basis_radius(r), basis_radius(s)\n",
    "\n",
    "    a = alpha + beta\n",
    "    b = gamma + delta\n",
    "\n",
    "    R[:,:,0:3] = a**(-0.5)*R[:,:,0:3]\n",
    "    R[:,:,3:6] = b**(-0.5)*R[:,:,3:6]\n",
    "\n",
    "    r1 = np.sqrt(R[:,:,0]**2 + R[:,:,1]**2 + R[:,:,2]**2)\n",
    "    r2 = np.sqrt(R[:,:,3]**2 + R[:,:,4]**2 + R[:,:,5]**2)\n",
    "    r12 = abs(r1-r2)\n",
    "    \n",
    "    print(\"Number of infinities in 1/r12: \", len(np.where(1/r12 == np.inf)[0]))\n",
    "    print(\"Number of infinities in 1/sampling_function: \", len(np.where(1/sampling_function_new(R[:,:,0:6]) == np.inf)[0]))\n",
    "    print(\"Number of infinities in wf(r1): \", len(np.where(HO_wf_3D_basis_new(R[:,:,0:3],p) == np.inf)[0]))\n",
    "    print(\"Number of infinities in wf(r2): \", len(np.where(HO_wf_3D_basis_new(R[:,:,3:6],p) == np.inf)[0]))\n",
    "    \n",
    "    I = (a*b)**(-1.5)*HO_wf_3D_basis_new(R[:,:,0:3],p)*HO_wf_3D_basis_new(R[:,:,3:6],r)*(1/r12)*HO_wf_3D_basis_new(R[:,:,0:3],q)*HO_wf_3D_basis_new(R[:,:,3:6],s)/sampling_function_new(R[:,:,0:6])\n",
    "\n",
    "    return I\n",
    "\n",
    "\n",
    "def HO_wf_3D_basis_new(R, k, omega_x=OMEGA_X, omega_y=OMEGA_X, omega_z=OMEGA_Z):\n",
    "\n",
    "    nx, ny, nz = index_to_q_numbers(k)\n",
    "\n",
    "    psi = HO_wf_3D(R[:,:,0], R[:,:,1], R[:,:,2], nx, ny, nz, omega_x, omega_y, omega_z)\n",
    "\n",
    "    return psi\n",
    "\n",
    "def sampling_function_new(R):\n",
    "\n",
    "    value = 1/(2*np.pi)**3*np.exp(-0.5*np.sum(R**2,axis=-1))\n",
    "\n",
    "    return value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "6b7b89bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of infinities in 1/r12:  0\n",
      "Number of infinities in 1/sampling_function:  0\n",
      "Number of infinities in wf(r1):  0\n",
      "Number of infinities in wf(r2):  0\n",
      "Number of infinities:  0\n"
     ]
    }
   ],
   "source": [
    "# CHECK IF THERE ARE DIVERGENCES IN E_LOCAL() WITH NEW INTEGRAND\n",
    "\n",
    "E_local = two_body_integrand_new(steps, indices)\n",
    "\n",
    "print(\"Number of infinities: \", len(np.where(E_local == np.inf)[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "d506fe51",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.523494833996234e-15\n"
     ]
    }
   ],
   "source": [
    "E_alpha_walkers = np.average(E_local, axis=0)\n",
    "print(np.average(E_alpha_walkers))"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}