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()[0]
# 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")
sad_basis_list = psi4.core.BasisSet.build(dft_wfn.molecule(), "ORBITAL",
psi4.core.get_global_option("BASIS"), puream=dft_wfn.basisset().has_puream(),
return_atomlist=True)
sad_fitting_list = psi4.core.BasisSet.build(dft_wfn.molecule(), "DF_BASIS_SAD",
psi4.core.get_option("SCF", "DF_BASIS_SAD"), puream=dft_wfn.basisset().has_puream(),
return_atomlist=True)
# Use Psi4 SADGuess object to build the SAD Guess
SAD = psi4.core.SADGuess.build_SAD(dft_wfn.basisset(), sad_basis_list)
SAD.set_atomic_fit_bases(sad_fitting_list)
SAD.compute_guess();
D_AO = SAD.Da()
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:
EHxc = Vpotential.quadrature_values()["FUNCTIONAL"]
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))