!pip install openfermion openfermioncirq pyscf openfermionpyscf

import openfermion as of

op = of.FermionOperator(((4, 1), (3, 1), (9, 0), (1, 0)), 1+2j) + of.FermionOperator(((3, 1), (1, 0)), -1.7)

print(op.terms)

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

print(op.terms)

print(op)

op = of.QubitOperator(((1, 'X'), (2, 'Y'), (3, 'Z')))
op += of.QubitOperator('X3 Z4', 3.0)

print(op)

op = of.FermionOperator('2^ 15')

print(of.jordan_wigner(op))
print()
print(of.bravyi_kitaev(op, n_qubits=16))

a2 = of.FermionOperator('2')
print(of.jordan_wigner(a2))
print()

a2dag = of.FermionOperator('2^')
print(of.jordan_wigner(a2dag*a2))
print()

a7 = of.FermionOperator('7')
a7dag = of.FermionOperator('7^')
print(of.jordan_wigner((1+2j)*(a2dag*a7) + (1-2j)*(a7dag*a2)))

a2 = of.FermionOperator('2')
a2dag = of.FermionOperator('2^')
a7 = of.FermionOperator('7')
a7dag = of.FermionOperator('7^')


print(of.jordan_wigner(a2dag))
print()
print(of.jordan_wigner(a2dag*a2))
print()

op = (2+3j)*a2dag*a7
op += of.hermitian_conjugated(op)
print(of.jordan_wigner(op))

a2_jw = of.jordan_wigner(a2)
a2dag_jw = of.jordan_wigner(a2dag)
a7_jw = of.jordan_wigner(a7)
a7dag_jw = of.jordan_wigner(a7dag)

print(a2_jw * a7_jw + a7_jw * a2_jw)
print(a2_jw * a7dag_jw + a7dag_jw * a2_jw)
print(a2_jw * a2dag_jw + a2dag_jw * a2_jw)

import openfermionpyscf as ofpyscf

# Set molecule parameters
geometry = [('H', (0.0, 0.0, 0.0)), ('H', (0.0, 0.0, 0.8))]
basis = 'sto-3g'
multiplicity = 1
charge = 0

# Perform electronic structure calculations and
# obtain Hamiltonian as an InteractionOperator
hamiltonian = ofpyscf.generate_molecular_hamiltonian(
  geometry, basis, multiplicity, charge)

# Convert to a FermionOperator
hamiltonian_ferm_op = of.get_fermion_operator(hamiltonian)

print(hamiltonian_ferm_op)

import scipy.sparse

# Map to QubitOperator using the JWT
hamiltonian_jw = of.jordan_wigner(hamiltonian_ferm_op)

# Convert to Scipy sparse matrix
hamiltonian_jw_sparse = of.get_sparse_operator(hamiltonian_jw)

# Compute ground energy
eigs, _ = scipy.sparse.linalg.eigsh(hamiltonian_jw_sparse, k=1, which='SA')
ground_energy = eigs[0]

print('Ground_energy: {}'.format(ground_energy))
print('JWT transformed Hamiltonian:')
print(hamiltonian_jw)

import scipy.sparse

# Map to QubitOperator using the JWT
hamiltonian_bk = of.bravyi_kitaev(hamiltonian_ferm_op)

# Convert to Scipy sparse matrix
hamiltonian_bk_sparse = of.get_sparse_operator(hamiltonian_bk)

# Compute ground energy
eigs, _ = scipy.sparse.linalg.eigsh(hamiltonian_bk_sparse, k=1, which='SA')
ground_energy = eigs[0]

print('Ground_energy: {}'.format(ground_energy))
print('BK transformed Hamiltonian:')
print(hamiltonian_bk)

# Map to QubitOperator using the BKT
hamiltonian_bk = of.bravyi_kitaev(hamiltonian_ferm_op)

# Convert to Scipy sparse matrix
hamiltonian_bk_sparse = of.get_sparse_operator(hamiltonian_bk)

# Compute ground state energy
eigs, _ = scipy.sparse.linalg.eigsh(hamiltonian_bk_sparse, k=1, which='SA')
ground_energy = eigs[0]

print('Ground_energy: {}'.format(ground_energy))
print('BKT transformed Hamiltonian:')
print(hamiltonian_bk)

# Create a random initial state
n_qubits = of.count_qubits(hamiltonian)
initial_state = of.haar_random_vector(2**n_qubits, seed=7)

# Set evolution time
time = 1.0

