# Task 1: Single-particle basis

We will here discuss how to set up the single-particle basis for the Hamiltonian used by {cite}`zanghellini_2004`. To achieve this we will solve the one-body Hamiltonian on a grid using finite differences. This will give us an approximation to the eigenstates and eigenenergies of the non-interacting single-particle problem.

## One-body harmonic oscillator Hamiltonian

In {cite}`zanghellini_2004` they look at two particles in a one-dimensional Hamiltonian coupled with a soft Coulomb interaction.
We will use a single-particle basis of one-dimensional harmonic oscillator eigenfunctions.
To construct this basis we diagonalize a one-dimensional harmonic oscillator Hamiltonian without the interaction term.
The Hamiltonian is then given by
\begin{align}
    \hat{h}
    = -\frac{1}{2} \frac{\text{d}^2}{\text{d}x^2}
    + \frac{1}{2} \omega^2 \hat{x}^2,
\end{align}
where we use natural units with $\hbar = m = 1$.
We proceed by solving the time-independent Schrödinger equation, viz.
\begin{align}
    \hat{h}|\chi_{\alpha}\rangle
    = \epsilon_{\alpha}|\chi_{\alpha}\rangle,
\end{align}
where $\epsilon_{\alpha}$ are the eigenenergies and $|\chi_{\alpha}\rangle$ the corresponding eigenstates.
We already know the answers to this equation from {cite}`griffiths2010introduction`.
The eigenenergies are $\epsilon_{\alpha} = \omega(\alpha + 1/2)$ for the energies and
\begin{align}
    \chi_{\alpha}(x) = \left(\frac{\omega}{\pi}\right)^{1 / 4}
    \frac{1}{\sqrt{2^{\alpha} \alpha!}} H_{\alpha}(\sqrt{\omega} x),
\end{align}
for the eigenstates projected onto $x$.
The $H_{\alpha}(x)$ are the physicists' Hermite polynomials.

## Setting up the basis using finite-differences

