Wyznaczanie geometrii molekularnej

W poprzedniej sekcji zaimplementowaliśmy VQE w celu wyznaczenia energii stanu podstawowego molecule. Jest to prawidłowe zastosowanie obliczeń kwantowych, ale jeszcze bardziej użyteczne byłoby wyznaczenie struktury molecule.

Krok 1: Odwzorowanie klasycznych danych wejściowych na problem kwantowy

Pozostając przy naszym podstawowym przykładzie dwuatomowego wodoru, jedynym parametrem geometry, który można zmieniać, jest długość wiązania. Aby to osiągnąć, postępujemy tak jak poprzednio, ale używając zmiennej w początkowej konstrukcji molecule (długość wiązania, x, jako argument). Jest to dość prosta zmiana, ale wymaga, aby zmienna była uwzględniona w funkcjach w całym procesie, ponieważ zaczyna się w konstrukcji fermionowego Hamiltonian i propaguje przez odwzorowanie, a w końcu do funkcji kosztu.

Najpierw ładujemy niektóre z pakietów, których używaliśmy wcześniej, i definiujemy funkcję Cholesky'ego.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy pyscf qiskit qiskit-aer qiskit-ibm-runtime scipy

from qiskit.quantum_info import SparsePauliOp
import matplotlib.pyplot as plt
import numpy as np

#!pip install pyscf==2.4.0
from pyscf import ao2mo, gto, mcscf, scf

def cholesky(V, eps):
    # see https://arxiv.org/pdf/1711.02242.pdf section B2
    # see https://arxiv.org/abs/1808.02625
    # see https://arxiv.org/abs/2104.08957
    no = V.shape[0]
    chmax, ng = 20 * no, 0
    W = V.reshape(no**2, no**2)
    L = np.zeros((no**2, chmax))
    Dmax = np.diagonal(W).copy()
    nu_max = np.argmax(Dmax)
    vmax = Dmax[nu_max]
    while vmax > eps:
        L[:, ng] = W[:, nu_max]
        if ng > 0:
            L[:, ng] -= np.dot(L[:, 0:ng], (L.T)[0:ng, nu_max])
        L[:, ng] /= np.sqrt(vmax)
        Dmax[: no**2] -= L[: no**2, ng] ** 2
        ng += 1
        nu_max = np.argmax(Dmax)
        vmax = Dmax[nu_max]
    L = L[:, :ng].reshape((no, no, ng))
    print(
        "accuracy of Cholesky decomposition ",
        np.abs(np.einsum("prg,qsg->prqs", L, L) - V).max(),
    )
    return L, ng

def identity(n):
    return SparsePauliOp.from_list([("I" * n, 1)])

def creators_destructors(n, mapping="jordan_wigner"):
    c_list = []
    if mapping == "jordan_wigner":
        for p in range(n):
            if p == 0:
                ell, r = "I" * (n - 1), ""
            elif p == n - 1:
                ell, r = "", "Z" * (n - 1)
            else:
                ell, r = "I" * (n - p - 1), "Z" * p
            cp = SparsePauliOp.from_list([(ell + "X" + r, 0.5), (ell + "Y" + r, -0.5j)])
            c_list.append(cp)
    else:
        raise ValueError("Unsupported mapping.")
    d_list = [cp.adjoint() for cp in c_list]
    return c_list, d_list

Aby teraz zdefiniować nasz Hamiltonian, użyjemy PySCF dokładnie tak samo jak w poprzednim przykładzie, ale teraz uwzględnimy zmienną x, która będzie pełnić rolę naszej odległości międzyatomowej. Funkcja ta zwróci, tak jak wcześniej, energię rdzenia, energię jednoelektronową oraz energie dwuelektronowe.

def ham_terms(x: float):
    distance = x
    a = distance / 2
    mol = gto.Mole()
    mol.build(
        verbose=0,
        atom=[
            ["H", (0, 0, -a)],
            ["H", (0, 0, a)],
        ],
        basis="sto-6g",
        spin=0,
        charge=0,
        symmetry="Dooh",
    )

    # mf = scf.RHF(mol)
    # mx = mcscf.CASCI(mf, ncas=2, nelecas=(1, 1))
    # mx.kernel()

    mf = scf.RHF(mol)
    mf.kernel()
    if not mf.converged:
        raise RuntimeError(f"SCF did not converge for distance {x}")

    mx = mcscf.CASCI(mf, ncas=2, nelecas=(1, 1))
    casci_energy = mx.kernel()
    if casci_energy is None:
        raise RuntimeError(f"CASCI failed for distance {x}")

    # Other variables that might come in handy:
    # active_space = range(mol.nelectron // 2 - 1, mol.nelectron // 2 + 1)
    #    E1 = mf.kernel()
    # mo = mx.sort_mo(active_space, base=0)
    #    E2 = mx.kernel(mo)[:2]

    h1e, ecore = mx.get_h1eff()
    h2e = ao2mo.restore(1, mx.get_h2eff(), mx.ncas)
    return ecore, h1e, h2e

