from openfermion.ops import FermionOperator

my_term = FermionOperator(((3,1), (1,0)))
print(my_term)

my_term = FermionOperator('3^ 1')
print(my_term)

"""
Explanation: Introduction to OpenFermion

A codealong of openfermion_tutorial.ipynb
Wayne H Nixalo – 2018/6/27

<div class="alert alert-info">
Note that all the examples below must be run sequentially within a section.
</div>

1. Initializing the FermionOperator data structure

Fermionic systems are often treated in second quantization where arbitrary operators can be expressed using the fermionic creation and annihilation operators $a_k^\dagger$ and $a_k$. The fermionic ladder operators play a similar role to their qubit ladder operator counterparts, $σ_k^\dagger$ and $σ_k^-$ but are distinguished by the canonical fermionic anticommutation relations, ${a_i^\dagger, a_j^\dagger} = {a_i, a_j} = 0$ and ${a_i, a_j^\dagger} = δ_{ij}$. Any weighted sums of products of these oeprators are represented with the FermionOperator data structure in OpenFermion. The following are examples of valid FermionOperators:

<img src="https://render.githubusercontent.com/render/math?math=%5Cbegin%7Balign%7D%0A%26amp%3B%20a_1%20%5Cnonumber%20%5C%5C%0A%26amp%3B%201.7%20a%5E%5Cdagger_3%20%5Cnonumber%20%5C%5C%0A%26amp%3B-1.7%20%5C%2C%20a%5E%5Cdagger_3%20a_1%20%5Cnonumber%20%5C%5C%0A%26amp%3B%281%20%2B%202i%29%20%5C%2C%20a%5E%5Cdagger_4%20a%5E%5Cdagger_3%20a_9%20a_1%20%5Cnonumber%20%5C%5C%0A%26amp%3B%281%20%2B%202i%29%20%5C%2C%20a%5E%5Cdagger_4%20a%5E%5Cdagger_3%20a_9%20a_1%20-%201.7%20%5C%2C%20a%5E%5Cdagger_3%20a_1%20%5Cnonumber%0A%5Cend%7Balign%7D&mode=display">

The FermionOperator class is contained in ops/_fermion_operators.py. In order to support fast addition of FermionOperator instances, the class is implemented as a hash table (python dictionary). The keys of the dictionary encoe the strings of ladder operators and values of the dictionary store the coefficients. The strings of ladder operators are encocded as a tuple of 2-tuples which we refer to as the "terms tuple". Each ladder operator is represented by a 2-tuple. The first element of the 2-tuple is an int indicating the tensor factor on which the ladder operator acts. The second element of the 2-tuple is Boole: 1 represents raising and 0 represents lowering. For instance $a_8^\dagger$ is represented in a 2-tuple as $(8,1)$. Note that indices start at 0 and the identity operator is an empty list. Below we give some examples of operators and their terms tuple:

<img src="https://render.githubusercontent.com/render/math?math=%5Cbegin%7Balign%7D%0AI%20%26amp%3B%20%5Cmapsto%20%28%29%20%5Cnonumber%20%5C%5C%0Aa_1%20%26amp%3B%20%5Cmapsto%20%28%281%2C%200%29%2C%29%20%5Cnonumber%20%5C%5C%0Aa%5E%5Cdagger_3%20%26amp%3B%20%5Cmapsto%20%28%283%2C%201%29%2C%29%20%5Cnonumber%20%5C%5C%0Aa%5E%5Cdagger_3%20a_1%20%26amp%3B%20%5Cmapsto%20%28%283%2C%201%29%2C%20%281%2C%200%29%29%20%5Cnonumber%20%5C%5C%0Aa%5E%5Cdagger_4%20a%5E%5Cdagger_3%20a_9%20a_1%20%26amp%3B%20%5Cmapsto%20%28%284%2C%201%29%2C%20%283%2C%201%29%2C%20%289%2C%200%29%2C%20%281%2C%200%29%29%20%5Cnonumber%0A%5Cend%7Balign%7D&mode=display">

Let's initialize our first term! We do it two different ways below.
End of explanation
"""


good_way_to_initialize = FermionOperator('3^ 1', -1.7)
print(good_way_to_initialize)

bad_way_to_initialize = -1.7 * FermionOperator('3^ 1')
print(bad_way_to_initialize)

