# DFT (from scratch) doesn't converge to the correct result

Hi,

I am currently trying to reproduce a DFT code using all features of PSI4 and the wavefunction object.
I believe I managed to reproduce the DFT self-consistent procedure, tested on Hydrogen chains (for instance H8 in the example here). For R<1.5, everything seems OK as far as I can see, but for larger interatomic distances it starts to deteriorate quite badly. It seems that the self-consistent procedure gets stuck in another minima where the KS orbitals get degenerate (which is however not the case with the psi4 reference calculation!).

Actually, one of my main concern is that my energy also goes below the one of psi4, thus violating the variational principle…

As a sanity check, I also computed the DFT energy by using the Hxc potential computed with the AO density matrix of the DFT wavefunction object of psi4, and I obtained the correct energy with an accuracy of 10^{-7} Hartree.

Do you have any clue of why I get this behaviour ? I want to get the exact same convergence as in PSI4. I enclose my code just below.

Thank you !

``````import numpy as np
import math
import psi4

functional   = "SVWN"
basis        = "sto-3g"
thrs_energy  = 1e-5
thrs_density = 1e-5
SCF_maxiter  = 100

# Define the system under study: (n_orbs hydrogens chain)
n_orbs       = 8
n_elec       = n_orbs
n_occ        = n_elec//2
R            = 2.5

psi4.core.set_output_file("H{}_R{}_{}_{}_Psi4.dat".format(n_orbs,R,basis,functional),True)
# Define the geometry
string_geo  = "0 1\n"
for d in range(n_occ):
string_geo += "H 0. 0. {}\n".format(- (R/2. + d*R))
string_geo += "H 0. 0. {}\n".format(+ (R/2. + d*R))
string_geo += "symmetry c1\n"
string_geo += "nocom\n"
string_geo += "noreorient\n"
psi4.geometry(string_geo)

psi4.set_options({'basis': basis,'save_jk':True})

# Perform the DFT calculation with PSI4 for comparison
dft_e, dft_wfn = psi4.energy(functional, return_wfn=True)

# Extract the important informations from the wavefunction
# Hcore and J (Coulomb) matrices:
Hcore = dft_wfn.H().clone()
E_nuc = dft_wfn.get_energies('Nuclear')
J_coulomb = dft_wfn.jk().J()
# Density matrix in the AO basis:
D_AO = dft_wfn.Da().clone()
D_AO_new = dft_wfn.Da().clone()
# Overlap matrix in the AO basis:
S_AO = dft_wfn.S().np
# KS-MO coefficient matrix in the AO basis:
C_AO = dft_wfn.Ca().np

# Transformation to the OAO basis
# Compute the inverse square root of the overlap matrix S
S_eigval, S_eigvec = np.linalg.eigh(S_AO)
S_sqrt_inv = S_eigvec @ np.diag((S_eigval)**(-1./2.)) @ S_eigvec.T
C_transformation = np.linalg.inv(S_sqrt_inv)
C_OAO = C_transformation @ C_AO

# Initialize the potential objects
V_xc = dft_wfn.Da().clone()
Vpotential = dft_wfn.V_potential()

# Construct the SAD Guess for the initial density D_AO
psi4.core.prepare_options_for_module("SCF")
psi4.core.get_global_option("BASIS"), puream=dft_wfn.basisset().has_puream(),
return_atomlist=True)
return_atomlist=True)

iteration   = 0
Etot        = 0
Delta_E     = 1e8
convergence = True

while ( Delta_E > thrs_energy or normocc > thrs_density) and convergence:

# Compute the KS potential with the new density matrix in the AO basis
Vpotential.set_D([D_AO])
Vpotential.compute_V([V_xc])

# Compute the Fock matrix in the OAO basis
F_AO = Hcore.np + 2*J_coulomb.np + V_xc.np
F_OAO = S_sqrt_inv @ F_AO @ S_sqrt_inv.T

