Second Quantization¶

Second quantization is a convenient formalism in quantum chemistry that allows one to represent quantum mechanical operators and wave functions in terms of operators that create or destroy particles (creation and annihilation operators, respectively). This project provides an overview of the properties of Fermionic creation and annihilation operators and introduces the second-quantized forms of one- and two-particle operators and many-particle wave functions. It also covers the concept of normal order, which is useful for evaluating integrals involving second-quantized operators and wave functions. Python code is provided for bringing strings of second-quantized operators to normal order, using the $p^\dagger q$ package.

Fermionic Creation and Annihilation Operators¶

Central to the second quantization formalism is the concept of Fock space, which is the space spanned by functions with different particle numbers (i.e., zero-, one-, two-particle functions, etc.). Given a set of $m$ spin-orbitals, $\{\phi\}$, the occupation number representation of an $N$-particle state $|\psi^N\rangle$ is

$$ |\psi^N\rangle = | n_0, n_1, n_2, ..., n_m \rangle $$

where $n_p$ is the occupation of the $p$th orbital (0 or 1), and $\sum_p n_p = N$. This state is just a Slater determinant of $N$ occupied spin-orbitals. Different particle-number states can be accessed via the application of creation and annihilation operators to $|\psi^N\rangle$. In this project, we are primarily concerned with electronic wave functions, so we focus our attention on Fermionic creation and annihilation operators. We define the operators $\hat{a}^\dagger_p$ and $\hat{a}_q$ such that they create an electron in spin-orbital $\phi_p$ or destroy an electron in spin-orbital $\phi_q$, respectively. Consider the action of these operators on an electronic state that involves only a single spin-orbital, $\phi_k$:

$$ \begin{align} \hat{a}_k | 0_k \rangle &= 0 \\ \hat{a}_k | 1_k \rangle &= |0_k\rangle \\ \hat{a}^\dagger_k | 0_k \rangle &= |1_k\rangle \\ \hat{a}^\dagger_k | 1_k \rangle &= 0 \\ \end{align} $$

The right-hand side of the first equation is zero because one cannot annihilate an electron in an empty spin-orbital. The right-hand side of the last equation is zero because the spin-orbitals can contain a maximum of one electron. Acting on a multielectron state is similar, except that we must be a little careful because we need to account for the antisymmetry of the state. For example, consider an annihilator acting on a two-electron / two-orbital state

$$ \hat{a}_1 = | 1_0, 1_1 \rangle = -\hat{a}_1| 1_1, 1_0 \rangle = -|0_1 1_0\rangle = -|1_0 0_1 \rangle$$

The minus sign resuts from the fact that we need to move orbital 1 to the front of the occupation number vector in order to annihilate the electron in that orbital, and we accumulate a minus sign every time we swap two occupied orbitals. For a general Fock state, we would have

$$ \begin{align} \hat{a}_k | n_0, n_1, ..., 1_k, n_{k+1}, ..., n_m \rangle &= (-1)^{\sum_{l < k} n_l} \hat{a}_k | 1_k, n_0, n_1, ..., n_{k+1}, ..., n_m \rangle \\ &= (-1)^{\sum_{l < k} n_l} | 0_k, n_0, n_1, ..., n_{k+1}, ..., n_m \rangle \\ & = (-1)^{\sum_{l < k} n_l} | n_0, n_1, ..., 0_k, n_{k+1}, ..., n_m \rangle \end{align} $$

Similar considerations should be made when adding electrons to a Fock state with creation operators

$$ \begin{align} \hat{a}^\dagger_k | n_0, n_1, ..., 0_k, n_{k+1} ... n_m \rangle &= \hat{a}^\dagger_k | 0_k, n_0, n_1, ..., n_{k+1}, ..., n_m \rangle \\ &= | 1_k, n_0, n_1, ..., n_{k+1}, ..., n_m \rangle \\ & = (-1)^{\sum_{l < k} n_l} | n_0, n_1, ..., 1_k, n_{k+1} ... n_m \rangle \end{align} $$

Given these definitions, it is not too difficult to see that the product of operators $\hat{a}^\dagger_k \hat{a}_k$ counts the occupation of the $k$th orbital, i.e.,

