{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### H2+ PDFT Inversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from CADMium import Pssolver\n",
    "from CADMium import Psgrid\n",
    "from CADMium import Partition\n",
    "from CADMium import Inverter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: If len(KS) > 1 Has not been migrated from matlab\n"
     ]
    }
   ],
   "source": [
    "a = 2.0/2\n",
    "#Nuclear charge for fragments A and B\n",
    "Za, Zb = 1,1\n",
    "#Set polarization 1-Unpolarized, 2-Polarized|\n",
    "pol = 2\n",
    "#Fragment a electrons [alpha, beta]\n",
    "Nmo_a = [[1  ,0]] #Number of molecular orbitals to calculate\n",
    "N_a   = [[0.5,0]]\n",
    "#Ensemble mix\n",
    "nu_a = 1\n",
    "#Fragment b electrons\n",
    "Nmo_b = [[1  ,0]]\n",
    "N_b   = [[0.5,0]]\n",
    "#Ensemble mix\n",
    "nu_b = 1\n",
    "\n",
    "#Molecular elctron configuration\n",
    "Nmo_m = [[1,1]]\n",
    "N_m   = [[1,1]]\n",
    "\n",
    "#Set up grid\n",
    "NP = 7\n",
    "NM = [4,4]\n",
    "L = np.arccosh(10/a)\n",
    "loc = np.array(range(-4,5)) #Stencil outline\n",
    "\n",
    "grid = Psgrid(NP, NM, a, L, loc)\n",
    "grid.initialize()\n",
    "\n",
    "part = Partition(grid, Za, Zb, pol, Nmo_a, N_a, nu_a, Nmo_b, N_b, nu_b, {  \"ab_sym\"            : True,\n",
    "                                                                           \"ens_spin_sym\"      : False,  \n",
    "                                                                           \"kinetic_part_type\" : \"libxcke\",\n",
    "                                                                           \"k_family\"          : \"gga\",\n",
    "                                                                           \"ke_func_id\"        : 500,\n",
    "                                                                           \"interaction_type\"  : \"ni\",\n",
    "                                                                           \"fractional\"        : True,\n",
    "                                                                            })\n",
    "\n",
    "#Setup inverter object\n",
    "mol_solver = Pssolver(grid, Nmo_m, N_m)\n",
    "part.inverter = Inverter(grid, mol_solver, {  \"AB_SYM\"         : True,\n",
    "                                              \"ENS_SPIN_SYM\"   : False,  \n",
    "                                              \"use_iterative\"  : False,\n",
    "                                              \"invert_type\"    : \"wuyang\",\n",
    "                                              \"Tol_lin_solver\" : 1e-3,\n",
    "                                              \"disp\"           : True,  \n",
    "                                            })\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                  Total Energy            \n",
      "\n",
      "                __________________        \n",
      "\n",
      "Iteration         A            B              res     \n",
      "\n",
      "_______________________________________________________\n",
      "\n",
      "    1           -0.25000     -0.25000       1.000e+00 \n",
      "    2           -0.25000     -0.25000       8.199e-10 \n"
     ]
    }
   ],
   "source": [
    "### Perform Isolated Fragment Calculation. \n",
    "part.optPartition.isolated = True\n",
    "part.scf({\"disp\"     : True,\n",
    "          \"e_tol\"    : 1e-7})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                  Total Energy            \n",
      "\n",
      "                __________________        \n",
      "\n",
      "Iteration         A            B              res     \n",
      "\n",
      "_______________________________________________________\n",
      "\n",
      "    1           -0.20949     -0.20949       1.000e+00 \n",
      "    2           -0.21997     -0.21997       7.238e-02 \n",
      "    3           -0.22295     -0.22295       1.886e-02 \n",
      "    4           -0.22339     -0.22339       2.846e-03 \n",
      "    5           -0.22347     -0.22347       5.259e-04 \n",
      "    6           -0.22348     -0.22348       9.523e-05 \n",
      "    7           -0.22349     -0.22349       1.723e-05 \n",
      "    8           -0.22349     -0.22349       3.114e-06 \n",
      "    9           -0.22349     -0.22349       5.623e-07 \n",
      "   10           -0.22349     -0.22349       1.015e-07 \n",
      "   11           -0.22349     -0.22349       1.831e-08 \n"
     ]
    }
   ],
   "source": [
    "### Perform PDFT Calculation. \n",
    "\n",
    "part.optPartition.isolated = False\n",
    "part.scf({'iterative' : False,\n",
    "          'disp'      : True,\n",
    "           'continuing' : True})"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}