Evaluate charge density at the atomic nuclei

Do an electronic structure calculation

[1]:
from pymolpro import Project

p = Project("evaluate_charge_density")
p.write_input('''
basis,cc-pVDZ
geometry={f;h,f,1.732}
set,sewprop=0
gexpec,delta,f
gexpec,delta,h
rhf
{casscf;closed,2;state,3;natorb,2251.2,state=1.1}
{casscf;closed,2;state,3;natorb,2252.2,state=2.1}
{casscf;closed,2;state,3;natorb,2253.2,state=3.1}
''')
p.run(wait=True)
assert p.status == 'completed'

Evaluate occupied orbitals then density on a grid at the nuclear positions

[2]:
import numpy as np

atoms = p.geometry()
print(str(len(atoms)) + ' atoms:')
for atom in atoms:
    print(atom['id'], atom['elementType'], atom['xyz'])

points = [atom['xyz'] for atom in atoms]

results = []
for state in range(len(p.xpath('//orbitals'))):
    orbitalsAtPoints = p.evaluateOrbitals(points, instance=state)
    results.append([0.0 for point in points])
    for orbital in orbitalsAtPoints:
        for ipoint in range(len(orbital['values'])):
            results[-1][ipoint] += orbital['values'][ipoint] ** 2 * orbital['occ']

print('Calculated densities at the nuclei\n', np.array(results))
2 atoms:
a1 F [0.0, 0.0, -0.08725992951335484]
a2 H [0.0, 0.0, 1.644740070486644]
Calculated densities at the nuclei
 [[4.28636637e+02 3.39348454e-01]
 [4.29484632e+02 2.80599128e-01]
 [4.28693250e+02 1.75420682e-01]
 [4.31473378e+02 3.58104652e-01]]

Check that the densities at the nuclei agree with values calculated directly in Molpro

[3]:
for state in range(len(p.xpath('//orbitals'))):
    for atom in atoms:
        search_ = '(//property[@name="DELTA(' + atom['elementType'].upper() + ')"])[' + str(state + 1) + ']'
        value = float(p.xpath_search(search_, 'value')[0])
        print(state+1, atom['elementType'],results[state][atoms.index(atom)],value,results[state][atoms.index(atom)]-value)
1 F 428.636637409104 428.63664514377 -7.734665985026368e-06
1 H 0.33934845419587656 0.339348457122595 -2.926718423168495e-09
2 F 429.4846324205379 429.484605946402 2.647413589329517e-05
2 H 0.2805991278930762 0.280599140980421 -1.308734481897389e-08
3 F 428.693249783628 428.693250939239 -1.1556110166566214e-06
3 H 0.17542068187713988 0.175420624927375 5.6949764876135234e-08
4 F 431.4733780183073 431.473358358035 1.9660272300825454e-05
4 H 0.35810465239098177 0.358104654722505 -2.3315232522413964e-09
[ ]: