V2: SCF optimization with VAMPyR#

Hydrogen atom#

We intend to solve the Schrödinger equation for the Hydrogen atom

\[\begin{equation*} \left[\hat{T} + \hat{V}\right]\phi = E\phi \end{equation*}\]

where the potential operator contains only the nuclear attraction

\[\begin{equation*} \hat{V}(r) = \frac{-Z}{|R-r|} \end{equation*}\]

where \(Z=1\) and \(R=(0,0,0)\) are the nuclear charge and position, respectively.

Solution using multiwavelets#

In order to solve this equation using multiwavelets we reformulate it in an integral form:

(1)#\[\begin{equation} \phi = -2\hat{G}_{\mu}\hat{V}\phi \end{equation}\]

and iterate until convergence. \(\hat{G}_\mu\) is the Bound-State Helmholtz (BSH) integral operator, whose kernel is defined as

(2)#\[\begin{equation} G_\mu(r - r') = \frac{\exp(-\mu |r - r'|)}{4\pi |r - r'|} \end{equation}\]

and \(\mu=\sqrt{-2E}\) takes the place of the energy in the original Schrödinger equation. For the Hydrogen atom this energy is known (E=-0.5 a.u.) so we can use the “exact” BSH operator when solving for the wavefunction, but in general the energy will have to be computed for each iteration and the BSH operator re-initialized with the updated \(\mu\) parameter. For simple systems with a handful of electrons, the straightforward power iteration of the BSH equation is sufficient to achieve convergence, but more complicated molecular systems require additional acceleration techniques (not discussed further).

Implementation in VAMPyR#

What we need for this exercise is:

  1. Define a suitable Multi-Resolution Analysis (MRA) for the problem

  2. Make analytic nuclear potential and project onto the MRA

  3. Make analytic initial guess for the orbital and project onto the MRA

  4. Create a Helmholtz operator \(\hat{G}\) using the exact value for the Hydrogen atom (\(\mu=\sqrt{-2E}\), \(E=-0.5 a.u.\))

  5. Compute new orbital through application of the Helmholtz operator \(\tilde{\phi}^{n+1} = -2\hat{G}\left[\hat{V}\phi^n\right]\)

  6. Check for convergence by computing the norm of the orbital update \(\Delta\phi^n = \tilde{\phi}^{n+1} - \phi^n\)

  7. Normalize the orbital \(\phi^{n+1} = \tilde{\phi}^{n+1}/||\tilde{\phi}^{n+1}||\)

  8. Update orbital \(\phi^{n} \leftarrow \phi^{n+1}\) for next iteration

  9. Repeat steps 5-8 until your orbital has converged

Preparing the notebook#

First we load all necessary Python modules for this exercise, and define a helper function for plotting that will be used later.

from vampyr import vampyr3d as vp
import numpy as np
import matplotlib.pyplot as plt

# Prepare simple plotting function for visualizations
def line_plot(f, x_min, x_max):
    # Evenly spaced points along the x-axis
    r_x = np.linspace(x_min, x_max, 1000)
    data = [f([x, 0.0, 0.0]) for x in r_x]
    plt.plot(r_x, data)
    return plt

Initializing the Multi-Resolution Analysis (MRA)#

We need to define a computational world that is large enough for the expected solution to vanish at the boundary. For the Hydrogen atom in atomic units, a cubic box of \(\pm 20\) bohrs in each dimension should be sufficient. Initially we choose a low precision of \(\epsilon = 10^{-3}\), but this can be increased later when we know that the algorithm is working. For this precision a suitable polynomial order is around \(k=5\), but this should be increased gradually when the precision is increased.

# Global parameters
k = 5                        # Polynomial order
L = [-20,20]                 # Simulation box size
epsilon = 1.0e-3             # Relative precision

# Define MRA and multiwavelet projector
print("FIXME: construct a MRA using polynomial order 'k' and box size 'L'")
# MRA = ...
print("FIXME: construct a ScalingProjector 'P_mra' for the MRA with precision 'epsilon'")
# P_mra = ...

Nuclear potential#

By placing the nucleus at the origin, the nuclear potential for Hydrogen (Z=1) reduces to \(\hat{V}(r) = -1/|r|\). In order to make a numerical MW representation of this potential we must first define an analytic function that takes a position coordinate \(r[3]\) and returns a scalar function value for that point. Once we have a Python function for this we can project it onto the MRA using the scaling projector. We can then plot the function to make sure it looks reasonable.

# Analytic nuclear potential (the argument r is a 3D coordinate)
def f_nuc(r):
    Z = 1.0
    R = np.sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2])
    return -Z / R
    
# Project nuclear potential onto MRA
print("FIXME: construct a MW representation 'V_nuc' of the analytic function 'f_nuc' using the 'P_mra' projector")
# V_nuc = ...

# Uncomment the following to plot function between x=[-1.0, 1.0]
# plt = line_plot(V_nuc, -1.0, 1.0)
# plt.show()

Initial guess for wavefunction#

In order to start the iteration process, we need an initial guess for the wavefunction solution. We choose a simple Gaussian function \(\phi^0(r) = e^{-r^2}\) centered at the origin, and project it in the same way as the nuclear potential above.

# Analytic initial guess (the argument r is a 3D coordinate)
def f_phi(r):
    print("FIXME: return the function value of the initial guess, e.g. exp(-r^2)")
    
# Project initial guess onto MRA and normalize
print("FIXME: construct a MW representation 'phi_n' of the analytic function 'f_phi' using the 'P_mra' projector")
print("FIXME: normalize the resulting orbital using 'phi_n.normalize()'")
# phi_n = ...


# Uncomment the following to plot function between x=[-10.0, 10.0]
# plt = line_plot(phi_n, -10.0, 10.0)
# plt.show()

Bound-State Helmholtz operator#

The BSH operator is implemented in VAMPyR, and you can construct it using any real-valued \(\mu\), along with the appropriate MRA and precision \(\epsilon\). Since we know the exact energy for the Hydrogen atom, we can construct a single operator and use it throughout the iteration procedure.

# Prepare Helmholtz operator using exact energy
E = -0.5
mu = np.sqrt(-2*E)
G = vp.HelmholtzOperator(MRA, exp=mu, prec=epsilon)

Iterative solution of the Schrödinger equation#

We now have all the componenents needed for the power iteration that will optimize the wavefunction solution. In the loop we need to

  1. Assemble Helmholtz argument : \(V_{nuc}*\phi^n\)

  2. Apply BSH operator : \(\tilde{\phi}^{n+1} = -2 \hat{G}\left[V_{nuc}*\phi^n\right]\)

  3. Compute norm of new orbital : \(||\tilde{\phi}^{n+1}||\)

  4. Compute norm of update : \(||\Delta\tilde{\phi}^n|| = ||\tilde{\phi}^{n+1} - \phi^n||\)

  5. Normalize wavefunction : \(\phi^{n+1} = \tilde{\phi}^{n+1}/||\tilde{\phi}^{n+1}||\)

  6. Prepare for next iteration : \(\phi^n \leftarrow \phi^{n+1}\)

  7. Repeat until : \(||\Delta\tilde{\phi}^n|| < \epsilon_{thrs}\)

i = 0            # Iteration counter
thrs = 1.0e-3    # Convergence threshold
norm = 1.0       # Dummy norm
update = 1.0     # Dummy update
while (i < 1):
    # Main optimization loop
    print("FIXME: assemble Helmholtz argument")
    print("FIXME: apply BSH operator to obtain next iteration 'phi_np1'")
    print("FIXME: compute the norm of new orbital : norm = ||phi_np1||")
    print("FIXME: compute the norm of update : update = ||phi_np1 - phi_n||")
    print("FIXME: use the 'line_plot' function to visualize the orbitals")
    print("FIXME: switch to (update < thrs) convergence check once the body of the loop is working")
    
    # Print the current convergence status
    print("iteration: {}    Norm: {}   Update: {}".format(i, norm, update))
    i += 1