# Diagonalize F_OAO to get the KS orbitals and orbital energies
KS_energies,KS_coefficients = np.linalg.eigh(F_OAO)

# Transform back to the AO basis:
C_occ = S_sqrt_inv @ KS_coefficients[:,:n_occ]

# Compute the density matrix
D_AO_new.np[:] = np.einsum('pi,qi->pq', C_occ, C_occ, optimize=True)

# Compute the norm of the change of density between the iterations.
normocc  = np.linalg.norm(D_AO.np - D_AO_new.np)

# Compute the Hxc potential contribution (old potential * new density !)
Hxcpot_energy = 2*np.einsum('pq,pq->', (J_coulomb.np + V_xc.np), D_AO_new.np) # factors 2 because D_AO is only alpha-D_AO.
Vpotential.set_D([D_AO_new])
Vpotential.compute_V([V_xc]) # otherwise it doesn't change the EHxc energy...
# Compute the Hxc energy with the new density:
Etot_new = 2*np.sum(KS_energies[:n_occ]) + EHxc - Hxcpot_energy + E_nuc
Delta_E  = abs(Etot_new - Etot)
Etot     = Etot_new

D_AO = D_AO_new

# Print results
print("*"*10 + " ITERATION {:3d} ".format(iteration) + "*"*10)
print("KS orb energies  : {}".format(KS_energies[:n_occ]))
print("Energy (hartree) : {:16.8f}".format(Etot))
print("Delta Energy     : {:16.8f}".format(Delta_E))
print("Norm Delta_occ   : {:16.8f}\n".format(normocc))

iteration += 1

# If maxiter is reached, stop the algorithm.
if iteration > SCF_maxiter: convergence=False

if (convergence):
print("*"*10 + " SUCCESS " + "*"*10)
print("Iteration        : {:16d}".format(iteration))
print("DFT energy       : {:16.8f}".format(Etot))
print("DFT energy exact : {:16.8f}".format(dft_e))
else:
print("*"*10 + " FAILURE " + "*"*10)
print("Iteration > {}, NO_CONVERGENCE_REACHED".format(SCF_maxiter))
``````

One thing I see is that you dont seem to recompute the Coulomb matrix.
The J integrals need to be contracted with the updated density.

``````# Initialize the JK object
jk = psi4.core.JK.build(dft_wfn.basisset())
jk.set_memory(int(1.25e8))  # 1GB
jk.initialize()
jk.set_do_K(False)
jk.set_do_J(True)
Cocc = psi4.core.Matrix(dft_wfn.nso(),dft_wfn.nalpha())
``````
``````# Compute JK
jk.compute()
jk.C_clear()
J_coulomb = jk.J()
``````

To get the initial Cocc from SAD psi4 does a Cholesky of the SAD density

``````from scipy.linalg import cholesky
``````

Alternatively you can request a `SADNO` guess where you should get a filled `SAD.Ca()`

Then, you can use the psi4 options

``````debug 1
print 5
``````

to print out Ca/Vpot/Da/Fa etc. matrices during the iterations and try to track deviations.

(code fragments in part taken from psi4numpy/Self-Consistent-Field at master · psi4/psi4numpy · GitHub)

Thanks ! I completely forgot about the Hartree term indeed…

With you code, apparently it should be

``````jk.C_left_add(SAD.Ca())
``````

instead of `Cocc`, otherwise I don’t get the correct Coulomb matrix when compared with the psi4 output obtained with the debug and print options.

I also noticed that I could use

``````J_coulomb = np.einsum('pqrs,rs->pq', I, D_AO.np)
``````

where

``````    mints = psi4.core.MintsHelper(dft_wfn.basisset())
I = np.asarray(mints.ao_eri())
``````

has been calculated before the SCF iterations.

With this new update of the Coulomb matrix, I couldn’t manage to converge at all at the begining, but my matrices were correct at the first iteration. I improved my code by just copy-pasting the DIIS algorithm (see RHF_DIIS.py in the link you shared), and it finally worked !

