{ "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-pVQZ\n", "geometry={f;h,f,1.732}\n", "set,sewprop=0\n", "gexpec,delta,f\n", "gexpec,delta,h\n", "{casscf;closed,2;state,3;natorb,2251.2,state=1.1};put,xml\n", "{casscf;closed,2;state,3;natorb,2252.2,state=2.1};put,xml\n", "{casscf;closed,2;state,3;natorb,2253.2,state=3.1};put,xml\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": { "collapsed": false, "jupyter": { "outputs_hidden": false }, "pycharm": { "name": "#%%\n" }, "ExecuteTime": { "end_time": "2023-11-16T20:17:35.664496Z", "start_time": "2023-11-16T20:17:35.643011Z" } }, "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.40214461e+02 3.42402191e-01]\n", " [4.39554369e+02 2.14828887e-01]\n", " [4.42399757e+02 4.06209932e-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" ], "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " 2023.1\n", " \n", " 2023-11-16T16:53:05\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " -1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1\n", " \n", " \n", " 1 0 0 0 0 -1 0 0 0 0 1 0 0 0 0 1\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "