State-preparation examples
Preparing the input register a primitive will consume. These examples build
chemistry-style states with the composable stateprep device kernels, each one
checked against a dense reference.
See also
State preparation — the ansatz kernels, operator pools, and injection contract behind these examples.
Hartree-Fock and fixed-parameter UCCSD state preparation
Prepares a 4-qubit / 2-electron UCCSD ansatz state at fixed amplitudes
(H2-style). It shows how the Hartree-Fock reference and the fixed-parameter
UCCSD singles and doubles compose into a single state_prep kernel.
"""Hartree-Fock + fixed-parameter UCCSD state preparation.
Prepares a 4-qubit / 2-electron UCCSD ansatz state at fixed amplitudes
(H2-style; in a real workflow the parameters would come from a classical
pre-optimization) and validates it against dense matrix exponentials of
the operator pool. Also demonstrates the open-shell Hartree-Fock
occupation, whose interleaved alpha/beta layout matches the UCCSD
excitation enumeration.
Run with: python3 hartree_fock_ucc.py
"""
import os
import numpy as np
from scipy.linalg import expm
import cudaq
from cudaq_algorithms import stateprep
NUM_QUBITS = 4
NUM_ELECTRONS = 2
PARAMETERS = [0.1129, -0.0421, 0.2839] # fixed UCCSD amplitudes
PAULI_MATRICES = {
"I": np.eye(2),
"X": np.array([[0.0, 1.0], [1.0, 0.0]]),
"Y": np.array([[0.0, -1.0j], [1.0j, 0.0]]),
"Z": np.array([[1.0, 0.0], [0.0, -1.0]]),
}
def dense_reference(pool, parameters, ket, num_qubits):
"""prod_g exp(i theta_g G_g) |ket> from the pool's Pauli terms."""
state = np.array(ket, dtype=np.complex128)
for theta, op in zip(parameters, pool):
generator = np.zeros((1 << num_qubits, 1 << num_qubits),
dtype=np.complex128)
for term in cudaq.SpinOperator(op):
word = str(term.get_pauli_word(num_qubits))
matrix = np.array([[1.0]], dtype=np.complex128)
# Qubit 0 is the least-significant statevector index.
for label in reversed(word):
matrix = np.kron(matrix, PAULI_MATRICES[label])
generator += term.evaluate_coefficient().real * matrix
state = expm(1.0j * theta * generator) @ state
return state
def main():
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
pool = stateprep.make_uccsd_operator_pool(NUM_QUBITS, NUM_ELECTRONS)
words, coeffs = stateprep.get_fixed_parameter_ucc_pauli_lists(
pool, NUM_QUBITS)
resources = stateprep.estimate_fixed_parameter_ucc_resources(
NUM_QUBITS, words)
prep = stateprep.hartree_fock_ucc_kernel(NUM_QUBITS,
PARAMETERS,
words,
coeffs,
num_electrons=NUM_ELECTRONS)
@cudaq.kernel
def entry():
qubits = cudaq.qvector(NUM_QUBITS)
prep(qubits)
state = np.array(cudaq.get_state(entry))
hf_ket = np.zeros(1 << NUM_QUBITS, dtype=np.complex128)
hf_ket[0b0011] = 1.0
expected = dense_reference(pool, PARAMETERS, hf_ket, NUM_QUBITS)
error = float(np.linalg.norm(state - expected))
print("Fixed-parameter UCCSD state preparation")
print("=" * 64)
print(f"qubits: {NUM_QUBITS}")
print(f"electrons: {NUM_ELECTRONS}")
print(f"excitation groups: {resources.num_excitations}")
print(f"pauli rotations: {resources.num_pauli_rotations}")
print(f"max rotations/group: "
f"{resources.max_pauli_rotations_per_excitation}")
print(f"L2 error vs dense: {error:.3e}")
print("dominant amplitudes:")
for index in np.argsort(np.abs(state))[::-1][:4]:
print(f" |{index:04b}> {state[index].real:+.6f}"
f"{state[index].imag:+.6f}j "
f"(P = {abs(state[index])**2:.4f})")
print()
print("Open-shell Hartree-Fock (8 qubits, 4 electrons, spin 2):")
occupation = stateprep.make_hartree_fock_occupation(8, 4, 2)
print(f" occupied spin orbitals: {occupation} (interleaved "
"alpha/beta, NOT the contiguous {0, 1, 2, 3})")
if error > 1e-6:
raise SystemExit("UCC state does not match the dense reference")
print("PASS")
return 0
if __name__ == "__main__":
raise SystemExit(main())
Givens-rotation Slater determinant preparation
Prepares a real 4-orbital / 2-electron determinant and a complex 5-orbital /
3-electron determinant with the composable Givens-rotation stateprep kernels.
Both are validated against their dense Slater-determinant references.
"""Givens-rotation Slater determinant preparation.
Prepares a real 4-orbital / 2-electron determinant and a complex
5-orbital / 3-electron determinant with the composable ``stateprep``
kernels, and checks both against the dense reference (all minors of the
orbital-coefficient matrix).
Run with: python3 givens_slater_determinant.py
"""
import os
import numpy as np
import cudaq
from cudaq_algorithms import stateprep
@cudaq.kernel
def prepare_real(num_spin_orbitals: int, orbital_indices: list[int],
angles: list[float], num_electrons: int):
qubits = cudaq.qvector(num_spin_orbitals)
stateprep.slater_determinant(qubits, orbital_indices, angles,
num_electrons)
@cudaq.kernel
def prepare_complex(num_spin_orbitals: int, orbital_indices: list[int],
angles: list[float], phases: list[float],
final_phases: list[float], num_electrons: int):
qubits = cudaq.qvector(num_spin_orbitals)
stateprep.complex_slater_determinant(qubits, orbital_indices, angles,
phases, final_phases, num_electrons)
def reference_slater_state(orbital_coefficients):
orbital_coefficients = np.asarray(orbital_coefficients, dtype=complex)
num_spin_orbitals, num_electrons = orbital_coefficients.shape
state = np.zeros(2**num_spin_orbitals, dtype=complex)
for basis_index in range(2**num_spin_orbitals):
occupied = [
orbital for orbital in range(num_spin_orbitals)
if (basis_index >> orbital) & 1
]
if len(occupied) != num_electrons:
continue
state[basis_index] = np.linalg.det(orbital_coefficients[np.ix_(
occupied, range(num_electrons))])
return state
def phase_aligned_l2(actual, expected):
actual = np.asarray(actual, dtype=complex)
expected = np.asarray(expected, dtype=complex)
pivot = int(np.argmax(np.abs(expected)))
phase = 1.0
if abs(expected[pivot]) > 1.0e-14:
phase = actual[pivot] / expected[pivot]
phase /= abs(phase)
return np.linalg.norm(actual - phase * expected)
def run_case(label, orbital_coefficients):
schedule = stateprep.make_givens_rotation_schedule(orbital_coefficients)
resources = stateprep.estimate_givens_resources(schedule)
indices = stateprep.get_givens_rotation_indices(schedule)
angles = stateprep.get_givens_rotation_angles(schedule)
if schedule.is_complex:
state = cudaq.get_state(prepare_complex, schedule.num_spin_orbitals,
indices, angles,
stateprep.get_givens_rotation_phases(schedule),
list(schedule.final_phases),
schedule.num_electrons)
else:
state = cudaq.get_state(prepare_real, schedule.num_spin_orbitals,
indices, angles, schedule.num_electrons)
error = phase_aligned_l2(np.asarray(state),
reference_slater_state(orbital_coefficients))
print(label)
print(f" spin orbitals: {schedule.num_spin_orbitals}")
print(f" electrons: {schedule.num_electrons}")
print(f" complex: {schedule.is_complex}")
print(f" Givens rotations: {resources.num_givens_rotations}")
print(f" exp_pauli calls: {resources.num_exp_pauli_calls}")
print(f" phase rotations: {resources.num_phase_rotations}")
print(f" phase-aligned L2 error: {error:.3e}")
if error > 1.0e-6:
raise SystemExit(f"{label}: prepared state does not match the "
"Slater determinant")
def main():
cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
rng = np.random.default_rng(11)
real_orbitals, _ = np.linalg.qr(rng.normal(size=(4, 2)))
run_case("real orbital-coefficient matrix (4 spin orbitals, 2 electrons)",
real_orbitals)
rng = np.random.default_rng(13)
raw = rng.normal(size=(5, 3)) + 1j * rng.normal(size=(5, 3))
complex_orbitals, _ = np.linalg.qr(raw)
run_case(
"complex orbital-coefficient matrix (5 spin orbitals, 3 electrons)",
complex_orbitals)
print("PASS")
return 0
if __name__ == "__main__":
raise SystemExit(main())