Quantum neural networks force fields generation

We adopt a supervised learning approach with training sets of the form $\mathcal{A} = \{(\vec{C}_\alpha,E_\alpha,\vec{F}_\alpha)\}$, namely collections of n-atom molecular configurations $\vec{C}_\alpha = (x_{a_1}^{\,\alpha},y_{a_1}^{\,\alpha},z_{a_1}^{\,\alpha},\ldots,x_{a_n}^{\,\alpha},y_{a_n}^{\,\alpha},z_{a_n}^{\,\alpha})$ with associated total energies E_α and forces $\vec{F}_\alpha = (F_{a_1,x}^{\,\alpha},F_{a_1,y}^{\,\alpha},F_{a_1,z}^{\,\alpha},\ldots,F_{a_n,x}^{\,\alpha},F_{a_n,y}^{\,\alpha},F_{a_n,z}^{\,\alpha})$. Here, the index α runs over the different elements of the training set, $c_{a_i}^{\,\alpha}$ for $c = x,y,z$ is the Cartesian coordinate of the ith atom and:

$$\begin{align} F_{a_i,c}^{\,\alpha} = - \frac{\partial E_\alpha}{\partial c^{\alpha}_{a_i}}, \end{align} \tag{ 1 }$$

is the force acting on atom a_i along the direction c. We will also denote with $|\mathcal{A}|$ the number of samples contained in $\mathcal{A}$.

The training data, derived from classical ab initio methods that will be specified in the following, are used to optimize the predictions made by quantum models, specifically QNNs. These are based on the general notion of PQCs [33, 35, 37], consisting of a reference initial state $|0\rangle$, an output observable $\mathcal{O}$ and a model unitary $M(\vec{x},\Theta)$, depending on both some input data $\vec{x}$ and a set of trainable parameters $\Theta = \{\vec{\theta}_0,\ldots,\vec{\theta}_D\}$. The function expressed by a QNN takes the general form:

$$\begin{align} f_\Theta(\vec{x}) = \langle0|M(\vec{x},\Theta)^\dagger \mathcal{O} M(\vec{x},\Theta)|0\rangle, \end{align} \tag{ 2 }$$

and can be represented schematically as shown in figure 1. For our application of QNNs to PES learning, the Cartesian coordinates $\vec{C}_\alpha$ will be playing the role of the input $\vec{x}$. In practical realizations, we allow for an additional classical preprocessing map $\vec{y} = W(\vec{x})$, which in our specific context can serve the purpose of converting between Cartesian coordinates and more suitable molecular descriptors, as well as enhancing the effective nonlinearity of the model [34].

Zoom In Zoom Out Reset image size

It is worth noticing that while the term QNNs is used to emphasize some significant similarities between parametrized quantum models and classical neural networks (NN)—most prominently the layered structure and the use of trainable parameters—some important differences also exists. Indeed, contrary to the usual feed-forward scheme of classical NN, QNNs are quite easily designed and interpreted in the form of re-uploading circuits [44, 45], where trainable layers $\mathcal{U}_\ell(\vec{\theta}_\ell)$, red in figure 1, are alternated with data encoding ones $\Phi_\ell(\vec{y})$, blue in figure 1, with the same classical input values appearing multiple times. The choice of this re-uploading architecture is actually related to a more fundamental difference between classical and quantum NNs, namely the fact that quantum unitary layers act linearly on quantum states, hence lacking an explicit implementation of the usual non-linear activation function as a separate step. Instead, the non-linear manipulation of the classical input data is essentially referred to the choice of the encoding map, the trainable quantum gates and the final measurements.

It has actually been shown that the re-uploading mechanism makes QNNs universal functional approximators [44–46]: more specifically, any QNN output function (equation (2)) can be recast into a truncated Fourier series with a set of available independent frequencies determined by the eigenvalues of the encoding map and growing with the number of re-uploading steps. While recent results suggest that, by expanding the size of the quantum registers, these models can actually be mapped back on simpler sequential ones in which all input operations appear at the beginning [47], one can still profit from the insights offered by the re-uploading picture as a guide for intuition in the actual design of application-specific QNNs. As an example, theory suggests that, by adjusting the number of input layers, the richness of the Fourier spectrum can be systematically increased.

