<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Copyright 2020 The OpenFermion Developers
Step2:  FQE vs OpenFermion vs Cirq
Step3:  The first example we will perform is diagonal Coulomb evolution on the Hartree-Fock state.  The diagonal Coulomb operator is defined as
\begin{align}
V = \sum_{\alpha, \beta \in {\uparrow, \downarrow}}\sum_{p,q} V_{pq,pq}n_{p,\alpha}n_{q,\beta}
\end{align}
The number of free parpameters are $\mathcal{O}(N^{2})$ where $N$ is the rank of the spatial basis. The DiagonalCoulomb Hamiltonian takes either a generic 4-index tensor or the $N \times N$ matrix defining $V$.  If the 4-index tensor is given the $N \times N$ matrix is constructed along with the diagonal correction.  If the goal is to just evolve under $V$ it is recommended the user input the $N \times N$ matrix directly.
All the terms in $V$ commute and thus we can evolve under $V$ exactly by counting the accumulated phase on each bitstring.
To start out let's define a Hartree-Fock wavefunction for 4-orbitals and 2-electrons $S_{z} =0$.
Step4:  Now we can define a random 2-electron operator $V$.  To define $V$ we need a $4 \times 4$ matrix.  We will generate this matrix by making a full random two-electron integral matrix and then just take the diagonal elements
Step5:  Evolution under $V$ can be computed by looking at each bitstring, seeing if $n_{p\alpha}n_{q\beta}$ is non-zero and then phasing that string by $V_{pq}$.  For the Hartree-Fock state we can easily calculate this phase accumulation.  The alpha and beta bitstrings are "0001" and "0001".
Step6:  We can now try this out for more than 2 electrons.  Let's reinitialize a wavefunction on 6-orbitals with 4-electrons $S_{z} = 0$ to a random state.
Step7:  We need to build our Diagoanl Coulomb operator For this bigger system.
Step8:  Now we can convert our wavefunction to a cirq wavefunction, evolve under the diagonal_coulomb operator we constructed and then compare the outputs.
Step9:  Finally, we can compare against evolving each term of $V$ individually.
<ASSISTANT_TASK:>
Python Code:
#@title Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
Explanation: Copyright 2020 The OpenFermion Developers
End of explanation
try:
    import fqe
except ImportError:
    !pip install fqe --quiet
from itertools import product
import fqe
from fqe.hamiltonians.diagonal_coulomb import DiagonalCoulomb
import numpy as np
import openfermion as of
from scipy.linalg import expm
#Utility function
def uncompress_tei(tei_mat, notation='chemistry'):
    
    uncompress chemist notation integrals
    tei_tensor[i, k, j, l] = tei_mat[(i, j), (k, l)]
    [1, 1, 2, 2] = [1, 1, 2, 2] = [1, 1, 2, 2]  = [1, 1, 2, 2]
    [i, j, k, l] = [k, l, i, j] = [j, i, l, k]* = [l, k, j, i]*
    For real we also have swap of i <> j and k <> l
    [j, i, k, l] = [l, k, i, j] = [i, j, l, k] = [k, l, j, i]
    tei_mat[(i, j), (k, l)] = int dr1 int dr2 phi_i(dr1) phi_j(dr1) O(r12) phi_k(dr1) phi_l(dr1)
    Physics notation is the notation that is used in FQE.
    Args:
        tei_mat: compressed two electron integral matrix
    Returns:
        uncompressed 4-electron integral tensor. No antisymmetry.
    
    if notation not in ['chemistry', 'physics']:
        return ValueError("notation can be [chemistry, physics]")
    norbs = int(0.5 * (np.sqrt(8 * tei_mat.shape[0] + 1) - 1))
    basis = {}
    cnt = 0
    for i, j in product(range(norbs), repeat=2):
        if i >= j:
            basis[(i, j)] = cnt
            cnt += 1
    tei_tensor = np.zeros((norbs, norbs, norbs, norbs))
    for i, j, k, l in product(range(norbs), repeat=4):
        if i >= j and k >= l:
            tei_tensor[i, j, k, l] = tei_mat[basis[(i, j)], basis[(k, l)]]
            tei_tensor[k, l, i, j] = tei_mat[basis[(i, j)], basis[(k, l)]]
            tei_tensor[j, i, l, k] = tei_mat[basis[(i, j)], basis[(k, l)]]
            tei_tensor[l, k, j, i] = tei_mat[basis[(i, j)], basis[(k, l)]]
            tei_tensor[j, i, k, l] = tei_mat[basis[(i, j)], basis[(k, l)]]
            tei_tensor[l, k, i, j] = tei_mat[basis[(i, j)], basis[(k, l)]]
            tei_tensor[i, j, l, k] = tei_mat[basis[(i, j)], basis[(k, l)]]
            tei_tensor[k, l, j, i] = tei_mat[basis[(i, j)], basis[(k, l)]]
    if notation == 'chemistry':
        return tei_tensor
    elif notation == 'physics':
        return np.asarray(tei_tensor.transpose(0, 2, 1, 3), order='C')
    return tei_tensor
