NVIDIA CUDA-Q チュートリアル＆H100ベンチマーク７：スピンハミルトニアンシミュレーション

このチュートリアルでは、1次元スピン鎖（スピンチェーン）に対する量子シミュレーションを CUDA-Q を使って実行する方法を解説します。特に、以下の2つの代表的なスピンハミルトニアンに焦点を当てます：

ハイゼンベルク模型（Heisenberg Hamiltonian）

H_{\text{Heisenberg}} = \sum_{i=1}^{N-1} \left( J_x X_i X_{i+1} + J_y Y_i Y_{i+1} + J_z Z_i Z_{i+1} \right)

横磁場イジング模型（Transverse Field Ising Model: TFIM）

H_{\text{TFIM}} = -J \sum_{i=1}^{N-1} Z_i Z_{i+1} - h \sum_{i=1}^{N} X_i

目的

これらのスピンモデルに基づく量子系の時間発展（Time Evolution）をシミュレートする
Trotter-Suzuki分解を用いて、複雑なハミルトニアンの指数関数的な演算を近似的に計算

実際に実行

早速行っていきます。まずはライブラリの読み込み

# Import Required Libraries
import cudaq
import numpy as np
import time
import sys
from cudaq import spin
import matplotlib.pyplot as plt
from typing import List

次に初期状態の準備。今回は01010101のような状態を想定して準備しています。

@cudaq.kernel
def get_initial_state(n_spins: int):
    """Create initial state |1010...>"""
    qubits = cudaq.qvector(n_spins)
    for i in range(0, n_spins, 2):
        x(qubits[i])

次の関数はハミルトニアンを1stepずつトロッター展開して時間発展を実行するものです。

@cudaq.kernel
def trotter_step(state: cudaq.State, dt: float, Jx: float, Jy: float, Jz: float,
                 h_x: list[float], h_y: list[float], h_z: list[float], _use_XXYYZZ_gate: bool,
                 coefficients: List[complex], words: List[cudaq.pauli_word]):
    """Perform single Trotter step"""
    qubits = cudaq.qvector(state)
    n_spins = len(qubits)

    # Apply two-qubit interaction terms
    if _use_XXYYZZ_gate:
        for j in range(2):
            for i in range(j % 2, n_spins - 1, 2):
                rx(-np.pi/2,qubits[i])
                rx(np.pi/2,qubits[i+1])
                x.ctrl(qubits[i], qubits[i+1])
                h(qubits[i])
                s(qubits[i])
                rz(-2*Jy*dt,qubits[i+1])
                x.ctrl(qubits[i], qubits[i+1])
                h(qubits[i])
                rx(2*Jx*dt,qubits[i])
                rz(-2*Jz*dt,qubits[i+1])
                x.ctrl(qubits[i], qubits[i+1])
    else:
        for i in range(len(coefficients)):
            exp_pauli(coefficients[i].real * dt, qubits, words[i])

こちらは初期状態とのオーバーラップを計算する関数

def compute_overlap_probability(initial_state: cudaq.State, evolved_state: cudaq.State):
    """Compute probability of the overlap with the initial state"""
    overlap = initial_state.overlap(evolved_state)
    return np.abs(overlap)**2

次にハイゼンベルグモデルのハミルトニアンを作るための関数。

def create_hamiltonian_heisenberg(n_spins: int, Jx: float, Jy: float, Jz: float, h_x: list[float], h_y: list[float], h_z: list[float]):
    """Create the Hamiltonian operator"""
    ham = cudaq.SpinOperator(num_qubits=n_spins)

    # Add two-qubit interaction terms for Heisenberg Hamiltonian
    for i in range(0, n_spins - 1):
        ham += Jx * spin.x(i) * spin.x(i + 1)
        ham += Jy * spin.y(i) * spin.y(i + 1)
        ham += Jz * spin.z(i) * spin.z(i + 1)

    return ham

こちらは横磁場イジングモデルを作るもの。

def create_hamiltonian_tfim(n_spins: int, h_field: float):
    """Create the Hamiltonian operator"""
    ham = cudaq.SpinOperator(num_qubits=n_spins)

    # Add single-qubit terms
    for i in range(0, n_spins):
        ham += -1 * h_field * spin.x(i)

    # Add two-qubit interaction terms for Ising Hamiltonian
    for i in range(0, n_spins-1):
        ham += -1 * spin.z(i) * spin.z(i + 1)

    return ham

トロッタ展開する時に使うパウリ演算子と相互作用を取り出すための関数

def extractCoefficients(hamiltonian: cudaq.SpinOperator):
    result = []
    hamiltonian.for_each_term(lambda term: result.append(term.get_coefficient()))
    return result