The structure of QNNs is completed by choosing a physical observable $\mathcal{O}$ and a suitable loss function $\mathcal{L}$, whose minimization drives the update of the trainable parameters via a classical optimization routine. In the following, we will make the simple choice:

$$\begin{align} \mathcal{O} = \sigma_z^1, \end{align} \tag{ 3 }$$

namely we will take the expectation value of the Pauli-z operator on the first qubit to construct the network output. Moreover, we will use a quadratic mean square error (MSE) loss function:

$$\begin{align} \mathcal{L}_\chi(\mathcal{A},\Theta) = \textrm{MSE(Energy)} + \chi\cdot \textrm{MSE(Forces)}, \end{align} \tag{ 4 }$$

with the hyperparameter χ weighting the contribution of energy and forces [48]. In practice, we directly associate the output of the QNN with the energy potential surface by defining:

$$\begin{align} \textrm{MSE(Energy)} = \frac{1}{|\mathcal{A}|}\sum_{\alpha \in \mathcal{A}} \Big(f_\Theta(\vec{C}_\alpha)-E_\alpha\Big)^2. \end{align} \tag{ 5 }$$

In parallel, consistent predictions for the forces are obtained by taking the derivative of the quantum circuit output with respect to the molecular coordinates:

$$\begin{align} \textrm{MSE(Forces)} = \frac{1}{3n\cdot|\mathcal{A}|}\sum_{\alpha \in \mathcal{A}} \Big\Vert \boldsymbol{\nabla}_c\ f_\Theta(\vec{C}_\alpha)+\vec{F}_\alpha\Big\Vert^2. \end{align} \tag{ 6 }$$

We leave for future works the investigation of alternative strategies, such as the use of two independent quantum circuits for the separate learning of energies and forces.

As mentioned in the previous section, the encoding operations $\Phi_\ell(\vec{y})$ are used to input the molecular configurations in the QNN model. Here, we assume for simplicity that these encoding unitaries are identical across different layers or re-uploading stages, i.e. $\Phi_\ell(\vec{y})\equiv \Phi(\vec{y})$. Moreover, let us denote by N the linear dimension of the feature vector $\vec{y}$, where in general N can differ from 3 n due to the preprocessing stage. Following standard practices [49], we will use N qubits to encode and manipulate N-dimensional features $\vec{y} = (y_1,\ldots,y_j,\ldots,y_N)$.

The quantum feature map $\Phi(\vec{y})$ is then constructed according to the expression:

$$\begin{align} \Phi(\vec{y}) = \mathcal{E}(\vec{y})\mathcal{S}(\vec{y}), \end{align} \tag{ 7 }$$

where $\mathcal{S}(\vec{y})$ is a collection of single-qubit Pauli-y rotations:

$$\begin{align} \mathcal{S}(\vec{y}) = \prod_{j = 1}^N \exp\left(-i\frac{\sigma_y^{(j)}}{2}y_j\right), \end{align} \tag{ 8 }$$

and $\mathcal{E}(\vec{y})$ is an entangling operation of the form:

$$\begin{align} \mathcal{E}(\vec{y}) = \prod_{(j,k)\in \mathcal{P}} \exp\left(-i(1/2)\sigma_z^{(j)}\sigma_z^{(k)}y_jy_k\right). \end{align} \tag{ 9 }$$

The latter directly resembles the so called ZZ feature map originally introduced in [49]. The set of qubit pairs $\mathcal{P}$ can be chosen in different ways, balancing ease of implementation on physical architectures with functional expressivity. Standard examples include linear entanglement:

$$\begin{align} \mathcal{P}_{\textrm{linear}} = \{(j,j+1)\, | \, j = 1,\ldots,N-1\}, \end{align} \tag{ 10 }$$

circular entanglement:

$$\begin{align} \mathcal{P}_{\textrm{circ}} = \{(j,j+1~\textrm{mod } N) \, | \, j = 1,\ldots,N\}, \end{align} \tag{ 11 }$$