Przypomnijmy, że powyższa konstrukcja tworzy fermionowy Hamiltonian na podstawie rodzajów atomów, geometry i orbital elektronowych. Poniżej odwzorowujemy ten fermionowy Hamiltonian na operatory Pauliego. Funkcja build_hamiltonian również będzie zawierać zmienną geometryczną jako argument.

def build_hamiltonian(distx: float) -> SparsePauliOp:
    ecore = ham_terms(distx)[0]
    h1e = ham_terms(distx)[1]
    h2e = ham_terms(distx)[2]

    ncas, _ = h1e.shape

    C, D = creators_destructors(2 * ncas, mapping="jordan_wigner")
    Exc = []
    for p in range(ncas):
        Excp = [C[p] @ D[p] + C[ncas + p] @ D[ncas + p]]
        for r in range(p + 1, ncas):
            Excp.append(
                C[p] @ D[r]
                + C[ncas + p] @ D[ncas + r]
                + C[r] @ D[p]
                + C[ncas + r] @ D[ncas + p]
            )
        Exc.append(Excp)

    # low-rank decomposition of the Hamiltonian
    Lop, ng = cholesky(h2e, 1e-6)
    t1e = h1e - 0.5 * np.einsum("pxxr->pr", h2e)

    H = ecore * identity(2 * ncas)
    # one-body term
    for p in range(ncas):
        for r in range(p, ncas):
            H += t1e[p, r] * Exc[p][r - p]
    # two-body term
    for g in range(ng):
        Lg = 0 * identity(2 * ncas)
        for p in range(ncas):
            for r in range(p, ncas):
                Lg += Lop[p, r, g] * Exc[p][r - p]
        H += 0.5 * Lg @ Lg

    return H.chop().simplify()

Załadujemy pozostałe pakiety potrzebne do uruchomienia samego VQE, takie jak ansatz efficient_su2 i minimalizatory SciPy:

# General imports

# Pre-defined ansatz circuit and operator class for Hamiltonian
from qiskit.circuit.library import efficient_su2
from qiskit.quantum_info import SparsePauliOp

# SciPy minimizer routine
from scipy.optimize import minimize

# Plotting functions

# Qiskit Runtime tools
from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()

Ponownie zdefiniujemy funkcję kosztu, ale zawsze przyjmowała ona jako argument w pełni zbudowany i odwzorowany Hamiltonian, więc w tej funkcji nic się nie zmienia.

def cost_func(params, ansatz, H, estimator):
    pub = (ansatz, [H], [params])
    result = estimator.run(pubs=[pub]).result()
    energy = result[0].data.evs[0]
    return energy

# def cost_func_sim(params, ansatz, H, estimator):
#    energy = estimator.run(ansatz, H, parameter_values=params).result().values[0]
#    return energy

Krok 2: Optymalizacja problemu pod wykonanie kwantowe

Ponieważ Hamiltonian będzie się zmieniał przy każdej nowej geometrii, transpilacja operatora będzie się zmieniać na każdym kroku. Możemy jednak zdefiniować ogólny pass manager, który będzie stosowany na każdym kroku, specyficzny dla sprzętu, którego chcemy użyć.

Tutaj użyjemy najmniej obciążonego dostępnego backendu. Wykorzystamy ten backend jako model dla naszego AerSimulator, pozwalając naszemu symulatorowi naśladować, na przykład, zachowanie szumowe rzeczywistego backendu. Te modele szumu nie są idealne, ale mogą pomóc zorientować się, czego można oczekiwać po rzeczywistym sprzęcie.

# Here, we select the least busy backend available:
backend = service.least_busy(operational=True, simulator=False)
print(backend)
# Or to select a specific real backend use the line below, and substitute 'ibm_strasbourg' for your chosen device.
# backend = service.get_backend('ibm_strasbourg')

