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 QubitHamiltonian objects 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

create_huckel_hamiltonian

Fermionic

Hamiltonian

Tight-binding with hopping only

create_hubbard_hamiltonian

Fermionic

Hamiltonian

Hopping + on-site Coulomb repulsion

create_ppp_hamiltonian

Fermionic

Hamiltonian

Hubbard + long-range Coulomb interactions

create_heisenberg_hamiltonian

Spin

QubitHamiltonian

Anisotropic spin coupling + external fields

create_ising_hamiltonian

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:

\[\hat{H}_\text{Hückel} = \sum_i \varepsilon_i\, \hat{n}_i - \sum_{\langle i,j \rangle} t_{ij}\, w_{ij}\, (\hat{a}_i^\dagger \hat{a}_j + \hat{a}_j^\dagger \hat{a}_i)\]

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:

\[\hat{H}_\text{Hubbard} = \hat{H}_\text{Hückel} + \sum_i U_i\, \hat{n}_{i\uparrow} \hat{n}_{i\downarrow}\]

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:

\[\hat{H}_\text{PPP} = \hat{H}_\text{Hubbard} + \frac{1}{2} \sum_{i \ne j} V_{ij}\, (\hat{n}_i - z_i)(\hat{n}_j - z_j)\]

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
\[V_{ij} = \frac{U_{ij}}{\sqrt{1 + \left(U_{ij}\,\varepsilon_r\,R_{ij}\right)^2}}\]

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
\[V_{ij} = \frac{U_{ij}}{1 + U_{ij}\,\varepsilon_r\,R_{ij}}\]
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:

\[\hat{H}_\text{Heisenberg} = \sum_{\langle i,j \rangle} w_{ij}\,\bigl( J_x^{ij}\,\hat{\sigma}_i^x \hat{\sigma}_j^x + J_y^{ij}\,\hat{\sigma}_i^y \hat{\sigma}_j^y + J_z^{ij}\,\hat{\sigma}_i^z \hat{\sigma}_j^z \bigr) + \sum_i \bigl( h_x^{i}\,\hat{\sigma}_i^x + h_y^{i}\,\hat{\sigma}_i^y + h_z^{i}\,\hat{\sigma}_i^z \bigr)\]

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:

\[\hat{H}_\text{Ising} = \sum_{\langle i,j \rangle} w_{ij}\,J^{ij}\,\hat{\sigma}_i^z \hat{\sigma}_j^z + \sum_i h^{i}\,\hat{\sigma}_i^x\]
# 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) or numpy.ndarray of 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:

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