%matplotlib inline
# Importing standard Qiskit libraries
from qiskit import QuantumCircuit, execute, Aer, IBMQ
from qiskit.compiler import transpile, assemble
from qiskit.tools.jupyter import *
from qiskit.visualization import *
from ibm_quantum_widgets import *
#from iqx import *

from qiskit.aqua.circuits import PhaseEstimationCircuit
from qiskit.chemistry.components.initial_states import HartreeFock
from qiskit.chemistry.core import Hamiltonian
from qiskit.circuit.library import PhaseEstimation
from qiskit import QuantumCircuit, execute, Aer
from qiskit.circuit import QuantumRegister, Qubit, Gate, ClassicalRegister
from qiskit.quantum_info import Statevector

# Loading your IBM Q account(s)
provider = IBMQ.load_account()

from qiskit.chemistry.drivers import PySCFDriver, UnitsType, Molecule, PSI4Driver
molecule = Molecule(geometry=[['H', [0., 0., 0.]],
                              ['H', [0., 0., 0.735]]],
                     charge=0, multiplicity=1)
driver = PySCFDriver(molecule = molecule, unit=UnitsType.ANGSTROM, basis='sto3g')
qmol = driver.run()

print(qmol.one_body_integrals)
print(qmol.two_body_integrals)

from qiskit.chemistry.transformations import (FermionicTransformation, 
                                              FermionicTransformationType, 
                                              FermionicQubitMappingType)

fermionic_transformation = FermionicTransformation(
            transformation=FermionicTransformationType.FULL,
            qubit_mapping=FermionicQubitMappingType.JORDAN_WIGNER,
            two_qubit_reduction=False,
            freeze_core=False)

qubit_op, _ = fermionic_transformation.transform(driver)
print('Qubit operator is', qubit_op)

num_orbitals = fermionic_transformation.molecule_info['num_orbitals']
num_particles = fermionic_transformation.molecule_info['num_particles']
qubit_mapping = fermionic_transformation.qubit_mapping
two_qubit_reduction = fermionic_transformation.molecule_info['two_qubit_reduction']
z2_symmetries = fermionic_transformation.molecule_info['z2_symmetries']

initial_state = HartreeFock(num_orbitals, num_particles, qubit_mapping,
                            two_qubit_reduction, z2_symmetries.sq_list)

initial_circ = initial_state.construct_circuit()

print(initial_state.bitstr)

initial_circ.draw()

qubit_op, _ = fermionic_transformation.transform(driver)
print(qubit_op)

num_evaluation_qubits = 5
normalization_factor = 1/2 # Since the min energy will be -1.8... we need to normalize it so that it is smaller than 1
unitary = (normalization_factor*qubit_op).exp_i().to_circuit()
phase_estimation = PhaseEstimation(num_evaluation_qubits = num_evaluation_qubits, unitary = unitary)

type(qubit_op.exp_i().to_circuit())

print(phase_estimation.qubits)
dir(phase_estimation)
phase_estimation.draw()

circ = initial_circ.combine(phase_estimation)
classical_eval = ClassicalRegister(num_evaluation_qubits, name = 'class_eval')
circ = circ + QuantumCircuit(classical_eval)
#circ.measure_all()
print(circ.qregs)
circ.measure(circ.qregs[1], classical_eval)
circ.draw()

# Use Aer's qasm_simulator
backend_sim = Aer.get_backend('qasm_simulator')

# Execute the circuit on the qasm simulator.
# We've set the number of repeats of the circuit
# to be 1024, which is the default.
job_sim = execute(circ, backend_sim, shots=1024)

# Grab the results from the job.
result_sim = job_sim.result()

counts = result_sim.get_counts(circ)
print(counts)

maxkey = '0'*num_evaluation_qubits
maxvalue = 0
for key, value in counts.items():
    if value > maxvalue:
        maxvalue = value
        maxkey = key

energy = 0
for i,s in zip(range(num_evaluation_qubits), maxkey[::-1]):
    energy -= (int(s)/2**(i))/normalization_factor
    
print('The energy is', energy, 'Hartrees')

36/32

26/16