# To run on a simulator:
# -----------
from qiskit_aer import AerSimulator

backend_sim = AerSimulator.from_backend(backend)

Importujemy pass manager oraz powiązane pakiety, aby pomóc nam zoptymalizować nasz obwód. Ten krok, oraz ten powyżej, są niezależne od Hamiltonianu, więc pozostają niezmienione w stosunku do poprzedniej lekcji.

from qiskit.transpiler import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
    ConstrainedReschedule,
)
from qiskit.circuit.library import XGate

target = backend.target
pm = generate_preset_pass_manager(target=target, optimization_level=3)
pm.scheduling = PassManager(
    [
        ALAPScheduleAnalysis(target=target),
        ConstrainedReschedule(
            acquire_alignment=target.acquire_alignment,
            pulse_alignment=target.pulse_alignment,
            target=target,
        ),
        PadDynamicalDecoupling(
            target=target,
            dd_sequence=[XGate(), XGate()],
            pulse_alignment=target.pulse_alignment,
        ),
    ]
)

Krok 3: Wykonanie przy użyciu prymitywów Qiskit.

W poniższym bloku kodu konfigurujemy tablicę do przechowywania wyników z każdego kroku w odległości międzyatomowej $x$. Wybraliśmy zakres $x$ na podstawie naszej wiedzy o eksperymentalnej wartości równowagowej długości wiązania: 0,74 Angstrema. Uruchomimy to najpierw na symulatorze, więc zaimportujemy nasz estimator (BackendEstimator) z qiskit.primitives. Dla każdego kroku geometrii budujemy Hamiltonian i pozwalamy na określoną liczbę kroków optymalizacji (tutaj 500) przy użyciu optymalizatora „cobyla". W każdym kroku geometrii przechowujemy zarówno całkowitą energię, jak i energię elektronową. Ze względu na dużą liczbę kroków optymalizatora może to zająć godzinę lub więcej. Możesz chcieć zmodyfikować poniższe dane wejściowe, aby skrócić wymagany czas.

from qiskit.primitives import BackendEstimatorV2

estimator = BackendEstimatorV2(backend=backend_sim)

distances_sim = np.arange(0.3, 1.3, 0.1)
vqe_energies_sim = []
vqe_elec_energies_sim = []

for dist in distances_sim:
    xx = dist

    # Random initial state and efficient_su2 ansatz
    H = build_hamiltonian(xx)
    ansatz = efficient_su2(H.num_qubits)
    ansatz_isa = pm.run(ansatz)
    x0 = 2 * np.pi * np.random.random(ansatz_isa.num_parameters)
    H_isa = H.apply_layout(ansatz_isa.layout)
    nuclear_repulsion = ham_terms(xx)[0]

    res = minimize(
        cost_func,
        x0,
        args=(ansatz_isa, H_isa, estimator),
        method="cobyla",
        options={"maxiter": 20, "disp": True},
    )

    # Note this returns the total energy, and we are often interested in the electronic energy
    tot_energy = getattr(res, "fun")
    electron_energy = getattr(res, "fun") - nuclear_repulsion
    print(electron_energy)
    vqe_energies_sim.append(tot_energy)
    vqe_elec_energies_sim.append(electron_energy)

    # Print all results
    print(res)

print("All energies have been calculated")

accuracy of Cholesky decomposition  1.1102230246251565e-15

/home/porter284/.pyenv/versions/3.11.12/lib/python3.11/site-packages/scipy/_lib/pyprima/common/preproc.py:68: UserWarning: COBYLA: Invalid MAXFUN; it should be at least num_vars + 2; it is set to 34
  warn(f'{solver}: Invalid MAXFUN; it should be at least {min_maxfun_str}; it is set to {maxfun}')

Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = 1.316011435623847
The corresponding X is:
[2.32948769 5.39918229 3.03787975 4.11789904 4.97130735 2.68662232
 1.76573151 2.48982571 5.40431972 3.65780829 1.33792786 5.48472494
 6.18738702 1.78741883 0.78195251 2.96658955 1.35827677 5.599321
 4.54850148 1.0939048  4.26158726 0.52100721 0.82318    4.76796961
 3.75795507 3.8526447  5.51100375 5.91023075 2.61494836 1.79908918
 2.65937756 5.53964148]

-0.44791260077615314
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: 1.316011435623847
       x: [ 2.329e+00  5.399e+00 ...  2.659e+00  5.540e+00]
    nfev: 34
   maxcv: 0.0