Explanation: FQE vs OpenFermion vs Cirq: Diagonal Coulomb Operators
<table class="tfo-notebook-buttons" align="left">
  <td>
    <a target="_blank" href="https://quantumai.google/openfermion/fqe/tutorials/diagonal_coulomb_evolution"><img src="https://quantumai.google/site-assets/images/buttons/quantumai_logo_1x.png" />View on QuantumAI</a>
  </td>
  <td>
    <a target="_blank" href="https://colab.research.google.com/github/quantumlib/OpenFermion/blob/master/docs/fqe/tutorials/diagonal_coulomb_evolution.ipynb"><img src="https://quantumai.google/site-assets/images/buttons/colab_logo_1x.png" />Run in Google Colab</a>
  </td>
  <td>
    <a target="_blank" href="https://github.com/quantumlib/OpenFermion/blob/master/docs/fqe/tutorials/diagonal_coulomb_evolution.ipynb"><img src="https://quantumai.google/site-assets/images/buttons/github_logo_1x.png" />View source on GitHub</a>
  </td>
  <td>
    <a href="https://storage.googleapis.com/tensorflow_docs/OpenFermion/docs/fqe/tutorials/diagonal_coulomb_evolution.ipynb"><img src="https://quantumai.google/site-assets/images/buttons/download_icon_1x.png" />Download notebook</a>
  </td>
