Use this file to discover all available pages before exploring further.
View on GitHub
Open this notebook in GitHub to run it yourself
Intermolecular dispersion forces can be determined by calculating the interaction between two molecules; however, the total force does not simply triple when a third molecule is introduced.This phenomenon is referred to as Many-Body Dispersion (MBD).While it is desirable to evaluate dispersion forces using first-principles calculation methods, attempting to do so with multiple specific molecules results in an exorbitant increase in computational complexity.The Quantum Drude Oscillator (QDO) model [1] is a coarse-grained model that describes “molecules interacting via dispersion forces” as “electrostatic interactions between QDOs.” In this model, a QDO is defined as a harmonic oscillator consisting of a pair of fictitious positively and negatively charged particles connected by a spring. QDOs behave as bosons, and in this tutorial, we see how we can treat such bosonic system using classiq [2].
!pip install -qq "classiq[chemistry]"
import timefrom functools import reducefrom operator import mulimport matplotlib.pyplot as pltimport numpy as npfrom classiq import *
In this tutorial, we consider a non-dimensional Hamiltonian for N one-dimensional QDOs which are interacting via quadratic coupling.H^=i=1∑N(x^i2+p^i2)+j>i∑γi,jx^ix^j,where \xi and \pi are the bosonic position and momentum operators for oscillator i and γi,j∈R is the coupling between oscillators i and j.For simplicity, we assume that all one-dimensional QDOs are identical and aligned parallel to each other and along the inter-oscillator axis, equally separated by distance R. In this case, the coupling constant is given by γi,j=−4α(R∣i−j∣)−3, where α is the dipole polarisability.We need to represent this Hamiltonian in terms of Pauli matrices to be able to implement on quantum circuit.Unlike fermions, bosonic Fock states cannot be represented by only zeros and ones. Therefore, when bosons occupy a certain quantum state, we need to encode this by representing the number n as a binary string in the qubit state. In this tutorial, we denote a Fock state expressed in decimal as ∣n⟩ (with an underline), and a Fock state in binary representation as ∣\bin(n)⟩ (without an underline). Here, \bin(n) represents the binary representation of n, and \bin(n)i denotes the i-th bit from the right of the bitstring.When adopting this binary encoding, there is an upper limit to the number of bosons that can occupy the same quantum state. If we allocate m qubits to the Fock state of each boson, the possible values range from 0 to 2m−1. We define the number of levels we can handle as d=2m. d must be large enough to make the result as reliable as possible.Although not used in this tutorial, Classiq’s QNum could potentially handle such bosonic Fock states more intuitively.The bosonic position and momentum operators can be written using ladder operators as:\x\p=21(\aD+\a),=−2i(\aD−\a),where the ladder operators are expressed using projection operators as follows:a^†a^=n=0∑d−1n+1∣n+1⟩⟨n∣,=n=0∑d−1n+1∣n⟩⟨n+1∣,What does this projection operator look like in binary representation? Let’s take ∣4⟩⟨3∣ with a bit length of m=4 as an example. In this case, the projection operator is:∣4⟩⟨3∣=∣0100⟩⟨0011∣=∣0⟩⟨0∣⊗∣1⟩⟨0∣⊗∣0⟩⟨1∣⊗2As we can see, the operator ∣n+1⟩⟨n∣ can be divided into three parts:
For the bits where ones continue from the 20 position, apply ∣0⟩⟨1∣.
For the rightmost zero, apply ∣1⟩⟨0∣.
For all other bits, there is no change, so apply ∣0⟩⟨0∣ or ∣1⟩⟨1∣.
Therefore, ∣n+1⟩⟨n∣ can be expressed using Pauli matrices as:∣n+1⟩⟨n∣=(i=k+1⨂m−1∣bin(n)i⟩⟨bin(n)i∣)⊗∣1⟩⟨0∣⊗(∣0⟩⟨1∣)⊗k=(i=k+1⨂m−12I+(−1)bin(n)iZ)⊗2X−iY⊗(2X+iY)⊗k,where k is the position of the rightmost 0 in \bin(n).By noting the relationship ∣n+1⟩⟨n∣=(∣n⟩⟨n+1∣)†, we can write the QDO Hamiltonian in terms of Pauli matrices.The following code generates the coupling terms for this Hamiltonian.
def get_ith_bit(n, i): """Get the i-th bit of an integer. Args: n (int): Input integer. i (int): Bit position to retrieve. Returns: int: The i-th bit (0 or 1). """ return (n >> i) & 1
def get_rightmost_zero_pos(n): """Get the position of the rightmost '0' in an integer. Args: n (int): Input integer. Returns: int: Position of the rightmost '0'. """ return (~n & (-~n)).bit_length() - 1
def create_ladder_projector(n: int, qubits: list[int], direction: str) -> SparsePauliOp: """Create ladder projector operator for given occupation number. 'up': |n+1><n| 'down': |n><n+1| Args: n (int): Occupation number. qubits (list[int]): List of qubit indices. direction (str): 'up' or 'down' for ladder direction. Returns: SparsePauliOp: Ladder projector operator. """ m = len(qubits) if not 0 <= n < 2**m: raise ValueError("Occupation number n must be in the range [0, 2^m).") paulis_coeffs = [] k = get_rightmost_zero_pos(n) if direction == "up": for i in range(k): paulis_coeffs.append((Pauli.X(qubits[i]) + 1j * Pauli.Y(qubits[i])) / 2) paulis_coeffs.append((Pauli.X(qubits[k]) - 1j * Pauli.Y(qubits[k])) / 2) elif direction == "down": # adjoint of "up" for i in range(k): paulis_coeffs.append((Pauli.X(qubits[i]) - 1j * Pauli.Y(qubits[i])) / 2) paulis_coeffs.append((Pauli.X(qubits[k]) + 1j * Pauli.Y(qubits[k])) / 2) for i in range(k + 1, m): paulis_coeffs.append( (Pauli.I(qubits[i]) + ((-1) ** get_ith_bit(n, i)) * Pauli.Z(qubits[i])) / 2 ) return reduce(mul, paulis_coeffs)
# Test the `create_ladder_projector` function# example: |3><2| = |11><10| for 2 qubitsexpected = (Pauli.I(1) - Pauli.Z(1)) / 2 * (Pauli.X(0) - 1j * Pauli.Y(0)) / 2actual = create_ladder_projector(2, [0, 1], "up")print("Is the created ladder projector correct?", expected == actual)print( "Are |3><2| and |2><3| Hermitian conjugates?", ( hamiltonian_to_matrix(create_ladder_projector(2, [0, 1], "up")).conj().T == hamiltonian_to_matrix(create_ladder_projector(2, [0, 1], "down")) ).all(),)
Output:
Is the created ladder projector correct? True Are |3><2| and |2><3| Hermitian conjugates? True
def create_x_operator(qubits): r"""Create $\hat{x}$ operator for given number of qubits. Args: qubits (list[int]): List of qubit indices. Returns: SparsePauliOp: $\hat{x}$ operator. """ ladder_ops = [] for n in range(2 ** len(qubits) - 1): ladder_ops.append( create_ladder_projector(n=n, qubits=qubits, direction="up") * np.sqrt((n + 1) / 2) ) ladder_ops.append( create_ladder_projector(n=n, qubits=qubits, direction="down") * np.sqrt((n + 1) / 2) ) return reduce(lambda x, y: x + y, ladder_ops)
Now, let’s verify that the coupling terms between oscillator 1 and oscillator 2 are expressed as follows when the number of qubits m assigned to each QDO’s Fock state is m=1 and m=2.x^1x^2∣m=1=21X1X2,x^1x^2∣m=2=43+2X1X3+\termplusX1X2X3+\termplusX1X3X4+41X1X2X3X4+\termplusY1Y2X3+41Y1Y2X3X4+\termplusX1Y3Y4+41X1X2Y3Y4−\termminusY1Y2X3Z4−\termminusX1Z2Y3Y4−41X1Z2X3−\termminusX1Z2X3X4−41X1X3Z4−\termminusX1X2X3Z4+42−3X1Z2X3Z4+41Y1Y2Y3Y4.
Are the numerically constructed and analytically derived x1x2 operators (m = 2) close? True
The non-interacting terms of the QDO Hamiltonian can be reduced to a simplified form as follows:x^2+p^2=2(a^†+a^)2−2(a^†−a^)2=a^†a^+a^a^†=2a^†a^+1=2(n=0∑d−1(n+1)∣n+1⟩⟨n+1∣)+1=2(j=0∑m−12j∣1⟩⟨1∣j)+1=j=0∑m−12j(Ij−Zj)+1=2mI⊗m−j=0∑m−12jZjNote that the number operator a^†a^, defined by its action a^†a^∣n⟩=n∣n⟩, is diagonal in the Fock basis.
def get_noninteractiog_qdo_hamiltonian( num_qdos: int, num_qubits_per_qdo: int,) -> SparsePauliOp: """Construct the non-interacting Quantum Drude Oscillator (QDO) Hamiltonian. Args: num_qdos (int): Number of QDOs. num_qubits_per_qdo (int): Number of qubits per QDO. Returns: SparsePauliOp: The non-interacting QDO Hamiltonian as a SparsePauliOp. """ hamiltonian = [] for qdo_index in range(num_qdos): qubits = np.arange( qdo_index * num_qubits_per_qdo, (qdo_index + 1) * num_qubits_per_qdo ) identity = reduce(mul, [Pauli.I(i) for i in qubits]) hamiltonian.append(2**num_qubits_per_qdo * identity) for i in range(num_qubits_per_qdo): hamiltonian.append(-(2**i) * Pauli.Z(qubits[i])) return reduce(lambda x, y: x + y, hamiltonian)
With the code above, we are now ready to generate the QDO Hamiltonian.
def get_hamiltonian( num_qdos: int, num_qubits_per_qdo: int, coupling_constants: list[list[float]]) -> SparsePauliOp: assert np.array(coupling_constants).shape == (num_qdos, num_qdos) hamiltonian = [ get_noninteractiog_qdo_hamiltonian( num_qdos=num_qdos, num_qubits_per_qdo=num_qubits_per_qdo, ) ] for j in range(num_qdos): for i in range(j): xi = create_x_operator( qubits=list(range(i * num_qubits_per_qdo, (i + 1) * num_qubits_per_qdo)) ) xj = create_x_operator( qubits=list(range(j * num_qubits_per_qdo, (j + 1) * num_qubits_per_qdo)) ) hamiltonian.append(coupling_constants[i][j] * xi * xj) return reduce(lambda x, y: x + y, hamiltonian)def calculate_coupling_constants( polarizability: float, distances: list[list[float]],) -> list[list[float]]: """Calculate coupling constants for given parameters. Args: polarizability (float): Polarizability. distances (list[list[float]]): Distances between QDOs. Returns: list[list[float]]: Coupling constants matrix. """ # polarizability = effective_charge**2 / (effective_mass * frequency**2) num_qdos = len(distances) coupling_constants = [[0.0 for _ in range(num_qdos)] for _ in range(num_qdos)] for i in range(num_qdos): for j in range(i): coupling_constants[i][j] = -4 * polarizability / distances[i][j] ** 3 coupling_constants[j][i] = coupling_constants[i][j] return coupling_constants
Quantum program link: https://platform.classiq.io/circuit/3AelNrcQSJdeHxmkmeD8D7kS0YD
# Test VQEwith ExecutionSession(qprog) as es: test_result = es.minimize( cost_function=vqe_hamiltonian_test, initial_params={"params": [1.0] * num_params}, max_iteration=500, )
print(f"VQE Energy: {test_result[-1][0]}")vqe_test_results = {k: np.real(test_result[k][0]) for k in range(len(test_result))}plt.plot(vqe_test_results.keys(), vqe_test_results.values(), "-")plt.ylabel(r"Ground state energy $E$ [$\frac{1}{2} \hbar \omega$]")plt.xlabel("Iteration")plt.tick_params(axis="both")plt.grid()plt.show()
Output:
VQE Energy: 1.5542974087915205
Let’s attempt to reproduce Fig. 1 from [2]. We will model I2 molecules as one-dimensional QDOs and calculate how the London dispersion energy ΔE changes according to the distance between the two parallel oscillators. We define the parameters such that the axial polarizability is α=14.5A˚3 and the frequency is ℏω/2=9.61eV.
The theoretical value of the London dispersion force in this case can be calculated as follows.The Hamiltonian is:=˝2μp^12+2μp^22+21μω2(x^12+x^22+γx^1x^2),where μ is the effective mass and γ=−4αR−3. By defining the transformed variablesp^±=21(p^1±p^2),x^±=21(x^1±x^2),the Hamiltonian can be rewritten in this new uncoupled basis as=˝(2μp^+2+21μω2(1+2γ)x^+2)+(2μp^−2+21μω2(1−2γ)x^−2).Therefore, the ground state energy is given byE=21ℏω(1+2γ+1−2γ).In the absence of the London dispersion force—specifically, when the distance between the oscillators is R=∞—we have γ=0, which yields E=ℏω.From this, the London dispersion energy can be calculated as:ΔE=21ℏω(1+2γ+1−2γ−2).Below, we will confirm that VQE can obtain energy values close to this exact solution regardless of the distance between oscillators.
From this graph, we can see that the VQE results yield energy values close to the analytical solution we calculated earlier.Next, let’s investigate how ΔE changes with respect to the number of qubits m (and the corresponding dimension d) assigned to each QDO.
def get_energies( interoscillator_distances: list[float], num_qubits_per_qdo: int, num_qdos: int = 2) -> list[float]: energies = [] for interoscillator_distance in interoscillator_distances: distances = np.zeros((num_qdos, num_qdos)) for i in range(num_qdos): for j in range(num_qdos): distances[i][j] = abs(i - j) * interoscillator_distance vqe_hamiltonian = get_hamiltonian( num_qdos=num_qdos, num_qubits_per_qdo=num_qubits_per_qdo, coupling_constants=calculate_coupling_constants( polarizability=polarizability, distances=distances, ), ) ham_matrix = hamiltonian_to_matrix(vqe_hamiltonian) energies.append(np.linalg.eigvalsh(ham_matrix).min()) return np.array(energies)energies_d_2 = get_energies(interoscillator_distances, num_qubits_per_qdo=1)energies_d_4 = get_energies(interoscillator_distances, num_qubits_per_qdo=2)energies_d_8 = get_energies(interoscillator_distances, num_qubits_per_qdo=3)energies_d_16 = get_energies(interoscillator_distances, num_qubits_per_qdo=4)
Looking at this graph, we can see that the accuracy of the calculated energy increases as more qubits are used for each QDO.The impact of m (or d) on accuracy becomes more apparent as we vary the coupling constant γ.
This graph makes it clearer that higher values of d lead to greater accuracy. We can see that d=4 is a good approximation to the exact solution until very large couplings γ.For reference, −γ=2 corresponds to a Hamiltonian that is no longer positive-semidefinite, at this point the charged Drude particles dissociate from the nuclei and the system breaks down.For physically relevant systems, we expect to require value of −γ much lower than 2, and we can therefore limit ourselves to small d [2].
Now that we know the exact London dispersion energy of two identical parallel QDOs along the inter-oscillator axis isΔE=21ℏω(1+2γ+1−2γ−2),where γ=−4αR−3, we might expect that adding a third QDO would result in a simple summation of these pairwise interactions.Under this assumption, the total dispersion energy for a linear chain of three QDOs (with a distance R between neighbors) would be given by:ΔE=ℏω(1+2γ1+1−2γ1−2)+21ℏω(1+2γ2+1−2γ2−2),where γ1=−4αR−3,γ2=−4α(2R)−3. Let’s compare this to the results obtained by diagonalizing the Hamiltonian for d=4.
From this plot, we can see that diagonalizing the Hamiltonian for d=4 yields a larger energy shift at short distances between QDOs than when considering only two-body pairwise interactions.This suggests that many-body effects are playing a role.
