What's wrong with my attempt to reproduce the dipole moment?

1

I want to compute the dipole moment after calculating the HF orbitals, but I get different results with psi4.core.variable. I tried to follow this tutorial and have taken into account that now the variable in SCF DIPOLE now and that it is in a.u. by default. I am using Psi4 version 1.6.1. Here is the Python code:

import time
import numpy as np
np.set_printoptions(precision=5, linewidth=200, suppress=True)
import psi4

# Memory for Psi4 in GB
psi4.core.set_output_file('output.dat', False)

# Memory for numpy in GB
numpy_memory = 2

mol = psi4.geometry("""
O
H 1 1.1
H 1 1.1 2 104
symmetry c1
""")


psi4.set_options({'basis': '3-21g',
                  'scf_type': 'pk',
                  'e_convergence': 1e-8,
                  'd_convergence': 1e-8})

# First compute SCF energy using Psi4
scf_e, wfn = psi4.energy('SCF', return_wfn=True)

print(f'Energy = {scf_e-mol.nuclear_repulsion_energy()}')

mints = psi4.core.MintsHelper(wfn.basisset())
Da = np.asarray(wfn.Da())
Db = np.asarray(wfn.Db())
D = Da + Db

Dx, Dy, Dz = mints.ao_dipole()

# Get nuclear dipole
nuc_dipole = mol.nuclear_dipole()

# Compute dipole moment components
mux = -np.einsum('pq, pq ->', D, Dx, optimize=True)
muy = -np.einsum('pq, pq ->', D, Dy, optimize=True)
muz = -np.einsum('pq, pq ->', D, Dz, optimize=True)

mux += nuc_dipole[0]
muy += nuc_dipole[1]
muz += nuc_dipole[2]

mu = np.array([mux, muy, muz])

mu_psi4 = psi4.core.variable('SCF DIPOLE')

print(mu_psi4) 
print(mu)
# output
# [ 0.     -0.      0.9497]
#[-0.       0.       1.30488]
2

There’s a sign error.

The tutorial you’re following constructed the integrals themselves. You’re getting your integrals from Psi4. Those integrals disagree by a factor of -1, so the corresponding term needs to disagree by a factor of -1. Get rid of the negative signs on your einsum, and you’ll match Psi4.