We will however be using the [quantum-systems library](https://github.com/Schoyen/quantum-systems) to set up our basis functions.
This library uses finite-differences to solve the one-body Hamiltonian as this allows for arbitrary potentials of the form $\hat{v}(x)$.


```{note}
Add theory on how the diagonalization is done.
```

We choose a basis of the $l$ first eigenstates.

(task_1_mel)=
### Matrix elements

The matrix elements we will be needing for this project are:
\begin{gather}
    h_{\alpha \beta} = \langle \chi_{\alpha} | \hat{h} | \chi_{\beta} \rangle
    = \epsilon_{\beta} \delta_{\alpha \beta}, \\
    u^{\alpha \beta}_{\gamma \delta}
    = \langle \chi_{\alpha} \chi_{\beta} | \hat{u} | \chi_{\gamma} \chi_{\delta} \rangle
    = \int \text{d}x_1 \text{d}x_2
    \frac{
        \chi^{*}_{\alpha}(x_1)
        \chi^{*}_{\beta}(x_2)
        \chi_{\gamma}(x_1)
        \chi_{\delta}(x_2)
    }{
        \sqrt{(x_1 - x_2)^2 + a^2}
    }, \\
    x_{\alpha \beta}
    = \langle \chi_{\alpha} | \hat{x} | \chi_{\beta} \rangle
    = \int \text{d} x
    \chi^{*}_{\alpha}(x)
    x
    \chi_{\beta}(x),
\end{gather}
where the position matrix elements $x_{\alpha \beta}$ will be used for the laser field.
We will also need a representation of the single-particle functions on a grid, i.e., $\chi_{\alpha}(x)$.
Luckily, quantum-systems will construct all of these matrix elements and single-particle functions for us.
Below we demonstrate how these states can be found from the library.


```python
%reload_ext nb_black
```


    <IPython.core.display.Javascript object>



```python
import numpy as np
import matplotlib.pyplot as plt

from quantum_systems import ODQD, GeneralOrbitalSystem
```


    <IPython.core.display.Javascript object>



```python
l = 10  # Number of eigenstates
grid_length = 10  # The width of the one-dimensional grid
num_grid_points = 2001  # The number of discretized points on the grid.
# More points give better results for the single-particle basis at the cost of slower setup.

alpha = 1  # The strength of the Coulomb interaction
a = 0.25  # The shielding parameter in the Coulomb interaction potential

omega = 0.25  # The frequency of the harmonic oscillator trap
```


    <IPython.core.display.Javascript object>



```python
potential = ODQD.HOPotential(omega=omega)  # Harmonic oscillator potential
```


    <IPython.core.display.Javascript object>


For examples on how to create your own potential see here: https://github.com/Schoyen/quantum-systems/blob/master/quantum_systems/quantum_dots/one_dim/one_dim_potentials.py

Below we construct the basis set.


```python
odqd = ODQD(
    l=l,
    grid_length=grid_length,
    num_grid_points=num_grid_points,
    alpha=alpha,
    a=a,
    potential=potential,
)
```


    <IPython.core.display.Javascript object>


The listed matrix elements can now be found from `odqd`.

1. The one-body Hamiltonian $h_{\alpha \beta}$ are found in `odqd.h`.
2. The two-body Coulomb interaction matrix elements $u^{\alpha \beta}_{\gamma \delta}$ are found in `odqd.u`.
3. The one-body position matrix elements $x_{\alpha \beta}$ are found in `odqd.position[0]`, where the index `0` is included in case of higher-dimensional systems. In our case there is only the one and `odqd.position.shape == (1, odqd.l, odqd.l)`.
4. The evaluated (evaluated at the grid) single-particle functions $\chi_{\alpha}(x)$ are found in `odqd.spf`, where `odqd.spf.shape == (odqd.l, num_grid_points)`. That is, the ground state harmonic oscillator single-particle function is found from `odqd.spf[0]`.
5. The discretized grid can be accessed from `odqd.grid`.

Note that all matrix elements and basis functions are inherently complex with the imaginary part being zero.
The ground-state solutions will be real, but when we include time-evolution everything will become complex and to avoid problems with data types when using NumPy we start of with complex matrix elements.

### Plotting the basis elements

Below we demonstrate how the single-particle functions can be plotted on the grid along with the harmonic oscillator potential.


```python
fig = plt.figure(figsize=(16, 10))

plt.plot(odqd.grid, potential(odqd.grid))

for i in range(l):
    plt.plot(
        odqd.grid,
        np.abs(odqd.spf[i]) ** 2 + odqd.h[i, i].real,
        label=r"$\chi_{" + f"{i}" + r"}$",
    )

plt.grid()
plt.legend()
plt.show()
```

## Types of orbital systems

In quantum chemistry it is common to call the initial single-particle functions for _atomic orbitals_ (AO).
We will also adhere to this convention, thus making the basis of harmonic oscillator single-particle functions the atomic orbitals of our system.
Now, as the harmonic oscillator Hamiltonian is independent of spin we have restricted our view to the spatial part of the atomic orbitals when using quantum-systems to construct our basis elements.
However, when moving to the full fermionic many-body picture (two = many) we have to include spin to get the full picture.
We often talk about three kinds of orbitals depending on the restrictions put on spin.

1. The general spin-orbital is on the form
    \begin{align}
        \phi(x, m_s)
        &= \varphi_1(x) \alpha(m_s) + \varphi_2(x) \beta(m_s),
    \end{align}
    where $\varphi_i(x)$ are the spatial oribtals, and $\sigma \in \{\alpha, \beta\}$ are spin basis functions with $m_s \in \{\uparrow, \downarrow\}$ the spin quantum number.
    Typical choices for the spin basis functions are
    \begin{gather}
        \alpha(\uparrow) = 1, \quad \beta(\uparrow) = 0, \\
        \alpha(\downarrow) = 0, \quad \beta(\downarrow) = 1.
    \end{gather}
2. The unrestricted spin-orbital is on the form
    \begin{align}
        \phi(x, m_s)
        = \varphi_{\sigma}(x) \sigma(m_s),
    \end{align}
    where an unrestricted spin-orbital has a definite spin-direction (either $\alpha$ or $\beta$), but allows for the spatial part to be different between spin-directions.
3. The restricted spin-orbital is on the form
    \begin{align}
        \phi(x, m_s)
        = \varphi(x) \sigma(m_s),
    \end{align}
    where the spin-orbitals have a definite spin-direction with the spatial part equal between the two spin-directions.

We will in this project create a general Hartree-Fock (GHF) solver which implies that we use general spin-orbitals.
This will turn out to be easiest mathematically (when deriving the equations for GHF) and the most general formulation.
However, the restricted and unrestricted formulations are more numerically efficient.

## Adding spin to the harmonic oscillator basis

We have constructed a spatial basis of harmonic oscillator eigenstates $\{\chi_{\alpha}\}_{\alpha = 1}^{l}$, with spin-independent matrix elements $h_{\alpha \beta}$, $u^{\alpha \beta}_{\gamma \delta}$, etc, as described [above](task_1_mel).
We now wish to include spin in the harmonic oscillator basis functions.
The way we do this is by doubling each spatial basis function (and hence each axis in the matrix elements) essentially creating duplicates.

We will denote this basis of general harmonic oscillator spin-orbitals by $\{\psi_{\mu}\}_{\mu = 1}^{2l}$.
The question is then how we can go from $h_{\alpha \beta} \to h_{\mu \nu}$ and similarly for the other matrix elements.
We let the general harmonic oscillator spin-orbitals be given by
\begin{align}
    |\psi_{\mu}\rangle
    &= |\chi_{\alpha}\rangle \otimes |\sigma\rangle,
\end{align}
where $\alpha = \mu // 2$ (integer division), and $\sigma = \alpha$ if $\mu$ is even and $\sigma = \beta$ if $\mu$ is odd.
The one-body Hamiltonian matrix elements are thus given by
\begin{align}
    h_{\mu \nu}
    &= \langle \psi_{\mu} | \hat{h} \otimes \hat{I} | \psi_{\nu} \rangle
    = \left(\langle \chi_{\alpha} | \otimes \langle \sigma | \right)
    \hat{h} \otimes \hat{I} \left( | \chi_{\beta} \rangle \otimes | \tau \rangle \right)
    \\
    &= \langle \chi_{\alpha} | \hat{h} | \chi_{\beta} \rangle
    \langle \sigma | \tau \rangle
    = h_{\alpha \beta} \delta_{\sigma \tau},
\end{align}
where $\hat{I}$ is the spin-basis identity operator, and we've chosen an orthonormal basis of spin states.
The full general spin-orbital matrix $\mathbf{h}_{GSO}$ can then be constructed from the (spatial) atomic orbital matrix $\mathbf{h}_{SpaS}$ and $\mathbf{I}_{2\times 2}$ the $2 \times 2$ identity matrix by
\begin{align}
    \mathbf{h}_{GSO} = \mathbf{h}_{SpaS} \otimes \mathbf{I}_{2 \times 2},
\end{align}
where $\otimes$ is the Kronecker product.
This ensures that all even (odd) indices correspond to $\alpha$ ($\beta$) spin-states.
From [quantum-systems](https://github.com/Schoyen/quantum-systems) this is done automatically when passing in a `BasisSet` to the class `GeneralOrbitalSystem`.
Below we demonstrate how this is done for harmonic oscillator basis set.


```python
print(f"Basis l: {odqd.l}")
print(f"Basis grid: {odqd.grid.shape}")
print(f"Basis h.shape: {odqd.h.shape}")
print(f"Basis u.shape: {odqd.u.shape}")
print(f"Basis x.shape: {odqd.position.shape}")
print(f"Basis spf.shape: {odqd.spf.shape}")
```

    Basis l: 10
    Basis grid: (2001,)
    Basis h.shape: (10, 10)
    Basis u.shape: (10, 10, 10, 10)
    Basis x.shape: (1, 10, 10)
    Basis spf.shape: (10, 2001)



    <IPython.core.display.Javascript object>



```python
n = 2  # Number of particles

# This adds spin to all matrix elements
# Note that gos.l = 2 * l
# Note also that odqd itself is transformed into a general orbital system.
# If you want to keep the pure spatial basis you can pass in a copy, viz.
# gos = GeneralOrbitalSystem(n, gos.copy_basis())
gos = GeneralOrbitalSystem(n, odqd)

print(f"GOS l: {gos.l}")
print(f"GOS grid: {gos._basis_set.grid.shape}")
print(f"GOS h.shape: {gos.h.shape}")
print(f"GOS u.shape: {gos.u.shape}")
print(f"GOS x.shape: {gos.position.shape}")
print(f"GOS spf.shape: {gos.spf.shape}")
```

    GOS l: 20
    GOS grid: (2001,)
    GOS h.shape: (20, 20)
    GOS u.shape: (20, 20, 20, 20)
    GOS x.shape: (1, 20, 20)
    GOS spf.shape: (20, 2001)



    <IPython.core.display.Javascript object>


Looking at the output above we see that the number of basis elements $l$ have been doubled.
However, the grid and the number of dimensions (the first axis in `gos.position`) are unchanged.
It is worth pointing out that `gos.spf[::2]` (all even indices) correspond to the spatial part of the $\alpha$ component of a general spin-orbital and `gos.spf[1::2]` (all odd indices) the $\beta$ component.

### Anti-symmetric two-body elements

When using general orbitals we often work with the anti-symmetric Coulomb elements.
These are defined to be
\begin{align}
    u^{\mu \nu}_{\kappa \lambda}
    = \langle \psi_{\mu} \psi_{\nu} | \hat{u} | \psi_{\kappa} \psi_{\lambda} \rangle_{AS}
    \equiv
    \langle \psi_{\mu} \psi_{\nu} | \hat{u} | \psi_{\kappa} \psi_{\lambda} \rangle
    -
    \langle \psi_{\mu} \psi_{\nu} | \hat{u} | \psi_{\lambda} \psi_{\kappa} \rangle,
\end{align}
where $u^{\mu \nu}_{\kappa \lambda}$ will be interpreted as the anti-symmetric two-body elements when using `GeneralOrbitalSystem` (this is the default behavior, however `gos = GeneralOrbitalSystem(n, odqd, anti_symmetrize=False)` will not make `gos.u` anti-symmetric).

b5fcfbd02cc72a8dbbbc7213985c2e70b9b6db0c
ipynb
Jupyter Notebook
docs/task-1-single-particle-basis.ipynb
Schoyen/tdhf-project-fys4411
b0231c0d759382c14257cc4572698aa80c1c94d0
MIT
2021-06-03T00:34:57.000Z
2021-06-03T00:34:57.000Z
docs/task-1-single-particle-basis.ipynb
Schoyen/tdhf-project-fys4411
b0231c0d759382c14257cc4572698aa80c1c94d0
MIT
docs/task-1-single-particle-basis.ipynb
Schoyen/tdhf-project-fys4411
b0231c0d759382c14257cc4572698aa80c1c94d0
MIT
Qwen/Qwen-72B
1. YES
2. YES