The aim of this project is to use the coupled cluster method method to evaluate the ground state energy of quantum dots with \( N=2 \), \( N=6 \), \( N=12 \) and \( N=20 \) electrons. These are so-called closed shell systems. These systems can however be changed to other ones. This is a project which favors people who have studied many-body physics at the level of FYS4480/9480.
We consider a system of electrons confined in a pure two-dimensional isotropic harmonic oscillator potential, with an idealized total Hamiltonian given by
$$ \begin{equation} \label{eq:finalH} \hat{H}=\sum_{i=1}^{N} \left( -\frac{1}{2} \nabla_i^2 + \frac{1}{2} \omega^2r_i^2 \right)+\sum_{i < j}^{N}\frac{1}{r_{ij}}, \end{equation} $$where natural units (\( \hbar=c=e=m_e=1 \)) are used and all energies are in so-called atomic units a.u. We will study systems of many electrons \( N \) as functions of the oscillator frequency \( \omega \) using the above Hamiltonian. The Hamiltonian includes a standard harmonic oscillator part
$$ \begin{equation*} \hat{H}_0=\sum_{i=1}^{N} \left( -\frac{1}{2} \nabla_i^2 + \frac{1}{2} \omega^2r_i^2 \right), \end{equation*} $$and the repulsive interaction between two electrons given by
$$ \begin{equation*} \hat{H}_I=\sum_{i < j}^{N}\frac{1}{r_{ij}}, \end{equation*} $$with the distance between electrons given by \( r_{ij}=\sqrt{\boldsymbol{r}_1-\boldsymbol{r}_2} \). We define the modulus of the positions of the electrons (for a given electron \( i \)) as \( r_i = \sqrt{r_{i_x}^2+r_{i_y}^2} \). There is no spin-orbit part in the two-body Hamiltonian and the different parts of the Hamiltonian have no spin operators. This means that the spin of a single-particle state or a many-particle state is a conserved quantity.
If only the harmonic oscillator part of the Hamiltonian, the so-called unperturbed part is given by
$$ \begin{equation*} \hat{H}_0=\sum_{i=1}^{N} \left( -\frac{1}{2} \nabla_i^2 + \frac{1}{2} \omega^2r_i^2 \right). \end{equation*} $$The wave function for one electron in an oscillator potential in two dimensions is
$$ \begin{equation*} \phi_{n_x,n_y}(x,y) = A H_{n_x}(\sqrt{\omega}x)H_{n_y}(\sqrt{\omega}y)\exp{(-\omega(x^2+y^2)/2}. \end{equation*} $$The functions \( H_{n_x}(\sqrt{\omega}x) \) are so-called Hermite polynomials, discussed in the appendix while \( A \) is a normalization constant. For the lowest-lying state we have \( n_x=n_y=0 \) and an energy \( \epsilon_{n_x,n_y}=\omega(n_x+n_y+1) = \omega \).
In our case the single-particle basis is provided by the harmonic oscillator functions. The first part of this project deals thus with the setup of the basis functions and the computation of matrix elements for the Coulomb interaction.
Make a program that sets up all the quantum numbers for the single-particle basis using a harmonic oscillator in two dimensions. The list should contain the quantum numbers \( n_x \) and \( n_y \), spin and its projection and the single particle energies in atomic units. Define a cutoff in the list according to the harmonic oscillator excitation energy. The table here lists the first four harmonic oscillator shells, with its pertinent degeneracies, and the total number of electrons which can be accomodated up to a given energy. These numbers are the so-called magic numbers for a two-dimensional quantum dot. The degeneracies take into account the two spin values an electron can take. Convince yourself about the correctness of this table and use it to check that your code is running correctly.
| Shell number | \( (n_x, n_y) \) | Energy | Degeneracy | \( N \) |
| 1 | \( (0,0) \) | \( \hbar\omega \) | 2 | 2 |
| 2 | \( (1,0) \), \( (0,1) \) | \( 2\hbar\omega \) | 4 | 6 |
| 3 | \( (2,0) \), \( (0,2) \), \( (1,1) \) | \( 3\hbar\omega \) | 6 | 12 |
| 3 | \( (3,0) \), \( (0,3) \), \( (2,1) \), \( (1,2) \) | \( 4\hbar\omega \) | 8 | 20 |
Convince yourself that the unperturbed lowest-lying energy for the two-electron system is simply \( 2\omega \). Similarly, find the corresponding energies for \( N=6 \), \( N=12 \), \( N=20 \) and \( N=30 \) electrons. These results will turn out to be very useful when you need to check your programs.
With the single-particle basis from the previous step, you should now write a class which encodes the information about the harmonic oscillator single-particle basis. This should include the single-particle energies, spin and its projections as well as the harmonic oscillator functions that depend on \( n_x \) and \( n_y \).
In order to set up the two-body matrix elements, we need to define the so-called direct and exchange matrix elements. This leads to what we call an anti-symmetrized matrix element.
We introduce the following shorthands for the integrals
$$ \langle pq|\hat{v}|rs\rangle = \int \psi_{p}^*(x_i)\psi_{q}^*(x_j)V(r_{ij})\psi_{r}(x_i)\psi_{s}(x_j) dx_idx_j, $$which defines the so-called direct matrix element and
$$ \langle pq|\hat{v}|sr\rangle = \int \psi_{p}^*(x_i)\psi_{q}^*(x_j) V(r_{ij})\psi_{s}(x_i)\psi_{r}(x_j) dx_idx_j, $$which defines the exchange element. The variables \( pqrs \) define all the single-particle quantum numbers. In our case these are \( n_x \), \( n_y \) and spin and its projections.
These matrix elements are defined in terms of two-body quantum numbers. Which quantum numbers are conserved? Which are the possible values the total spin of a two-body state can have?
The direct and exchange matrix elements can be brought together if we define the anti-symmetrized matrix element
$$ \langle pq|\hat{v}|rs\rangle_{\mathrm{AS}}= \langle pq|\hat{v}|rs\rangle-\langle pq|\hat{v}|sr\rangle. $$It has the symmetry property
$$ \langle pq|\hat{v}|rs\rangle_{\mathrm{AS}}= -\langle pq|\hat{v}|sr\rangle_{\mathrm{AS}}=-\langle qp|\hat{v}|rs\rangle_{\mathrm{AS}}, $$and
$$ \langle pq|\hat{v}|rs\rangle_{\mathrm{AS}}= \langle qp|\hat{v}|sr\rangle_{\mathrm{AS}}. $$Your task is to write a function which calculates the integral
$$ \langle pq|\hat{v}|rs\rangle = \int \psi_{p}^*(x_i)\psi_{q}^*(x_j)V(r_{ij})\psi_{r}(x_i)\psi_{s}(x_j) dx_idx_j, $$and then assemble the direct and exchange terms in order to construct the anti-symmetrized matrix elements to be used as inputs to the coupled cluster code. The integral needs to be parallelized. Chapter five of the lecture notes discusses Gaussian quadrature. You may find it convenient to use Hermite polynomials in order to set up the integration points and weights.
The aim of this project is to develop a coupled cluster doubles (CCD) code, where \( 2p-2h \) excitations are included only.
We will start with a two-electron problem and compare our results to those of Taut, see reference [1] below.
The ansatz for the ground state is given by
$$ \begin{equation*} \vert \Psi_0\rangle = \vert \Psi_{CC}\rangle = e^{\hat{T}} \vert \Phi_0\rangle = \left( \sum_{n=1}^{N} \frac{1}{n!} \hat{T}^n \right) \vert \Phi_0\rangle, \end{equation*} $$where \( N \) represents the maximum number of particle-hole excitations and \( \hat{T} \) is the cluster operator defined as
$$ \begin{align*} \hat{T} &= \hat{T}_1 + \hat{T}_2 + \ldots + \hat{T}_N \\ \hat{T}_n &= \left(\frac{1}{n!}\right)^2 \sum_{\substack{ i_1,i_2,\ldots i_n \\ a_1,a_2,\ldots a_n}} t_{i_1i_2\ldots i_n}^{a_1a_2\ldots a_n} a_{a_1}^\dagger a_{a_2}^\dagger \ldots a_{a_n}^\dagger a_{i_n} \ldots a_{i_2} a_{i_1}. \end{align*} $$The energy is given by
$$ \begin{equation*} E_{\mathrm{CC}} = \langle\Phi_0\vert \overline{H}\vert \Phi_0\rangle, \end{equation*} $$where \( \overline{H} \) is a similarity transformed Hamiltonian
$$ \begin{align*} \overline{H}&= e^{-\hat{T}} \hat{H}_N e^{\hat{T}} \\ \hat{H}_N &= \hat{H} - \langle\Phi_0\vert \hat{H} \vert \Phi_0\rangle. \end{align*} $$The coupled cluster energy is a function of the unknown cluster amplitudes \( t_{i_1i_2\ldots i_n}^{a_1a_2\ldots a_n} \), given by the solutions to the amplitude equations
$$ \begin{equation}\label{eq:amplitudeeq} 0 = \langle\Phi_{i_1 \ldots i_n}^{a_1 \ldots a_n}\vert \overline{H}\vert \Phi_0\rangle. \end{equation} $$In order to set up the above equations, the similarity transformed Hamiltonian \( \overline{H} \) is expanded using the Baker-Campbell-Hausdorff expression,
$$ \begin{equation}\label{eq:bch} \overline{H}= \hat{H}_N + \left[ \hat{H}_N, \hat{T} \right] + \frac{1}{2} \left[\left[ \hat{H}_N, \hat{T} \right], \hat{T}\right] + \ldots + \frac{1}{n!} \left[ \ldots \left[ \hat{H}_N, \hat{T} \right], \ldots \hat{T} \right] +\dots \end{equation} $$and simplified using the connected cluster theorem
$$ \begin{equation*} \overline{H}= \hat{H}_N + \left( \hat{H}_N \hat{T}\right)_c + \frac{1}{2} \left( \hat{H}_N \hat{T}^2\right)_c + \dots + \frac{1}{n!} \left( \hat{H}_N \hat{T}^n\right)_c +\dots \end{equation*} $$We will discuss parts of the the derivation below.
We will now approximate the cluster operator \( \hat{T} \) to include only \( 2p-2h \) correlations. This leads to the so-called CCD approximation, that is
$$ \hat{T}\approx \hat{T}_2=\frac{1}{4}\sum_{abij}t_{ij}^{ab}a^{\dagger}_aa^{\dagger}_ba_ja_i, $$meaning that we have
$$ \vert \Psi_0 \rangle \approx \vert \Psi_{CCD} \rangle = \exp{\left(\hat{T}_2\right)}\vert \Phi_0\rangle. $$Inserting these equations in the expression for the computation of the energy we have, with a Hamiltonian defined with respect to a general reference vacuum
$$ \hat{H}=\hat{H}_N+E_{\mathrm{ref}}, $$with
$$ \hat{H}_N=\sum_{pq}\langle p \vert \hat{f} \vert q \rangle a^{\dagger}_pa_q + \frac{1}{4}\sum_{pqrs}\langle pq \vert \hat{v} \vert rs \rangle a^{\dagger}_pa^{\dagger}_qa_sa_r, $$we obtain that the energy can be written as
$$ \langle \Phi_0 \vert \exp{\left(-\hat{T}_2\right)}\hat{H}_N\exp{\left(\hat{T}_2\right)}\vert \Phi_0\rangle = \langle \Phi_0 \vert \hat{H}_N(1+\hat{T}_2)\vert \Phi_0\rangle = E_{CCD}. $$This quantity becomes
$$ E_{CCD}=E_{\mathrm{ref}}+\frac{1}{4}\sum_{abij}\langle ij \vert \hat{v} \vert ab \rangle t_{ij}^{ab}, $$where the latter is the correlation energy from this level of approximation of coupled cluster theory. Similarly, the expression for the amplitudes reads
$$ \langle \Phi_{ij}^{ab} \vert \exp{\left(-\hat{T}_2\right)}\hat{H}_N\exp{\left(\hat{T}_2\right)}\vert \Phi_0\rangle = 0. $$These equations can be reduced to (after several applications of Wick's theorem), for all \( i > j \) and all \( a > b \),
$$ \begin{align} 0 = \langle ab \vert \hat{v} \vert ij \rangle + \left(\epsilon_a+\epsilon_b-\epsilon_i-\epsilon_j\right)t_{ij}^{ab}+\frac{1}{2}\sum_{cd} \langle ab \vert \hat{v} \vert cd \rangle t_{ij}^{cd}+\frac{1}{2}\sum_{kl} \langle kl \vert \hat{v} \vert ij \rangle t_{kl}^{ab}+\hat{P}(ij\vert ab)\sum_{kc} \langle kb \vert \hat{v} \vert cj \rangle t_{ik}^{ac} & \nonumber \\ +\frac{1}{4}\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ij}^{cd}t_{kl}^{ab}+\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{ac}t_{jl}^{bd}-\frac{1}{2}\hat{P}(ij)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{ik}^{dc}t_{lj}^{ab}-\frac{1}{2}\hat{P}(ab)\sum_{klcd} \langle kl \vert \hat{v} \vert cd \rangle t_{lk}^{ac}t_{ij}^{db},& \label{eq:ccd} \end{align} $$where we have defined
$$ \hat{P}\left(ab\right)= 1-\hat{P}_{ab}, $$where \( \hat{P}_{ab} \) interchanges two particles occupying the quantum numbers \( a \) and \( b \). The operator \( \hat{P}(ij\vert ab) \) is defined as
$$ \hat{P}(ij\vert ab) = (1-\hat{P}_{ij})(1-\hat{P}_{ab}). $$The single-particle energies \( \epsilon_p \) are normally taken to be either plain harmonic oscillator ones or Hartree-Fock single-particle energies. Recall also that the unknown amplitudes \( t_{ij}^{ab} \) represent anti-symmetrized matrix elements, meaning that they obey the same symmetry relations as the two-body interaction, that is
$$ t_{ij}^{ab}=-t_{ji}^{ab}=-t_{ij}^{ba}=t_{ji}^{ba}. $$The two-body matrix elements are also anti-symmetrized, meaning that
$$ \langle ab \vert \hat{v} \vert ij \rangle = -\langle ab \vert \hat{v} \vert ji \rangle= -\langle ba \vert \hat{v} \vert ij \rangle=\langle ba \vert \hat{v} \vert ji \rangle. $$The non-linear equations for the unknown amplitudes \( t_{ij}^{ab} \) are solved iteratively.
In order to develop a program, chapter 8 of the recent Lecture Notes in Physics (volume 936) is highly recommended as literature. All material is available from the source site. Example of CCD codes are available from the program site. These can be used to benchmark your own program.
Here you should feel free to use either a plain harmonic oscillator basis or Hartree-Fock basis. If you have performed Hartree-Fock calculations and are familiar with these, the Hartree-Fock basis defines the so-called reference energy
$$ \begin{equation} E_{\mathrm{ref}} = \sum_{i\le F} \sum_{\alpha\beta} C^*_{i\alpha}C_{i\beta}\langle \alpha | h | \beta \rangle + \frac{1}{2}\sum_{ij\le F}\sum_{{\alpha\beta\gamma\delta}} C^*_{i\alpha}C^*_{j\beta}C_{i\gamma}C_{j\delta}\langle \alpha\beta|\hat{v}|\gamma\delta\rangle. \label{_auto1} \end{equation} $$If you plan to use Hartree-Fock based matrix elements, you will need to transform the matrix elements from the harmonic oscillator basis to the Hartree-Fock basis. The first step is to program
$$ \begin{equation} \langle pq \vert \hat{v} \vert rs\rangle_{AS}= \sum_{{\alpha\beta\gamma\delta}} C^*_{p\alpha}C^*_{q\beta}C_{r\gamma}C_{s\delta}\langle \alpha\beta|\hat{v}|\gamma\delta\rangle_{AS}, \label{_auto2} \end{equation} $$where the coefficients are those from the last Hartree-Fock iteration and the matrix elements are all anti-symmetrized. You can extend your Hartree-Fock program to write out these matrix elements after the last Hartree-Fock iteration. Make sure that your matrix elements are structured according to conserved quantum numbers, avoiding thereby the write out of many zeros.
To test that your matrix elements are set up correctly, when you read in these matrix elements in the CCD code, make sure that the reference energy from your Hartree-Fock calculations are reproduced. Alternatively, you can just use the standard harmonic oscillator one-body and two-body matrix elements.
Set up a code which solves the CCD equation by encoding the equations as they stand, that is follow the mathematical expressions and perform the sums over all single-particle states. Compute the energy of the two-electron systems using all single-particle states. Compare these with Taut's results for \( \omega=1 \) a.u. Since you do not include singles you will not get the exact result. If you wish to include singles, you will able to obtain the exact results in a basis with at least ten major oscillator shells. Perform also calculations with \( N=6 \), \( N=12 \) and \( N=20 \) electrons and compare with reference [2] of Pedersen et al below.
The next step consists in rewriting the equations in terms of matrix-matrix multiplications and subdividing the matrix elements and operations in terms of two-particle configuration that conserve total spin projection and projection of the orbital momentum. Rewrite also the equations in terms of so-called intermediates, as detailed in section 8.7 of Lietz et al. This section gives a detailed description on how to build a coupled cluster code and is highly recommended.
Rerun your calculations for \( =2 \), \( N=6 \), \( N=12 \) and \( N=20 \) electrons using your optimal Hartree-Fock basis. Make sure your results from 2b) stay the same.
Calculate as well ground state energies for \( \omega=0.5 \) and \( \omega=0.1 \). Try to compare with eventual variational Monte Carlo results from other students, if possible.
The final step is to parallelize your CCD code using either OpenMP or MPI and do a performance analysis. Use the \( N=6 \) case. Make a performance analysis by timing your serial code with and without vectorization. Perform several runs and compute an average timing analysis with and without vectorization. Comment your results.
Compare thereafter your serial code(s) with the speedup you get by parallelizing your code, running either OpenMP or MPI or both. Do you get a near \( 100\% \) speedup with the parallel version? Comment again your results and perform timing benchmarks several times in order to extract an average performance time.
The Hermite polynomials are the solutions of the following differential equation
$$ \begin{equation} \frac{d^2H(x)}{dx^2}-2x\frac{dH(x)}{dx}+ (\lambda-1)H(x)=0. \label{eq:hermite} \end{equation} $$The first few polynomials are
$$ \begin{equation*} H_0(x)=1, \end{equation*} $$ $$ \begin{equation*} H_1(x)=2x, \end{equation*} $$ $$ \begin{equation*} H_2(x)=4x^2-2, \end{equation*} $$ $$ \begin{equation*} H_3(x)=8x^3-12x, \end{equation*} $$and
$$ \begin{equation*} H_4(x)=16x^4-48x^2+12. \end{equation*} $$They fulfil the orthogonality relation
$$ \begin{equation*} \int_{-\infty}^{\infty}e^{-x^2}H_n(x)^2dx=2^nn!\sqrt{\pi}, \end{equation*} $$and the recursion relation
$$ \begin{equation*} H_{n+1}(x)=2xH_{n}(x)-2nH_{n-1}(x). \end{equation*} $$Here follows a brief recipe and recommendation on how to write a report for each project.
The preferred format for the report is a PDF file. You can also use DOC or postscript formats or as an ipython notebook file. As programming language we prefer that you choose between C/C++, Fortran2008 or Python. The following prescription should be followed when preparing the report:
Finally, we encourage you to work two and two together. Optimal working groups consist of 2-3 students. You can then hand in a common report.