identity = FermionOperator('')
print(identity)

zero_operator = FermionOperator()
print(zero_operator)

"""
Explanation: Ther preffered way to specify the coefficient in openfermion is to provide an optional coefficient argument. If not provided, the coefficient defaults to 1. In the code below, the first method is preferred. The multiplication in the second method actually creates a copy of the term, which introduces some additional cost. All inplace operants (such as +=) modify classes wheras binary operands (such as +) create copies. Important caveats are that the empty tuple FermionOperator(()) and the empty string FermionOperator('') initializes identically. The empty initializer FermonOperator() initializes the zero operator.
End of explanation
"""


my_operator = FermionOperator('4^ 1^ 3 9', 1. + 2.j)
print(my_operator)
print(my_operator.terms)

"""
Explanation: <div class="alert alert-info">
Note that FermionOperator has only 1 attribute: .terms. This attribute is the dictionary which stores the term tuples.
</div>
End of explanation
"""


from openfermion.ops import FermionOperator

term_1 = FermionOperator('4^ 3^ 9 1', 1. + 2.j)
term_2 = FermionOperator('3^ 1', -1.7)
my_operator = term_1 + term_2
print(my_operator)

my_operator = FermionOperator('4^ 3^ 9 1', 1. + 2.j)
term_2 = FermionOperator('3^ 1', -1.7)
my_operator += term_2
print(); print(my_operator)

"""
Explanation: 2. Manipulating the FermionOperator data structure

So far we have explained how to initialize a single FermionOperator such as $-1.7a_3^\dagger a_1$. However, in general we will want to represent sums of these operators such as $(1 + 2i) \,a_4^\dagger a_3^\dagger a_9 a_1 - 1.7 \,a_3^\dagger a_1$. To do this, just add together two FermionOperators! We demonstrate below.
End of explanation
"""


term_1 = FermionOperator('4^ 3^ 9 1', 1. + 2.j)
term_2 = FermionOperator('3^ 1', -1.7)

my_operator = term_1 -33. * term_2
print(my_operator)

my_operator *= 3.17 * (term_2 + term_1) ** 2
print(); print(my_operator)

print(); print(term_2 ** 3)

print(); print(term_1 == 2. * term_1 - term_1)
print(term_1 == my_operator)

"""
Explanation: The print function prints each term in the operator on a different line. Note that the line my_operator = term_1 + term_2 creates a new object, which involves a copy of term_1 and term_2. The second block of code uses the inplace method +=, which is more efficient. This is especially import when trying to construct a very large FermionOperator. FermionOperators also support a wide range of builtins including str(), repr(), ==, !=, =, /, /=, +, +=, -, -=, and **. Note that since FermionOperators involve floats, == and != check from (in)equality up to numerical precision. We demonstrate some of these methods below.
End of explanation
"""


from openfermion.utils import commutator, count_qubits, hermitian_conjugated, normal_ordered

# Get the Hermitian conjugate of a FermionOperator, 
# count its qubits check if it's normal-ordered
term_1 = FermionOperator('4^ 3 3^', 1. + 2.j)
print(hermitian_conjugated(term_1))
print(term_1.is_normal_ordered())
print(count_qubits(term_1))

# Normal order the term
term_2 = normal_ordered(term_1)
print(); print(term_2); print(term_2.is_normal_ordered())

# Compute a commutator of the terms
print(); print(commutator(term_1, term_2))

"""
Explanation: Additionally, there are a variety of methods that act on the FermionOperator data structure. We demonstrate a small subset of those methods here.
End of explanation
"""


from openfermion.ops import QubitOperator

my_first_qubit_oeprator = QubitOperator('X1 Y2 Z3')
print(my_first_qubit_oeprator)
print(my_first_qubit_oeprator.terms)

operator_2 = QubitOperator('X3 Z4', 3.17)
operator_2 -= 77. * my_first_qubit_oeprator
print(f'\n{operator_2}')

"""
Explanation: 3. The QubitOperator data structure

The QubitOperator data structure is another essential part of OpenFermion. As the name suggests, QubitOperator is used to store qubit operators in almost exactly the same way that FermionOperators is used to store fermion operators. For instance, $X_0Z_3Y_4$ is a QubitOperator. The inernal representation of this as a terms tuple would be $((0,"X"),(3,"Z"),(4,"Y"))$. Note that one important difference between QubitOperator and FermionOperator is that the terms in QubitOperator are always sorted in order of tensor factor. In some cases, this enables faster manipulation. We initialize some QubitOperators below.
End of explanation
"""


