In this project, you will learn how to use the pdaggerq project to bring complicated strings of second-quantized operators to normal order, relative to both the true vacuum state and the Fermi vacuum. You will generate singles and doubles amplitudes equations for coupled cluster with single and double excitations (CCSD). You can also follow a video version of this tutorial, here, but note that much of the syntax used in pdaggerq has changed since that video was made:





Background:

Many-body quantum chemistry methods are often formulated using second-quantized creation and annihiliation operators. Second-quantized operators give us compact representations of many-body wave functions, and the antisymmetry properties of fermionic wave functions are cooked into the anti commutation relations for the operators. We have





and



where and represent creation and annihilation operators, respectively. Oftentimes, we encounter products of creation and annihilation operators that are more conveniently dealt with when in so-called "normal order," that is, when all of the creation operators lie to the left of the annihilation operators (note that this definition of normal order is relative to the true vacuum state, whereas normal order relative to the Fermi vacuum is slightly different and will be discussed below). Bringing a string of 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 allows one to simply define and rearrange complicated strings of creation and annihilation operators within a Python interface. For example, the following Python code will evaluate the expression

    import pdaggerq
   
    pq = pdaggerq.pq_helper("true")

    pq.set_print_level(1)

    pq.set_string(['i*','j','l*','k'])
    pq.add_new_string()

    pq.simplify()
    pq.print()

which when executed prints

    // starting string:
    //     + 1.00000 i* j l* k 

    // normal-ordered strings:
    //     + 1.00000 i* k d(jl) 
    //     - 1.00000 i* l* j k 

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


Coupled-cluster theory and the Fermi vacuum:

In correlated methods like coupled cluster (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. The HF state can be generated from the true vacuum state as



Now, given this reference state, rearranging operator strings 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.

For the vacuum state, normal order is defined when the annihilation operators are all to the right of the creation operators. In this way, any expectation value of a string of normal-ordered operators will be zero. For the Fermi vacuum, normal order is defined in a different, yet similar way. We would like to arrange our operators such that any expectation value involving operators is zero, and we are left only with "fully contracted" terms 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, and , are treated as annihilation operators (i.e., they annihilate the HF reference), while and are treated as creation operators. To bring a string of operators to normal order relative to the Fermi vacuum, all and must lie to the right of all and .

Now, we consider the coupled cluster with single and double excitations (CCSD) energy expression, which should be evaluated with normal-order defined relative to the Fermi vacuum:



Here, T=T1+T2 and H represent the cluster operator and Hamiltonian, respectively, and the similarity transformation should be carried out using commutators according to the Baker Campbell Hausdorff (BCH) expansion. The following code makes use of the built-in operator string types in pdaggerq (i.e., t1, t2, etc.), as well as commutator and double commutator functions, to evaluate the CCSD energy. Note that the BCH expansion truncates after only double commutators in this case. For specific definitions of these operators and functions, see here.

    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()

    pq.print_fully_contracted()

    pq.clear()

The output of this code will be

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


    // fully-contracted strings:
    //     +     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)

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:

    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()
    pq.print_fully_contracted()
    pq.clear()

and this code should produce the same output as above.

Now, your task is to evaluate the CCSD singles and doubles equations according to



and



You should initialize pdaggerq relative to the Fermi vacuum and use the built-in "set_left_operators" function to indicate that your bra state is or . For example, for the singles equations, you would begin with

    import pdaggerq

    pq = pdaggerq.pq_helper('fermi')

    pq.set_left_operators([['e1(m,e)']])

Working python scripts for this exercise can be found here: singles, doubles.