from_cirq for non-singlets

The following code is used to test the conversion for a doublet state. It fails to generate an FQE wavefunction with the correct symmetry or any amplitudes.

import openfermion as of
import numpy as np
from functools import reduce
import fqe

def lowdin_projector(op, target_eig, other_eigs):
    identity = np.eye(op.shape[0])
    projector_components = [(op - identity * oe) / (target_eig - oe) for oe in other_eigs]
    return reduce(np.dot, projector_components)

if __name__ == "__main__":
    sites = 5
    na, nb = 3, 2
    U = 4
    hubbard = of.hamiltonians.fermi_hubbard(1, sites, tunneling=1, coulomb=U,
                                            chemical_potential=0,
                                            magnetic_field=0,
                                            periodic=True,
                                            spinless=False)

    # get operators on the full Fock space
    hamiltonian = of.get_interaction_operator(hubbard)
    h_mat = of.get_sparse_operator(hamiltonian).toarray().real
    s2 = of.get_sparse_operator(of.s_squared_operator(sites))
    sz = of.get_sparse_operator(of.sz_operator(sites))
    num_op = of.get_sparse_operator(of.number_operator(2 * sites))
    # get all possible quantum numbers supported on the Fock space
    particle_num_vals = list(range(2 * sites + 1))
    del particle_num_vals[particle_num_vals.index(sites)]
    print(particle_num_vals)  # should be all particles except 5
    sz_num_vals = list(range(-sites, sites + 1, 1))
    del sz_num_vals[sz_num_vals.index(1)]
    sz_num_vals = [xx / 2 for xx in sz_num_vals]
    print(sz_num_vals) # should be all Sz except 0.5
    s_s2_vals = np.arange(sites + 1) / 2
    s2_num_vals = [s * (s + 1) for s in s_s2_vals]
    del s2_num_vals[s2_num_vals.index(0.75)]
    print(s2_num_vals)  # should [0, 2, 3.75, 6, 8.75]

    # get projectors with the Lowdin projection formula..
    # aka lagrange interpolation
    proj_n = lowdin_projector(num_op, sites, particle_num_vals)
    proj_sz = lowdin_projector(sz, 0.5, sz_num_vals)
    proj_s2 = lowdin_projector(s2, 0.75, s2_num_vals)

    # get operators for comparing local expectation value
    # and RDM elements
    local_na = [of.number_operator(2 * sites, 2 * x) for x in range(sites)]
    local_nb = [of.number_operator(2 * sites, 2 * x + 1) for x in range(sites)]
    na_site = [of.get_sparse_operator(op, n_qubits=2 * sites) for op in local_na]
    nb_site = [of.get_sparse_operator(op, n_qubits=2 * sites) for op in local_nb]

    # solve the Schrodinger equation
    # projected into the sector I care about Sz = 0.5, S = 0.5 <S^2> = 0.75
    w_full, v_full = np.linalg.eigh(h_mat)
    w_n, v_n = np.linalg.eigh(proj_s2 @ proj_sz @ proj_n @ h_mat @ proj_n @ proj_sz @ proj_sz)

    print("N exp ", v_n[:, [0]].conj().T @ num_op @ v_n[:, [0]])
    print("Sz exp ", v_n[:, [0]].conj().T @ sz @ v_n[:, [0]])
    print("Sz exp ", v_n[:, [0]].conj().T @ s2 @ v_n[:, [0]])
    wfn = fqe.from_cirq(np.array(v_n[:, [0]]).flatten(), 1.0E-12)
    wfn.print_wfn()  # print's empty

The FQE.Wavefunction thinks this sector is Sz=-1 which is opposite of OpenFermion. It also doesn't seem to think there are any coefficients in that sector.

quantumlib / OpenFermion-FQE

from_cirq for non-singlets #40