</table>
Special routines are available for evolving under a diagonal Coulomb operator.  This notebook describes how to use these built in routines and how they work.
End of explanation
norbs = 4
tedim = norbs * (norbs + 1) // 2
if (norbs // 2) % 2 == 0:
    n_elec = norbs // 2
else:
    n_elec = (norbs // 2) + 1
sz = 0
fqe_wfn = fqe.Wavefunction([[n_elec, sz, norbs]])
fci_data = fqe_wfn.sector((n_elec, sz))
fci_graph = fci_data.get_fcigraph()
hf_wf = np.zeros((fci_data.lena(), fci_data.lenb()), dtype=np.complex128)
hf_wf[0, 0] = 1  # right most bit is zero orbital.
fqe_wfn.set_wfn(strategy='from_data',
                raw_data={(n_elec, sz): hf_wf})
fqe_wfn.print_wfn()
Explanation: The first example we will perform is diagonal Coulomb evolution on the Hartree-Fock state.  The diagonal Coulomb operator is defined as
\begin{align}
V = \sum_{\alpha, \beta \in {\uparrow, \downarrow}}\sum_{p,q} V_{pq,pq}n_{p,\alpha}n_{q,\beta}
\end{align}
The number of free parpameters are $\mathcal{O}(N^{2})$ where $N$ is the rank of the spatial basis. The DiagonalCoulomb Hamiltonian takes either a generic 4-index tensor or the $N \times N$ matrix defining $V$.  If the 4-index tensor is given the $N \times N$ matrix is constructed along with the diagonal correction.  If the goal is to just evolve under $V$ it is recommended the user input the $N \times N$ matrix directly.
All the terms in $V$ commute and thus we can evolve under $V$ exactly by counting the accumulated phase on each bitstring.
To start out let's define a Hartree-Fock wavefunction for 4-orbitals and 2-electrons $S_{z} =0$.
End of explanation
tei_compressed = np.random.randn(tedim**2).reshape((tedim, tedim))
tei_compressed = 0.5 * (tei_compressed + tei_compressed.T)
tei_tensor = uncompress_tei(tei_compressed, notation='physics')
diagonal_coulomb = of.FermionOperator()
diagonal_coulomb_mat = np.zeros((norbs, norbs))
for i, j in product(range(norbs), repeat=2):
    diagonal_coulomb_mat[i, j] = tei_tensor[i, j, i, j]
    for sigma, tau in product(range(2), repeat=2):
        diagonal_coulomb += of.FermionOperator(
            ((2 * i + sigma, 1), (2 * i + sigma, 0), (2 * j + tau, 1),
             (2 * j + tau, 0)), coefficient=diagonal_coulomb_mat[i, j])
dc_ham = DiagonalCoulomb(diagonal_coulomb_mat)
Explanation: Now we can define a random 2-electron operator $V$.  To define $V$ we need a $4 \times 4$ matrix.  We will generate this matrix by making a full random two-electron integral matrix and then just take the diagonal elements
End of explanation
alpha_occs = [list(range(fci_graph.nalpha()))]
beta_occs = [list(range(fci_graph.nbeta()))]
occs = alpha_occs[0] + beta_occs[0]
diag_ele = 0.
for ind in occs:
    for jnd in occs:
        diag_ele += diagonal_coulomb_mat[ind, jnd]
evolved_phase = np.exp(-1j * diag_ele)
print(evolved_phase)
# evolve FQE wavefunction
evolved_hf_wfn = fqe_wfn.time_evolve(1, dc_ham)
# check they the accumulated phase is equivalent!
assert np.isclose(evolved_hf_wfn.get_coeff((n_elec, sz))[0, 0], evolved_phase)
Explanation: Evolution under $V$ can be computed by looking at each bitstring, seeing if $n_{p\alpha}n_{q\beta}$ is non-zero and then phasing that string by $V_{pq}$.  For the Hartree-Fock state we can easily calculate this phase accumulation.  The alpha and beta bitstrings are "0001" and "0001".
End of explanation
norbs = 6
tedim = norbs * (norbs + 1) // 2
if (norbs // 2) % 2 == 0:
    n_elec = norbs // 2
else:
    n_elec = (norbs // 2) + 1
sz = 0
fqe_wfn = fqe.Wavefunction([[n_elec, sz, norbs]])
fqe_wfn.set_wfn(strategy='random')
inital_coeffs = fqe_wfn.get_coeff((n_elec, sz)).copy()
print("Random initial wavefunction")
fqe_wfn.print_wfn()
Explanation: We can now try this out for more than 2 electrons.  Let's reinitialize a wavefunction on 6-orbitals with 4-electrons $S_{z} = 0$ to a random state.
End of explanation
tei_compressed = np.random.randn(tedim**2).reshape((tedim, tedim))
tei_compressed = 0.5 * (tei_compressed + tei_compressed.T)
tei_tensor = uncompress_tei(tei_compressed, notation='physics')
diagonal_coulomb = of.FermionOperator()
diagonal_coulomb_mat = np.zeros((norbs, norbs))
for i, j in product(range(norbs), repeat=2):
    diagonal_coulomb_mat[i, j] = tei_tensor[i, j, i, j]
    for sigma, tau in product(range(2), repeat=2):
        diagonal_coulomb += of.FermionOperator(
            ((2 * i + sigma, 1), (2 * i + sigma, 0), (2 * j + tau, 1),
             (2 * j + tau, 0)), coefficient=diagonal_coulomb_mat[i, j])
dc_ham = DiagonalCoulomb(diagonal_coulomb_mat)
Explanation: We need to build our Diagoanl Coulomb operator For this bigger system.
End of explanation
cirq_wfn = fqe.to_cirq(fqe_wfn).reshape((-1, 1))
final_cirq_wfn = expm(-1j * of.get_sparse_operator(diagonal_coulomb)) @ cirq_wfn
# recover a fqe wavefunction
from_cirq_wfn = fqe.from_cirq(final_cirq_wfn.flatten(), 1.0E-8)
fqe_wfn = fqe_wfn.time_evolve(1, dc_ham)
print("Evolved wavefunction")
fqe_wfn.print_wfn()
print("From Cirq Evolution")
from_cirq_wfn.print_wfn()
assert np.allclose(from_cirq_wfn.get_coeff((n_elec, sz)),
                   fqe_wfn.get_coeff((n_elec, sz)))
print("Wavefunctions are equivalent")
Explanation: Now we can convert our wavefunction to a cirq wavefunction, evolve under the diagonal_coulomb operator we constructed and then compare the outputs.
End of explanation
fqe_wfn = fqe.Wavefunction([[n_elec, sz, norbs]])
fqe_wfn.set_wfn(strategy='from_data',
                raw_data={(n_elec, sz): inital_coeffs})
for term, coeff in diagonal_coulomb.terms.items():
    op = of.FermionOperator(term, coefficient=coeff)
    fqe_wfn = fqe_wfn.time_evolve(1, op)
assert np.allclose(from_cirq_wfn.get_coeff((n_elec, sz)),
               fqe_wfn.get_coeff((n_elec, sz)))
print("Individual term evolution is equivalent")
Explanation: Finally, we can compare against evolving each term of $V$ individually.
End of explanation
<END_TASK>