from qiskit_nature.second_q.drivers import PySCFDriver
driver = PySCFDriver(
    atom="H 0 0 0; H 0 0 0.72"  # Two Hydrogen atoms, 0.72 Angstrom apart
)
molecule = driver.run()

from qiskit_nature.second_q.mappers import QubitConverter, ParityMapper
qubit_converter = QubitConverter(ParityMapper())
hamiltonian = qubit_converter.convert(molecule.second_q_ops()[0])

from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver
sol = NumPyMinimumEigensolver().compute_minimum_eigenvalue(hamiltonian)
real_solution = molecule.interpret(sol)

real_solution.groundenergy

from qiskit_ibm_runtime import QiskitRuntimeService, Estimator, Session, Options

service = QiskitRuntimeService()

backend = "ibmq_qasm_simulator"

from qiskit.algorithms.minimum_eigensolvers import VQE

# Use RealAmplitudes circuit to create trial states
from qiskit.circuit.library import RealAmplitudes
ansatz = RealAmplitudes(num_qubits=2, reps=2)

# Search for better states using SPSA algorithm
from qiskit.algorithms.optimizers import SPSA
optimizer = SPSA(150)

# Set a starting point for reproduceability
import numpy as np
np.random.seed(6)
initial_point = np.random.uniform(-np.pi, np.pi, 12)

# Create an object to store intermediate results
from dataclasses import dataclass
@dataclass
class VQELog:
    values: list
    parameters: list
    def update(self, count, parameters, mean, _metadata):
        self.values.append(mean)
        self.parameters.append(parameters)
        print(f"Running circuit {count} of ~350", end="\r", flush=True)

log = VQELog([],[])

# Main calculation
with Session(service=service, backend=backend) as session:
    options = Options()
    options.optimization_level = 3

    vqe = VQE(Estimator(session=session, options=options),
              ansatz, optimizer, callback=log.update, initial_point=initial_point)
    result = vqe.compute_minimum_eigenvalue(hamiltonian)
    print("Experiment complete.".ljust(30))
    print(f"Raw result: {result.optimal_value}")
    
    if 'simulator' not in backend:
        # Run once with ZNE error mitigation
        options.resilience_level = 2
        vqe = VQE(Estimator(session=session, options=options),
                  ansatz, SPSA(1), initial_point=result.optimal_point)
        result = vqe.compute_minimum_eigenvalue(hamiltonian)
        print(f"Mitigated result: {result.optimal_value}")

import matplotlib.pyplot as plt

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

# Plot energy and reference value
plt.figure(figsize=(12, 6))
plt.plot(log.values, label="Estimator VQE")
plt.axhline(y=real_solution.groundenergy, color="tab:red", ls="--", label="Target")

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

import qiskit_ibm_runtime
qiskit_ibm_runtime.version.get_version_info()

import qiskit.tools.jupyter

%qiskit_version_table
%qiskit_copyright