Chemistry¶
In chemistry applications, the electronic Hamiltonian written as a linear combination of Pauli strings is the principal object.
This a setting in which the TermSum types in Zixy are very useful.
The TermSum is instantiable from a sparse string representation of the Pauli strings and coefficients.
import numpy as np
from zixy.qubit.pauli import RealTermSum as RealPauliOperator
H2_STO3G_HAM_JW_INPUT = "\
(-0.05962058276034765, ), (0.1757594291831968, Z0), (0.17575942918319679, Z1), (0.17001546439603182, Z0 Z1), \
(0.044917169890753894, X0 Y1 Y2 X3), (-0.044917169890753894, X0 X1 Y2 Y3), (-0.044917169890753894, Y0 Y1 X2 X3), \
(0.044917169890753894, Y0 X1 X2 Y3), (-0.23667117678035537, Z2), (0.12222714936261826, Z0 Z2), (0.16714431925337217, Z1 Z2), \
(-0.23667117678035543, Z3), (0.16714431925337217, Z0 Z3), (0.12222714936261826, Z1 Z3), (0.1757033833190701, Z2 Z3) \
"
H2_STO3G_HF_ENERGY = -1.1175058842043306
H2_STO3G_FCI_ENERGY = -1.136846575472054
ham_op = RealPauliOperator.from_str(H2_STO3G_HAM_JW_INPUT)
Alternatively we can use fermionic encodings to construct the qubit Pauli operators
from zixy.fermion.mappings import JordanWignerMapper
jw = JordanWignerMapper(4, mode_ordering=None)
tmp = RealPauliOperator(4)
tmp += jw.encode_ccaa(2, 3, 0, 1) * 8
tmp += jw.encode_caca(0, 1, 2, 3) * 1.2
tmp -= jw.encode_ca(2, 1) * 2.3
tmp -= jw.encode_n(2) * 2
tmp += jw.encode_nn(2, 3) * 2
tmp
assert tmp.conserves_hamming_weight()
The Rust implementation constructs temporary number operators to detect symmetries via commutation:
ham_op.conserves_hamming_weight(), ham_op.conserves_odd_bit_hamming_weight()
(True, True)
It’s convenient to be able to view the terms in a pandas DataFrame
ham_op.to_dataframe()
| Component | Coefficient | |
|---|---|---|
| 0 | -0.05962058276034765 | |
| 1 | Z0 | 0.1757594291831968 |
| 2 | Z1 | 0.17575942918319679 |
| 3 | Z0 Z1 | 0.17001546439603182 |
| 4 | X0 Y1 Y2 X3 | 0.044917169890753894 |
| 5 | X0 X1 Y2 Y3 | -0.044917169890753894 |
| 6 | Y0 Y1 X2 X3 | -0.044917169890753894 |
| 7 | Y0 X1 X2 Y3 | 0.044917169890753894 |
| 8 | Z2 | -0.23667117678035537 |
| 9 | Z0 Z2 | 0.12222714936261826 |
| 10 | Z1 Z2 | 0.16714431925337217 |
| 11 | Z3 | -0.23667117678035543 |
| 12 | Z0 Z3 | 0.16714431925337217 |
| 13 | Z1 Z3 | 0.12222714936261826 |
| 14 | Z2 Z3 | 0.1757033833190701 |
The terms of a list can be sorted in various ways. Either by Pauli string (integer value of the X part, then that of the Z part):
terms = ham_op.to_terms()
terms.lexicographic_sort(ascending=False)
terms.to_dataframe()
| Component | Coefficient | |
|---|---|---|
| 0 | X0 X1 Y2 Y3 | -0.044917169890753894 |
| 1 | Y0 X1 X2 Y3 | 0.044917169890753894 |
| 2 | X0 Y1 Y2 X3 | 0.044917169890753894 |
| 3 | Y0 Y1 X2 X3 | -0.044917169890753894 |
| 4 | Z2 Z3 | 0.1757033833190701 |
| 5 | Z1 Z3 | 0.12222714936261826 |
| 6 | Z0 Z3 | 0.16714431925337217 |
| 7 | Z3 | -0.23667117678035543 |
| 8 | Z1 Z2 | 0.16714431925337217 |
| 9 | Z0 Z2 | 0.12222714936261826 |
| 10 | Z2 | -0.23667117678035537 |
| 11 | Z0 Z1 | 0.17001546439603182 |
| 12 | Z1 | 0.17575942918319679 |
| 13 | Z0 | 0.1757594291831968 |
| 14 | -0.05962058276034765 |
Or numerically:
terms.numeric_sort(ascending=False, by_magnitude=True)
terms.to_dataframe()
| Component | Coefficient | |
|---|---|---|
| 0 | Z3 | -0.23667117678035543 |
| 1 | Z2 | -0.23667117678035537 |
| 2 | Z0 | 0.1757594291831968 |
| 3 | Z1 | 0.17575942918319679 |
| 4 | Z2 Z3 | 0.1757033833190701 |
| 5 | Z0 Z1 | 0.17001546439603182 |
| 6 | Z1 Z2 | 0.16714431925337217 |
| 7 | Z0 Z3 | 0.16714431925337217 |
| 8 | Z0 Z2 | 0.12222714936261826 |
| 9 | Z1 Z3 | 0.12222714936261826 |
| 10 | -0.05962058276034765 | |
| 11 | Y0 Y1 X2 X3 | -0.044917169890753894 |
| 12 | X0 Y1 Y2 X3 | 0.044917169890753894 |
| 13 | Y0 X1 X2 Y3 | 0.044917169890753894 |
| 14 | X0 X1 Y2 Y3 | -0.044917169890753894 |
Linear combinations of Pauli strings can be converted to sparse matrix representation in either little- or big-endian ordering.
ham_mat = ham_op.to_sparse_matrix(False).toarray()
assert np.isclose(np.linalg.eigh(ham_mat)[0][0], H2_STO3G_FCI_ENERGY)
Generally, the full computational basis is too large for enumeration and energy estimation, so it is useful to have access to subspace-projected matrices.
In this example, only the basis vector corresponding to the Hartree-Fock state and its double excitation are non-zero in the ground state wavefunction.
from zixy.qubit.state import Strings as StateStrings
subspace = StateStrings.from_iterable(([1, 1, 0, 0], [0, 0, 1, 1]), ham_op.qubits)
ham_mat = ham_op.project_into_ortho_subspace(subspace)
ham_mat
array([[-1.11750588, 0.17966868],
[ 0.17966868, 0.53221654]])
As such, the [1, 1, 0, 0], [0, 0, 1, 1] subspace is sufficient to obtain the correct correlation energy
assert np.isclose(ham_mat[0, 0], H2_STO3G_HF_ENERGY)
assert np.isclose(np.linalg.eigh(ham_mat)[0][0], H2_STO3G_FCI_ENERGY)
The subspace is checked for orthonormality before projection of the operator into it is performed.
subspace.append([1, 1, 0, 0])
ham_mat = ham_op.project_into_ortho_subspace(subspace)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[11], line 2
1 subspace.append([1, 1, 0, 0])
----> 2 ham_mat = ham_op.project_into_ortho_subspace(subspace)
File ~/Desktop/zixy/zixy-py/docs/.venv/lib/python3.14/site-packages/zixy/qubit/pauli/_terms.py:655, in RealTermSum.project_into_ortho_subspace(self, subspace)
646 """Project the term into an orthogonal subspace defined by an ordered sequence of states.
647
648 Args:
(...) 652 The resulting vector of coefficients in the subspace basis.
653 """
654 assert isinstance(self._impl._coeffs, RealCoeffs)
--> 655 return QubitPauliArray.lincomb_project_into_ortho_subspace_real(
656 self._impl._cmpnts._impl,
657 self._impl._coeffs._impl,
658 subspace._impl,
659 )
ValueError: The subspace is non-orthogonal (it contains repeated basis vectors)
Sometimes, non-orthonormality of the basis is a desirable property, and Zixy also has functionality to handle these cases.
from zixy.qubit.state import RealTermSum as RealState
subspace = [
RealState.from_iterable([
([1, 1, 0, 0], 0.4), ([0, 0, 1, 1], -0.5)
], ham_op.qubits),
RealState.from_iterable([
([1, 1, 0, 0], 0.4), ([0, 0, 1, 1], 0.9)
], ham_op.qubits)
]
ham_mat, ovlp_mat = ham_op.project_into_nonortho_subspace(subspace)
The metric of the subspace is returned along with the projected operator.
ovlp_mat
array([[ 0.41, -0.29],
[-0.29, 0.97]])
Note that the two linear combinations of basis vectors have the same span as the orthonormal subspace does, so we expect the correct correlation energy here too.
import scipy
assert np.allclose(scipy.linalg.eigh(ham_mat, ovlp_mat)[0][0], H2_STO3G_FCI_ENERGY)
State TermSum instances can be added together just as those of the Pauli operator TermSum type can.
state_sum = sum(subspace, start = RealState(ham_op.qubits))
assert state_sum == RealState.from_iterable([([1, 1, 0, 0], 0.8), ([0, 0, 1, 1], 0.4)], ham_op.qubits)
Of course, it is standard to unit-normalize the states
state_sum.l2_normalize()
assert np.isclose(state_sum.l2_norm, 1.0)
assert np.isclose(state_sum.l2_norm_square, 1.0)
Additionally, we have functionality to handle products of operators.
Note that even pairs of hermitian Pauli operators (those with real coefficients) multiply to give complex results (when the two are non-commuting). Rather than going to the additional cost of checking for commutivity to deduce the return type, we leave it to the user to check and discard the imaginary part where appropriate.
ham2_op = ham_op * ham_op
assert not ham2_op.imag_part.to_terms().coeffs.any_significant(0.0)
ham2_op = ham2_op.real_part
assert ham_op.commutes_with(ham2_op)
ham2_mat = ham2_op.to_sparse_matrix(False)
assert np.isclose(scipy.sparse.linalg.eigsh(ham2_mat)[0][0], H2_STO3G_FCI_ENERGY**2)
ham2_mat = ham2_mat.toarray()
ham_mat = ham_op.to_sparse_matrix(False).toarray()
assert np.allclose(ham2_mat, ham_mat @ ham_mat)
States are able to be created from dense vectors in little- or big-endian ordering
evals, evecs = np.linalg.eigh(ham_mat)
ground_state = RealState.from_dense(ham_op.qubits, evecs[:, 0].real, False)
ground_state
(-0.9942559956062316, [1, 1, 0, 0]), (0.10702810472516666, [0, 0, 1, 1])
There are methods to apply operators on states and take the inner product of states.
ham_times_ground_state = ham_op.apply(ground_state)
assert not ham_times_ground_state.imag_part.to_terms().coeffs.any_significant(0.0)
ham_times_ground_state = ham_times_ground_state.real_part
assert np.isclose(ground_state.vdot(ham_times_ground_state), H2_STO3G_FCI_ENERGY)
Expectation values and matrix elements are also available.
assert np.isclose(ham_op.exp_val(ground_state), H2_STO3G_FCI_ENERGY)
assert np.isclose(ham_op.mat_elem(ground_state, ground_state), H2_STO3G_FCI_ENERGY)