accuracy of Cholesky decomposition  5.551115123125783e-16
Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = 0.7235003672327549
The corresponding X is:
[2.56282915 5.63369524 5.58059887 4.049643   4.2021266  3.06866011
 6.01619635 1.52520776 4.35403161 0.33673958 0.32623161 1.2179545
 2.84001371 3.98956684 4.89632562 1.38303588 1.96194695 2.13182089
 0.29739166 1.77895165 3.29151585 3.54355374 4.49626674 0.95756626
 0.87103927 4.53068385 1.31051302 0.37103108 1.02961355 3.13342311
 5.65815319 2.24770604]

-0.5994426600672451
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: 0.7235003672327549
       x: [ 2.563e+00  5.634e+00 ...  5.658e+00  2.248e+00]
    nfev: 34
   maxcv: 0.0
accuracy of Cholesky decomposition  5.551115123125783e-16
Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = 0.34960914928810116
The corresponding X is:
[5.44143165 6.75955835 1.56836472 3.09522093 4.67873235 1.67071481
 0.3056494  0.65998337 1.02197668 5.21162959 0.43690354 3.56522934
 4.56033119 1.90736037 0.40863891 2.87007312 3.2516952  5.90360196
 1.99057799 5.20726456 0.74710237 6.03179202 3.80685028 0.03844391
 5.88580196 3.62233258 3.98723567 2.50591888 5.44020267 2.2792993
 5.57102303 4.46548617]

-0.7087452725518989
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: 0.34960914928810116
       x: [ 5.441e+00  6.760e+00 ...  5.571e+00  4.465e+00]
    nfev: 34
   maxcv: 0.0
accuracy of Cholesky decomposition  2.220446049250313e-16
Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = 0.10594558882184543
The corresponding X is:
[5.35675483 2.26629567 1.45430546 5.56758296 5.76309509 0.73239338
 5.1216998  3.03258872 4.33624828 1.93197674 0.5292902  3.32274987
 3.43247633 0.81490741 0.48060245 1.9944799  5.67519646 5.12534057
 0.06510627 2.52989834 6.1699519  0.94828957 5.91634548 1.5994961
 4.27902164 2.3129213  1.82353095 2.10634209 1.43740426 4.06988733
 0.59624074 4.93925418]

-0.7760164293781545
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: 0.10594558882184543
       x: [ 5.357e+00  2.266e+00 ...  5.962e-01  4.939e+00]
    nfev: 34
   maxcv: 0.0
accuracy of Cholesky decomposition  1.1102230246251565e-16
Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = -0.06473600797229297
The corresponding X is:
[6.07735568 0.18019501 0.20743128 4.15445985 3.59388894 5.10047555
 6.09938474 6.54707528 3.36251167 2.05475223 3.67078456 5.96010605
 2.58589996 5.2723619  3.26352977 2.47432334 3.50289983 2.06620525
 6.0946056  1.22751903 0.97320057 2.19564095 5.73174941 2.05127682
 5.73805165 3.84046105 1.84816963 2.1247504  3.11106736 2.44136052
 3.39002685 0.81596991]

-0.8207034521437214
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: -0.06473600797229297
       x: [ 6.077e+00  1.802e-01 ...  3.390e+00  8.160e-01]
    nfev: 34
   maxcv: 0.0
accuracy of Cholesky decomposition  5.551115123125783e-17
Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = -0.19562982094782935
The corresponding X is:
[-0.02184462  3.67041038  7.25918653  5.89799546  0.63583624  1.84214506
  2.84059837  5.31485182  1.6053784   0.04556618  0.32018993 -0.03884066
  0.69131496  0.24203727  1.97397262  3.59723495  0.43355775  2.30131056
  4.63482292  3.9857415   4.32320753  4.55388437  2.18753433  5.99034987
  2.50489913  0.90650534  4.82518088  2.32954849  2.29901832  5.33658863
  5.91246716  3.2405013 ]

-0.8571013345978292
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: -0.19562982094782935
       x: [-2.184e-02  3.670e+00 ...  5.912e+00  3.241e+00]
    nfev: 34
   maxcv: 0.0