and full entanglement:

$$\begin{align} \mathcal{P}_{\textrm{full}} = \{(j,k) \, | \, j = 1,\ldots,N-1, k \gt j\}. \end{align} \tag{ 12 }$$

For further generality, we also define a natural extension of the ZZ feature map to degree l interactions, adding factors of the form:

$$\begin{align} \exp\left(-i(1/2)\prod_{k = 1}^ly_{j_k}\sigma_z^{(j_k)}\right), \end{align} \tag{ 13 }$$

where, in principle, the l qubits $(j_1,\ldots,j_l)$ can be arbitrarily chosen among the N available.

In figure 2 we explicitly show a 3-qubit example with a full 2-qubit entangling map and an additional l = 3 operation. All multiple qubit operations are already decomposed into a standard universal set made of single-qubit rotations and CNOTs [50], which is typical of superconducting quantum computing architectures.

Zoom In Zoom Out Reset image size

Before moving to the description of the trainable part of the QNN models, let us also remark a few points about the classical preprocessing step W (see figure 1). As suggested above, this classical manipulation is generally used to enhance the nonlinear behaviour of the network, for example by taking inverse trigonometric functions of the original inputs [34]. At the same time, we can use this initial step to embed the relevant set of physical symmetries into the abstract representation of the target molecular systems seen by the QNN. This is known to be crucial already for classical ML methods, where the role of symmetry preserving features is played, e.g. by the so called symmetry functions [8, 10].

To limit the complexity of the quantum models, in the following we will use a set of internal coordinates, namely bond distances and angles, which by design respect translational and rotational symmetries. The integration of more advanced techniques, including fragmentation of large systems into local atomic environments, the use of other classes of molecular fingerprints and, possibly, the realization of symmetry adapted quantum circuits all represent natural future extensions of the present work.

We explicitly notice that the conversion between Cartesian and internal coordinates ($\vec{I}_\alpha = W(\vec{C}_\alpha)$) must be taken into account when computing the quantum forces predictions, as these are defined with respect to the former class (equation (1)). As a result, equations (5) and (6) are more properly rewritten as:

$$\begin{align} \textrm{MSE(Energy)} = \frac{1}{|\mathcal{A}|}\sum_{\alpha \in \mathcal{A}} \Big(f_\Theta\left(W(\vec{C}_\alpha)\right)-E_\alpha\Big)^2, \end{align} \tag{ 14 }$$

and

$$\begin{align} \textrm{MSE(Forces)} = \frac{1}{3n\cdot|\mathcal{A}|}\sum_{\alpha \in \mathcal{A}} \sum_{c,a_i}\Big( \boldsymbol{\nabla}_I\ f_\Theta(\vec{I}_\alpha)\cdot \frac{\partial W}{\partial c^{\alpha}_{a_i}}+F^{\,\alpha}_{a_i,c}\Big)^2.\end{align} \tag{ 15 }$$

The parametrized operations $\mathcal{U}_\ell(\vec{\theta}_\ell)$ represent the adjustable components of the QNN model. If suitably optimized during the learning phase, these select the most appropriate mapping from inputs (molecular coordinates) to outputs (energy and forces).

The precise structure of the $\mathcal{U}_\ell(\vec{\theta}_\ell)$ can vary significantly across different models and applications, with several popular choices known in the literature. For our specific functional regression problem, we follow the steps of Mitarai et al [34] and make use of a physics-inspired ansatz which is known to generate highly entangled states. In particular, we consider the following fully connected transverse field Ising Hamiltonian:

$$\begin{align} H = \sum_{j = 1}^N a_j \sigma_y^{(j)} + \sum_{j = 1}^N\sum_{k = j}^NJ_{jk}\sigma_z^{(j)}\sigma_z^{(k)}, \end{align} \tag{ 16 }$$

and we use the induced time evolution operator $U(t) = e^{-iHt}$ as a template for the trainable layers. By using the approximate product formula:

$$\begin{align} e^{-itH} \approx \prod_{j = 1}^{N}e^{-ita_j\sigma_y^{(j)}}\cdot \prod_{j = 1}^N \prod_{k = j}^N e^{-it\,\,J_{jk}\sigma_z^{(j)}\sigma_z^{(k)}}, \end{align} \tag{ 17 }$$

and by replacing a_j t and $J_{jk}t$ with free parameters, we obtain a unitary operation of the form:

$$\begin{align} \mathcal{U}_\ell(\vec{\theta}_\ell) = \prod_{j = 1}^{N}e^{-i\frac{\theta_\ell^j}{2}\sigma_y^{(j)}}\cdot \prod_{j = 1}^N \prod_{k = j}^N e^{-i\frac{\theta_\ell^{jk}}{2}\sigma_z^{(j)}\sigma_z^{(k)}}, \end{align} \tag{ 18 }$$

where now $\vec{\theta}_\ell$ is the collection of all $\theta_\ell^{\,j}$ and $\theta_\ell^{\,jk}$. It is easy to see that $\mathcal{U}_\ell(\vec{\theta}_\ell)$ can be implemented with single-qubit Pauli-y rotations and two-qubit ZZ operations, similarly to what happens for the encoding map Φ defined in the previous section. We stress however that, while the encoding layers are identical across different reuploading stages, as they repeatedly input the same classical data $\vec{y}$, the trainable parameters $\vec{\theta}_\ell$ are allowed to vary across different layers $\ell = 0,\ldots,D$.

For $\ell = 0$, we make the special choice $\theta_\ell^{jk} = 0~\textrm{, } \forall j,k$, namely we only use single-qubit rotations at the beginning of the computation (see figure 1). Moreover, we empirically find that including generalized l-qubit interactions up to l = 3, i.e. terms of the form $e^{-i(1/2)\theta_{jkm}\sigma_z^{(j)}\sigma_z^{(k)}\sigma_z^{(m)}}$ improves the overall performances of the model at the cost of only a modest increase of circuit complexity. The resulting $M(\vec{x},\Theta)$, as seen in figure 1, can thus be written as:

$$\begin{align} M(\vec{x},\Theta) = \left[ \prod_{\ell = D}^1 \mathcal{U}_\ell(\vec{\theta}_\ell)\cdot \Phi(W(\vec{x})) \right] \cdot \mathcal{U}_0(\vec{\theta}_0). \end{align} \tag{ 19 }$$

Notice that the product runs from $\ell = D$ to $\ell = 1$ from left to right.

We train QNN models by minimizing an average quadratic error that in principle contains both energy and forces labels, see equation (4). In most instances, we make use of an update rule which follows the negative gradient direction of the loss function:

$$\begin{align} \Theta^{t+1} = \Theta^t-\eta\nabla \mathcal{L}_\chi(\mathcal{A},\Theta^t), \end{align} \tag{ 20 }$$

where η is the learning rate and $\Theta^t$ the set of trainable parameters at the optimization step t. For energy predictions derived from a quantum circuit, derivatives with respect to any given trainable parameter µ can be easily computed with the parameter shift rule [51]:

$$\begin{align} \frac{\partial E_\alpha^\textrm{QNN}}{\partial \mu} = \frac{1}{2}\left[f_{\mu+\frac{\pi}{2};\Theta_\mu}(W(\vec{C}_\alpha))-f_{\mu-\frac{\pi}{2};\Theta_\mu}(W(\vec{C}_\alpha))\right], \end{align} \tag{ 21 }$$

where $E_\alpha^\textrm{QNN} = f_\Theta\left(W(\vec{C}_\alpha)\right)$ and $\Theta_\mu$ is the set of all trainable parameters without µ. The corresponding result for the forces predictions is obtained with the iterative parameter shift rule [52], and schematically reads:

$$\begin{align} \frac{\partial F_{a_i,c}^{\alpha,\textrm{QNN}}}{\partial \mu} & = -\frac{\partial^2 E^{\textrm{QNN}}_\alpha}{\partial c_{a_i} \partial \mu} = - \frac{1}{4}\sum_j\Big[f_{\mu+\frac{\pi}{2};\Theta_\mu}\left(\vec{I}_\alpha +\frac{\pi}{2}\vec{e}_j\right) \nonumber\\ &\quad -f_{\mu-\frac{\pi}{2};\Theta_\mu}\left(\vec{I}_\alpha +\frac{\pi}{2}\vec{e}_j\right) -f_{\mu+\frac{\pi}{2};\Theta_\mu}\left(\vec{I}_\alpha -\frac{\pi}{2}\vec{e}_j\right) \nonumber\\ &\quad +f_{\mu-\frac{\pi}{2};\Theta_\mu}\left(\vec{I}_\alpha -\frac{\pi}{2}\vec{e}_j\right) \Big] \cdot \frac{\partial I_j^{\,\alpha}}{\partial c_{a_i}}, \end{align} \tag{ 22 }$$

with $(\vec{e}_j)_k = \delta_{jk}$ the standard basis vectors.

Notice that the parameter shift rule of equation (21) can also be used to retrieve the actual QNN forces predictions, namely:

$$\begin{align} F^{\alpha,\textrm{QNN}}_{a_i,c} = -\boldsymbol{\nabla}_I\,\, f_\Theta(\vec{I}_\alpha)\cdot \frac{\partial W}{\partial c_{a_i}} = -\sum_j \frac{\partial f_\Theta(\vec{I}_\alpha)}{\partial I^{\,\alpha}_j}\frac{\partial I^{\,\alpha}_j}{\partial c_{a_i}}. \end{align} \tag{ 23 }$$

In this case, factors of the form $\partial f/\partial I_j$ must be computed from the quantum circuit, while the chain rule factors taking into account the classical preprocessing are known analytically and are determined solely by the mapping function W.

At the beginning of the training, all free parameters are initialized at zero, so that each hidden layer acts as the identity operator. This essentially follows the recommendations given in [53] to promote effective optimization steps in the early training phase.

In the remaining part of this work, we will provide a series of examples to demonstrate how QNN models are able to learn and give consistent predictions of PES and FFs for individual molecules. The overall performances will be assessed primarily through the evaluation of the average root MSE on suitable test sets.

In addition to that, we will also make use of the concept of effective model dimension—originally introduced in [39]—to inform the comparison with classical counterparts and to investigate potential advantages. As a brief summary, the effective dimension quantifies the capacity of both classical and quantum parametrized ML models [54] by measuring the portion of model space that they occupy. In other words, it estimates the capability of a model in covering the functional space defined by a particular model class by making a productive use of all its parameters, going beyond naive parameter-counting arguments. A high effective dimension is therefore related to a richer set of expressible functions and better trainability, as it can lead to a more favourable landscape for gradient based methods [39]. Moreover, the effective dimension can sensibly bound the model generalization error [39, 54].

For a given statistical model $y = f_\Theta(x)$, the effective dimension is defined as:

$$\begin{align} d_n(f_\Theta) = \frac{2\log{\left(\frac{1}{V_\Omega}\int_{\Omega}\sqrt{\det\left(\unicode{x1D7D9}+\frac{n}{2\pi\log{n}}F(\Theta)\right)}\,d\Theta\right)}}{\log{\left(\frac{n}{2\pi \log{n}}\right)}}, \end{align} \tag{ 24 }$$

where $\Omega\subset \mathbb{R}^d$ is a d-dimensional parameter space of volume $V_\Omega$, $n\in\mathbb{N}$ is the number of data samples. Here, we have also introduced the Fisher information matrix:

$$\begin{align} F & = \mathbb{E}\left[\nabla_\Theta \log{(p(x,y;\Theta)}\nabla_\Theta \log{(p(x,y;\Theta)}^T\right] \end{align} \tag{ 25 }$$

$$\begin{align} &\approx \frac{1}{K}\sum_{k = 1}^K\nabla_\Theta f_\Theta(x_k)\nabla_\Theta f_\Theta(x_k)^T, \end{align} \tag{ 26 }$$

with $p(x,y;\Theta) = p(y|x;\Theta)p(x)$ and $p(y|x;\Theta) = \exp(-\Vert{y-f_\Theta(x)}\Vert^2/2)/\sqrt{2\pi}$ being the probability distribution mass of the model. The effective dimension is bounded by the rank of the Fisher information matrix and is usually normalized with the total number of parameters d.

3.1.1. Lithium hydride

As a first proof-of-concept, we consider a single LiH molecule and we design a QNN model that learns its dissociation curve, and the corresponding FF, as a function of the bond length r.

The presence of a single internal coordinate, which is obtained from the Cartesian positions of the Li and H atoms with a mapping $W(\vec{C}) = |\vec{C}_\mathrm{Li}-\vec{C}_\mathrm{H}|$, makes the problem effectively one dimensional. We construct a classical dataset by numerically diagonalizing the second quantized Hamiltonian expressed in the STO3G basis set [41] for different bond lengths r in the range $[0.9, 4.5]$ Å. Forces are computed via finite differences over the exact PES.

To make better use of the Fourier series structure of the QNN, we make the exact PES symmetric by mirroring it around the r = 4.5 Å and selecting the data set over the extended range, see appendix A. Furthermore, we use an additional preprocessing step—on top of the Cartesian to internal coordinate transformation—to construct 3-dimensional features from the original 1D problem. First, we scale all inputs in $[-1,1]$ with the scikit-learn MinMaxScaler, then we apply the map:

$$\begin{align} W : [-1,1] &\rightarrow [-\pi,\pi]^3 \nonumber\\[6pt] r & \mapsto \begin{pmatrix} \pi r \\ \arcsin{(r)} \\ \arccos{(r)}\end{pmatrix}.\end{align} \tag{ 27 }$$

We employ a 3-qubit QNN model with D = 10 trainable layers, alternated with the same amount of input stages. Notice that, even if a single qubit could in principle be used to encode the 1-dimensional problem under study, we instead choose to engineer 3-dimensional feature vectors from the single independent internal coordinate available (interatomic distance) to enrich by design the set of expressible functions and to enable the use of genuine multi-qubit quantum resources such as entanglement. More in general, a very common choice when using a parametric type of encoding is to use one qubit per feature. However, similarly to the depth of the model and the overall number of trainable parameters, this is by no means a strict requirement and could be freely modified, subject to appropriate convergence tests, to improve performances or in response to additional constraints (e.g. limited number of qubits or finite quantum coherence). We use a full feature map, as introduced in section 2.2, with l = 3 (see also figure 2), and similar trainable layers with full entanglement and degree l = 3 interactions. The input of the feature map is $\vec{y} = W(r)$, with W as in equation (27) and r the Li-H inter-atomic distance.

We benchmark our model against a fully connected classical NN, which is given access to the same set of seven engineered inputs and products thereof shown in figure 2. The classical NN is composed of five layers with $[7,4,5,2,1]$ units using the hyperbolic tangent activation function. Such network topology is chosen through a random search over all possible configurations with the same number of parameters (identical to the QNN), selecting the one that achieves the best validation loss.

Both the classical and the QNN models are trained with 4000 steps of the ADAM [55] algorithm on 50 data points using the $\mathcal{L}_0$ (χ = 0) loss. We remark that the small size of the classical NN makes it quite sensitive to parameters initialization, even when using the Glorot and Bengio [56] scheme. Hence, the model underfits in $65\%$ of the cases and needs 1600 epochs to reach convergence. On the other hand, the QNN appears to be robust against weights initialization and converges in 230 epochs.

The validation root mean square errors (RMSE) for both models, computed on a test set with 120 data points, are given in table 1 and the respective predictions are shown in figure 3. We notice that, in this setting, the QNN outperforms the best classical counterpart with the same number of parameters in terms of stability and quality of the predictions, and also exhibits a larger effective dimension. Although the classical model can match the QNN results in absolute terms when its size grows or if the forces are explicitly included in the loss function (χ ≠ 0), the present comparison, supported by the effective dimension analysis, certifies the competitiveness of the proposed quantum models. In appendix B we also report an example of MD trajectories driven by the exact and learned FFs.

Zoom In Zoom Out Reset image size

Table 1. Validation RMSE for LiH energy (eV) and forces (eV Å⁻¹), together with the normalized effective dimension ($d_{50}/d$), the total number of parameters (d) and the number of training epochs for both the QNN and the classical NN models.

LiH	RMSE (E)	RMSE (F)	$d_{50}/d$	d	Epochs
QNN	0.003	0.044	0.9	73	230
NN	0.034	0.21	0.38	73	1600

3.1.2. Argon dimer

As a second experiment, we consider two Argon atoms separated by a bond length r. This example covers a richer set of intermolecular interactions, as Ar₂ is loosely bounded via Van der Waals forces. The dataset is constructed using density functional theory (DFT/B3LYP with Grimme D3 corrections as implemented in Gaussian [57]), and smoothened by fitting to a $C_1/r^{12} - C_2/r^{6}$ function.

The training strategy is almost identical to the LiH case, including the QNN architecture (with D = 9) and the dataset mirroring. One notable difference is the use of the COBYLA optimizer, which empirically showed better performance than the ADAM one. The benchmark is performed against a fully connected classical NN, taking as input the bond length, composed of ten layers with $[1, 7, 4, 3, 2, 5, 3, 4, 2, 1]$ hidden units. The architecture is chosen with a random search through the models with the same number of parameters as the QNN and variable number of inputs, selecting the one that achieves the best validation loss.

The validation RMSE for both models, computed on a test set with 120 data points, are given in table 2 and the respective predictions are shown in figure 4. As for the LiH case, we observe that the QNN outperforms the best classical counterpart with the same number of parameters in terms of stability and quality of the predictions, and also exhibits a larger effective dimension. Notably, the region around the PES minimum is better reconstructed with the QNN than the classical NN.

Zoom In Zoom Out Reset image size

Table 2. Validation RMSE for Ar₂ energy (eV) and forces (eV Å⁻¹), together with the normalized effective dimension ($d_{165}/d$), the total number of parameters (d) for both the QNN and the classical NN models.

Ar2	RMSE (E)	RMSE (F)	$d_{165}/d$	d
QNN	$2.6\times10^{-5^{}}$	0.058	0.77	66
NN	$1\times 10^{-4}$	0.6	0.11	66

We move towards the multi-dimensional case by considering an individual H₂O molecule. For this more challenging test, we choose 3 internal coordinates, namely the two O–H bond lengths and the H–H planar angle, which we scale again in $[-1,1]$ and preprocess to create the following 3-dimensional feature vector:

$$\begin{align} \vec{y} = \begin{pmatrix} \arcsin{(r_{\mathrm{OH},1})} \\ \arcsin{(r_{\mathrm{OH},2})} \\ \arcsin{(\phi_{\mathrm{HH}})}\end{pmatrix}, \end{align} \tag{ 28 }$$

which is then inputed in the model via the feature map. The classical configurations and the corresponding energy and forces labels, computed with DFT, are retrieved from [58]. To simplify the problem, we concentrate on configurations around the equilibrium position, discarding those with bond length outside $[1.6,2.1]$ Å.

Our QNN model is constructed similarly to the LiH case, with three qubits and a depth of D = 12 encoding and trainable layers. We compare again its results with those of a classical NN model with the same number of parameters and designed with the same conditions applied in the LiH case presented above. Furthermore, we also use as a reference the state-of-the-art specialized n2p2 [59] classical package, which makes full use of symmetry functions and of the sub-network fragmentation idea [8], and from which we expect peak performance.

More specifically, the simplified classical NN model is a fully-connected 6-layer network with $[7, 4, 6, 2, 2, 1]$ units and hyperbolic tangent activation function. Instead, the n2p2 network uses three sub-networks with two hidden layers, each with 15 units, and hyperbolic tangent activation. This model takes as input a set of 15 symmetry functions for the Oxygen atom and 20 for the two Hydrogen ones. The models are trained by minimizing a $\mathcal{L}_1$ loss function on 300 data points (1000 for n2p2). Notice that here we set χ = 1: indeed, as pointed out in [60], it is important to incorporate the forces in the training if we wish to predict them accurately. However, this makes a full gradient descent computationally very intensive for the simulation of the quantum model, as it requires the calculation of the circuit Hessian matrix with recursive parameter shifts (see equation (22)). For this reason, we choose the gradient free optimizer COBYLA for the present demonstration, while the classical models are trained with the ADAM method.

The validation loss results (computed on a test set with 650 data points) are presented in table 3, while figure 5 shows the predicted energy and forces against the reference data points. The QNN is once again competitive with respect to classical counterparts, outperforming the non-specialized classical NN and reporting the largest normalized effective dimension.

Zoom In Zoom Out Reset image size

Table 3. Validation RMSE for H₂O energy (eV) and forces (eV Å⁻¹), together with the normalized effective dimension ($d_{300}/d$) and the total number of parameters (d) for the QNN and the classical NN and n2p2 models.

H2O	RMSE (E)	RMSE (F)	$d_{300}/d$	d
QNN	0.005	0.06	0.72	87
NN	0.006	0.1	0.25	87
n2p2	$7 \times 10^{-4}$	0.01	0.04	1642

In our last test, we consider a single hydronium (H₃O⁺) molecule. Following [61], we prepare a training set by sampling configurations traversed along the inversion pathway shown in figure 6 using DFT-based MD simulations at 400 K and collecting the corresponding energies and forces. In order to sample the full profile along the dihedral angle ‘HHHO’ from −0.78 to 0.78 rad, we also applied a dynamical constraint with increments of 0.005 rad at each MD time step. All calculations were performed with the plane-wave (PW) code CPMD [62] using unrestricted Kohn–Sham DFT with the PBE functional [63], a PW cutoff of 70 Ry, a cubic simulation box of edge 14 Å, and Trouiller–Martins pseudo-potentials [64]. For the MD, a time step 10 a.u. (0.242 fs) was used.

Zoom In Zoom Out Reset image size

We describe each configuration with 6 internal degrees of freedom, namely three O–H bond lengths, two O–H angles and the dihedral angle formed by the four atoms (see appendix C for the formal definitions). The internal coordinates are preprocessed with the function $W(\cdot) = \arcsin{(\cdot)}$ and loaded in the model via the feature map, see appendix D. We also make use of the explicit formulas for derivatives with respect to Cartesian coordinates provided in [65].

A 6-qubit QNN model is constructed with D = 10 repetitions of the encoding map and the trainable layers, both with linear entanglement and order l = 3 interactions. The latter are also placed with linear couplings across all neighbouring 3-qubit groups. Due to the complexity of the system and the high computational cost, we only include energy labels in the training set, leaving a more refined study of the forces for future investigations. Hence, we train the model by minimizing a $\mathcal{L}_0$ loss on 500 data points (9000 for n2p2) with 5000 steps of the ADAM algorithm.

As in the previous examples, we compare the QNN model with a classical NN, whose topology is optimized under the constraint that the number of trainable parameters be the same of the QNN, and with a n2p2-build model. In this case, the classical NN has [6, 14, 2, 1] units, while the n2p2 NN is composed of four sub-networks, each one with a structure similar to the case if the H₂O molecule in section 3.2.

The RMSE results for a test set with 500 data points are reported in table 4, and energy predictions are also depicted in figure 6. Here, the QNN model is still competitive with respect to the generic classical NN, although both of then are clearly outperformed by the specialized and much larger n2p2 model. It is also worth noticing that forces predictions are much worse than energy ones, which confirms the necessity of including them in the loss function to achieve good results in non-trivial systems.

Table 4. Validation RMSE for H₃O⁺ energy (eV) and forces (eV Å⁻¹), together with the normalized effective dimension ($d_{500}/d$) and the total number of parameters (d) for the QNN and the classical NN and n2p2 models.

H3O+	RMSE (E)	RMSE (F)	$d_{500}/d$	d
QNN	$3.7\times 10^{-3}$	0.26	0.67	135
NN	$4.2\times 10^{-3}$	0.207	0.19	135
n2p2	$5\times 10^{-4}$	0.19	0.03	2214