from openfermion.ops import FermionOperator
from openfermion.transforms import jordan_wigner, bravyi_kitaev
from openfermion.utils import eigenspectrum, hermitian_conjugated

# Initialize an operator
fermion_operator = FermionOperator('2^ 0', 3.17)
fermion_operator += hermitian_conjugated(fermion_operator)
print(fermion_operator)

# Transform to qubits under the Jordan-Wigner transformation and print its spectrum
jw_operator = jordan_wigner(fermion_operator)
print(f'\n{jw_operator}')
jw_spectrum = eigenspectrum(jw_operator)
print(jw_spectrum)

# Transform to qubits under the Bravyi-Kitaev transformation and print its spectrum
bk_operator = bravyi_kitaev(fermion_operator)
print(f'\n{bk_operator}')
bk_spectrum = eigenspectrum(bk_operator)
print(f'{bk_spectrum}')

"""
Explanation: 4. Jordan-Wigner and Barvyi-Kitaev

openfermion provides functions for mapping FermionOpeartors to QubitOperators.
End of explanation
"""


from openfermion.transforms import reverse_jordan_wigner

# Initialize the QubitOperator
my_operator = QubitOperator('X0 Y1 Z2', 88.)
my_operator += QubitOperator('Z1 Z4', 3.17)
print(my_operator)

# Map QubitOperator to a FermionOperator
mapped_operator = reverse_jordan_wigner(my_operator)
print(f'\n{mapped_operator}')

# Map the operator back to qubits and make sure it's the same
back_to_normal = jordan_wigner(mapped_operator)
back_to_normal.compress()
print(f'\n{back_to_normal}')

"""
Explanation: We see that despite the different representation, these operators are iso-spectral. We can also apply the Jordan-Wigner transform in reverse to map arbitrary QubitOperators to FermionOperators. Note that we also demonstrate the .compress() method (a method on both FermionOperators and QubitOperators) which removes zero entries.
End of explanation
"""


from openfermion.hamiltonians import fermi_hubbard
from openfermion.transforms import get_sparse_operator, jordan_wigner
from openfermion.utils import get_ground_state

# Set model
x_dimension = 2
y_dimension = 2
tunneling = 2.
coulomb = 1.
magnetic_field = 0.5
chemical_potential = 0.25
periodic = 1
spinless = 1

# Get fermion operator
hubbard_model = fermi_hubbard(
    x_dimension, y_dimension, tunneling, coulomb, chemical_potential,
    magnetic_field, periodic, spinless)
print(hubbard_model)

# Get qubit operator under Jordan-Wigner
jw_hamiltonian = jordan_wigner(hubbard_model)
jw_hamiltonian.compress()
print(f'\n{jw_hamiltonian}')

# Get the scipy.sparse.csc representation
sparse_operator = get_sparse_operator(hubbard_model)
print(f'\n{sparse_operator}\nEnergy of the model is {get_ground_state(sparse_operator)[0]} in units of T and J.')

"""
Explanation: 5. Sparse matrices and the Hubbard model

Often, one would like to obtain a sparse matrix representation of an operator which can be analyzed numerically. There's code in both openfermion.transforms and openfermion.utils which facilitates this. The function get_sparse_operator converts either a FermionOperator, a QubitOperator, or other more advanced classes such as InteractionOperator to a scipy.sparse.csc matrix. There're numerous functions in openfermion.utils which one can call on the sparse operators such as "get_gap", "get_hartree_fock_state", "get_ground_state", etc. We show this off by computing the ground state energy of the Hubbard model. To do that, we use code from the openfermion.hamiltonians module which constructs lattice models of fermions such as Hubbard models.
End of explanation
"""


from openfermion.hamiltonians import jellium_model
from openfermion.utils import eigenspectrum, fourier_transform, Grid
from openfermion.transforms import jordan_wigner

# Let's look at a very small model of jellium in 1D
grid = Grid(dimensions=1, length=3, scale=1.0)
spinless = True

# Get the momentum Hamiltonian
momentum_hamiltonian = jellium_model(grid, spinless)
momentum_qubit_operator = jordan_wigner(momentum_hamiltonian)
momentum_qubit_operator.compress()
print(momentum_qubit_operator)

