{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3484c3bd",
   "metadata": {},
   "source": [
    "# Evaluate charge density at the atomic nuclei"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10f4a3b5",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Do an electronic structure calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "22b8911d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-16T20:17:35.571994Z",
     "start_time": "2023-11-16T20:17:35.317558Z"
    }
   },
   "outputs": [],
   "source": [
    "from pymolpro import Project\n",
    "\n",
    "p = Project(\"evaluate_charge_density\")\n",
    "p.write_input('''\n",
    "basis,cc-pVDZ\n",
    "geometry={f;h,f,1.732}\n",
    "set,sewprop=0\n",
    "gexpec,delta,f\n",
    "gexpec,delta,h\n",
    "rhf\n",
    "{casscf;closed,2;state,3;natorb,2251.2,state=1.1}\n",
    "{casscf;closed,2;state,3;natorb,2252.2,state=2.1}\n",
    "{casscf;closed,2;state,3;natorb,2253.2,state=3.1}\n",
    "''')\n",
    "p.run(wait=True)\n",
    "assert p.status == 'completed'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c25f4fb1",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Evaluate occupied orbitals then density on a grid at the nuclear positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8a1dce6a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-16T20:17:35.664496Z",
     "start_time": "2023-11-16T20:17:35.643011Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 atoms:\n",
      "a1 F [0.0, 0.0, -0.08725992951335484]\n",
      "a2 H [0.0, 0.0, 1.644740070486644]\n",
      "Calculated densities at the nuclei\n",
      " [[4.28636637e+02 3.39348454e-01]\n",
      " [4.29484632e+02 2.80599128e-01]\n",
      " [4.28693250e+02 1.75420682e-01]\n",
      " [4.31473378e+02 3.58104652e-01]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "atoms = p.geometry()\n",
    "print(str(len(atoms)) + ' atoms:')\n",
    "for atom in atoms:\n",
    "    print(atom['id'], atom['elementType'], atom['xyz'])\n",
    "\n",
    "points = [atom['xyz'] for atom in atoms]\n",
    "\n",
    "results = []\n",
    "for state in range(len(p.xpath('//orbitals'))):\n",
    "    orbitalsAtPoints = p.evaluateOrbitals(points, instance=state)\n",
    "    results.append([0.0 for point in points])\n",
    "    for orbital in orbitalsAtPoints:\n",
    "        for ipoint in range(len(orbital['values'])):\n",
    "            results[-1][ipoint] += orbital['values'][ipoint] ** 2 * orbital['occ']\n",
    "\n",
    "print('Calculated densities at the nuclei\\n', np.array(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1021bbc1",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Check that the densities at the nuclei agree with values calculated directly in Molpro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ba6c7c79",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-16T20:17:35.702185Z",
     "start_time": "2023-11-16T20:17:35.666809Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 F 428.636637409104 428.63664514377 -7.734665985026368e-06\n",
      "1 H 0.33934845419587656 0.339348457122595 -2.926718423168495e-09\n",
      "2 F 429.4846324205379 429.484605946402 2.647413589329517e-05\n",
      "2 H 0.2805991278930762 0.280599140980421 -1.308734481897389e-08\n",
      "3 F 428.693249783628 428.693250939239 -1.1556110166566214e-06\n",
      "3 H 0.17542068187713988 0.175420624927375 5.6949764876135234e-08\n",
      "4 F 431.4733780183073 431.473358358035 1.9660272300825454e-05\n",
      "4 H 0.35810465239098177 0.358104654722505 -2.3315232522413964e-09\n"
     ]
    }
   ],
   "source": [
    "for state in range(len(p.xpath('//orbitals'))):\n",
    "    for atom in atoms:\n",
    "        search_ = '(//property[@name=\"DELTA(' + atom['elementType'].upper() + ')\"])[' + str(state + 1) + ']'\n",
    "        value = float(p.xpath_search(search_, 'value')[0])\n",
    "        print(state+1, atom['elementType'],results[state][atoms.index(atom)],value,results[state][atoms.index(atom)]-value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f975a13-1e60-4dec-9da0-ccdf0f66faf4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.13.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}