accuracy of Cholesky decomposition  1.1102230246251565e-16
Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = -0.2833766309947055
The corresponding X is:
[ 3.1700088   5.05055456  1.2545611   4.28751811  0.6255103   1.67526577
  5.48201473  4.83820497  7.34880059  5.99705431  4.2502643   0.32066274
  0.41001404  0.27271241  4.15682546  4.22393693  4.35148115  0.64538137
  5.26288622  5.03810489  4.62426621  4.74997689  1.09603919  0.34752466
  1.8116275   0.7474807   5.31754143  4.11181763  1.58797998  5.6299796
  3.0109383  -0.19062772]

-0.8713513097947054
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: -0.2833766309947055
       x: [ 3.170e+00  5.051e+00 ...  3.011e+00 -1.906e-01]
    nfev: 34
   maxcv: 0.0
accuracy of Cholesky decomposition  1.1102230246251565e-16
Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = -0.3527503628484244
The corresponding X is:
[3.90513622 4.61398739 5.92552705 1.99953405 4.82157369 1.35702441
 2.77701782 5.73612247 4.22710527 1.83463189 0.45796297 4.62509318
 0.98998668 0.11666217 3.0234641  4.54298546 0.14034033 4.15635797
 1.41257357 4.48719602 2.39365535 0.19672041 5.0763044  1.86357581
 3.657757   4.60298344 2.49769577 1.88086199 3.00108725 1.84475841
 5.24047385 4.91142914]

-0.8819275737684243
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: -0.3527503628484244
       x: [ 3.905e+00  4.614e+00 ...  5.240e+00  4.911e+00]
    nfev: 34
   maxcv: 0.0
accuracy of Cholesky decomposition  2.7755575615628914e-17
Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = -0.4022181851996095
The corresponding X is:
[6.09453981 3.5109422  3.37216019 4.94732621 1.25662002 5.89645164
 5.06403334 2.68073141 4.40385083 1.13638366 1.73347762 6.82932871
 1.15265014 2.07145964 4.36520459 1.14960341 1.62288871 4.32315915
 5.45622821 0.93554005 3.17418483 0.47230243 1.31535502 5.77698726
 2.04927925 2.50663538 5.9706002  5.4984681  2.9421232  1.56636313
 1.09394523 4.62582   ]

-0.8832883769450639
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: -0.4022181851996095
       x: [ 6.095e+00  3.511e+00 ...  1.094e+00  4.626e+00]
    nfev: 34
   maxcv: 0.0
accuracy of Cholesky decomposition  1.1102230246251565e-16
Return from COBYLA because the objective function has been evaluated MAXFUN times.
Number of function values = 34   Least value of F = -0.44423031870708934
The corresponding X is:
[4.05765050e+00 3.99144950e+00 3.13287593e+00 3.28855137e+00
 4.32613515e+00 4.91104512e+00 1.86521867e+00 2.18822879e+00
 6.01336171e+00 1.82501276e+00 2.64830637e+00 5.53045823e+00
 2.36110093e+00 3.98821703e+00 4.69013438e-01 4.38996815e+00
 7.78103801e-04 1.72994378e+00 2.24970934e+00 1.11978200e+00
 2.24846445e+00 4.90745512e+00 5.38474921e+00 5.03587994e+00
 3.54297277e+00 4.78147533e+00 1.25990218e+00 1.99168068e+00
 5.89203503e+00 1.77673987e+00 5.37848357e+00 5.60245198e-01]

-0.8852113278070892
 message: Return from COBYLA because the objective function has been evaluated MAXFUN times.
 success: False
  status: 3
     fun: -0.44423031870708934
       x: [ 4.058e+00  3.991e+00 ...  5.378e+00  5.602e-01]
    nfev: 34
   maxcv: 0.0
All energies have been calculated

xx

np.float64(1.2000000000000004)

Wyniki tego wyjścia zostały omówione poniżej w sekcji dotyczącej postprocessingu; na razie po prostu zauważ, że symulacja zakończyła się pomyślnie. Teraz jesteś gotowy do uruchomienia na prawdziwym sprzęcie. Ustawimy poziom odporności (resilience) na 1, co oznacza, że zostanie użyta mitygacja błędów TREX. Ponieważ pracujemy teraz z prawdziwym sprzętem, użyjemy Qiskit Runtime oraz prymitywów Runtime. Zauważ, że zarówno pętla for związana z geometrią, jak i wielokrotne próby wariacyjne znajdują się wewnątrz sesji.

Ponieważ z uruchomieniami na prawdziwym sprzęcie wiążą się koszty i limity czasowe, zmniejszyliśmy poniżej liczbę kroków geometrii oraz kroków optymalizatora. Dopasuj te kroki do swoich celów dokładnościowych i limitów czasowych.