# Fourier transform the Hamiltonian to the position basis
position_hamiltonian = fourier_transform(momentum_hamiltonian, grid, spinless)
position_qubit_operator = jordan_wigner(position_hamiltonian)
position_qubit_operator.compress()
print(f'\n{position_qubit_operator}')

# Check the spectra to make sure these representations are iso-spectral
spectral_difference = eigenspectrum(momentum_qubit_operator) - eigenspectrum(position_qubit_operator)
print(f'\n{spectral_difference}')

"""
Explanation: 6. Hamiltonians in the plane wave basis

A user can write plugins to openfermion which allow for the use of, eg: 3rd-party electronic structure packages to compute molecular orbitals, Hamiltonians, energies, reduced density matrices, coupled cluster amplitudes, etc, using Gaussian basis sets. We may provide scripts which interface between such packages and openfermion in the future but do not discuss them in this tutorial.
When using simpler basis sets such as plane waves, these packages are not needed. OpenFermion comes with code which computes Hamiltonians in the plane wave basis. Note that when using plane waves, one is working with the periodized Coulomb operator, best suited for condensed phase calculations such as studying the electronic structure of a solid. To obtain these Hamiltonians one must choose to study the system without a spin degree of freedom (spinless), one must specify the dimension in which the calculation is performed (n_dimensions, usually 3), one must specify how many plane waves are in each dimension (grid_length) and one must specify the length scale of the plane wave harmonics in each dimension (length_scale) and also the locations and charges of the nuclei. One can generate these models with plane_wave_hamiltonian() found in openfermion.hamiltonians. For simplicity, below we compute the Hamiltonian in the case of zero external charge (corresponding to the uniform electron gas, aka Jellium). We also demonstrate that one can transform the plane wave Hamiltonian using a Fourier transform without effecting the spectrum of the operator.
End of explanation
"""


from openfermion.hamiltonians import MolecularData

# Set parameters to make a simple molecule
diatomic_bond_length = .7414
geometry = [('H', (0.,0.,0.)), ('H', (0.,0.,diatomic_bond_length))]
basis = 'sto-3g'
multiplicity = 1
charge = 0
description = str(diatomic_bond_length)

# Make molecule and print out a few interesting facts about it
molecule = MolecularData(geometry, basis, multiplicity,
                         charge, description)
print(f'Molecule has automatically generated name {molecule.name}')
print(f'Information about this molecule would be saved at:\n{molecule.filename}\n')
print(f'This molecule has {molecule.n_atoms} atoms and {molecule.n_electrons} electrons.')
for atom,atomic_number in zip(molecule.atoms, molecule.protons):
    print(f'Contains {atom} atom, which has {atomic_number} protons.')

"""
Explanation: 7. Basics of MolecularData class

Data from electronic structure calculations can be saved in an OpenFermion data structure  called MolecularData, which makes it easy to access within our library. Often, one would like to analyze a chemical series or look at many different Hamiltonians and sometimes the elecronic structure calculations are either expensive to compute or difficult to converge (eg: one needs to mess around with different types of SCF routines to make things converge). Accordingly, we anticipate that users will want some way to automatically database the results of their electronic structure calculations so that important data (such as the SCF integrals) can be looked up on-the-fly if the user has computed them in the past. OpenFermion supports a data provenance strategy which saves key results of the electronic structure calculatin (including pointers to files containing large amounts of data, such as the molecular integrals) in an HDF5 container.
The MolecularData class stores information about molecules. One initializes a MolecularData object by specifying parameters of a molecule such as its geometry, basis, multiplicity, charge, and an optional string describing it. One can also initialize MolecularData simply by providing a string giving a filename where a previous MolecularData object was saved in an HDF5 container. One can save a MolecularData instance by calling the class' .save() method. This automatically saves the instance in a data folder specified during OpenFermion installation. The name of the file is generated automatically from the instance attributes and optionally provided description. Alternatively, a filename can also be provided as an optional input if one wishes to manually name the file.
When electronic structure calculations are run, the data files for the molecule can be automatically updated. If one wishes to later use the data they either initialize MolecularData with the instance filename or initialize the instance and then later call the .load() method.
Basis functions are provided to initialization using a string such as "6-31g". Geometries can be specified using a simple txt input file (see the geometry_from_file function in molecular_data.py) or can be passed using a simple python list format demonstrated below. Atoms are specified using a string for their atomic symbol. Distances should be provided in ångstrom. Below we initialize a simple instance of MolecularData without performing any electronic structure calculations.
End of explanation
"""


