Show code cell content
"""Beryllium atom"""
__author__ = "Magnar Bjørgve"
__credit__ = ["Magnar Bjørgve", "Stig Rune Jensen"]
__date__ = "2021-02-16"
Beryllium atom#
We solve the Hartree-Fock equations for a closed-shell many-electron system
\[\begin{equation*}
\left[\hat{T} + \hat{V}\right]\Phi = \Phi F
\end{equation*}\]
In the equation above \(\Phi\) is a row vector of doubly occupied orbitals, \(F = \langle\Phi|\hat{T}+\hat{V}|\Phi\rangle\) is the Fock matrix, and \(\hat{V} = \hat{V}_{nuc} + 2\hat{J} - \hat{K}\) is the potential operator. We here define the nuclear attraction potential
\[\begin{align*}
\hat{V}_{nuc}(r) = \sum_I \frac{-Z_I}{|R_I - r|}
\end{align*}\]
with \(Z_I\) and \(R_I\) being the nuclear charges and positions, and the Coulomb and Exchange operators
\[\begin{align*}
\hat{J} \varphi_i &= \sum_j \hat{P}[\varphi_j \varphi_j^\dagger] \varphi_i \\
\hat{K} \varphi_i &= \sum_j \hat{P}[\varphi_i \varphi_j^\dagger] \varphi_j
\end{align*}\]
where \(\hat{P}\) is the convolution operator with the Poisson kernel \(P=1/r\). We iterate the Hartree-Fock equations in their integral form
\[\begin{equation*}
\tilde{\Phi}^{n+1} = -2 \hat{G}^{\vec{\mu}^n}\left[\hat{V} \Phi^n + \Phi^n\left(\Lambda^n - F^n\right)\right]
\end{equation*}\]
where \(\hat{G}^{\vec{\mu}}\) is a vector of convolution operators with the bound-state Helmholtz kernel \(G^{\mu} = e^{-\mu r}/r\) and \(\mu_i^n = \sqrt{-2\lambda_i^n}\) is a positive real number. \(\Lambda_{ij}^n = \lambda_i^n\delta_{ij}\) is a diagonal matrix containing the same parameters as used in the construction of the Helmholz vector. In practice we choose \(\lambda_i^n = F_{ii}^n\) such that diagonal elements vanish in \(\left(\Lambda^n - F^n\right)\).
from vampyr import vampyr3d as vp
import numpy as np
import numpy.linalg as LA
import matplotlib.pyplot as plt
class HelmholtzOperator():
"""
Vectorized Helmholtz operator
Parameters
----------
mra : The multiresolution analysis we work on
lamb : vector of lambda parameters, mu_i = sqrt(-2*lambda_i)
prec : Precision requirement
Attributes
----------
operators : list containing HelmholtzOperators for each orbital
"""
def __init__(self, mra, lamb, prec):
self.mra = mra
self.lamb = lamb
self.prec = prec
self.operators = []
self.setup()
def setup(self):
mu = [np.sqrt(-2.0*l) if l < 0 else 1.0 for l in self.lamb]
for m in mu:
self.operators.append(vp.HelmholtzOperator(mra=self.mra, exp=m, prec=self.prec))
def __call__(self, Psi):
"""Operate the Helmholtz operator onto an orbital vector"""
return np.array([self.operators[i](Psi[i]) for i in range(len(Psi))])
class ExchangeOperator():
"""
Vectorized Exchange operator
K(phi) = \sum_i P[phi * psi_i] * psi_i
Parameters
----------
mra : The multiresolution analysis we work on
Psi : Orbital vector defining the operator
prec : Precision requirement
Attributes
----------
poisson : Poisson Operator
"""
def __init__(self, mra, Psi, prec):
self.mra = mra
self.Psi = Psi
self.prec = prec
self.poisson = vp.PoissonOperator(mra=mra, prec=self.prec)
def __call__(self, Phi):
"""Apply the exchange operator onto an orbital vector Phi"""
Phi_out = []
for j in range(len(Phi)):
V_j0 = self.poisson(Phi[j] * self.Psi[0])
tmp = (self.Psi[0] * V_j0).crop(self.prec)
for i in range(1, len(self.Psi)):
V_ji = self.poisson(Phi[j] * self.Psi[i])
tmp += (self.Psi[i] * V_ji).crop(self.prec)
tmp *= 4.0*np.pi
Phi_out.append(tmp)
return np.array([phi.crop(self.prec) for phi in Phi_out])
class KineticEnergyOperator():
def __init__(self, mra, prec):
self.prec = prec
self.mra = mra
self.derivative = vp.ABGVDerivative(mra, 0.5, 0.5)
def laplace_operator(self, phi):
return self.derivative(self.derivative(phi, 0), 0) + self.derivative(self.derivative(phi, 1), 1) + self.derivative(self.derivative(phi, 2), 2)
def __call__(self, Phi):
return np.array([(-0.5*self.laplace_operator(phi)).crop(self.prec) for phi in Phi])
class CouloumbOperator():
"""
Vectorized Couloumb operator
J(phi) = \sum_i P[psi_i * psi_i] * phi
Parameters
----------
mra : The multiresolution analysis we work on
Psi : Orbital vector defining the operator
prec : Precision requirement
Attributes
----------
poisson : Poisson operator
potential : Coulomb potential
"""
def __init__(self, mra, Psi, prec):
self.mra = mra
self.Psi = Psi
self.prec = prec
self.poisson = vp.PoissonOperator(mra=mra, prec=self.prec)
self.potential = None
self.setup()
def setup(self):
rho = self.Psi[0]**2
for i in range(1, len(self.Psi)):
rho += self.Psi[i]**2
rho.crop(self.prec)
self.potential = (4.0*np.pi)*self.poisson(rho).crop(self.prec)
def __call__(self, Phi):
"""Apply Couloumb operator onto an orbital vector Phi"""
return np.array([(self.potential*phi).crop(self.prec) for phi in Phi])
class NuclearOperator():
"""
Vectorized Nuclear potential operator
Parameters
----------
mra : The multiresolution analysis we work on
atoms : List of dicts containing charge and coordinates of the atoms
prec : Precision requirement
Attributes
----------
potential : Nuclear potential
"""
def __init__(self, mra, atoms, prec):
self.mra = mra
self.prec= prec
self.atoms = atoms
self.potential = None
self.setup()
def setup(self):
# Project analytic function onto MRA basis
P_mra = vp.ScalingProjector(mra=self.mra, prec=self.prec)
f_nuc = NuclearFunction(self.atoms)
self.potential = P_mra(f_nuc)
def __call__(self, Phi):
"""Apply nuclear potential operator onto an orbital vector"""
return np.array([(self.potential*phi).crop(self.prec) for phi in Phi])
class NuclearFunction():
"""
Vectorized Nuclear potential operator
Parameters
----------
atoms : List of dicts containing charge and coordinates of the atoms
"""
def __init__(self, atoms):
self.atoms = atoms
def __call__(self, r):
"Returns the nuclear potential value in R"
tmp = 0.0
for atom in self.atoms:
Z = atom["Z"]
R = atom["R"]
R2 = (R[0]-r[0])**2 + (R[1]-r[1])**2 + (R[2]-r[2])**2
tmp += -Z / np.sqrt(R2)
return tmp
def scf_solver(mra, atoms, Phi_n, F_n, precision, threshold, max_iter=30):
"""Kinetric-free Hartree-Fock SCF solver
Parameters:
mra : The Multiresolution analysis to work on
atoms : List of dicts containing charge and coordinates of the atoms
Phi_n : Starting guess orbitals
F_n : Starting guess for Fock matrix
precision : Precision requirement
threshold : Usually set to the same as precision, set to -1 to limit iterations by max_iter
max_iter : Set maximum iterations
Returns:
Updates : Vector of orbital residual norms at each iteration
Energies : List of energy contributions at each iteration
Phi_n : Converged orbital vector
"""
# Setup nuclear potential
V_nuc = NuclearOperator(mra, atoms, precision)
# Loop parameters
iteration = 0 # Iteration counter
update = np.ones(len(Phi_n)) # Initialize error measure (norm of orbital updates)
updates = [] # Will capture wavefunction updates for visusualization
energies = [] # Will capture energies for visualization
# SCF loop
while (max(update) > threshold):
if iteration > max_iter-1:
break
# Initialize operators for first iteration
J_n = CouloumbOperator(mra, Phi_n, precision)
K_n = ExchangeOperator(mra, Phi_n, precision)
# Initialize vector of Helmholtz operators based on Fock matrix diagonal
Lambda_n = np.diag(np.diag(F_n))
G = HelmholtzOperator(mra, np.diag(Lambda_n), precision)
# Apply potential operator to all orbitals
VPhi = V_nuc(Phi_n) + 2*J_n(Phi_n) - K_n(Phi_n)
KinOpr = KineticEnergyOperator(mra, precision)
# Apply Helmholtz operators to all orbitals
Phi_np1 = -2*G(VPhi + (Lambda_n - F_n) @ Phi_n)
dPhi_n = Phi_np1 - Phi_n
update = np.array([phi.norm() for phi in dPhi_n])
# Compute overlap matrices
S_tilde = calc_overlap(Phi_np1, Phi_np1)
# Löwdin orthonormalization S^{-1/2} = U * Sigma^{-1/2} * U^T
sigma, U = LA.eig(S_tilde)
Sm5 = U @ np.diag(sigma**(-0.5)) @ U.transpose()
Phi_bar = Sm5 @ Phi_np1
# Initialize n+1 operators
J_np1 = CouloumbOperator(mra, Phi_bar, precision)
K_np1 = ExchangeOperator(mra, Phi_bar, precision)
# Prepare for next iteration
Phi_n = Phi_bar
VPhikin = V_nuc(Phi_n) + 2*J_np1(Phi_n) - K_np1(Phi_n) + KinOpr(Phi_n)
F_n = calc_overlap(Phi_n, VPhikin)
# Compute energy contributions
energy = calc_energies(F_n, Phi_n, V_nuc, J_np1, K_np1)
# Collect output
updates.append(update)
energies.append(energy)
print(iteration, " | E_tot:", energy["$E_{tot}$"], " | dPhi:", max(update))
iteration += 1
return np.array(updates), energies, Phi_n
def calc_energies(F_mat, Phi, V, J, K):
""""Calcuate all energy contributions"""
V_mat = calc_overlap(Phi, V(Phi))
J_mat = calc_overlap(Phi, J(Phi))
K_mat = calc_overlap(Phi, K(Phi))
E_orb = 2.0*F_mat.trace()
E_en = 2.0*V_mat.trace()
E_coul = 2.0*J_mat.trace()
E_ex = -K_mat.trace()
E_tot = E_orb - E_coul - E_ex
E_kin = E_tot - E_en - E_coul - E_ex
return {
"$E_{orb}$": E_orb,
"$E_{en}$": E_en,
"$E_{coul}$": E_coul,
"$E_{ex}$": E_ex,
"$E_{kin}$": E_kin,
"$E_{tot}$": E_tot
}
def calc_overlap(Bra, Ket):
"""Calculate the overlap matrix between the orbitals <Bra| and |Ket>
Parameters:
Bra : bra vector <Phi|
Ket : ket vector |Phi>
Returns:
Overlap matrix
"""
S = np.empty((len(Bra), len(Ket)))
for i in range(len(Bra)):
for j in range(len(Ket)):
S[i, j] = vp.dot(Bra[i], Ket[j])
return S
def starting_guess(mra, atom, n_orbs, prec):
"""Primitive starting guess, works for Be"""
# Define projector onto the MRA basis
P_mra = vp.ScalingProjector(mra=mra, prec=prec)
Phi = []
for i in range(1, n_orbs+1):
R0 = atom["R"]
def f_gauss(r):
R2 = (r[0]-R0[0])**2 + (r[1]-R0[1])**2 + (r[2]-R0[2])**2
return np.exp(-R2/i)
phi = P_mra(f_gauss)
phi.normalize()
Phi.append(phi)
Phi = np.array(Phi)
# Löwdin orthonormalization S^{-1/2} = U * Sigma^{-1/2} * U^T
S = calc_overlap(Phi, Phi)
sigma, U = LA.eig(S)
Sm5 = U @ np.diag(sigma**(-0.5)) @ U.transpose()
return Sm5 @ Phi
# Define molecule
mol = [
{"Z": 4.0, "R": [0.0, 0.0, 0.0]}
]
# Set target precision
precision = 1.0e-3
# Define MRA basis and computational domain
k = 5 # Polynomial order of MRA basis
world = [-20, 20] # Computational domain [-L,L]^3 (a.u.)
MRA = vp.MultiResolutionAnalysis(order=k, box=world)
# Setup guess for the orbitals and the Fock matrix
Phi_n = starting_guess(mra=MRA, atom=mol[0], n_orbs=2, prec=precision)
F_n = np.array([
[-1.0, 0.0],
[ 0.0,-0.5]
])
# Run SCF solver to optimize the orbitals
orbital_updates, energies, Phi_n = scf_solver(MRA, mol, Phi_n, F_n, precision, precision, max_iter=30)
0 | E_tot: -12.158806160497663 | dPhi: 0.9554261103233852
1 | E_tot: -13.910545962062471 | dPhi: 0.6030032083595198
2 | E_tot: -14.425154911621195 | dPhi: 0.11248971009686566
3 | E_tot: -14.550559502300937 | dPhi: 0.056367744794410625
4 | E_tot: -14.573091352701868 | dPhi: 0.024464066721782796
5 | E_tot: -14.576656010668405 | dPhi: 0.009607118621517596
6 | E_tot: -14.577184437878596 | dPhi: 0.0036909439005734903
7 | E_tot: -14.577255599761454 | dPhi: 0.0017015854922230098
8 | E_tot: -14.577262672068981 | dPhi: 0.0009093846224262461