The coupled-cluster method is an efficient tool to compute properties of many-body systems with an effort that grows polynomial with system size. While this might still be expensive, it is now possible to compute for example atomic nuclei with mass numbers about \( A\approx 208 \) (which corresponds to lead) with this method. Recall that full configuration interaction (FCI) exhibits an exponential cost and is therefore limited to systems with few degrees of freedom.
For some recent applications see
We start from the reference state
$$ \begin{equation} \vert\Phi_0\rangle = \prod_{i=1}^N a^\dagger_i \vert 0\rangle \label{_auto1} \end{equation} $$for the description of a system with \( N \) particles. Usually, this reference is the Hartree-Fock state, but that is not necessary as discussed throughout this course.
After we introduced the particle-hole formalism, we have opted for a convention where the indices \( i,j,k,\ldots \) run over hole states, i.e. orbitals occupied in the reference state, while \( a,b,c,\ldots \) run over particle states, i.e. unoccupied orbitals. Indices \( p,q,r,s \) can identify any orbital. Let \( n_u \) be the number of unoccupied states, and \( N \) is of course the number of particles, which is also the number of occupied states used to define the ansatz for the ground state.
We consider the Hamiltonian
$$ \begin{equation} \label{Ham} H = \sum_{pq} \langle p\vert h_0\vert q \rangle a^\dagger_p a_q + \frac{1}{4}\sum_{pqrs}\langle pq\vert v\vert rs\rangle a^\dagger_pa^\dagger_q a_sa_r \end{equation} $$The reference state is a non-trivial vacuum of our theory. We normal order this Hamiltonian with respect to the nontrivial vacuum state given by the Hartree-Fock reference and obtain the normal-ordered Hamiltonian
$$ \begin{equation} \label{HN} H_N = \sum_{pq} \langle p\vert f\vert q \rangle \left\{a^\dagger_p a_q\right\} + \frac{1}{4}\sum_{pqrs}\langle pq\vert V\vert rs\rangle \left\{a^\dagger_pa^\dagger_q a_sa_r\right\}. \end{equation} $$Here,
$$ \begin{equation} \label{Fock} \langle p\vert f\vert q \rangle = \langle p\vert h_0\vert q \rangle + \sum_i \langle pi\vert V\vert qi\rangle \end{equation} $$is the Fock matrix. We note that the Fock matrix is diagonal in the Hartree-Fock basis. The brackets \( \{\cdots\} \) denote normal ordering, i.e. all operators that annihilate the nontrivial vaccum are to the right of those operators that create with respect to that vaccum. Normal ordering implies that \( \langle \Phi_0\vert H_N\vert \Phi_0\rangle = 0 \).
The coupled-cluster method is a very efficient tool to compute properties of many-particle systems when a "good" reference state is available. Let us assume that the reference state results from a Hartree-Fock calculation.
How do you know whether a Hartree-Fock state is a "good" reference? Which results of the Hartree-Fock computation will inform you?
Once the Hartree-Fock equations are solved, the Fock matrix becomes diagonal, and its diagonal elements can be viewed as single-particle energies. Hopefully, there is a clear gap in the single-particle spectrum at the Fermi surface, i.e. after \( N \) orbitals are filled.
If symmetry-restricted Hartree-Fock is used, one is limited to compute systems with closed subshells for neutrons and for protons. On a first view, this might seem as a severe limitation. But is it?
If one limits oneself to nuclei with mass number up to mass number \( A=60 \), how many nuclei can potentially be described with the coupled-cluster method? Which of these nuclei are potentially interesting? Why?
Nuclear shell closures are at \( N,Z=2,8,20,28,50,82,126 \), and subshell closures at \( N,Z=2,6,8,14,16,20,28,32,34,40,50,\ldots \).
In the physics of nuclei, the evolution of nuclear structure as neutrons are added (or removed) from an isotope is a key
interest. Examples are the rare isotopes of helium (He-8,10) oxygen (O-22,24,28), calcium (Ca-52,54,60), nickel (Ni-78) and tin (Sn-100,132). The coupled-cluster method has the potential to address questions regarding these nuclei, and in a several cases was used to make predictions before experimental data was available. In addition, the method can be used to compute neighbors of nuclei with closed subshells.
We had
$$ (\hat{H} -E)|\Psi_0\rangle = (\hat{H} -E)\sum_{P'H'}C_{H'}^{P'}|\Phi_{H'}^{P'} \rangle=0. $$As stated previously, one solves this equation normally by diagonalization. If we are able to solve this equation exactly (that is numerically exactly) in a large Hilbert space (it will be truncated in terms of the number of single-particle states included in the definition of Slater determinants), it can then serve as a benchmark for other many-body methods which approximate the correlation operator \( \hat{C} \).
For reasons to come (link with Coupled-Cluster theory and Many-Body perturbation theory), we will rewrite the full FCI equation as a set of coupled non-linear equations in terms of the unknown coefficients \( C_H^P \).
To see this, we look at $ \langle \Phi_H^P | = \langle \Phi_0 |$, that is we multiply with \( \langle \Phi_0 | \) from the left in
$$ (\hat{H} -E)\sum_{P'H'}C_{H'}^{P'}|\Phi_{H'}^{P'} \rangle=0, $$and we assume that we have a two-body operator at most. Using the Condon-Slater rules gives then and equation for the correlation energy in terms of \( C_i^a \) and \( C_{ij}^{ab} \). We get then
$$ \langle \Phi_0 | \hat{H} -E| \Phi_0\rangle + \sum_{ai}\langle \Phi_0 | \hat{H} -E|\Phi_{i}^{a} \rangle C_{i}^{a}+ \sum_{abij}\langle \Phi_0 | \hat{H} -E|\Phi_{ij}^{ab} \rangle C_{ij}^{ab}=0, $$or
$$ E-E_0 =\Delta E=\sum_{ai}\langle \Phi_0 | \hat{H}|\Phi_{i}^{a} \rangle C_{i}^{a}+ \sum_{abij}\langle \Phi_0 | \hat{H}|\Phi_{ij}^{ab} \rangle C_{ij}^{ab}, $$where the \( E_0 \) is the reference energy and \( \Delta E \) becomes the correlation energy. We have already computed the expectation values \( \langle \Phi_0 | \hat{H}|\Phi_{i}^{a} \) and \( \langle \Phi_0 | \hat{H}|\Phi_{ij}^{ab}\rangle \).
We can rewrite
$$ E-E_0 =\Delta E=\sum_{ai}\langle \Phi_0 | \hat{H}|\Phi_{i}^{a} \rangle C_{i}^{a}+ \sum_{abij}\langle \Phi_0 | \hat{H}|\Phi_{ij}^{ab} \rangle C_{ij}^{ab}, $$as
$$ \Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle C_{i}^{a}+ \sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}. $$This equation determines the correlation energy but not the coefficients \( C \). We need more equations. Our next step is to mulitply from the left with \( \Phi_i^a \)
$$ \langle \Phi_i^a | \hat{H} -E| \Phi_0\rangle + \sum_{bj}\langle \Phi_i^a | \hat{H} -E|\Phi_{j}^{b} \rangle C_{j}^{b}+ \sum_{bcjk}\langle \Phi_i^a | \hat{H} -E|\Phi_{jk}^{bc} \rangle C_{jk}^{bc}+ \sum_{bcdjkl}\langle \Phi_i^a | \hat{H} -E|\Phi_{jkl}^{bcd} \rangle C_{jkl}^{bcd}=0, $$as this equation will allow us to find an expression for the coefficents \( C_i^a \) since we can rewrite this equation as
$$ \langle i | \hat{f}| a\rangle +\langle \Phi_i^a | \hat{H}-E|\Phi_{i}^{a} \rangle C_{i}^{a}+ \sum_{bj\ne ai}\langle \Phi_i^a | \hat{H}|\Phi_{j}^{b} \rangle C_{j}^{b}+ \sum_{bcjk}\langle \Phi_i^a | \hat{H}|\Phi_{jk}^{bc} \rangle C_{jk}^{bc}+ \sum_{bcdjkl}\langle \Phi_i^a | \hat{H}|\Phi_{jkl}^{bcd} \rangle C_{jkl}^{bcd}=0. $$We rewrite this equation as
$$ C_{i}^{a}=-(\langle \Phi_i^a | \hat{H}-E|\Phi_{i}^{a})^{-1}\left(\langle i | \hat{f}| a\rangle+ \sum_{bj\ne ai}\langle \Phi_i^a | \hat{H}|\Phi_{j}^{b} \rangle C_{j}^{b}+\sum_{bcjk}\langle \Phi_i^a | \hat{H}|\Phi_{jk}^{bc} \rangle C_{jk}^{bc}+ \sum_{bcdjkl}\langle \Phi_i^a | \hat{H}|\Phi_{jkl}^{bcd} \rangle C_{jkl}^{bcd}\right). $$Since these equations are solved iteratively ( that is we can start with a guess for the coefficients \( C_i^a \)), it is common to start the iteration by setting
$$ C_{i}^{a}=-\frac{\langle i | \hat{f}| a\rangle}{\langle \Phi_i^a | \hat{H}-E|\Phi_{i}^{a}\rangle}, $$and the denominator can be written as
$$ C_{i}^{a}=\frac{\langle i | \hat{f}| a\rangle}{\langle i | \hat{f}| i\rangle-\langle a | \hat{f}| a\rangle+\langle ai | \hat{v}| ai\rangle-E}. $$The observant reader will however see that we need an equation for \( C_{jk}^{bc} \) and \( C_{jkl}^{bcd} \) as well. To find equations for these coefficients we need then to continue our multiplications from the left with the various \( \Phi_{H}^P \) terms.
For \( C_{jk}^{bc} \) we need then
$$ \langle \Phi_{ij}^{ab} | \hat{H} -E| \Phi_0\rangle + \sum_{kc}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{k}^{c} \rangle C_{k}^{c}+ \sum_{cdkl}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{kl}^{cd} \rangle C_{kl}^{cd}+ $$ $$ \sum_{cdeklm}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{klm}^{cde} \rangle C_{klm}^{cde}+\sum_{cdefklmn}\langle \Phi_{ij}^{ab} | \hat{H} -E|\Phi_{klmn}^{cdef} \rangle C_{klmn}^{cdef}=0, $$and we can isolate the coefficients \( C_{kl}^{cd} \) in a similar way as we did for the coefficients \( C_{i}^{a} \).
At the end we can rewrite our solution of the Schr\"odinger equation in terms of \( n \) coupled equations for the coefficients \( C_H^P \). This is a very cumbersome way of solving the equation. However, by using this iterative scheme we can illustrate how we can compute the various terms in the wave operator or correlation operator \( \hat{C} \). We will later identify the calculation of the various terms \( C_H^P \) as parts of different many-body approximations to full CI. In particular, we will relate this non-linear scheme with Coupled Cluster theory and many-body perturbation theory.
If we use a Hartree-Fock basis, how can one simplify the equation
$$ \Delta E=\sum_{ai}\langle i| \hat{f}|a \rangle C_{i}^{a}+ \sum_{abij}\langle ij | \hat{v}| ab \rangle C_{ij}^{ab}? $$There are several ways to view and understand the coupled-cluster method. A first simple view of coupled-cluster theory is that the method induces correlations into the reference state by expressing a correlated state as
$$ \begin{equation} \label{psi} \vert\Psi\rangle = e^T \vert\Phi_0\rangle , \end{equation} $$Here, \( T \) is an operator that induces correlations. We can now demand that the correlate state \eqref{psi} becomes and eigenstate of the Hamiltonian \( H_N \), i.e. \( H\vert \Psi\rangle = E\vert \Psi\rangle \). This view, while correct, is not the most productive one. Instead, we left-multiply the Schroedinger equation with \( e^{-T} \) and find
$$ \begin{equation} \label{Schroedinger} \overline{H_N}\vert \Phi_0\rangle = \Delta E \vert \Phi_0\rangle . \end{equation} $$Here, \( \Delta E \) is the correlation energy, and the total energy is \( E=\Delta E+E_{0}^{\mathrm{ref}} \).
The similarity-transformed Hamiltonian is defined as
$$ \begin{equation} \label{Hsim} \overline{H_N} \equiv e^{-T} H_N e^T . \end{equation} $$A more productive view on coupled-cluster theory thus emerges: This method seeks a similarity transformation such that the uncorrelated reference state becomes an exact eigenstate of the similarity-transformed Hamiltonian \eqref{Hsim}.
What are the conditions on \( T \) such that \( \overline{H_N} \) is Hermitian?
For a Hermitian \( \overline{H_N} \), we need a unitary \( e^T \), i.e. an anti-Hermitian \( T \) with \( T = -T^\dagger \)
As we will see below, coupld-cluster theory employs a non-Hermitian Hamiltonian.
Show that \( \overline{H_N} \) has the same eigenvalues as \( H_N \) for arbitrary \( T \). What is the spectral decomposition of a non-Hermitian \( \overline{H_N} \) ?
Let \( H_N\vert E\rangle = E\vert E\rangle \). Thus
$$ \begin{align*} H_N e^{T} e^{-T} \vert E\rangle &= E\vert E\rangle , \\ \left(e^{-T} H_N e^T\right) e^{-T} \vert E\rangle &= Ee^{-T} \vert E\rangle , \\ \overline{H_N} e^{-T} \vert E\rangle &= E e^{-T}\vert E\rangle . \end{align*} $$Thus, if \( \vert E\rangle \) is an eigenstate of \( H_N \) with eigenvalue \( E \), then \( e^{-T}\vert E\rangle \) is eigenstate of \( \overline{H_N} \) with the same eigenvalue.
A non-Hermitian \( \overline{H_N} \) has eigenvalues \( E_\alpha \) corresponding to left \( \langle L_\alpha\vert \) and right \( \vert R_\alpha \rangle \) eigenstates. Thus
$$ \begin{align} \overline{H_N} = \sum_\alpha \vert R_\alpha\rangle E_\alpha \langle L_\alpha \vert \label{_auto2} \end{align} $$with bi-orthonormal \( \langle L_\alpha\vert R_\beta\rangle = \delta_\alpha^\beta \).
To make progress, we have to specify the cluster operator \( T \). In coupled cluster theory, this operator is
$$ \begin{equation} \label{Top} T \equiv \sum_{ia} t_i^a a^\dagger_a a_i + \frac{1}{4}\sum_{ijab}t_{ij}^{ab} a^\dagger_aa^\dagger_ba_ja_i + \cdots + \frac{1}{(N!)^2}\sum_{i_1\ldots i_N a_1 \ldots a_N} t_{i_1\ldots i_N}^{a_1\ldots a_N} a^\dagger_{a_1}\cdots a^\dagger_{a_N} a_{i_N}\cdots a_{i_1} . \end{equation} $$Thus, the operator \eqref{Top} induces particle-hole (p-h) excitations with respect to the reference. In general, \( T \) generates up to \( Ap-Ah \) excitations, and the unknown parameters are the cluster amplitides \( t_i^a \), \( t_{ij}^{ab} \), ..., \( t_{i_1,\ldots,i_A}^{a_1,\ldots,a_A} \).
Thus, the coupled-cluster method with the full cluster operator \eqref{Top} is exponentially expensive, just as FCI. To make progress, we need to make an approximation by truncating the operator. Here, we will use the CCSD (coupled clusters singles doubles) approximation, where
$$ \begin{equation} \label{Tccsd} T \equiv \sum_{ia} t_i^a a^\dagger_a a_i + \frac{1}{4}\sum_{ijab}t_{ij}^{ab} a^\dagger_aa^\dagger_ba_ja_i . \end{equation} $$We need to determine the unknown cluster amplitudes that enter in CCSD. Let
$$ \begin{align} \vert\Phi_i^a\rangle &= a^\dagger_a a_i \vert \Phi_0\rangle , \label{_auto3}\\ \vert\Phi_{ij}^{ab}\rangle &= a^\dagger_a a^\dagger_b a_j a_i \vert \Phi_0\rangle \label{_auto4} \end{align} $$be 1p-1h and 2p-2h excitations of the reference. Computing matrix elements of the Schroedinger Equation \eqref{Schroedinger} yields
$$ \begin{align} \label{ccsd} \langle \Phi_0\vert \overline{H_N}\vert \Phi_0\rangle &= E_c , \\ \langle \Phi_i^a\vert \overline{H_N}\vert \Phi_0\rangle &= 0 , \label{_auto5}\\ \langle \Phi_{ij}^{ab}\vert \overline{H_N}\vert \Phi_0\rangle &= 0 . \label{_auto6} \end{align} $$The first equation states that the coupled-cluster correlation energy is an expectation value of the similarity-transformed Hamiltonian. The second and third equations state that the similarity-transformed Hamiltonian exhibits no 1p-1h and no 2p-2h excitations. These equations have to be solved to find the unknown amplitudes \( t_i^a \) and \( t_{ij}^{ab} \). Then one can use these amplitudes and compute the correlation energy from the first line of Eq. \eqref{ccsd}.
We note that in the CCSD approximation the reference state is not an exact eigenstates. Rather, it is decoupled from simple states but \( \overline{H} \) still connects this state to 3p-3h, and 4p-4h states etc.
At this point, it is important to recall that we assumed starting from a "good" reference state. In such a case, we might reasonably expect that the inclusion of 1p-1h and 2p-2h excitations could result in an accurate approximation. Indeed, empirically one finds that CCSD accounts for about 90% of the corelation energy, i.e. of the difference between the exact energy and the Hartree-Fock energy. The inclusion of triples (3p-3h excitations) typically yields 99% of the correlation energy.
We see that the coupled-cluster method in its CCSD approximation yields a similarity-transformed Hamiltonian that is of a two-body structure with respect to a non-trivial vacuum. When viewed in this light, the coupled-cluster method "transforms" an \( A \)-body problem (in CCSD) into a two-body problem, albeit with respect to a nontrivial vacuum.
Above we argued that a similarity transformation preserves all eigenvalues. Nevertheless, the CCD correlation energy is not the exact correlation energy. Explain!
The CCD approximation does not make \( \vert\Phi_0\rangle \) an exact eigenstate of \( \overline{H_N} \); it is only an eigenstate when the similarity-transformed Hamiltonian is truncated to at most 2p-2h states. The full \( \overline{H_N} \), with \( T=T_2 \), would involve six-body terms (do you understand this?), and this full Hamiltonian would reproduce the exact correlation energy. Thus CCD is a similarity transformation plus a truncation, which decouples the ground state only from 2p-2h states.
The solution of the CCSD equations, i.e. the second and third line of Eq. \eqref{ccsd}, and the computation of the correlation energy requires us to compute matrix elements of the similarity-transformed Hamiltonian \eqref{Hsim}. This can be done with the Baker-Campbell-Hausdorff expansion
$$ \begin{align} \overline{H_N} &= e^{-T} H_N e^T \label{_auto7}\\ &=H_N + \left[ H_N, T\right]+ \frac{1}{2!}\left[ \left[ H_N, T\right], T\right] + \frac{1}{3!}\left[\left[ \left[ H_N, T\right], T\right], T\right] +\ldots . \label{_auto8} \end{align} $$We now come to a key element of coupled-cluster theory: the cluster operator \eqref{Top} consists of sums of terms that consist of particle creation and hole annihilation operators (but no particle annihilation or hole creation operators). Thus, all terms that enter \( T \) commute with each other. This means that the commutators in the Baker-Campbell-Hausdorff expansion can only be non-zero because each \( T \) must connect to \( H_N \) (but no \( T \) with another \( T \)). Thus, the expansion is finite.
We see that the (disadvantage of having a) non-Hermitian Hamiltonian \( \overline{H_N} \) leads to the advantage that the Baker-Campbell-Hausdorff expansion is finite, thus leading to the possibility to compute \( \overline{H_N} \) exactly.
We write the similarity-transformed Hamiltonian as
$$ \begin{align} \overline{H_N}=\sum_{pq} \overline{H}^p_q a^\dagger_q a_p + {1\over 4} \sum_{pqrs} \overline{H}^{pq}_{rs} a^\dagger_p a^\dagger_q a_s a_r + \ldots \label{_auto9} \end{align} $$with
$$ \begin{align} \overline{H}^p_q &\equiv \langle p\vert \overline{H_N}\vert q\rangle , \label{_auto10}\\ \overline{H}^{pq}_{rs} &\equiv \langle pq\vert \overline{H_N}\vert rs\rangle . \label{_auto11} \end{align} $$If we truncate the operator \( T \) at the level of singles and doubles, the CCSD Eqs. \eqref{ccsd} for the amplitudes can be written as \( \overline{H}_i^a = 0 \) and \( \overline{H}_{ij}^{ab}=0 \).
In what follows, we will consider the coupled cluster doubles (CCD) approximation. This approximation is valid in cases where the system cannot exhibit any particle-hole excitations (such as nuclear matter when formulated on a momentum-space grid) or for the pairing model (as the pairing interactions only excites pairs of particles). In this case \( t_i^a=0 \) for all \( i, a \), and \( \overline{H}_i^a=0 \). The CCD approximation is also of some sort of leading order approximation in the Hartree-Fock basis (as the Hartree-Fock Hamiltonian exhibits no particle-hole excitations).
Let us consider the matrix element \( \overline{H}_{ij}^{ab} \). Clearly, it consists of all diagrams (i.e. all combinations of \( T_2 \), and a single \( F \) or \( V \) that have two incoming hole lines and two outgoing particle lines. Write down all these diagrams.
We start systematically and consider all combinations of \( F \) and \( V \) diagrams with 0, 1, and 2 cluster amplitudes \( T_2 \).
You learned about the pairing Hamiltonian earlier in the second midterm. Convince yourself that this Hamiltonian does not induce any 1p-1h excitations. Let us solve the CCD equations for this problem. This consists of the following steps
The Hamiltonian is
$$ \begin{align} H = \delta \sum_{p=1}^\Omega (p-1)\left(a^\dagger_{p+}a_{p+} + a^\dagger_{p-}a_{p-}\right) -{g \over 2} \sum_{p, q=1}^\Omega a^\dagger_{p+}a^\dagger_{p-} a_{q-} a_{q+} . \label{_auto12} \end{align} $$#!/usr/bin/python
from sympy import *
from pylab import *
import matplotlib.pyplot as plt
ga = linspace(-1,1,20)
e1 = []
for g_val in ga:
H1 = matrix([[2-g_val , -g_val/2., -g_val/2., -g_val/2., -g_val/2., 0],
[-g_val/2., 4-g_val, -g_val/2., -g_val/2., 0., -g_val/2.],
[-g_val/2., -g_val/2., 6-g_val, 0, -g_val/2., -g_val/2.],
[-g_val/2., -g_val/2., 0, 6-g_val, -g_val/2., -g_val/2.],
[-g_val/2., 0, -g_val/2., -g_val/2., 8-g_val, -g_val/2.],
[0 , -g_val/2., -g_val/2., -g_val/2., -g_val/2., 10-g_val]])
u1, v1 = linalg.eig(H1)
e1.append(min(u1))
exact = e1 - (2-ga)
print(exact)
plt.axis([-1,1,-0.5,0.05])
plt.xlabel(r'Interaction strength, $g$', fontsize=16)
plt.ylabel(r'Correlation energy', fontsize=16)
exact = plt.plot(ga, exact,'b-*',linewidth = 2.0, label = 'Exact')
plt.legend()
plt.savefig('pairing.pdf', format='pdf')
plt.show()
## Coupled clusters in CCD approximation
## Implemented for the pairing model of Lecture Notes in Physics 936, Chapter 8.
import numpy as np
def init_pairing_v(g,pnum,hnum):
"""
returns potential matrices of the pairing model in three relevant channels
param g: strength of the pairing interaction, as in Eq. (8.42)
param pnum: number of particle states
param hnum: number of hole states
return v_pppp, v_pphh, v_hhhh: np.array(pnum,pnum,pnum,pnum),
np.array(pnum,pnum,hnum,hnum),
np.array(hnum,hnum,hnum,hnum),
The interaction as a 4-indexed tensor in three channels.
"""
v_pppp=np.zeros((pnum,pnum,pnum,pnum))
v_pphh=np.zeros((pnum,pnum,hnum,hnum))
v_hhhh=np.zeros((hnum,hnum,hnum,hnum))
gval=-0.5*g
for a in range(0,pnum,2):
for b in range(0,pnum,2):
v_pppp[a,a+1,b,b+1]=gval
v_pppp[a+1,a,b,b+1]=-gval
v_pppp[a,a+1,b+1,b]=-gval
v_pppp[a+1,a,b+1,b]=gval
for a in range(0,pnum,2):
for i in range(0,hnum,2):
v_pphh[a,a+1,i,i+1]=gval
v_pphh[a+1,a,i,i+1]=-gval
v_pphh[a,a+1,i+1,i]=-gval
v_pphh[a+1,a,i+1,i]=gval
for j in range(0,hnum,2):
for i in range(0,hnum,2):
v_hhhh[j,j+1,i,i+1]=gval
v_hhhh[j+1,j,i,i+1]=-gval
v_hhhh[j,j+1,i+1,i]=-gval
v_hhhh[j+1,j,i+1,i]=gval
return v_pppp, v_pphh, v_hhhh
def init_pairing_fock(delta,g,pnum,hnum):
"""
initializes the Fock matrix of the pairing model
param delta: Single-particle spacing, as in Eq. (8.41)
param g: pairing strength, as in eq. (8.42)
param pnum: number of particle states
param hnum: number of hole states
return f_pp, f_hh: The Fock matrix in two channels as numpy arrays np.array(pnum,pnum), np.array(hnum,hnum).
"""
# the Fock matrix for the pairing model. No f_ph needed, because we are in Hartree-Fock basis
deltaval=0.5*delta
gval=-0.5*g
f_pp = np.zeros((pnum,pnum))
f_hh = np.zeros((hnum,hnum))
for i in range(0,hnum,2):
f_hh[i ,i ] = deltaval*i+gval
f_hh[i+1,i+1] = deltaval*i+gval
for a in range(0,pnum,2):
f_pp[a ,a ] = deltaval*(hnum+a)
f_pp[a+1,a+1] = deltaval*(hnum+a)
return f_pp, f_hh
def init_t2(v_pphh,f_pp,f_hh):
"""
Initializes t2 amlitudes as in MBPT2, see first equation on page 345
param v_pphh: pairing tensor in pphh channel
param f_pp: Fock matrix in pp channel
param f_hh: Fock matrix in hh channel
return t2: numpy array in pphh format, 4-indices tensor
"""
pnum = len(f_pp)
hnum = len(f_hh)
t2_new = np.zeros((pnum,pnum,hnum,hnum))
for i in range(hnum):
for j in range(hnum):
for a in range(pnum):
for b in range(pnum):
t2_new[a,b,i,j] = v_pphh[a,b,i,j] / (f_hh[i,i]+f_hh[j,j]-f_pp[a,a]-f_pp[b,b])
return t2_new
# CCD equations. Note that the "->abij" assignment is redundant, because indices are ordered alphabetically.
# Nevertheless, we retain it for transparency.
def ccd_iter(v_pppp,v_pphh,v_hhhh,f_pp,f_hh,t2):
"""
Performs one iteration of the CCD equations (8.34), using also intermediates for the nonliniar terms
param v_pppp: pppp-channel pairing tensor, numpy array
param v_pphh: pphh-channel pairing tensor, numpy array
param v_hhhh: hhhh-channel pairing tensor, numpy array
param f_pp: Fock matrix in pp channel
param f_hh: Fock matrix in hh channel
param t2: Initial t2 amplitude, tensor in form of pphh channel
return t2_new: new t2 amplitude, tensor in form of pphh channel
"""
pnum = len(f_pp)
hnum = len(f_hh)
Hbar_pphh = ( v_pphh
+ np.einsum('bc,acij->abij',f_pp,t2)
- np.einsum('ac,bcij->abij',f_pp,t2)
- np.einsum('abik,kj->abij',t2,f_hh)
+ np.einsum('abjk,ki->abij',t2,f_hh)
+ 0.5*np.einsum('abcd,cdij->abij',v_pppp,t2)
+ 0.5*np.einsum('abkl,klij->abij',t2,v_hhhh)
)
# hh intermediate, see (8.47)
chi_hh = 0.5* np.einsum('cdkl,cdjl->kj',v_pphh,t2)
Hbar_pphh = Hbar_pphh - ( np.einsum('abik,kj->abij',t2,chi_hh)
- np.einsum('abik,kj->abji',t2,chi_hh) )
# pp intermediate, see (8.46)
chi_pp = -0.5* np.einsum('cdkl,bdkl->cb',v_pphh,t2)
Hbar_pphh = Hbar_pphh + ( np.einsum('acij,cb->abij',t2,chi_pp)
- np.einsum('acij,cb->baij',t2,chi_pp) )
# hhhh intermediate, see (8.48)
chi_hhhh = 0.5 * np.einsum('cdkl,cdij->klij',v_pphh,t2)
Hbar_pphh = Hbar_pphh + 0.5 * np.einsum('abkl,klij->abij',t2,chi_hhhh)
# phph intermediate, see (8.49)
chi_phph= + 0.5 * np.einsum('cdkl,dblj->bkcj',v_pphh,t2)
Hbar_pphh = Hbar_pphh + ( np.einsum('bkcj,acik->abij',chi_phph,t2)
- np.einsum('bkcj,acik->baij',chi_phph,t2)
- np.einsum('bkcj,acik->abji',chi_phph,t2)
+ np.einsum('bkcj,acik->baji',chi_phph,t2) )
t2_new=np.zeros((pnum,pnum,hnum,hnum))
for i in range(hnum):
for j in range(hnum):
for a in range(pnum):
for b in range(pnum):
t2_new[a,b,i,j] = ( t2[a,b,i,j]
+ Hbar_pphh[a,b,i,j] / (f_hh[i,i]+f_hh[j,j]-f_pp[a,a]-f_pp[b,b]) )
return t2_new
def ccd_energy(v_pphh,t2):
"""
Computes CCD energy. Call as
energy = ccd_energy(v_pphh,t2)
param v_pphh: pphh-channel pairing tensor, numpy array
param t2: t2 amplitude, tensor in form of pphh channel
return energy: CCD correlation energy
"""
erg = 0.25*np.einsum('abij,abij',v_pphh,t2)
return erg
###############################
######## Main Program
# set parameters as for model
pnum = 4 # number of particle states
hnum = 4 # number of hole states
delta = 1.0
g = 1.0
print("parameters")
print("delta =", delta, ", g =", g)
# Initialize pairing matrix elements and Fock matrix
v_pppp, v_pphh, v_hhhh = init_pairing_v(g,pnum,hnum)
f_pp, f_hh = init_pairing_fock(delta,g,pnum,hnum)
# Initialize T2 amplitudes from MBPT2
t2 = init_t2(v_pphh,f_pp,f_hh)
erg = ccd_energy(v_pphh,t2)
# Exact MBPT2 for comparison, see last equation on page 365
exact_mbpt2 = -0.25*g**2*(1.0/(2.0+g) + 2.0/(4.0+g) + 1.0/(6.0+g))
print("MBPT2 energy =", erg, ", compared to exact MBPT(2):", exact_mbpt2)
# iterate CCD equations niter times
niter=60
for iter in range(niter):
t2_new = ccd_iter(v_pppp,v_pphh,v_hhhh,f_pp,f_hh,t2)
erg = ccd_energy(v_pphh,t2_new)
print("iter=", iter, "erg=", erg)
t2 = 0.5 * (t2_new + t2)