%matplotlib inline

# Set molecule parameters
basis = 'sto-3g'
multiplicity = 1
bond_length_interval = 0.1
n_points = 25

# Generate molecule at different bond lengths
hf_energies = []
fci_energies = []
bond_lengths = []
for point in range(3, n_points + 1):
    bond_length = bond_length_interval * point
    bond_lengths += [bond_length]
    description = str(round(bond_length,2))
    print(f'\n{description}')
    geometry = [('H', (0.,0.,0.)), ('H', (0.,0., bond_length))]
    molecule = MolecularData(
        geometry, basis, multiplicity, description=description)
    
    # Load data
    molecule.load()
    
    # Print out some calculation results
    print(f'At bond length of {bond_length} ångstrom, molecular hydrogen has:')
    print(f'Hartree-Fock energy of {molecule.hf_energy} Hartree.')
    print(f'MP2 energy of {molecule.mp2_energy} Hartree.')
    print(f'FCI energy of {molecule.fci_energy} Hartree.')
    print(f'Nuclear repulsion energy between protons is {molecule.nuclear_repulsion} Hartree.')
    for orbital in range(molecule.n_orbitals):
        print(f'Spatial orbital {orbital} has energy of {molecule.orbital_energies[orbital]} Hartree.')
    hf_energies += [molecule.hf_energy]
    fci_energies += [molecule.fci_energy]
    
# Plot
import matplotlib.pyplot as plt

plt.figure(0)
plt.plot(bond_lengths, fci_energies, 'x-'); plt.plot(bond_lengths, hf_energies, 'o-');
plt.ylabel('Energy in Hartree'); plt.xlabel('Bond length in ångstrom');

"""
Explanation: If we had previously computed this molecule using an electronic structure package, we can call molecule.load() to populate all sorts of interesting fields in the data structure. Though we make no assumptions about what electronic structure packages users might install, we assume that the calculations are saved in OpenFermion's MolecularData objects. Currently plugins are availabel for Psi4, (OpenFermion-Psi4), and PySCF (OpenFermion-PrSCF), and there may be more in the future. For the purposes of this example, we'll load data that ships with OpenFermion to make a plot of the energy surface of hydrogen. Note that helper functions to initialize some interesting chemical benchmarks are found in openfermion.utils.
End of explanation
"""


from openfermion.utils import geometry_from_pubchem

methane_geometry = geometry_from_pubchem('methane')
print(methane_geometry)

"""
Explanation: The geometry data needed to generate MolecularData can aplso be retrived from the PubChem online database by inputting the molecule's name.
End of explanation
"""


from openfermion.hamiltonians import MolecularData
from openfermion.transforms import get_fermion_operator, get_sparse_operator, jordan_wigner
from openfermion.utils import get_ground_state
import numpy as np
import scipy
import scipy.linalg

# Load saved file for LiH
diatimoc_bond_length = 1.45
geometry = [('Li', (0.,0.,0.,)), ('H', (0.,0.,diatomic_bond_length))]
basis = 'sto-3g'
multiplicity = 1

# Set Hamiltonian parameters
active_space_start = 1
active_space_stop = 3

# Generate and populate instance of MolecularData
molecule = MolecularData(geometry, basis, multiplicity, description='1.45')
molecule.load()

# Get the Hamiltonian in an active space
molecular_hamiltonian = molecule.get_molecular_hamiltonian(
    occupied_indices=range(active_space_start),
    active_indices=range(active_space_start, active_space_stop))

# Map operator to fermions and qubits
fermion_hamiltonian = get_fermion_operator(molecular_hamiltonian)
qubit_hamiltonian = jordan_wigner(fermion_hamiltonian)
qubit_hamiltonian.compress()
print(f'The Jordan-Wigner Hamiltonian in canonical basis follows:\n{qubit_hamiltonian}')

# Get sparse operator and ground state energy
sparse_hamiltonian = get_sparse_operator(qubit_hamiltonian)
energy,state = get_ground_state(sparse_hamiltonian)
print(f'Ground state energy before rotation is {energy} Hartree.\n')