$$ \hat{a}^\dagger_k \hat{a}_k | n_0, n_1, ..., n_k, ..., n_m \rangle = n_k | n_0, n_1, ..., n_k, ..., n_{k+1}, ..., n_m \rangle $$

which suggests that we can define a number operator to count the total number of electrons in the state as

$$ \hat{N} = \sum_p \hat{a}^\dagger_p \hat{a}_p $$

The existence of this operator is our first clue that we should be able build useful quantum mechanical operators from sums / products of creation and annihilation operators.

The algebraic properties of Fermionic creation and annihilation operators encode the antisymmetry properties that should be satisfied by many-electron wave functions. Specifically, Fermionic creation and annihilation operators satisfy the following anticommutation relations

$$ \{\hat{a}_p, \hat{a}_q \} = \hat{a}_p \hat{a}_q + \hat{a}_q \hat{a}_p = 0$$

$$ \{\hat{a}^\dagger_p, \hat{a}^\dagger_q \} = \hat{a}^\dagger_p \hat{a}^\dagger_q + \hat{a}^\dagger_q \hat{a}^\dagger_p = 0$$

and

$$ \{\hat{a}^\dagger_p, \hat{a}_q \} = \hat{a}^\dagger_p \hat{a}_q + \hat{a}_q \hat{a}^\dagger_p = \delta_{pq}$$

where $\delta_{pq}$ is the Kronecker delta function.

Second-Quantized Operators¶

As suggested above, sums / products of creation and annihilation operators can be used to represent quantum mechanical operators.

One-Body Operators¶

Consider a one-body operator

$$ \mathcal{O}_1 = \sum_i^N \hat{h}({\bf x}_i) $$

where $i$ sums over all the $N$ electrons in the system and ${\bf x}_i$ represents the coordinates (spatial plus spin) of electron $i$. The second-quantized form of this operator, expanded within a basis of $m$ spin orbitals, $\{ \phi \}$, is

$$ \mathcal{O}_1 = \sum_{p}^m \sum_q^m \langle \phi_p({\bf x}_1) | \hat{h}({\bf x}_1) | \phi_q({\bf x}_1) \rangle \hat{a}^\dagger_p \hat{a}_q = \sum_{pq} h_{pq} \hat{a}^\dagger_p \hat{a}_q $$

Note that the sum in the original expression runs over the particles in the system, whereas the sums in the second-quantized form of the operator run over the spin-orbitals. Moving forward, it is assumed that the sums have these limits, and this notation is suppressed.

Two-Body Operators¶

Consider a general two-body operator

$$\mathcal{O}_2 = \frac{1}{2} \sum_{i \neq j} \hat{g}({\bf x}_i, {\bf x}_j) $$

The second-quantized form of this operator, expanded within a basis of $m$ spin orbitals, is

$$ \mathcal{O}_2 = \frac{1}{2} \sum_{pqrs} \langle \phi_p({\bf x}_1)\phi_q({\bf x}_2) | \hat{g}({\bf x}_i, {\bf x}_j) | \phi_r({\bf x}_1) \phi_s({\bf x}_2) \rangle \hat{a}_p^\dagger \hat{a}_q^\dagger \hat{a}_s \hat{a}_r = \frac{1}{2} \sum_{pqrs} g_{pqrs} \hat{a}_p^\dagger \hat{a}_q^\dagger \hat{a}_s \hat{a}_r $$

Note that the double sum in the original operator is restricted, whereas no restrictions appear in the second-quantized form of the operator.

Products of Operators¶

Occasionally one may encounter products of operators, and it is important to recognize that the second-quantized form of the product operator depends on whether the one takes the product of the original operators or their second-quantized forms. Consider two one-body operators, $\hat{A}$ and $\hat{B}$, defined as

$$ \begin{align} \hat{A} &= \sum_i \hat{A}({\bf x}_i) \\ &= \sum_{pq} A_{pq} \hat{a}^\dagger_p \hat{a}_q \end{align} $$

and

$$ \begin{align} \hat{B} &= \sum_i \hat{B}({\bf x}_i) \\ &= \sum_{pq} B_{pq} \hat{a}^\dagger_p \hat{a}_q \end{align} $$

