I would like to do compute myself the dipole moment with a finite-field method. I have a working version with psi4. In my code, I only have DIIS to accelerate the convergence of my code. The dipole moment perturbation is included in the Fock operator and is diagonalized in the SCF loop. I get small differences between the obtained energy in my implementation and the energy obtained via psi4. It makes the finite-field approximation of the dipole energy wrong. For comparison with a sto-3g basis:
# Psi4 dipole moment
[ 1.39698386e-09 -1.39698386e-09 6.03521254e-01]
# My dipole moment
[ 1.86264515e-09 4.65661287e-10 -2.61884902e-01]
. My questions are:
- what is the energy criteria for convergence? Absolute/relative error with respect to previous iterate?
- How is treated the dipole moment perturbation in the computation of the energy? Should the perturbation be treated with MP theory?
So the outline of my code to compute the energy for a given field strength is:
import time
import numpy as np
from utils import *
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
""")
# parts of code not need for the explanation...
# SCF DIIS procedure
"""compute RHF energy with a DIIS SCF loop"""
# PERTURBATION OPERATORS
# h is a list of length 3
# DFx, DFy, DFz is the dipole moment in atomic basis
DF = h[0] * DFx + h[1] * DFy + h[2] * DFz
D = np.zeros_like(H_core)
E_old = 0 #np.inf
# INITIAL CORE GUESS
F_total = A.dot(H_core).dot(A)
_, C_p = np.linalg.eigh(F_total)
C = A.dot(C_p)
C_occ = C[:, :ndocc]
D = np.einsum('pi,qi->pq', C_occ, C_occ, optimize=True)
# trial + residuals vector list
F_list = []
DIIS_res_list = []
# SCF LOOP
for i in range(MAX_ITER):
# total Hamiltonian
J = np.einsum('pqrs,rs->pq', I, D, optimize=True)
K = np.einsum('prqs,rs->pq', I, D, optimize=True)
F = H_core + 2*J - K
F_total = F + DF
# save residuals and trial of the Fock matrix
DIIS_res = A.dot(F_total.dot(D).dot(S) - S.dot(D).dot(F_total)).dot(A)
F_list.append(F_total)
DIIS_res_list.append(DIIS_res)
# check convergence
E = np.einsum('pq,pq->', (H_core + F_total), D, optimize=True)
dE = E - E_old
dRMS = np.mean(DIIS_res**2)**0.5
if verbose:
print('SCF Iteration %3d: Energy = %4.16f dE = % 1.5E dRMS = %1.5E' % (i, E+8.0023552092614807, abs(dE), dRMS))
# E_eps = D_eps = 1e-8
if (abs(dE) < E_eps) and (dRMS < D_eps):
if verbose:
print(f'Converged in {i} iterations!')
return E
E_old = E
#"""
if i>0:
# Build B matrix
B_dim = len(F_list) + 1
B = np.empty((B_dim, B_dim))
B[-1, :] = -1
B[:, -1] = -1
B[-1, -1] = 0
for i in range(len(F_list)):
for j in range(len(F_list)):
B[i, j] = np.einsum('ij,ij->', DIIS_res_list[i], DIIS_res_list[j], optimize=True)
# normalize
B[abs(B) < 1.e-14] = 1.e-14
B[:-1, :-1] /= np.abs(B[:-1, :-1]).max()
# Build RHS of Pulay equation
rhs = np.zeros((B_dim))
rhs[-1] = -1
# Solve Pulay equation for c_i's with NumPy
coeff = np.linalg.solve(B, rhs)
# Build DIIS Fock matrix
F_total = np.zeros_like(F)
for x in range(coeff.shape[0] - 1):
F_total += coeff[x] * F_list[x]
#"""
# diagonalize Hamiltonian with the perturbation
F_total = A.dot(F_total).dot(A)
_, C_p = np.linalg.eigh(F_total)
C = A.dot(C_p)
# new density matrix
C_occ = C[:, :ndocc]
D = np.einsum('pi,qi->pq', C_occ, C_occ, optimize=True)
print("Not converged...")
return E
Psi4 working finite code field:
import time
import numpy as np
np.set_printoptions(precision=8, 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
""")
basis_list = ['sto-3g']
n_pts = 2
h = 1e-5
pts = np.array([h, -h])
weights = 1/2/h*np.array([1, -1])
for basis in basis_list:
psi4.set_options({'basis': basis,
'scf_type': 'pk',
'e_convergence': 1e-8,
'd_convergence': 1e-8,
})
dipole_moment = np.zeros(3)
for p in range(n_pts):
for i in range(3):
Dh = [0, 0, 0]
Dh[i] = pts[p]
psi4.set_options({
'perturb_h': True,
'perturb_with': 'dipole',
'perturb_dipole':Dh,
'guess' : 'CORE',
'SCF_INITIAL_ACCELERATOR': 'none',
'SOSCF_MAX_ITER': 0,
})
E_perturbed = (psi4.energy('scf') - mol.nuclear_repulsion_energy())
dipole_moment[i] += weights[p] * E_perturbed
print("basis = ", basis)
print(dipole_moment)
print('-'*5)
I deactivated the soscf iterations to have the closest implementation to my code. For h = [0, 0, 5e-4], the respective output of the DIIS SCF iterations for psi4 and mine are:
Total Energy Delta E RMS |[F,P]|
@RHF iter 1: -73.28582353882528 -7.32858e+01 1.20729e-01 DIIS
@RHF iter 2: -74.82814071610819 -1.54232e+00 5.08273e-02 DIIS
@RHF iter 3: -74.93872767506176 -1.10587e-01 7.31197e-03 DIIS
@RHF iter 4: -74.94185673694147 -3.12906e-03 1.73504e-03 DIIS
@RHF iter 5: -74.94206033343640 -2.03596e-04 5.26032e-04 DIIS
@RHF iter 6: -74.94208580610406 -2.54727e-05 4.36830e-05 DIIS
@RHF iter 7: -74.94208593404171 -1.27938e-07 2.95119e-07 DIIS
@RHF iter 8: -74.94208593404690 -5.18696e-12 8.44139e-09 DIIS
and
SCF Iteration 0: Energy = -73.2858158156854529 dE = 8.12882E+01 dRMS = 1.20729E-01
SCF Iteration 1: Energy = -74.8281477990109067 dE = 1.54233E+00 dRMS = 5.08272E-02
SCF Iteration 2: Energy = -74.9387306821756596 dE = 1.10583E-01 dRMS = 7.31198E-03
SCF Iteration 3: Energy = -74.9418593505401844 dE = 3.12867E-03 dRMS = 1.73504E-03
SCF Iteration 4: Energy = -74.9420629619417298 dE = 2.03611E-04 dRMS = 5.26032E-04
SCF Iteration 5: Energy = -74.9420884273501002 dE = 2.54654E-05 dRMS = 4.36830E-05
SCF Iteration 6: Energy = -74.9420885527631100 dE = 1.25413E-07 dRMS = 2.95118E-07
SCF Iteration 7: Energy = -74.9420885527437406 dE = 1.93694E-11 dRMS = 8.44139E-09
Converged in 7 iterations!