Automatically adding the "missing" amplitudes from Psi4.

[bsenjean]

Hi,

When I tried to recover the CCSD energy based on the CCSD amplitudes extracted from OpenFermion-Psi4, I realized that some were missing: the same-spin amplitudes.

See here for a question I asked on the Psi4 forum: http://forum.psicode.org/t/understanding-mp2-ccsd-amplitudes-for-h4/1866

Also this is not an issue as they can be determined later on, I think it would actually be better to extract those directly from the OpenFermion-Psi4 interface, because I'm pretty sure most user will not know that they are missing. And CCSD amplitudes could be used in some UCC ansatz.

[ncrubin]

I would assume this problem persists with ROHF right? If so then yes we should extract the amplitudes for the users!

[bsenjean]

The problem is solved when using ROHF (which in OpenFermion-Psi4 is used when multiplicity is not equal to 1).
But this should still be fixed when using RHF right ?

[ncrubin]

Hey @bsenjean, sorry getting back to this. I too ran into something similar where I learned the CCSD closed shell amplitudes are S2 and Sz spin adapted (from a very old paper in the 80's which everyone cites). Therefore, as psi4 forum mentions same-spin t2 are defined as T_IJAB = T_IjAb - T_IjBa = T_ijab. Is that what you mean?

[bsenjean]

Yes, typically these amplitudes are required and computed inside Psi4 but they do not print them as they can be determined from the others by the expression T_IJAB = T_IjAb - T_IjBa = T_ijab. So one should also extract them in OpenFermion otherwise you'll see that you don't get the correct CCSD energy.

One of our student made a subroutine to add them to the ccsd operator in OpenFermion, see below (it could certainly be more compacted):

def add_same_spin_amplitudes(
        ccsd_op: openfermion.FermionOperator,
        n_electrons: int,
        n_qubits: int) -> openfermion.FermionOperator:
    """
    Add same-spin amplitudes to a normal ordered ccsd operator.

    Largest TIJAB Amplitudes:
          1   0   1   0        -0.0060909463
    Largest Tijab Amplitudes:
          1   0   1   0        -0.0060909463
    they are equal to T_IJAB = T_ijab = T_IjAb - T_IjBa
    (see
    http://forum.psicode.org/t/understanding-mp2-ccsd-amplitudes-for-h4/1866/3)
    In normal-ordered form: (if not normal order, add the symmetric term with
    an appropriate factor 1/2 in the amplitudes)

    Args:
        ucc_operator (FermionOperator): Coupled-cluster operator.
        n_electrons (int): Number of electrons.
        n_qubits (int): Number of qubits.

    Returns:
        ccsd_op(openfermin.FermionOperator): normal ordered ccsd operator.
    """
    if not isinstance(ccsd_op, openfermion.FermionOperator):
        raise TypeError('CC operator must be FermionOperator.')

    ccsd_op = openfermion.normal_ordered(ccsd_op)
    for i in range(0, n_electrons, 2):
        for j in range(i + 2, n_electrons, 2):
            for a in range(n_electrons, n_qubits, 2):
                for b in range(a + 2, n_qubits, 2):
                    # Make sure that T_IjAb is present in ccsd_op
                    if ((b, 1), (a + 1, 1), (j, 0),
                            (i + 1, 0)) not in list(ccsd_op.terms.keys()):
                        # if one is not present then both are not
                        # (ie if T_IjAb=0 then T_IjBa=0)
                        if (((b, 1), (a + 1, 1), (j, 0), (i + 1, 0)) or
                                ((b + 1, 1), (a, 1), (j, 0),
                                    (i + 1, 0))) in list(ccsd_op.terms.keys()):
                            print(
                                'CCSD operator incomplete, missing: ', ((
                                    b + 1, 1), (a, 1), (j, 0), (i + 1, 0)))
                    else:
                        add_op = openfermion.FermionOperator(
                            '{}^ {}^ {} {}'.format(b, a, j, i), ccsd_op.terms[(
                                (b, 1), (a + 1, 1), (j, 0), (i + 1, 0))] +
                            ccsd_op.terms[
                                ((b + 1, 1), (a, 1), (j, 0), (i + 1, 0))])
                        if list(add_op.terms.keys())[0] not in list(
                                ccsd_op.terms.keys()):
                            ccsd_op += add_op
                        add_op = openfermion.FermionOperator(
                            '{}^ {}^ {} {}'.format(b + 1, a + 1, j + 1, i + 1),
                            ccsd_op.terms[(
                                (b, 1), (a + 1, 1), (j, 0), (i + 1, 0))] +
                            ccsd_op.terms[
                                ((b + 1, 1), (a, 1), (j, 0), (i + 1, 0))])
                        if list(add_op.terms.keys())[0] not in list(
                                ccsd_op.terms.keys()):
                            ccsd_op += add_op
    return ccsd_op

[ncrubin]

Great. Do you or the student want to submit this as a PR? It would be nice to have so others know how to "unspin-adapt" the amplitudes.

[bsenjean]

Sure, I guess a good place would be openfermionpsi4/_psi4_conversion_functions.py ?

I could either modify directly parse_psi4_ccsd_amplitudes and add the same-spin amplitudes if restricted == True, or just create another function. In my opinion, it is more a correction than an addition and should be in parse_psi4_ccsd_amplitudes.

[bsenjean]

Alright, I added the pull request. I tried to minimize the coding as much as possible, i.e. by only adding one dictionary to store the amplitudes that I needed to compute the same-spin amplitudes for restricted/closed-shell calculations.
I realized I needed an additional factor 1/2, I'm not entirely sure where it comes from.

I checked the CCSD energy resulting from these new amplitudes for H4 and LiH, and it worked well (I recover the correct energy up to single precision).