However, I still have a very slight difference with psi4. For instance, at R=2.5, I get:

``````E DFT         :      -3.46244834
E DFT psi4 :      -3.46353008
``````

Do you know what could cause this ? Actually, I saw that my Coulomb matrix was a bit different from psi4 (in the order of 10^{-3}), while I use the D_AO matrix constructed from the SAG guess to build it. So it might come from that, but I don’t know what this is different.

I enclose my complete code below:

``````import os,sys
import numpy as np
import math
sys.path.insert(1, os.path.abspath('/usr/local/psi4/lib/'))
import psi4

functional  = "SVWN"
basis       = "sto-3g"
E_conv      = 1e-6
D_conv      = 1e-6
SCF_maxiter = 100

n_orbs = 8
n_elec = n_orbs
n_occ  = n_elec//2
R      = 2.5

psi4.core.clean()
# Define the geometry
psi4.core.set_output_file("H{}_R{}_{}_{}_Psi4.dat".format(n_orbs,R,basis,functional),True)
geometry = "0 1\n"
for d in range(n_orbs//2):
geometry += "H 0. 0. {}\n".format(- (R/2. + d*R))
geometry += "H 0. 0. {}\n".format(+ (R/2. + d*R))
geometry += "symmetry c1\n"
geometry += "nocom\n"
geometry += "noreorient\n"

psi4.geometry(geometry)

psi4.set_options({'basis': basis,'save_jk':True, 'debug':1, 'print':5})
dft_e, dft_wfn = psi4.energy(functional, return_wfn=True)

Hcore = dft_wfn.H().clone()
E_nuc = dft_wfn.get_energies('Nuclear')
mints = psi4.core.MintsHelper(dft_wfn.basisset())
I = np.asarray(mints.ao_eri())

# Overlap matrix in the AO basis:
S_AO = dft_wfn.S().np
# KS-MO coefficient matrix in the AO basis:
C_AO = dft_wfn.Ca().np
# Compute the inverse square root of the overlap matrix S
S_eigval, S_eigvec = np.linalg.eigh(S_AO)
S_sqrt_inv = S_eigvec @ np.diag((S_eigval)**(-1./2.)) @ S_eigvec.T
C_transformation = np.linalg.inv(S_sqrt_inv)

# Construct the SAD Guess for the initial density D_AO
psi4.core.prepare_options_for_module("SCF")
psi4.core.get_global_option("BASIS"), puream=dft_wfn.basisset().has_puream(),
return_atomlist=True)
return_atomlist=True)

# Initialize the potential object
V_xc = dft_wfn.Da().clone()
Vpotential = dft_wfn.V_potential()

Etot        = 0
Delta_E     = 1e8
convergence = True
dRMS        = 1e8
Fock_list   = []
DIIS_error  = []

for SCF_ITER in range(1, SCF_maxiter+1):

# Compute the Coulomb potential
J_coulomb = np.einsum('pqrs,rs->pq', I, D_AO.np)
# Compute the XC potential with the new density matrix in the AO basis
Vpotential.set_D([D_AO])
Vpotential.compute_V([V_xc])

# Compute the Fock matrix in the AO basis:
F_AO = Hcore.np + 2*J_coulomb + V_xc.np

# DIIS
diis_e = np.einsum('ij,jk,kl->il', F_AO, D_AO.np, S_AO) - np.einsum('ij,jk,kl->il', S_AO, D_AO.np, F_AO)
diis_e = S_sqrt_inv @ diis_e @ S_sqrt_inv
Fock_list.append(F_AO)
DIIS_error.append(diis_e)
dRMS = np.mean(diis_e**2)**0.5

if SCF_ITER >= 2:

# Limit size of DIIS vector
diis_count = len(Fock_list)
if diis_count > 6:
# Remove oldest vector
del Fock_list
del DIIS_error
diis_count -= 1

