Model Hamiltonians
QDK/Chemistry provides functionality to construct and manipulate model Hamiltonians used in quantum chemistry and condensed matter physics. These model Hamiltonians serve as simplified representations of complex quantum systems to study their properties and behaviors using quantum computing techniques.
Unlike molecular Hamiltonians, model Hamiltonians do not require a molecular structure or precomputed integrals. They are defined directly in terms of their parameters and a LatticeGraph that specifies the site connectivity.
Overview
QDK/Chemistry supports two families of model Hamiltonians:
- Fermionic models
Operate on fermionic degrees of freedom (creation and annihilation operators) and produce Hamiltonian objects that are compatible with all QDK/Chemistry algorithms.
Hückel — tight-binding model with one-body hopping only
Hubbard — extends Hückel with on-site electron-electron repulsion
Pariser-Parr-Pople (PPP) — extends Hubbard with long-range intersite Coulomb interactions
- Spin models
Operate on spin-½ degrees of freedom and produce
QubitHamiltonianobjects expressed as sums of Pauli operators.Heisenberg — anisotropic spin-spin coupling with external magnetic fields
Ising — special case of Heisenberg with ZZ coupling and transverse X field
All model Hamiltonian builders take a LatticeGraph as their first argument, which defines the site connectivity and hopping structure. For a brief description of the available model Hamiltonian builders, see the table below. For a more detailed description of each model Hamiltonian and their parameters, see the following sections.
Builder function |
Type |
Output |
Description |
|---|---|---|---|
|
Fermionic |
Hamiltonian |
Tight-binding with hopping only |
|
Fermionic |
Hamiltonian |
Hopping + on-site Coulomb repulsion |
|
Fermionic |
Hamiltonian |
Hubbard + long-range Coulomb interactions |
|
Spin |
QubitHamiltonian |
Anisotropic spin coupling + external fields |
|
Spin |
QubitHamiltonian |
ZZ coupling + transverse X field |
Fermionic models
Hückel model
The Hückel (tight-binding) model describes non-interacting electrons hopping on a lattice:
where \(\hat{a}_i^\dagger\) and \(\hat{a}_i\) are the fermionic creation and annihilation operators for site i, \(\hat{n}_i = \sum_\sigma \hat{a}_{i,\sigma}^\dagger \hat{a}_{i,\sigma}\) is the number operator, \(\varepsilon_i\) are on-site energies, \(t_{ij}\) are hopping integrals, \(w_{ij}\) is the edge weight from the lattice adjacency matrix, and the sum runs over connected site pairs. This model produces a Hamiltonian with one-body integrals only.
// Create a 6-site chain for the Hückel model
auto lattice = LatticeGraph::chain(6);
// Uniform parameters: all sites have the same on-site energy and hopping
auto hamiltonian =
create_huckel_hamiltonian(lattice, /*epsilon=*/0.0, /*t=*/1.0);
bool has_h1 = hamiltonian.has_one_body_integrals(); // true
bool has_h2 = hamiltonian.has_two_body_integrals(); // false
# Create a 6-site chain for the Hückel model
lattice = LatticeGraph.chain(6)
# Uniform parameters: all sites have the same on-site energy and hopping
hamiltonian = create_huckel_hamiltonian(lattice, epsilon=0.0, t=1.0)
# Print the one-body integrals
print(f"Has one-body integrals: {hamiltonian.has_one_body_integrals()}")
print(f"Has two-body integrals: {hamiltonian.has_two_body_integrals()}")
Hubbard model
The Hubbard model extends the Hückel model with on-site Coulomb repulsion:
where \(U_i\) is the on-site repulsion strength. This model produces a Hamiltonian with both one-body and two-body integrals.
// Create a 4-site Hubbard chain
auto hub_lattice = LatticeGraph::chain(4);
auto hub_hamiltonian = create_hubbard_hamiltonian(
hub_lattice, /*epsilon=*/0.0, /*t=*/1.0, /*U=*/4.0);
bool hub_has_h1 = hub_hamiltonian.has_one_body_integrals(); // true
bool hub_has_h2 = hub_hamiltonian.has_two_body_integrals(); // true
// Verify the on-site repulsion
for (int i = 0; i < hub_lattice.num_sites(); ++i) {
double Uii = hub_hamiltonian.get_two_body_element(i, i, i, i); // 4.0
}
# Create a 4-site Hubbard chain
lattice = LatticeGraph.chain(4)
hamiltonian = create_hubbard_hamiltonian(lattice, epsilon=0.0, t=1.0, U=4.0)
print(f"Has one-body integrals: {hamiltonian.has_one_body_integrals()}")
print(f"Has two-body integrals: {hamiltonian.has_two_body_integrals()}")
# Verify the on-site repulsion
for i in range(lattice.num_sites):
print(f" U({i},{i},{i},{i}) = {hamiltonian.get_two_body_element(i, i, i, i):.1f}")
The Hubbard model naturally extends to 2D lattices for studying strongly correlated materials:
// Create a 4x4 square lattice Hubbard model with periodic boundaries
auto square_lattice =
LatticeGraph::square(4, 4, /*periodic_x=*/true, /*periodic_y=*/true);
auto square_hamiltonian = create_hubbard_hamiltonian(
square_lattice, /*epsilon=*/0.0, /*t=*/1.0, /*U=*/4.0);
bool sq_has_h2 = square_hamiltonian.has_two_body_integrals(); // true
# Create a 4x4 square lattice Hubbard model with periodic boundaries
lattice = LatticeGraph.square(4, 4, periodic_x=True, periodic_y=True)
hamiltonian = create_hubbard_hamiltonian(lattice, epsilon=0.0, t=1.0, U=4.0)
print(f"Has two-body integrals: {hamiltonian.has_two_body_integrals()}")
Pariser-Parr-Pople (PPP) model
The PPP model extends the Hubbard model with long-range intersite Coulomb interactions:
where \(V_{ij}\) is the intersite Coulomb repulsion and \(z_i\) are effective core charges. The intersite potential \(V_{ij}\) is typically computed using the Ohno or Mataga-Nishimoto parametrizations (see Intersite potentials below).
Note
The stored two-body integrals do not include the \(\frac{1}{2}\) prefactor. This follows the standard quantum chemistry convention where the factor is applied at contraction time rather than stored in the integrals.
// Create a 6-site ring for the PPP model (benzene-like)
auto ring = LatticeGraph::chain(6, /*periodic=*/true);
// Compute interpair Coulomb repulsion with the Ohno potential
auto V = ohno_potential(ring, /*U=*/0.414, /*R=*/2.65);
// Create the PPP Hamiltonian
auto ppp_hamiltonian = create_ppp_hamiltonian(ring,
/*epsilon=*/0.0,
/*t=*/0.088,
/*U=*/0.414,
/*V=*/V,
/*z=*/1.0);
bool ppp_has_h1 = ppp_hamiltonian.has_one_body_integrals(); // true
bool ppp_has_h2 = ppp_hamiltonian.has_two_body_integrals(); // true
double core_e = ppp_hamiltonian.get_core_energy(); // > 0
# Create a 6-site ring for the PPP model (benzene-like)
lattice = LatticeGraph.chain(6, periodic=True)
# Compute interpair Coulomb repulsion with the Ohno potential
V = ohno_potential(lattice, U=0.414, R=2.65)
# Create the PPP Hamiltonian
hamiltonian = create_ppp_hamiltonian(
lattice,
epsilon=0.0,
t=0.088,
U=0.414,
V=V,
z=1.0,
)
print(f"Has one-body integrals: {hamiltonian.has_one_body_integrals()}")
print(f"Has two-body integrals: {hamiltonian.has_two_body_integrals()}")
print(f"Core energy: {hamiltonian.get_core_energy():.6f}")
Intersite potentials
For the PPP model, the intersite Coulomb interaction \(V_{ij}\) is typically computed from a distance-dependent parametrization. QDK/Chemistry provides two standard potentials and a custom potential interface.
By default, all potential functions compute \(V_{ij}\) for every pair of sites, not just lattice-connected neighbours.
This is consistent with the PPP Hamiltonian, which sums the Coulomb term over all pairs \(i \ne j\).
All three potential functions accept an optional nearest_neighbor_only flag (default false) that restricts the evaluation to lattice-connected pairs only, setting \(V_{ij} = 0\) for non-adjacent sites.
Ohno potential
where \(U_{ij} = \sqrt{U_i U_j}\) is the geometric mean of on-site parameters, \(R_{ij}\) is the intersite distance, and \(\varepsilon_r\) is the relative permittivity.
Mataga-Nishimoto potential
Custom pairwise potential
The pairwise_potential function accepts a user-defined callable func(i, j, U_ij, R_ij) -> V_ij for arbitrary distance-dependent potentials.
auto pot_lattice = LatticeGraph::chain(4);
// Ohno potential: V_ij = U_ij / sqrt(1 + (U_ij * epsilon_r * R_ij)^2)
auto V_ohno =
ohno_potential(pot_lattice, /*U=*/0.414, /*R=*/2.65, /*epsilon_r=*/1.0);
// Mataga-Nishimoto potential: V_ij = U_ij / (1 + U_ij * epsilon_r * R_ij)
auto V_mn = mataga_nishimoto_potential(pot_lattice, /*U=*/0.414, /*R=*/2.65,
/*epsilon_r=*/1.0);
// Custom pairwise potential using a user-defined function
auto V_custom = pairwise_potential(
pot_lattice,
/*U=*/0.414,
/*R=*/2.65,
[](int, int, double Uij, double Rij) { return Uij / (1.0 + Rij); });
lattice = LatticeGraph.chain(4)
# Ohno potential: V_ij = U_ij / sqrt(1 + (U_ij * epsilon_r * R_ij)^2)
V_ohno = ohno_potential(lattice, U=0.414, R=2.65, epsilon_r=1.0)
# Mataga-Nishimoto potential: V_ij = U_ij / (1 + U_ij * epsilon_r * R_ij)
V_mn = mataga_nishimoto_potential(lattice, U=0.414, R=2.65, epsilon_r=1.0)
# Custom pairwise potential using a user-defined function
V_custom = pairwise_potential(
lattice,
U=0.414,
R=2.65,
func=lambda i, j, Uij, Rij: Uij / (1.0 + Rij),
)
print(f"Ohno V(0,1) = {V_ohno[0, 1]:.4f}")
print(f"Mataga-Nishimoto V(0,1) = {V_mn[0, 1]:.4f}")
print(f"Custom V(0,1) = {V_custom[0, 1]:.4f}")
Spin models
Heisenberg model
The anisotropic Heisenberg model describes spin-½ particles interacting on a lattice with external magnetic fields:
where \(J_x, J_y, J_z\) are the spin-spin coupling constants, \(h_x, h_y, h_z\) are external magnetic field components, and \(w_{ij}\) is the edge weight from the lattice adjacency matrix. Each qubit corresponds to a lattice site.
# Create an anisotropic Heisenberg model on a square lattice
lattice = LatticeGraph.square(3, 3)
qubit_hamiltonian = create_heisenberg_hamiltonian(
lattice,
jx=1.0,
jy=1.0,
jz=1.0,
hz=0.5, # External longitudinal field
)
print(f"Heisenberg Hamiltonian ({lattice.num_sites} qubits):")
print(f" Number of Pauli terms: {len(qubit_hamiltonian.pauli_strings)}")
print(f" Is Hermitian: {qubit_hamiltonian.is_hermitian()}")
Special cases of the Heisenberg model include:
Isotropic (XXX): \(J_x = J_y = J_z\)
XXZ: \(J_x = J_y \ne J_z\)
XY: \(J_z = 0\)
Ising model
The transverse-field Ising model is a special case of the Heisenberg model with ZZ coupling and a transverse X field:
# Create a transverse-field Ising model on a chain
lattice = LatticeGraph.chain(6)
qubit_hamiltonian = create_ising_hamiltonian(lattice, j=1.0, h=0.5)
print(f"Ising Hamiltonian ({lattice.num_sites} qubits):")
print(f" Number of Pauli terms: {len(qubit_hamiltonian.pauli_strings)}")
print(f" Is Hermitian: {qubit_hamiltonian.is_hermitian()}")
Parameter flexibility
All model Hamiltonian builders accept parameters as either scalars (applied uniformly to all sites or pairs) or arrays (specifying per-site or per-pair values). This allows modelling inhomogeneous systems such as impurities or spatially varying fields.
- Per-site parameters
Scalar
float(broadcast to all sites) ornumpy.ndarrayof length n (one value per site). Used for: on-site energy (\(\varepsilon\)), on-site repulsion (\(U\)), core charges (\(z\)), magnetic fields (\(h_x, h_y, h_z\)).- Per-pair parameters
Scalar
float(broadcast to all pairs) or(n, n)numpy.ndarray(one value per pair). Used for: hopping (\(t\)), intersite potential (\(V\)), spin couplings (\(J_x, J_y, J_z\)).
// Site-dependent Hubbard model: different U on each site
auto chain = LatticeGraph::chain(4);
Eigen::VectorXd U_values(4);
U_values << 2.0, 4.0, 4.0, 2.0; // Weaker repulsion at edges
auto site_dep_hamiltonian = create_hubbard_hamiltonian(
chain, /*epsilon=*/0.0, /*t=*/1.0, /*U=*/U_values);
# Site-dependent Hubbard model: different U on each site
lattice = LatticeGraph.chain(4)
U_values = np.array([2.0, 4.0, 4.0, 2.0]) # Weaker repulsion at edges
hamiltonian = create_hubbard_hamiltonian(lattice, epsilon=0.0, t=1.0, U=U_values)
# Bond-dependent Ising model: different coupling on each bond
j_matrix = np.ones((4, 4))
j_matrix[0, 1] = 2.0 # Stronger coupling on first bond
h_fields = np.array([0.5, 0.3, 0.3, 0.5]) # Inhomogeneous field
qh = create_ising_hamiltonian(lattice, j=j_matrix, h=h_fields)
Using model Hamiltonians with algorithms
Fermionic model Hamiltonians produce Hamiltonian objects that are fully compatible with all QDK/Chemistry algorithms, including:
Multi-configuration calculators (FCI, ASCI, etc.)
Qubit mapping (Jordan-Wigner, Bravyi-Kitaev, etc.)
Phase estimation (IQPE, standard QPE)
Spin model Hamiltonians produce QubitHamiltonian objects directly, which can be used with quantum algorithms without an intermediate qubit mapping step.
Example: exact diagonalization of the Hubbard model
// Create a 4-site half-filled Hubbard chain
auto solve_lattice = LatticeGraph::chain(4);
auto solve_hamiltonian = std::make_shared<Hamiltonian>(
create_hubbard_hamiltonian(solve_lattice,
/*epsilon=*/0.0, /*t=*/1.0, /*U=*/4.0));
// Run exact diagonalization (CASCI) with half filling (2 alpha + 2 beta)
auto mc_calculator =
qdk::chemistry::algorithms::MultiConfigurationCalculatorFactory::create(
"macis_cas");
auto [energy, wavefunction] = mc_calculator->run(solve_hamiltonian, 2, 2);
# Create a 4-site half-filled Hubbard chain
lattice = LatticeGraph.chain(4)
hamiltonian = create_hubbard_hamiltonian(lattice, epsilon=0.0, t=1.0, U=4.0)
# Run exact diagonalization (CASCI) with half filling (2 alpha + 2 beta electrons)
mc_calculator = algorithms.create("multi_configuration_calculator", "macis_cas")
energy, wavefunction = mc_calculator.run(hamiltonian, 2, 2)
print(f"Ground state energy: {energy:.6f} a.u.")
Example: exact diagonalization of the Ising model
# Create a transverse-field Ising model on a 6-site chain
lattice = LatticeGraph.chain(6)
qubit_hamiltonian = create_ising_hamiltonian(lattice, j=1.0, h=0.5)
# Solve for the ground state energy using exact diagonalization
q_solver = algorithms.create("qubit_hamiltonian_solver", "qdk_sparse_matrix_solver")
energy, ground_state = q_solver.run(qubit_hamiltonian)
print(f"Ising ground state energy: {energy:.6f}")
print(f"Number of qubits: {qubit_hamiltonian.num_qubits}")
Further reading
The above examples can be downloaded as complete C++ and Python scripts.
Multi-configuration calculations — Solving fermionic model Hamiltonians with exact diagonalization
QubitHamiltonianSolver— Solving spin model Hamiltonians with exact diagonalization