Automatically adding the "missing" amplitudes from Psi4. about openfermion-psi4 HOT 7 OPEN

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