# Randomly rotate
n_orbitals = molecular_hamiltonian.n_qubits // 2
n_variables = int(n_orbitals * (n_orbitals - 1) / 2)
np.random.seed(1)
random_angles = np.pi * (1. - 2. * np.random.rand(n_variables))
k = np.zeros((n_orbitals, n_orbitals))
index = 0
for p in range(n_orbitals):
    for q in range(p + 1, n_orbitals):
        k[p, q] = random_angles[index]
        k[q, p] = -np.conjugate(random_angles[index])
        index += 1
        
    # Build the unitary rotation matrix
    difference_matrix = k + k.transpose()
    rotation_matrix = scipy.linalg.expm(k)
    
    # Apply the unitary
    molecular_hamiltonian.rotate_basis(rotation_matrix)
    
# Get qubit Hamiltonian in rotate basis
qubit_hamiltonian = jordan_wigner(molecular_hamiltonian)
qubit_hamiltonian.compress()
print(f'The Jordan-Wigner Hamiltonian in rotate basis follows:\n{qubit_hamiltonian}')

# Get sparse Hamiltonian and energy in rotated basis
sparse_hamiltonian = get_sparse_operator(qubit_hamiltonian)
energy,state = get_ground_state(sparse_hamiltonian)
print(f'Ground state energy after rotation is {energy} Hartree.')

"""
Explanation: 8. InteractionOperator and InteractionRDM for efficient numerical representations

Fermion Hamiltonians can be expressed as $H = h_0 + \sum_{pq}h_{pq}a^\dagger_p a^\dagger_q a_r a_s$, where $h_0$ is a constant shift due to nuclear repulsion and $h_{pq}$ and $h_{pqrs}$ are the famous molecular integrals. Since fermions interact pairwise, their energy is thus a unique function of the one-particle and two-particle reduced densitty matrices which are expressed in second quantization as $ρ_{pq} = \big\langle p\lvert a^\dagger_p a_q \rvert q\big\rangle$ and $ρ_{pqrs} = \big\langle pq \lvert a_p^\dagger a_q^\dagger a_r a_s \rvert rs \big\rangle$, respectively.
Because the RDMs and molecular Hamiltonians are both compactly represented and manipulated as 2- and 4- index tensors, we can represent them in a particularly efficient form using similar data structures. The InteractionOperator data structure can be initialized for a Hamiltonian by passing the constant $h_0$ (or 0), as well as numpy arrays representing $h_{pq}$ (or $ρ_{pq}$) and $h_{pqrs}$ (or $ρ_{pqrs}$). Importantly, InteractionOperators can also be obtained by calling MolecularData.get_molecular_hamiltonian() or by calling the function get_interaction_operator() found in openfermion.transforms) on a FermionOperator. The InteractionRDM data structure is similar but represents RDMs. For instance, one can get a molecular RDM by calling MolecularData.get_molecular_rdm(). When generating Hamiltonians from the MolecularData class, one can choose to restrict the system to an active space.
These classes inherit from the same base class, PolynomialTensor. This data structure overloads the slice oeprator [] so that one can get or set the key attributes of the InteractionOperator: .constant, .one_body_coefficients, and .two_body_coefficients. For instance, InteractionOperator[(p,1), (q,1), (r,0), (s,0)] would return $h_{pqrs}$ and InteractionRDM would return $ρ_{pqrs}$. Importantly, the class supports fast basis transformations using the method PolynomialTensor.rotate_basis(rotation_matrix). But perhaps most importantly, one can map the InteractionOperator to any of the other data structures we've described here.
Below, we load MolecularData from a saved calculation of LiH. We then obtain an InteractionOperator representation of this system in an active space. We then map that operator to qubits. We then demonstrate that one can rotate the orbital basis of the InteractionOperator using random angles to obtain a totally different operator that is still iso-spectral.
End of explanation
"""


from openfermion.hamiltonians import mean_field_dwave
from openfermion.transforms import get_quadratic_hamiltonian

# Set model
x_dim = 2
y_dim = 2
tunenling = 2.
sc_gap = 1.
periodic = True

# Get FermionOperator
mean_field_model = mean_field_dwave(x_dim, y_dim,
                        tunneling, sc_gap, periodic=periodic)

# Covnert to Quadratic Hamiltonian
quadratic_hamiltonian = get_quadratic_hamiltonian(mean_field_model)

