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-pVQZ
geometry={f;h,f,1.732}
set,sewprop=0
gexpec,delta,f
gexpec,delta,h
{casscf;closed,2;state,3;natorb,2251.2,state=1.1};put,xml
{casscf;closed,2;state,3;natorb,2252.2,state=2.1};put,xml
{casscf;closed,2;state,3;natorb,2253.2,state=3.1};put,xml
''')
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.40214461e+02 3.42402191e-01]
 [4.39554369e+02 2.14828887e-01]
 [4.42399757e+02 4.06209932e-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:
        value = float(
            p.xpath('//property[@name="DELTA(' + atom['elementType'] + ')" and @stateNumber="' + str(state + 1) + '"]')[
                0].xpath('@value')[0])
        print(state+1, atom['elementType'],results[state][atoms.index(atom)],value,results[state][atoms.index(atom)]-value)

1 F 440.2144610686323 440.214461088477 -1.9844662801915547e-08
1 H 0.34240219096060126 0.342402151500937 3.945966425833447e-08
2 F 439.5543685301112 439.554369380541 -8.50429785259621e-07
2 H 0.21482888674378783 0.2148289453077 -5.856391216418899e-08
3 F 442.3997566224122 442.399756622169 2.432329893053975e-10
3 H 0.4062099317051013 0.406209931733528 -2.842670543401482e-11