The product of these operators contains both one-body and two-body operators

$$ \begin{align} \hat{A}\hat{B} &= \sum_{ij} \hat{A}({\bf x}_i) \hat{B}({\bf x}_j) \\ &= \sum_{i\neq j} \hat{A}({\bf x}_i) \hat{B}({\bf x}_j) + \sum_{i} \hat{A}({\bf x}_i) \hat{B}({\bf x}_i) \end{align} $$

Based on the definitions above, the second-quantized form of this operator would be

$$ \begin{align} \hat{A}\hat{B} &= \sum_{pqrs} \langle \phi_p({\bf x}_1)\phi_q({\bf x}_2) | \hat{A}({\bf x}_1)\hat{B}({\bf x}_2) | \phi_r({\bf x}_1) \phi_s({\bf x}_2) \rangle \hat{a}_p^\dagger \hat{a}_q^\dagger \hat{a}_s \hat{a}_r \\ &+ \sum_{pq} \langle \phi_p({\bf x}_1) | \hat{A}({\bf x}_1)\hat{B}({\bf x}_1)|\phi_q({\bf x}_1)\rangle \hat{a}^\dagger_p \hat{a}_q \\ &= \sum_{pqrs} \langle \phi_p({\bf x}_1) | \hat{A}({\bf x}_1) | \phi_r({\bf x}_1) \rangle \langle \phi_q({\bf x}_2) | \hat{B}({\bf x}_2) | \phi_s({\bf x}_2) \rangle \hat{a}_p^\dagger \hat{a}_q^\dagger \hat{a}_s \hat{a}_r \\ &+ \sum_{pq} (AB)_{pq} \hat{a}^\dagger_p \hat{a}_q \\ &= \sum_{pqrs} A_{pr} B_{qs} \hat{a}_p^\dagger \hat{a}_q^\dagger \hat{a}_s \hat{a}_r + \sum_{pq} (AB)_{pq} \hat{a}^\dagger_p \hat{a}_q \end{align} $$

where we have used the fact that the operator $\hat{A}({\bf x}_1)\hat{B}({\bf x}_2)$ is separable so we can perform the integrals over electron coordinates ${\bf x}_1$ and ${\bf x}_2$ separately.

If instead we take the product of the second-quantized form of these operators, we obtain a slightly different result

$$ \begin{align} \left ( \sum_{pr} A_{pr} \hat{a}^\dagger_p \hat{a}_r \right ) \left ( \sum_{qs} B_{qs} \hat{a}^\dagger_q \hat{a}_s \right ) &= \sum_{pqrs} A_{pr} B_{qs} \hat{a}^\dagger_p \hat{a}_r \hat{a}^\dagger_q \hat{a}_s \\ &= \sum_{pqrs} A_{pr} B_{qs} \hat{a}^\dagger_p (\delta_{qr} - \hat{a}^\dagger_q \hat{a}_r ) \hat{a}_s \\ &= -\sum_{pqrs} A_{pr} B_{qs} \hat{a}^\dagger_p \hat{a}^\dagger_q \hat{a}_r \hat{a}_s + \sum_{pqrs} \delta_{pr} A_{pr} B_{qs} \hat{a}^\dagger_p \hat{a}_s \\ &= \sum_{pqrs} A_{pr} B_{qs} \hat{a}^\dagger_p \hat{a}^\dagger_q \hat{a}_s \hat{a}_r + \sum_{pq} \left ( \sum_r A_{pr} B_{rq} \right ) \hat{a}^\dagger_p \hat{a}_q \end{align} $$

where, in the second and fourth lines, we have made use of two of the anticommutator identies given above. We can see that these two representations of $\hat{A}\hat{B}$ differ in the one-body part, and it turns out that these representations are only equivalent in the limit that the basis set $\{\phi\}$ is complete. This result makes sense because we can relate the two forms of the one-body part through the resolution of the identity, i.e.,