# Build error matrix B, [Pulay:1980:393], Eqn. 6, LHS
B = np.empty((diis_count + 1, diis_count + 1))
B[-1, :] = -1
B[:, -1] = -1
B[-1, -1] = 0
for num1, e1 in enumerate(DIIS_error):
for num2, e2 in enumerate(DIIS_error):
if num2 > num1: continue
val = np.einsum('ij,ij->', e1, e2)
B[num1, num2] = val
B[num2, num1] = val

# normalize
B[:-1, :-1] /= np.abs(B[:-1, :-1]).max()

# Build residual vector, [Pulay:1980:393], Eqn. 6, RHS
resid = np.zeros(diis_count + 1)
resid[-1] = -1

# Solve Pulay equations, [Pulay:1980:393], Eqn. 6
ci = np.linalg.solve(B, resid)

# Calculate new fock matrix as linear
# combination of previous fock matrices
F_AO = np.zeros_like(F_AO)
for num, c in enumerate(ci[:-1]):
F_AO += c * Fock_list[num]

# Build the Fock matrix in the OAO basis:
F_OAO = S_sqrt_inv @ F_AO @ S_sqrt_inv.T
eigvals,eigvecs = np.linalg.eigh(F_OAO)

# Transform back to the AO basis:
C_occ = S_sqrt_inv @ eigvecs[:,:n_occ] # AO --> OAO = C_transformation = S^{-1/2}, OAO --> AO = S^{1/2}

# Compute the density matrix
D_AO.np[:] = np.einsum('pi,qi->pq', C_occ, C_occ, optimize=True)

# Compute the Hxc energy with the new density:
Hxcpot_energy = 2*np.einsum('pq,pq->', (J_coulomb + V_xc.np), D_AO.np) # factors 2 because D_AO is only alpha-D_AO.
Vpotential.set_D([D_AO])
Vpotential.compute_V([V_xc]) # otherwise it doesn't change the EHxc energy...
# Compute the Hxc potential contribution (old potential * new density !)
Etot_new = 2*np.sum(eigvals[:n_occ]) + EHxc - Hxcpot_energy + E_nuc
Delta_E  = abs(Etot_new - Etot)
Etot     = Etot_new

# Print results
print("*"*10 + " ITERATION {:3d} ".format(SCF_ITER) + "*"*10)
print("J matrix         : {}".format(J_coulomb))
#print("F matrix         : {}".format(F_AO))
#print("V(xc) matrix     : {}".format(V_xc.np))
#print("C_occ matrix     : {}".format(C_occ))
#print("D_AO matrix      : {}".format(D_AO.np))
#print("KS energies      : {}".format(eigvals))
#print("KS energy orbs   : {:16.8f}".format(2*np.sum(eigvals[:n_occ])))
#print("Hxcpot energy    : {:16.8f}".format(Hxcpot_energy))
print("Functional energy: {}".format(EHxc))
print("Energy (hartree) : {:16.8f}".format(Etot))
print("Delta Energy     : {:16.8f}".format(Delta_E))
print("dRMS             : {:16.8f}\n".format(dRMS))

if (Delta_E < E_conv) and (dRMS < D_conv):
print("*"*10 + " SUCCESS " + "*"*10)
print("R = ",R)
print("Iteration        : {:16d}".format(SCF_ITER))
print("DFT energy       : {:16.8f}".format(Etot))
print("DFT energy psi4  : {:16.8f}".format(dft_e))
break

if SCF_ITER == SCF_maxiter:
psi4.core.clean()
raise Exception("Maximum number of SCF cycles exceeded.")

``````

Oh sorry yes, my mistake. I didn’t really test it but it needs the full orbital matrix for sure.

The SCF convergence without any help from tricks like damping, shifting or DIIS is surprisingly bad yeah.

One remaining error might be the RI approximation in the “psi4 reference”. Need `scf_type pk` or `direct` to request conventional integrals.

1 Like

Perfect, I get the exact same energy now with numerical precision using `scf_type pk`!

1 Like

This topic was automatically closed 60 days after the last reply. New replies are no longer allowed.