# Apply exp(-i H t) to the state
exact_state = scipy.sparse.linalg.expm_multiply(-1j*hamiltonian_jw_sparse*time, initial_state)

import cirq
import openfermioncirq as ofc
import numpy as np

# Initialize qubits
qubits = cirq.LineQubit.range(n_qubits)

# Create circuit
circuit = cirq.Circuit.from_ops(
    ofc.simulate_trotter(
        qubits, hamiltonian, time,
        n_steps=10,
        order=0,
        algorithm=ofc.trotter.LOW_RANK)
)

# Apply the circuit to the initial state
result = circuit.apply_unitary_effect_to_state(initial_state)

# Compute the fidelity with the final state from exact evolution
fidelity = abs(np.dot(exact_state, result.conj()))**2

print(fidelity)

print(circuit.to_text_diagram(transpose=True))

import cirq
import openfermioncirq as ofc

ansatz = ofc.LowRankTrotterAnsatz(hamiltonian)

cirq.DropNegligible().optimize_circuit(ansatz.circuit)
print(ansatz.circuit.to_text_diagram(transpose=True))

import scipy.optimize

def energy_from_params(x):
  param_resolver = ansatz.param_resolver(x)
  circuit = cirq.resolve_parameters(ansatz.circuit, param_resolver)
  final_state = circuit.apply_unitary_effect_to_state(0b1100)
  return of.expectation(hamiltonian_jw_sparse, final_state).real

initial_guess = ansatz.default_initial_params()
result = scipy.optimize.minimize(energy_from_params, initial_guess)

print('Initial energy: {}'.format(energy_from_params(initial_guess)))
print('Optimized energy: {}'.format(result.fun))

n_qubits = 5
quad_ham = of.random_quadratic_hamiltonian(
    n_qubits, conserves_particle_number=True, seed=7)

print(of.get_fermion_operator(quad_ham))

_, basis_change_matrix, _ = quad_ham.diagonalizing_bogoliubov_transform()

qubits = cirq.LineQubit.range(n_qubits)
circuit = cirq.Circuit.from_ops(
    ofc.bogoliubov_transform(
        qubits,
        basis_change_matrix))

print(circuit.to_text_diagram(transpose=True))

orbital_energies, constant = quad_ham.orbital_energies()

print(orbital_energies)
print(constant)

# Apply the circuit with initial state having the first two modes occupied.
result = circuit.apply_unitary_effect_to_state(0b11000)

# Compute the expectation value of the final state with the Hamiltonian
quad_ham_sparse = of.get_sparse_operator(quad_ham)
print(of.expectation(quad_ham_sparse, result))

# Print out the ground state energy; it should match
print(quad_ham.ground_energy())

from scipy.sparse.linalg import expm_multiply

# Create a random initial state
initial_state = of.haar_random_vector(2**n_qubits)

# Set evolution time
time = 1.0

# Apply exp(-i H t) to the state
final_state = expm_multiply(-1j*quad_ham_sparse*time, initial_state)

# Initialize qubits
qubits = cirq.LineQubit.range(n_qubits)

# Write code below to create the circuit
# You should define the `circuit` variable here
# ---------------------------------------------


# ---------------------------------------------

# Apply the circuit to the initial state
result = circuit.apply_unitary_effect_to_state(initial_state)

# Compute the fidelity with the correct final state
fidelity = abs(np.dot(final_state, result.conj()))**2

# Print fidelity; it should be 1
print(fidelity)

# Initialize qubits
qubits = cirq.LineQubit.range(n_qubits)

# Write code below to create the circuit
# You should define the `circuit` variable here
# ---------------------------------------------
def exponentiate_quad_ham(qubits, quad_ham):
    _, basis_change_matrix, _ = quad_ham.diagonalizing_bogoliubov_transform()
    orbital_energies, _ = quad_ham.orbital_energies()
    
    yield cirq.inverse(
        ofc.bogoliubov_transform(qubits, basis_change_matrix))
    for i in range(len(qubits)):
        yield cirq.Rz(rads=-orbital_energies[i]).on(qubits[i])
    yield ofc.bogoliubov_transform(qubits, basis_change_matrix)

circuit = cirq.Circuit.from_ops(exponentiate_quad_ham(qubits, quad_ham))
# ---------------------------------------------

# Apply the circuit to the initial state
result = circuit.apply_unitary_effect_to_state(initial_state)

# Compute the fidelity with the correct final state
fidelity = abs(np.dot(final_state, result.conj()))**2

# Print fidelity; it should be 1
print(fidelity)