Extension to Helium#

A few things change when you go from Hydrogen to Helium:

  1. The energy is no longer known exactly, and thus will have to be computed from the wavefunction

  2. The Helmholtz operator which depends on the energy through \(\mu = \sqrt{-2E}\) needs to be updated in every iteration

  3. The potential operator \(V\) depends on the wavefunction and must be updated in every iteration

In this example we will use the Hartree-Fock model, which for a single-orbital system like Helium, reduces to the following potential operator:

(3)#\[\begin{align} \hat{V} &= \hat{V}_{nuc} + 2\hat{J} - \hat{K}\\ &= \hat{V}_{nuc} + \hat{J} \end{align}\]

since \(\hat{K} = \hat{J}\) for a doubly occupied single orbital.

The Coulomb potential \(\hat{J}\) can be computed by application of the Poisson operator \(P\):

(4)#\[\begin{equation} \hat{J}(r) = P\left[4\pi\rho\right] \end{equation}\]

Where \(\rho\) is the square of the orbital

(5)#\[\begin{equation} \rho = \phi*\phi \end{equation}\]

Pen and paper exercise:#

Use the fact that

(6)#\[\begin{equation} \tilde{\phi}^{n+1} = -\Big[\hat{T} - E^n\Big]^{-1} V^n\phi^n \end{equation}\]

to show that

(7)#\[\begin{equation} E^{n+1} = \frac{\langle\tilde{\phi}^{n+1}|\hat{T} + \hat{V}^{n+1}|\tilde{\phi}^{n+1}\rangle} {||\tilde{\phi}^{n+1}||^2} \end{equation}\]

can be written as a pure update \(\Delta E^n\) involving only the potentials \(\hat{V}^{n+1}\), \(\hat{V}^n\) as well as the orbitals \(\tilde{\phi}^{n+1}\) and \(\phi^n\)

(8)#\[\begin{equation} E^{n+1} = E^{n} + \Delta E^n \end{equation}\]

Implementation exercise:#

  1. Make a nuclear potential function f_nuc(r) for the Helium atom

  2. Make an initial guess for the orbital as a python function f_phi(r)

  3. Project both nuclear potential (\(V_{nuc}\)) and orbital (\(\phi_n\)) to the MW basis using a vp.ScalingProjector with precision \(\epsilon=10^{-3}\)

  4. Create a Helmholtz operator \(G^n\) with \(\mu^n = \sqrt{-2E^n}\) using the current energy (start with some reasonable guess)

  5. Compute total potential \(V^n = V_{nuc} + J^n\), where the Coulomb potential is computed using the vp.PoissonOperator on the current squared orbital \(J = P\left[4\pi\phi^2\right]\)

  6. Compute new orbital through application of the Helmholtz operator on \(\tilde{\phi}^{n+1} = -2\hat{G}^nV^n\phi^n\)

  7. Compute the update to the orbital energy (\(\Delta E^n\)) using the expression from the pen and paper exercise

  8. Compute the norm of the orbital update \(||\Delta\tilde{\phi}^n|| = ||\tilde{\phi}^{n+1} - \phi^n||\)

  9. Normalize the orbital \(\phi^{n+1} = \tilde{\phi}^{n+1}/||\tilde{\phi}^{n+1}||\)

  10. Repeat steps 4-9 until your orbital has converged

You should expect the orbital energy to converge towards \(E_n \approx -0.918\).

Bonus exercise:#

The total energy can be computed after convergence as \(E_{tot} = 2E_n - \langle\rho|J\rangle\), should be around \(E_{tot} \approx -2.86\).