from qiskit import BasicAer, Aer, IBMQ
from qiskit.aqua import QuantumInstance, aqua_globals
from qiskit.aqua.algorithms import VQE, ExactEigensolver, NumPyEigensolver
from qiskit.aqua.components.initial_states import Zero
from qiskit.aqua.components.optimizers import COBYLA, L_BFGS_B, SLSQP, SPSA, CRS
from qiskit.aqua.components.variational_forms import RY, RYRZ, SwapRZ
from qiskit.aqua.operators import WeightedPauliOperator, Z2Symmetries
from qiskit.chemistry import FermionicOperator
from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from qiskit.chemistry.components.variational_forms import UCCSD
from qiskit.chemistry.components.initial_states import HartreeFock

from qiskit.providers.aer import QasmSimulator
from qiskit.providers.aer.noise import NoiseModel
from qiskit.providers.aer.noise.errors import QuantumError, ReadoutError
from qiskit.providers.aer.noise.errors import pauli_error
from qiskit.providers.aer.noise.errors import depolarizing_error
from qiskit.providers.aer.noise.errors import thermal_relaxation_error

from qiskit.providers.aer import noise

from qiskit.tools.monitor import job_monitor

from qiskit.providers.ibmq import least_busy

import warnings

IBMQ.load_account()
provider = IBMQ.get_provider(group='open')
backend = least_busy(provider.backends(simulator=False, operational=True))

import numpy as np
import matplotlib.pyplot as plt
from functools import partial

inter_dist = 1.5108585
driver = PySCFDriver(atom='Li .0 .0 .0; H .0 .0' + str(inter_dist), unit=UnitsType.ANGSTROM, charge=0, spin=0, basis='sto3g')

# please be aware that the idx here with respective to original idx
freeze_list = [0]
remove_list = [-3, -2] # negative number denotes the reverse order

#Classically solve for the lowest eigenvalue
def exact_solver(qubitOp):
    ee = ExactEigensolver(qubitOp)
    result = ee.run()
    ref = result['energy']
    print('Reference value: {}'.format(ref))
    return ref

backend = BasicAer.get_backend('statevector_simulator')

#Define our noise model based on the ibmq_essex chip
chip_name = 'ibmq_essex'
device = provider.get_backend(chip_name)
coupling_map = device.configuration().coupling_map
noise_model = noise.device.basic_device_noise_model(device.properties())
basis_gates = noise_model.basis_gates

# Classically solve for the lowest eigenvalue
# This is used just to compare how well you VQE approximation is performing
def exact_solver(qubitOp):
    ee = ExactEigensolver(qubitOp)
    result = ee.run()
    ref = result['energy']
    print('Reference value: {}'.format(ref))
    return ref

counts = []
values = []
params = []
deviation = []
def store_intermediate_result(eval_count, parameters, mean, std):
    counts.append(eval_count)
    values.append(mean)
    params.append(parameters)
    deviation.append(std)

# Define your function for computing the qubit operations of LiH
def compute_LiH_qubitOp(map_type, inter_dist, basis='sto3g'):
    
    # Specify details of our molecule
    driver = PySCFDriver(atom='Li 0 0 0; H 0 0 ' + str(inter_dist), unit=UnitsType.ANGSTROM, charge=0, spin=0, basis=basis)

    # Compute relevant 1 and 2 body integrals.
    molecule = driver.run()
    h1 = molecule.one_body_integrals
    h2 = molecule.two_body_integrals
    
    num_particles = molecule.num_alpha + molecule.num_beta
    num_spin_orbitals = molecule.num_orbitals * 2
    nuclear_repulsion_energy = molecule.nuclear_repulsion_energy
    
    print("# of electrons: {}".format(num_particles))
    print("# of spin orbitals: {}".format(num_spin_orbitals))

    # Please be aware that the idx here with respective to original idx
    freeze_list = [0,1,6,7]
    remove_list = [0,1,4,5] # negative number denotes the reverse order
    
    # Prepare full idx of freeze_list and remove_list
    # Convert all negative idx to positive
    
    num_spin_orbitals -= len(remove_list)

    # Prepare fermionic hamiltonian with orbital freezing and eliminating, and then map to qubit hamiltonian
    # and if PARITY mapping is selected, reduction qubits
    energy_shift = 0.0
    qubit_reduction = True if map_type == 'parity' else False
    ferOp = FermionicOperator(h1=h1, h2=h2)
    if len(freeze_list) > 0:
        ferOp, energy_shift = ferOp.fermion_mode_freezing(freeze_list)
        num_spin_orbitals -= len(freeze_list)
        num_particles -= len(freeze_list)//2
    if len(freeze_list) > 0:
        ferOp = ferOp.fermion_mode_elimination(remove_list)
    qubitOp = ferOp.mapping(map_type)
    qubitOp = Z2Symmetries.two_qubit_reduction(qubitOp, num_particles)
    return qubitOp, num_spin_orbitals, num_particles, qubit_reduction, molecule, energy_shift, nuclear_repulsion_energy

map_type = "parity"
qubitOp, num_spin_orbitals, num_particles, qubit_reduction, molecule, energy_shift, nuclear_repulsion_energy = compute_LiH_qubitOp(map_type, inter_dist)

# Classically solve for the exact solution and use that as your reference value
ref = exact_solver(qubitOp)

# Specify your initial state
init_state = HartreeFock(num_spin_orbitals, num_particles, qubit_mapping=map_type)

# Select a state preparation ansatz
# Equivalently, choose a parameterization for our trial wave function.
UCCSD_var_form = UCCSD(num_orbitals=num_spin_orbitals, num_particles=num_particles, qubit_mapping=map_type,
                       initial_state=init_state, two_qubit_reduction=qubit_reduction, reps=2)

# Choose where to run/simulate our circuit
quantum_instance = backend

max_eval = 10000
# Choose the classical optimizer
optimizer = SPSA(max_eval)

warnings.filterwarnings("ignore")
vqe = VQE(qubitOp, UCCSD_var_form, optimizer, callback=store_intermediate_result)
vqe_eigenvalue = np.real(vqe.run(backend)['eigenvalue'])
# Now compare the results of different compositions of your VQE algorithm!
print("VQE value: ", vqe_eigenvalue)
vqe_result = np.real(energy_shift + vqe_eigenvalue)
print("VQE energy: ", vqe_result)
print("HF energy: {}".format(molecule.hf_energy))
print("Reference Value Percent Error: " + str(abs((vqe_eigenvalue-ref)/ref)*100) + "%")
print("Energy Value Percent Error: " + str(np.real(abs((vqe_result-molecule.hf_energy)/molecule.hf_energy))*100) + "%")

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