# To continue running on real hardware use
from qiskit_ibm_runtime import Session
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import EstimatorOptions

estimator_options = EstimatorOptions(resilience_level=1, default_shots=2000)

distances = np.arange(0.5, 0.9, 0.1)
vqe_energies = []
vqe_elec_energies = []

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    for dist in distances:
        xx = dist

        # Random initial state and efficient_su2 ansatz

        H = build_hamiltonian(xx)
        ansatz = efficient_su2(H.num_qubits)
        ansatz_isa = pm.run(ansatz)
        H_isa = H.apply_layout(ansatz_isa.layout)
        nuclear_repulsion = ham_terms(xx)[0]
        x0 = 2 * np.pi * np.random.random(ansatz_isa.num_parameters)

        res = minimize(
            cost_func,
            x0,
            args=(ansatz_isa, H_isa, estimator),
            method="cobyla",
            options={"maxiter": 50, "disp": True},
        )

        # Note this returns the total energy, and we are often interested in the electronic energy
        tot_energy = getattr(res, "fun")
        electron_energy = getattr(res, "fun") - nuclear_repulsion
        print(electron_energy)
        vqe_energies.append(tot_energy)
        vqe_elec_energies.append(electron_energy)

        # Print all results
        print(res)

print("All energies have been calculated")

Krok 4: Przetwarzanie końcowe

Zarówno w przypadku symulatora, jak i rzeczywistego sprzętu, możemy wykreślić energie stanu podstawowego obliczone dla każdej odległości międzyatomowej i zobaczyć, gdzie osiągnięta jest najniższa energia. Powinna to być odległość międzyatomowa występująca w naturze, i rzeczywiście jest ona zbliżona. Gładszą krzywą można uzyskać, próbując innych ansatzów, optymalizatorów oraz wielokrotnie uruchamiając obliczenia dla każdego kroku geometrii i uśredniając wyniki dla kilku losowych warunków początkowych.

# Here we can plot the results from this simulation.
plt.plot(distances_sim, vqe_energies_sim, label="VQE Energy")
plt.xlabel("Atomic distance (Angstrom)")
plt.ylabel("Energy")
plt.legend()
plt.show()

Warto zauważyć, że samo zwiększenie liczby kroków optymalizacji prawdopodobnie nie poprawi wyników z symulatora, ponieważ wszystkie optymalizacje faktycznie osiągnęły zbieżność z wymaganą tolerancją w mniejszej liczbie iteracji niż maksymalna.

Wyniki z rzeczywistego sprzętu są porównywalne, pomijając nieco inny zakres próbkowanych wartości.

plt.plot(distances, vqe_energies, label="VQE Energy")
plt.xlabel("Atomic distance (Angstrom)")
plt.ylabel("Energy")
plt.legend()
plt.show()

Oprócz oczekiwanej długości wiązania H2 wynoszącej 0,74 Angstrema, całkowita energia powinna wynosić -1,17 Hartree. Widzimy, że wyniki z rzeczywistego sprzętu były bliższe tym wartościom niż te z symulatora. Wynika to prawdopodobnie z tego, że w obu przypadkach obecny (lub symulowany) był szum, ale tylko w przypadku rzeczywistego sprzętu zastosowano łagodzenie błędów.

Zakończenie

To kończy nasz kurs dotyczący VQE dla chemii kwantowej. Jeśli jesteś zainteresowany zrozumieniem niektórych fundamentów teorii informacji wykorzystywanej w obliczeniach kwantowych, zapoznaj się z kursem Johna Watrousa Podstawy informacji kwantowej. Dodatkowy krótki przykład przepływu pracy VQE znajdziesz w naszym samouczku dotyczącym estymacji energii stanu podstawowego łańcucha Heisenberga z użyciem VQE. Możesz też przeglądać samouczki oraz kursy, aby znaleźć więcej materiałów edukacyjnych na temat najnowszych technologii w obliczeniach kwantowych.

Nie zapomnij przystąpić do egzaminu kończącego ten kurs. Wynik 80% lub wyższy zapewni Ci odznakę Credly, która zostanie automatycznie wysłana do Ciebie e-mailem. Dziękujemy za bycie częścią sieci IBM Quantum® Network!

import qiskit
import qiskit_ibm_runtime

print(qiskit.version.get_version_info())
print(qiskit_ibm_runtime.version.get_version_info())

1.3.2
0.35.0