def extractWords(hamiltonian: cudaq.SpinOperator):
    result = []
    hamiltonian.for_each_term(lambda term: result.append(term.to_string(False)))
    return result

こちらがメインの実行で、パラメータを決めて先ほどのものを実行します。

# Import Required Libraries
import cudaq
import numpy as np
import time
import sys
from cudaq import spin
import matplotlib.pyplot as plt
from typing import List

# Parameters
n_spins = 4  # Number of spins in the chain
ham_type = "heisenberg"  # Choose between "heisenberg" and "tfim"
Jx, Jy, Jz = 1.0, 1.0, 1.0  # Coupling coefficients for Heisenberg Hamiltonian
h_field = 1.0  # Transverse field strength for TFIM
K = 100  # Number of Trotter steps
t = np.pi  # Total evolution time
dt = t / K  # Time step size

# Optimized XXYYZZ exponentiation. Works only for Heisenberg Hamiltonian
_use_XXYYZZ_gate = False
if _use_XXYYZZ_gate == True and ham_type == "tfim":
    print ("XXYYZZ exponentiation works only for Heisenberg")
    sys.exit(0)

# Create Hamiltonian
if ham_type == "heisenberg":
    # Initialize field for Heisenberg Hamiltonian
    h_x = np.ones(n_spins)
    h_y = np.ones(n_spins)
    h_z = np.ones(n_spins)
    hamiltonian = create_hamiltonian_heisenberg(n_spins, Jx, Jy, Jz, h_x, h_y,h_z)
elif ham_type == "tfim":
    hamiltonian = create_hamiltonian_tfim(n_spins, h_field)
else:
    raise ValueError("Invalid Hamiltonian type. Choose 'heisenberg' or 'tfim'.")

# Extract coefficients and words
coefficients = extractCoefficients(hamiltonian)
words = extractWords(hamiltonian)

# Initialize and save the initial state
print ("Initialize state")
initial_state = cudaq.get_state(get_initial_state, n_spins)
state = initial_state

# Store probabilities over time
probabilities = []
probabilities.append(1.0)

# Time evolution
start_time = time.time()
for k in range(1, K):

    # Apply single Trotter step
    state = cudaq.get_state(trotter_step, state, dt, Jx, Jy, Jz, h_x, h_y, h_z, _use_XXYYZZ_gate, coefficients, words)

    # Calculate probability between initial and current states
    probability = compute_overlap_probability(initial_state, state)
    probabilities.append(probability)

total_time = time.time() - start_time
print(f"Circuit execution time: {total_time:.3f} seconds")

H100での実行時間です。
Circuit execution time: 0.074 seconds

# Plot probability over time
plt.plot(range(K), probabilities, marker='o')
plt.xlabel('Step (k)')
plt.ylabel('Probability of initial state')
plt.title('Probability vs. Step in Trotter Evolution')
plt.grid()
plt.show()

そして、次は初期状態との重なりではなく、期待値を求めるもの。

# Initialize list to store expectation values
exp_values = []

# Time evolution
start_time = time.time()
for k in range(1, K):
    # Apply single Trotter step
    state = cudaq.get_state(trotter_step, state, dt, Jx, Jy, Jz, h_x, h_y, h_z, _use_XXYYZZ_gate, coefficients, words)

    # Calculate expectation value
    exp_val = cudaq.observe(trotter_step, hamiltonian, state, dt, Jx, Jy, Jz, h_x, h_y, h_z, _use_XXYYZZ_gate, coefficients, words).expectation()
    exp_values.append(exp_val.real)
    #print(f"Step {k}, Energy: {exp_val.real:.6f}")

# Plot expectation value over time
x = np.arange(1, K)
plt.plot(x, exp_values, marker='o')
plt.xlabel('Step (k)')
plt.ylabel('Expectation Value (Energy)')
plt.title('Expectation Value vs. Step in Trotter Evolution')
plt.grid()
plt.show()

ベンチマーク

CPUの方は時間がかかりすぎて全然計算できないのでGPUの方で初期状態とのオーバーラップのベンチマークを取ります。

トロッターステップは100に固定して、ハイゼンベルグモデルとTFIMの両方をスピンの長さを4から28スピンまで4スピンずつ増やして実行しています。

NVIDIA CUDA-Q チュートリアル＆H100ベンチマーク７： スピンハミルトニアンシミュレーション

Yuichiro Minato

目的

実際に実行

ベンチマーク

NVIDIA CUDA-Q チュートリアル＆H100ベンチマーク７：スピンハミルトニアンシミュレーション