from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import ParityMapper
from qiskit_nature.second_q.properties import ParticleNumber
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer

bond_distance = 2.5  # in Angstrom

# specify driver
driver = PySCFDriver(
    atom=f"Li 0 0 0; H 0 0 {bond_distance}",
    basis="sto3g",
    charge=0,
    spin=0,
    unit=DistanceUnit.ANGSTROM,
)
problem = driver.run()

# specify active space transformation
active_space_trafo = ActiveSpaceTransformer(
    num_electrons=problem.num_particles, num_spatial_orbitals=3
)

# transform the electronic structure problem
problem = active_space_trafo.transform(problem)

# construct the parity mapper with 2-qubit reduction
qubit_mapper = ParityMapper(num_particles=problem.num_particles)

from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver
from qiskit_nature.second_q.algorithms.ground_state_solvers import GroundStateEigensolver

np_solver = NumPyMinimumEigensolver()
np_groundstate_solver = GroundStateEigensolver(qubit_mapper, np_solver)

np_result = np_groundstate_solver.solve(problem)

target_energy = np_result.total_energies[0]
print(np_result)

from qiskit.circuit.library import EfficientSU2

ansatz = EfficientSU2(num_qubits=4, reps=1, entanglement="linear", insert_barriers=True)
ansatz.decompose().draw("mpl", style="iqx")

import numpy as np
from qiskit.utils import algorithm_globals

# fix random seeds for reproducibility

np.random.seed(5)
algorithm_globals.random_seed = 5

from qiskit.algorithms.optimizers import SPSA

optimizer = SPSA(maxiter=100)

initial_point = np.random.random(ansatz.num_parameters)

from qiskit.algorithms.minimum_eigensolvers import VQE
from qiskit.primitives import Estimator

estimator = Estimator()

local_vqe = VQE(
    estimator,
    ansatz,
    optimizer,
    initial_point=initial_point,
)

local_vqe_groundstate_solver = GroundStateEigensolver(qubit_mapper, local_vqe)

local_vqe_result = local_vqe_groundstate_solver.solve(problem)

print(local_vqe_result)

from qiskit import IBMQ

IBMQ.load_account()
provider = IBMQ.get_provider(group="open")  # replace by your runtime provider

backend = provider.get_backend("ibmq_qasm_simulator")  # select a backend that supports the runtime

from qiskit_nature.runtime import VQEClient

runtime_vqe = VQEClient(
    ansatz=ansatz,
    optimizer=optimizer,
    initial_point=initial_point,
    provider=provider,
    backend=backend,
    shots=1024,
    measurement_error_mitigation=True,
)  # use a complete measurement fitter for error mitigation

runtime_vqe_groundstate_solver = GroundStateEigensolver(qubit_mapper, runtime_vqe)
runtime_vqe_result = runtime_vqe_groundstate_solver.solve(problem)

print(runtime_vqe_result)

runtime_result = runtime_vqe_result.raw_result
history = runtime_result.optimizer_history
loss = history["energy"]

import matplotlib.pyplot as plt

plt.rcParams["font.size"] = 14

# plot loss and reference value
plt.figure(figsize=(12, 6))
plt.plot(loss + runtime_vqe_result.nuclear_repulsion_energy, label="Runtime VQE")
plt.axhline(y=target_energy + 0.2, color="tab:red", ls=":", label="Target + 200mH")
plt.axhline(y=target_energy, color="tab:red", ls="--", label="Target")

plt.legend(loc="best")
plt.xlabel("Iteration")
plt.ylabel("Energy [H]")
plt.title("VQE energy");

import qiskit.tools.jupyter

%qiskit_version_table
%qiskit_copyright