# Compute the ground energy
ground_energy = quadratic_hamiltonian.ground_energy()
print(ground_energy)

"""
Explanation: 9. Quadratic Hamiltonians and Slater determinants

The general electronic structure Hamiltonian $H = h_0 + \sum_{pq} h_{pq} a^\dagger_p a_q + \frac{1}{2}\sum_{pqrs} h_{pqrs} a^\dagger_p a^\dagger_q a_r a_s$ contains terms that act on up to 4 sites, or is quartic in the fermionic creation and annihilation operators. However, in many situations we may fruitfully approximate these Hamiltonians by replacing these quartic terms with terms that act on at most 2 fermionic sites, or quadratic terms, as in mean-field approximation theory.
These Hamiltonians have a number of special properties one can exploit for efficient simulation and manipulation of the Hamiltonian, thus warranting a special data structure. We refer to Hamiltonians which only contain terms that are quadratic in the fermionic creation and annihilation operators as quadtratic Hamiltonians, and include the general case of non-particle conserving terms as in a general Боголюбов transformation. Eigenstates of quadtratic Hamiltonians can be prepared efficiently on both a quantum and classical computer, making them amenable to initial guesses for many more challenging problems.
A general quadratic Hamiltonian takes the form
<img src="https://render.githubusercontent.com/render/math?math=H%20%3D%20%5Csum_%7Bp%2C%20q%7D%20%28M_%7Bpq%7D%20-%20%5Cmu%20%5Cdelta_%7Bpq%7D%29%20a%5E%5Cdagger_p%20a_q%20%2B%20%5Cfrac%7B1%7D%7B2%7D%20%5Csum_%7Bp%2C%20q%7D%20%28%5CDelta_%7Bpq%7D%20a%5E%5Cdagger_p%20a%5E%5Cdagger_q%20%2B%20%5CDelta_%7Bpq%7D%5E%2A%20a_q%20a_p%29%20%2B%20%5Ctext%7Bconstant%7D%2C&mode=display">
where $M$ is a Hermitian matrix, $Δ$ is an antisymmetric matrix, $δ_{pq}$ is the Kronecker delta symbol, and $μ$ is a chemical potential term which we keep spearate from $M$ so that we can use it to adjust the expectation of the total number of particles. In OpenFermion, quadtratic Hamiltonians are conveniently represented and manipulated using the QuadraticHamiltonian class, which stores $M$, $Δ$, $μ$, and the constant. It's specialized to exploit the properties unique to quadratic Hamiltonians. Like InteractionOperator and InteractionRDM, it inherits from the PolynomialTensor class.
The BCS mean-field model of superconductivity is a quadratic Hamiltonian. The following code constructs an instance of this model as a FermionOperator, converts it to a QuadraticHamiltonian, and then computes its ground energy:
End of explanation
"""


orbital_energies,constant = quadratic_hamiltonian.orbital_energies()
print(f'{orbital_energies}\n{constant}')

"""
Explanation: Any quadtratic Hamiltonian may be rewritten in the form
<img src="https://render.githubusercontent.com/render/math?math=H%20%3D%20%5Csum_p%20%5Cvarepsilon_p%20b%5E%5Cdagger_p%20b_p%20%2B%20%5Ctext%7Bconstant%7D%2C&mode=display">
where the $b_p$ are new annihilation operators that satisfy the fermionic anticommutation relations, and which are linear combinations of the old creation and annihilation operators. This form of $H$ makes it easy to deduce its eigenvalues; they are sums of subsets of the $ε_p$, which we call the orbital energies of $H$. The following code computes the orbital energies and the constant:
End of explanation
"""


from openfermion.utils import gaussian_state_preparation_circuit

circuit_description,start_orbitals = gaussian_state_preparation_circuit(quadratic_hamiltonian)
for parallel_ops in circuit_description:
    print(parallel_ops)
print(f'\n{start_orbitals}')

"""
Explanation: Eigenstates of quadratic hamiltonians are known as fermionic Gaussian states, and they can be prepared efficiently on a quantum computer. One can use OpenFermion to obtain circuits for preparing these states. The following code obtains the description of a circuit which prepares the ground state (operations that can be performed in parallel are grouped together), along with a description of the starting state to which the circuit should be applied:
End of explanation
"""