$$ \begin{align} \sum_r A_{pr} B_{rq} &= \sum_r \langle \phi_p({\bf x}_1)| \hat{A}({\bf x}_1) |\phi_r({\bf x}_1\rangle \langle \phi_r({\bf x}_1)| \hat{B}({\bf x}_1) | \phi_q({\bf x}_1)\rangle \\ & = \langle \phi_p({\bf x}_1)| \hat{A}({\bf x}_1) \hat{B}({\bf x}_1) | \phi_q({\bf x}_1)\rangle \\ & = (AB)_{pq} \end{align} $$

where we have used the fact that

$$\sum_r |\phi_r({\bf x}_1\rangle \langle \phi_r({\bf x}_1)| = \hat{1}$$

if the basis set is complete.

Second-Quantized Wave Functions¶

A Slater determinant for an $N$-electron system can be generated from the vacuum state by the repeated application of $N$ creation operators

$$|\Phi_0\rangle = \prod_i^N \hat{a}_i^{\dagger} |\rangle$$

More complicated many-electron wave functions can then be created by the application of additional creation and annihilation operators to this state. For example, consider the full configuration interaction (CI) wave function, which is a linear combination of all possible $N$-electron Slater determinants. We can generate the full CI wave function from the reference Slater determinant, $|\Phi_0\rangle$, by applying one-body, two-body, etc. transition operators to this state

$$ |\Psi_{\rm CI}\rangle = (1 + \hat{C}_1 + \hat{C}_2 + ... + \hat{C}_N)|\Phi_0\rangle$$

with

$$ \begin{align} \hat{C}_1 &= \sum_{ia} c_i^a \hat{a}^\dagger_a \hat{a}_i \\ \hat{C}_2 &= \frac{1}{(2!)^2} \sum_{ijab} c_{ij}^{ab} \hat{a}^\dagger_a \hat{a}^\dagger_b \hat{a}_j \hat{a}_i\\ \hat{C}_3 &= \frac{1}{(3!)^2} \sum_{ijkabc} c_{ijk}^{abc} \hat{a}^\dagger_a \hat{a}^\dagger_b \hat{a}^\dagger_c \hat{a}_k \hat{a}_j \hat{a}_i\\ ... \end{align} $$

Here, the labels $i$, $j$, $k$, ... refer to orbitals that are occupied in the reference determinant, $|\Phi_0\rangle$, and $a$, $b$, $c$, ... refer to unoccupied orbitals. So, the operator $\hat{C}_1$ is a single excitation operator, which generates Slater determinants where an electron in occupied orbital $\phi_i$ is moved to virtual orbital $\phi_a$. Similarly, $\hat{C}_2$ is a double excitation operator that generates configurations that are doubly substituted with respect to $|\Phi_0\rangle$. The coefficients, $c_i^a$, $c_{ij}^{ab}$, etc. are variable parameters that can be determined according to the variation theorem (i.e., by minimizing the electronic energy with respect to variations in the coefficients). Additional details regarding the properties of CI wave functions and how the CI coefficients can be determined in practice will be provided in another tutorial project.

Other types of wave functions can be defined in a similar way. For example, the coupled-cluster family of wave functions has the form

$$ |\Psi_{\rm CC}\rangle = \text{exp}(\hat{T})|\Phi_0\rangle$$

where $\hat{T} = \hat{T}_1 + \hat{T}_2 + ... \hat{T}_N$ is the cluster operator and

$$ \begin{align} \hat{T}_1 &= \sum_{ia} t_i^a \hat{a}^\dagger_a \hat{a}_i \\ \hat{T}_2 &= \frac{1}{(2!)^2} \sum_{ijab} t_{ij}^{ab} \hat{a}^\dagger_a \hat{a}^\dagger_b \hat{a}_j \hat{a}_i\\ \hat{T}_3 &= \frac{1}{(3!)^2} \sum_{ijkabc} t_{ijk}^{abc} \hat{a}^\dagger_a \hat{a}^\dagger_b \hat{a}^\dagger_c \hat{a}_k \hat{a}_j \hat{a}_i\\ ... \end{align} $$

The amplitudes $t_i^a$, $t_{ij}^{ab}$, etc. are the cluster amplitudes, which can be determined using a projective solution to the Schrödinger equation. Additional details regarding the properties of CC wave functions and how the cluster amplitudes can be determined in practice will be provided in another tutorial project.

Normal Order¶

The True Vacuum¶

One of the benefits of second quantization is that it simplifies the way in which we evaluate integrals involving quantum mechanical operators and many-body wave functions. It turns out that integrals over these operators are most conveniently dealt with when the operators are in "normal order," that is, when all of the creation operators lie to the left of the annihilation operators. Technically, this this definition of normal order is relative to the true vacuum state ($|\rangle$), whereas many electronic structure methods rely on a normal order relative to a different state, which is usually the Fermi vacuum. This case will be discussed in the next section.

For the true vacuum state, normal order is defined when the creation operators are all to the left of the annihilation operators. In this way, any expectation value of a string of normal-ordered operators (with respect to the vacuum state) will be zero. For example,

$$ \begin{align} \langle | \hat{a}^\dagger_p \hat{a}_q | \rangle &= 0 \\ \langle | \hat{a}_p \hat{a}_q | \rangle &= 0 \\ \langle | \hat{a}^\dagger_p \hat{a}^\dagger_q | \rangle &= 0 \end{align} $$

Note that the last integral vanishes, despite $\hat{a}^\dagger_p \hat{a}^\dagger_q$ not annihilating the ket state, $|\rangle$. The reason this term vanishes is that these operators will annihilate the bra state, $\langle |$, because

$$ \begin{align} \langle | \hat{a}^\dagger_p \hat{a}^\dagger_q = \left (\hat{a}_q \hat{a}_p|\rangle\right )^\dagger = 0 \\ \end{align} $$

Now, consider an operator that is not normal ordered, $\hat{a}_p \hat{a}^\dagger_q$. Briging this operator to normal order is easy; we simply apply the Fermionic anticommutation relations given above. We then find

$$ \langle | \hat{a}_p \hat{a}^\dagger_q | \rangle = \langle | \delta_{pq} | \rangle - \langle | \hat{a}^\dagger_p \hat{a}_q | \rangle = \delta_{pq} $$

The term that does not vanish is the one that does not include any operators. This is the power of the concept of normal order. If we would like to evaluate the expectation value of some second-quantized object, we simply need to bring that object to normal order. The only non-zero terms will be the "fully-contracted" ones that do not involve any operators.

Bringing more general strings of creation and annihilation operators to normal order is straightforward enough, given the anti-commutator relations obeyed by the operators. The challenge is really just that the application of these rules can become tedious. Thankfully, computers are well suited for tedious tasks. We have developed a tool (pdaggerq) that uses a simple Python interface to define and rearrange complicated strings of creation and annihilation operators. For example, the following Python code will bring the expression $\hat{a}_i^{\dagger}\hat{a}_j \hat{a}_l^{\dagger} \hat{a}_k$ to normal order.

In [1]:

import pdaggerq

pq = pdaggerq.pq_helper("true")

pq.add_operator_product(1.0, ['a*(i)','a(j)','a*(l)','a(k)'])
pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

['+1.00000000000000', 'i*', 'k', 'd(j,l)']
['-1.00000000000000', 'i*', 'l*', 'j', 'k']

This simple example can be verified by hand. Other more complicated strings can be much more annoying to rearrange.

The Fermi Vacuum¶

In correlated methods like CI or CC theory, normal order is not defined relative to the true vacuum state, but, rather, relative to the "Fermi vacuum," which is the Hartree-Fock (HF) state (a single Slater determinant). Recall that such a state can be generated by adding electrons to the true vacuum state

$$|\Phi_0\rangle = \prod_i^N \hat{a}_i^{\dagger} |\rangle$$

Now, given this reference state, rearranging operator strings for normal order relative to the true vacuum will become far too complicated because of the need to deal with all of the creation operators associated with constructing the HF state. The solution is to define normal order in a different way.

Our goal is to arrange our operators such that the expectation value of any term involving operators is zero, and the only non-zero terms are the fully-contracted ones involving no creation or annihilation operators. According to this rule, the usual creation or annihilation operators can act as either creators or annihilators on the Fermi vacuum, depending on the space (particle or hole) in which they act. So, $\hat{a}^{\dagger}_i$ and $\hat{a}_a$, should be treated as annihilation operators with respect to the Fermi vacuum because

$$ \begin{align} \hat{a}^{\dagger}_i |\Phi_0\rangle &= 0 \\ \hat{a}_a |\Phi_0\rangle &= 0 \end{align} $$

while $\hat{a}_j$ and $\hat{a}^{\dagger}_b$ are treated as creation operators with respect to $|\Phi_0\rangle$ because

$$ \begin{align} \hat{a}_j |\Phi_0\rangle &\neq 0\\ \hat{a}^{\dagger}_b |\Phi_0\rangle &\neq 0 \end{align} $$

To bring a string of operators to normal order relative to the Fermi vacuum, all $\hat{a}^{\dagger}_i$ and $\hat{a}_a$ must lie to the right of all $\hat{a}_j$ and $\hat{a}^{\dagger}_b$. In this way, the expectation value of any normal-ordered string of operators will be zero. For example

$$ \begin{align} \langle \Phi_0 | \hat{a}^{\dagger}_i |\Phi_0\rangle &= 0 \\ \langle \Phi_0 | \hat{a}_a |\Phi_0\rangle &= 0 \\ \langle \Phi_0 | \hat{a}_j |\Phi_0\rangle &= 0 \\ \langle \Phi_0 | \hat{a}^{\dagger}_b |\Phi_0\rangle &= 0 \end{align} $$

Note that the last two integrals vanish, despite $\hat{a}_j$ and $\hat{a}^{\dagger}_b$ not annihilating the ket state, $|\Phi_0\rangle$. The reason is that these operators will annihilate the bra state, $\langle \Phi_0 |$, because

$$ \begin{align} \langle \Phi_0 | \hat{a}_j = \left (\hat{a}^\dagger_j|\Phi_0\rangle\right )^\dagger = 0 \\ \langle \Phi_0 | \hat{a}^{\dagger}_b = \left (\hat{a}_b|\Phi_0\rangle\right )^\dagger = 0 \end{align} $$

These results are generalizable to longer strings of operators. For example, the following normal-ordered products of two operators also vanish:

$$ \begin{align} \langle \Phi_0 | \hat{a}_j \hat{a}^\dagger_i |\Phi_0\rangle = 0 \\ \langle \Phi_0 | \hat{a}^\dagger_i \hat{a}_a |\Phi_0\rangle = 0 \\ \langle \Phi_0 | \hat{a}^\dagger_a \hat{a}_i |\Phi_0\rangle = 0 \\ \langle \Phi_0 | \hat{a}^\dagger_a \hat{a}_b |\Phi_0\rangle = 0 \\ \end{align} $$

Again, an expectation value over any normal-ordered string of creation and annihilation operators will be zero, regardless of how many operators appear in the string.

Let us use what we have learned so far to evaluate the expectation value of the one-body operator

$$\hat{\mathcal{O}}_1 = \sum_{pq} h_{pq} \hat{a}^\dagger_p \hat{a}_q$$

with respect to the Hartree-Fock state, $|\Phi_0\rangle$. In order to bring this operator to normal order with respect to the Fermi vacuum, we need to split each sum over general orbitals into separate sums over occupied and virtual orbitals. In this way, we can more clearly identify which operators annihilate the Fermi vacuum. We have

$$\hat{\mathcal{O}}_1 = \sum_{ij} h_{ij} \hat{a}^\dagger_i \hat{a}_j + \sum_{ia} h_{ia} \hat{a}^\dagger_i \hat{a}_a + \sum_{ai} h_{ai} \hat{a}^\dagger_a \hat{a}_i + \sum_{ab} h_{ia} \hat{a}^\dagger_a \hat{a}_b $$

Here, only the first sum involves an operator that is not in normal order. Bringing it to normal order $(\hat{a}^\dagger_i \hat{a}_j = \delta_{ij} - \hat{a}_j \hat{a}^\dagger_i)$, we have

$$ \hat{\mathcal{O}}_1 = \sum_{i} h_{ii} - \sum_{ij} h_{ij} \hat{a}_j \hat{a}^\dagger_i + \sum_{ia} h_{ia} \hat{a}^\dagger_i \hat{a}_a + \sum_{ai} h_{ai} \hat{a}^\dagger_a \hat{a}_i + \sum_{ab} h_{ia} \hat{a}^\dagger_a \hat{a}_b $$

and the expectation value of $\hat{\mathcal{O}_1}$ is

$$ \langle \Phi_0 | \hat{\mathcal{O}}_1 |\Phi_0\rangle = \sum_{i} h_{ii} $$

because, as shown above, all of the expectation values involving the normal-ordered operators vanish. Note that this result is consistent with the Slater-Condon rules, as we should expect! Just to be sure, let's check if we can get the same result from the pdaggerq package. We can use the built-in operator 'h' to represent the one-body operator here.

In [2]:

pq = pdaggerq.pq_helper("fermi")

pq.add_operator_product(1.0, ['h'])
pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

['+1.00000000000000', 'h(i,i)']

Let us try a slightly more complicated example. We will use pdaggerq to evaluate the Hartree-Fock energy, which is the expectation value of the Hamiltonian with respect to the Hartree-Fock state. In second-quantization, the Hamiltonian has the form

$$ \hat{H} = \sum_{pq} h_{pq} \hat{a}^\dagger_p \hat{a}_q + \frac{1}{2} \sum_{pqrs} (pr|qs) \hat{a}^\dagger_p \hat{a}^\dagger_q \hat{a}_s \hat{a}_r $$

where $h_{pq}$ is the matrix representation of the core Hamiltonian (kinetic energy plus electron-nucleus potential energy operators), and $(pr|qs)$ is an electron repulsion integral in chemists' notation. Two-body operators in pdaggerq are assumed to be antisymmetrized, so we should instead use

$$ \hat{H} = \sum_{pq} h_{pq} \hat{a}^\dagger_p \hat{a}_q + \frac{1}{4} \sum_{pqrs} \langle pq || rs \rangle \hat{a}^\dagger_p \hat{a}^\dagger_q \hat{a}_s \hat{a}_r $$

Here, $\langle pq || rs \rangle = \langle pq | rs \rangle - \langle pq | sr \rangle$, where $\langle pq | rs \rangle = (pr|qs)$ is an electron repulsion integral in physicists' notation. The Hartree-Fock energy (evaluated in the basis of orthonormal molecular orbitals) is

$$ E_0 = \langle \Phi_0 | \hat{H} | \Phi_0 \rangle = \sum_{i} h_{ii} + \frac{1}{2} \sum_{ij} \langle ij || ij \rangle $$

In pdaggerq, we can use the built-in operators 'h' and 'g' to represent the core Hamiltonian operator and an antisymmetrized two-electron operator.

In [3]:

pq = pdaggerq.pq_helper("fermi")

pq.add_operator_product(1.0, ['h'])
pq.add_operator_product(0.25, ['g'])

pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

['+1.00000000000000', 'h(i,i)']
['+0.50000000000000', 'g(j,i,j,i)']

Indeed, this is the correct form for the Hartree-Fock energy! We can also use a slightly different representation of the Hamiltonian,

$$ \hat{H} = \hat{f} + \hat{v} $$

where $\hat{f}$ is the Fock operator and $\hat{v} = \hat{H} - \hat{f}$ is the fluctuation potential operator.

In [4]:

pq = pdaggerq.pq_helper("fermi")

pq.add_operator_product(1.0, ['f'])
pq.add_operator_product(1.0, ['v'])

pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

['+1.00000000000000', 'f(i,i)']
['-0.50000000000000', '<j,i||j,i>']

This result is equivalent to the one above because, in the MO basis,

$$ f_{ii} = \sum_{ii} h_{ii} + \sum_{ij} \langle ji || ji \rangle$$

$$ \sum_{i} f_{ii} - \frac{1}{2} \sum_{ij} \langle ji||ji \rangle = \sum_{i} h_{ii} + \frac{1}{2} \sum_{ij} \langle ji||ji \rangle$$

Now, let's consider one last, slightly more complicated example. We will use pdaggerq to evaluate the CC with single and double excitations (CCSD) energy expression:

$$E_{\rm {CC}} = \langle \Phi_0|\text{exp}(-\hat{T}) \hat{H} \text{exp}(\hat{T})|\Phi_0 \rangle$$

Here, $\hat{T}=\hat{T}_1+\hat{T}_2$ represents the cluster operator, and $\text{exp}(-\hat{T}) \hat{H} \text{exp}(\hat{T})$ is a similarity transformation of the Hamiltonian. This similarity transformation can be carried out using commutators according to the Baker Campbell Hausdorff (BCH) expansion of the Hamiltonian, which truncates after four nested commutators

$$ \hat{H} = \hat{H} + [\hat{H}, \hat{T}] + \frac{1}{2!} [[\hat{H}, \hat{T}], \hat{T}] + \frac{1}{3!} [[[\hat{H}, \hat{T}], \hat{T}], \hat{T}] + \frac{1}{4!} [[[[\hat{H}, \hat{T}], \hat{T}], \hat{T}], \hat{T}]$$

The following code makes use of the built-in operator string types in pdaggerq ('t1', 't2', etc.), as well as commutator and double commutator functions, to evaluate the CCSD energy. Note that, for the energy expression, we only need to consider the BCH expansion up to double commutators. For specific definitions of these operators and functions, see here.

In [5]:

import pdaggerq

pq = pdaggerq.pq_helper('fermi')

pq.set_print_level(0)

print('')
print('# < 0 | e(-T) H e(T) | 0> :')
print('')

# H = f + v (fock operator plus fluctuation potential)

# one-electron part: 

# f
pq.add_operator_product(1.0,['f'])

# [f, T1]
pq.add_commutator(1.0,['f'],['t1'])

# [f, T2]
pq.add_commutator(1.0,['f'],['t2'])

# two-electron part: 

# v
pq.add_operator_product(1.0,['v'])

# [v, T1]
pq.add_commutator(1.0,['v'],['t1'])

# [v, T2]
pq.add_commutator(1.0,['v'],['t2'])

# [[v, T1], T1]]
pq.add_double_commutator(0.5, ['v'],['t1'],['t1'])

pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

pq.clear()

# < 0 | e(-T) H e(T) | 0> :

['+1.00000000000000', 'f(i,i)']
['+1.00000000000000', 'f(i,a)', 't1(a,i)']
['-0.50000000000000', '<j,i||j,i>']
['+0.25000000000000', '<j,i||a,b>', 't2(a,b,j,i)']
['-0.50000000000000', '<j,i||a,b>', 't1(a,i)', 't1(b,j)']

We can see that the expectation value of the similarity-transformed Hamiltonian matrix can be broken down as

$$ E_{\rm CC} = E_0 + E_{\rm c}$$

where $E_0$ is the Hartree-Fock energy discussed above, and $E_{\rm c}$ is a correlation contribution defined as

$$ E_{\rm c} = \sum_{ia} f_{ia} t^a_i + \frac{1}{4}\sum_{ijab} \langle ij||ab\rangle \left ( t_{ij}^{ab} + 2 t^a_i t^b_j \right ) $$

Note that pdaggerq also has a built-in "add_st_operator" function that will automatically generate the nested commutators that arise in the BCH expansion (up to four nested commutators), which greatly simplifies the code given above:

In [6]:

import pdaggerq

pq = pdaggerq.pq_helper('fermi')

pq.set_print_level(0)

print('')
print('# < 0 | e(-T) H e(T) | 0> :')
print('')

# H = f + v (fock operator plus fluctuation potential)

pq.add_st_operator(1.0,['f'],['t1','t2'])
pq.add_st_operator(1.0,['v'],['t1','t2'])

pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)
    
pq.clear()

# < 0 | e(-T) H e(T) | 0> :

['+1.00000000000000', 'f(i,i)']
['+1.00000000000000', 'f(i,a)', 't1(a,i)']
['-0.50000000000000', '<j,i||j,i>']
['+0.25000000000000', '<j,i||a,b>', 't2(a,b,j,i)']
['-0.50000000000000', '<j,i||a,b>', 't1(a,i)', 't1(b,j)']

We will explore more of the functionality in pdaggerq in the tutorials on CI and CC theory that follow this one.