A Variable Neighbourhood Descent Heuristic for Conformational Search Using a Quantum Annealer

Discovering the low-energy conformations of a molecule is of great interest to computational chemists, with applications in in silico materials design and drug discovery. In this paper, we propose a variable neighbourhood search heuristic for the conformational search problem. Using the structure of a molecule, neighbourhoods are chosen to allow for the efficient use of a binary quadratic optimizer for conformational search. The method is flexible with respect to the choice of molecular force field and the number of discretization levels in the search space, and can be further generalized to take advantage of higher-order binary polynomial optimizers. It is well-suited for the use of devices such as quantum annealers. After carefully defining neighbourhoods, the method easily adapts to the size and topology of these devices, allowing for seamless scaling alongside their future improvements.

To address the computational complexity presented by larger molecules, many metaheuristic approaches have been studied. Examples of such approaches include genetic algorithms 16,17 , conformational space annealing 18,19 , tabu search 20,21 , molecular dynamics (MD) 22,23 , and basin/funnel hopping 24 . Variations of the Monte Carlo (MC) method have also been widely used 25,26 as a less computationally expensive alternative to MD 8 . In addition, parallel tempering (PT), also known as replica exchange 27,28 , can be applied to both MC and MD to further improve their sampling performance of the conformational search space.
A comparatively recent approach to address the growing computational complexity of optimization problems relies on the putative future advantage of specialized hardware like quantum annealers to solve binary quadratic optimization problems (see Supplementary Information for more details). One challenge lies in reformulating the optimization problems, a task that often requires approximations or simplifications. Our motivation was to develop such a formulation for the conformational search problem that avoids drastic compromises while producing good conformers.
For this purpose, we propose an iterative heuristic method for the conformational search problem based on variable neighbourhood descent (VND). In each iteration of the method, we use the molecular structure to choose specialized conformational neighbourhoods that can be minimized efficiently. More specifically, using the structural graph of a given molecule, subsets of rotatable bonds are selected at each iteration. Fixing the values of other torsion angles, the problem of minimizing the molecular energy with respect to the selected torsion angles is then formulated as a binary program with an objective function that is a polynomial of a chosen degree. This allows the method to be adapted to the specifics of the optimizer by limiting the degree of the binary program. The values of the selected torsion angles are then set to the solution of this binary program before starting each subsequent iteration during which a new subset of rotatable bonds are optimized. The process continues until some stopping criteria are met.
Although the method can be readily extended to any chosen degree, we assume in this paper that a binary quadratic program is desirable as it is well-suited for optimization using quantum annealing 29,30 . Furthermore, by changing the parameters of the neighbourhood selection procedure, the method can be easily adapted to the size of the conformational search problem in terms of the number of rotatable bonds, as well as the size and connectivity of available quantum annealers. The flexibility of solving any conformational search problem using currently available quantum annealers, without imposing restrictions on the granularity of the conformational space, differentiates our work from a previous study on protein folding using quantum annealing 31 .
We evaluate the performance of our proposed algorithm over three families of molecules relevant to industry, using an algorithm that returns an exact optimal solution and the D-Wave 2000Q quantum annealer [32][33][34] . The latter provides an assessment using the latest available hardware at the time of writing of this work, whereas the former can be seen as a limiting ideal case. For each molecule, we compare the lowest-energy conformations found by our algorithm with those found by both parallel tempering MC (PTMC) 27,35 and a simple local search method.

preliminaries
We give some preliminaries before presenting the details of our proposed conformational search method. Problem definition. As discussed above, we consider the conformational search problem as a special case of molecular structure analysis, where the structure is kept fixed except for rotations around selected bonds. Each of these torsional degrees of freedoms we hereafter refer to as a torsion for simplicity. We denote the i-th rotatable bond by T i and assign its rotational angle a variable t i , with i representing the torsion index. It is convenient to identify a conformation of a molecule with M torsions by a torsion vector , for all i, knowing that the method remains unchanged if each torsion has its own range chosen based on prior knowledge, experimental data, or known symmetries. For simplicity, let us assume all torsion angle values are chosen from the same set of d values The theoretical precision of this discretization scheme increases with d, while the size of the search space d M grows exponentially with the number of torsions.
Although it is natural to describe a molecule using a molecular graph, where the atoms and their bonds are represented by vertices and edges, respectively, we find it helpful to use the torsions to partition the molecule into + M 1 subsets called rigid bodies. The partitioning is performed such that all atoms within a rigid body are interconnected through non-torsion bonds. As a result, the relative positions of the atoms within a rigid body, denoted by R a , remain invariant under rotation and are therefore independent of t. This simplified representation of the molecule is now easily described by a rigid-body graph G = R T ( , ), where  is the set of + M 1 vertices and  is the set of M edges. In G, each vertex represents a rigid body and each edge represents a rotatable bond. Two vertices are connected by an edge if their associated rigid bodies are connected by the rotatable bond that the edge represents. We will therefore use T i to refer to both torsion i and its associated edge in the rigid-body graph. We further assume that each torsion is free to rotate independently of others, thus restricting the presence of ring systems or other cycles in the molecular graph to individual rigid bodies. Under this assumption, the rigid-body graph has no cycles and is a tree. An example of a simple molecule and its rigid-body graph is shown in Fig. 1.
The search space of the conformational search problem is a hypersurface described by an energy model (function)  Θ → U: M . For a given t, U t ( ) is the molecular energy consisting of the sum of all interatomic potentials in the molecule (e.g., van der Waals, torsional, and bending), which are dependent on the relative coordinates of the atoms. Various energy models or effective force fields can be used for our purpose, such as the widely used "Universal force field" (UFF) 36 . The conformational search problem, with the objective of finding the global minimum-energy conformer of a given molecule, can then be formulated as www.nature.com/scientificreports www.nature.com/scientificreports/ Upon changing the torsion angles of the molecule, some of the interatomic potential contributions will remain unchanged, while other contributions will change depending on the torsion angle vector t. To be more specific, let us denote the values of the torsion angles on the path connecting R a to R b on the rigid-body graph by a vector t ab . The length of this path is represented by m ab , meaning that t ab has m ab elements. The relative positions of the atoms in R a with respect to the atoms in R b depend only on the torsions on this path. Now, where U a is the sum of the interatomic potentials of the atoms within rigid body R a , which is invariant under rotation, while  Θ → U : ab m ab such that for a given t ab , U t ( ) ab ab is the sum of the interatomic potentials of all pairs of atoms where one atom is in R a and the other is in R b .
Binary optimization formulation for the conformational search problem. In order to use a quantum annealer to solve the conformational search problem, one needs to reformulate it as a quadratic unconstrained binary optimization (QUBO) problem (see the Supplementary Information for more details). To this end, we start by applying a one-hot encoding to the discrete values of the torsion angles, establishing a mapping between the torsion angle vector space and a binary solution space. That is, for each t i , we assign a binary variable As a result, t i can be expressed as to ensure t i takes one and only one value at a time. The constraint (5) is commonly referred to as a one-hot encoding constraint. Note that after applying the one-hot encoding, any arbitrary function f t ( ) i can be written as Similar to Eq. (6), a binary representation for U t ( ) ab can be found. For simplicity of presentation, let us assume that the torsion angles in t ab are indexed sequentially from 1 to m, that is,   Figure 1. Illustration of the rigid-body graph for a simple organic molecule (pregabalin). The molecular structure, with seven rotatable bonds highlighted in yellow, is shown on the left. The rotatable bonds connected to the two methyl groups have been discarded for visual clarity. The associated rigid-body graph is depicted on the right. www.nature.com/scientificreports www.nature.com/scientificreports/ Substituting U t ( ) ab ab from Eq. (7) into Eq. (2) results in a representation of the molecular energy U t ( ) in terms of the binary vector x. We denote this representation of the energy function by To solve the above binary optimization problem using a quantum annealer, one faces three challenges. First, the objective function in formulation (8) is not necessarily quadratic as U t ( ) ab ab may depend on more than two torsions, that is, > m 2 ab . Second, it is a constrained binary optimization problem. These two challenges indicate that the problem cannot be solved directly on a quantum annealer. The third challenge is that if is not much smaller than M, constructing an instance of formulation (8) becomes very computationally expensive due to the pre-evaluation of the coefficients in Eq. (7). In the following section, we propose a method that addresses these challenges in order to be able to use a quantum annealer for solving the conformational search problem.

Variable neighbourhood Descent for the conformational Search problem
Neighbourhood search, or local search (LS), is known to be an effective heuristic algorithm for solving a large number of combinatorial optimization problems. In defining a neighbourhood relation between solutions of a problem, local search begins from an initial solution and iteratively explores the neighbourhood of the current solution for improvement. It has been shown that a solution produced by a local search algorithm will often not be globally optimal, but will be suboptimal with respect to another neighbourhood relation 37 . When multiple neighbourhood relations are considered, the algorithm is often referred to as variable neighbourhood search 38 . In the context of conformational search, a solution refers to a vector of torsion angles t.
Let N k , for ∈ … k K {1, , }, denote a finite set of neighbourhood structures and N t ( ) k be the set of all solutions in the k-th neighbourhood of t. Starting from an initial solution and a neighbourhood structure, in each iteration, variable neighbourhood descent finds the best solution in the neighbourhood of the current solution. It then updates the current solution with the best solution found, and changes the neighbourhood structure before proceeding with the next iteration. The VND method is summarized in Table 1.
A simple neighbourhood structure is obtained by considering two solutions as neighbours if and only if they differ by exactly one torsion angle value. Such a VND heuristic is exactly the LS heuristic described above. While computationally inexpensive, the performance of LS can suffer in cases where a decrease in the molecular energy cannot be achieved by changing only a single torsion angle value in an iteration. Our proposed VND method improves upon LS by exploring more-complex neighbourhoods. In the following, we describe the components of the method. initial Solution. The initial solution in Table 1 can be selected in a variety of ways. One may simply choose a randomly generated torsion angle vector for the given molecule as the initial solution. Alternatively, one can use a greedy construction method. Another approach is to start from a known high-quality solution. This applies when using our VND method in conjunction with another conformational search method or by exploiting some prior knowledge about a given molecule Table 2. neighbourhood Structures. We now describe a more powerful neighbourhood structure which, to our knowledge, has not been previously studied. Let = G ( , ) R T be a rigid-body graph ′ ⊆   , and G′ be the graph resulting from contracting all edges in ′ \   (see Fig. 2 for an example). If G′ is a star graph (a tree graph with at most one vertex of degree >1), then we say that ′  has the property of 2-torsion dependency. The motivation for using this terminology is that any two vertices of G′ are connected with at most two edges (torsions). We label the maximal 2-torsion-dependent subsets of  as … , , K 1   and their associated star graphs as G 1 , …, G K . Then, the neighbourhood structure N k defines neighbourhoods containing all solutions which differ only in torsion angle values corresponding to edges in k  , for = … k K 1, , . Solutions in neighbourhoods defined by an arbitrary N k are also called neighbours under N k .   . The neighbourhood structure N 1 defines neighbourhoods containing all torsion vectors that differ only in torsion angle values for T 1 and T 3 . For example, =°°°t [5 , 10 , 20 ] and ′ =°°°t [0 , 10 , 90 ] are neighbours under N 1 while =°°°t [5 , 10 , 20 ] and ′ =°°°t [5 , 15 , 20 ] are not because they differ in the torsion angle value for T 2 .
neighbourhood Search. Based on the discussion in Section 2.2, the problem of finding the best solution in N t ( ) k can be formulated as a QUBO problem by restricting the optimization problem in (8) to the binary variables corresponding to the torsions in  k and moving the one-hot encoding constraints to the objective function using the quadratic penalty method as follows: terms represent the interaction energy of the two vertices of G k that are connected by T i and T j when terms represent the interaction energy of the two vertices connected by T i on G k when θ = t i k i and p is a sufficiently large penalty coefficient that enforces the one-hot encoding constraints. The above QUBO problem can be solved using various methods [39][40][41] , as well as specialized hardware devices such as quantum annealers. neighbourhood change. At each iteration, the neighbourhood is selected based on a random ordering, , of the torsions. The pseudocode for the neighbourhood change function is given below. www.nature.com/scientificreports www.nature.com/scientificreports/ practical considerations. Here, we detail the practical considerations of the proposed VND method to accommodate the use of existing and future quantum annealers.
The formulated QUBO problem in (9) is fully connected, meaning that for all i and j ( ≠ i j), the term x x ik jk i j appears in the objective function. On the other hand, the connectivity of the qubits on the D-Wave 2000Q follows a "Chimera graph"; thus, the problem in (9) must be embedded onto the hardware graph using an embedding strategy 42 . There is a limit on the number of variables a fully connected QUBO problem that can be embedded onto the graph can have. We take this limitation into consideration by imposing a limit on the number of variables in the formulated QUBO problem for each selected subset of torsion angles k  . We denote this parameter by s. In the following, we explain how this limit is imposed on the formulated QUBO problem.
For the selected neighbourhood structure at each iteration, we (randomly) select a total of s discrete values. That is, for each T i in  k , we randomly choose a set of Θ ⊆ Θ We have already defined the N t ( ) k neighbourhood as the set of torsion angle vectors t′ that are different from t only in the angle values of the torsions in k  . In addition, for any T i in  k , ′ t i (i.e., the value associated with T i in t′) takes on values only from Θ i .
With the above choice of neighbours, for an arbitrary t, N t ( ) solutions. However, to find the best solution in N t ( ) k , we need to pre-evaluate only i j T T i j i j , : , i j k energy terms to formulate the QUBO problem. As seen above, S k grows linearly with the product of |Θ | i , for all i, whereas the growth of (11) is quadratic. This means that the number of energy pre-evaluations grows more slowly than the size of the neighbourhood as  | | k increases. Another practical consideration for our proposed VND method is the stopping criteria. The first stopping criterion sets a limit on the computational effort of the method by introducing a maximum number of iterations, denoted by B. The second stopping criterion aims to terminate the method early if it becomes stuck at a local minimum or finds the global minimum of the problem. For this purpose, we introduce a parameter called A that represents the maximum number of consecutive iterations to can be performed without decreasing the energy.
With the above-mentioned practical considerations, the implemented VND method is summarized in Table 3. Two important features of the VND method is its scalability and ease of adaptation to the size and connectivity of the quantum annealer. More specifically, one can simply increase s in the described method to take advantage of improvements in the number of qubits and their connectivity. The rest of the method remains intact. It is worth mentioning that although the focus here has been on quantum annealers, other QUBO problem solvers could be used.
Effect of the molecular structure. As previously discussed, due to the limitations of current quantum annealers, we are able to jointly optimize only a subset of torsions that are 2-torsion dependent at each iteration of VND. The potential speedup of the quantum annealer over a naïve exhaustive QUBO problem solver depends on the cardinality of the selected subset of torsions. Since is fixed, the number of solutions in N t ( ) k , S k , is maximized when the cardinality of each Θ i is the same. That is, when

resulting in the maximum number of solutions
, where e is the base of the natural logarithm. As the potential speedup of the quantum annealer over an exhaustive QUBO solver is dependent on S k , it is favourable to select each  k with a large cardinality. The speedup diminishes when  | | = 2 k , which is the case for molecules with linear rigid-body graphs. The above discussion suggests that star-like molecules that have rigid-body graphs with high-degree nodes stand to benefit more (than those that do not) from using a quantum annealer to solve the QUBO problem at each VND iteration.

Experimental Results
In this section, we evaluate the performance of the proposed VND method. We first provide details about the molecules used in our experiments and then present the experimental results. trial molecules. The performance of VND was evaluated using a testbed containing the nine molecules depicted in Fig. 4. These model systems include three organometallic molecules useful for catalyzing reactions relevant to industry (labelled "A" 43,44 , "B" 45,46 , and "C" 47,48 ), a set of three n-alkanes whose basic structural motif appears in fuels, lubricants, solvents 49 , and resins 50 (labelled by D, E, and F), as well as a set of three ortho-phenylene oligomers that are of interest as electronic materials and nanomaterials (labelled "G", "H", and "I") 51,52 .
The choice of model systems A-I was motivated by several considerations. First, it is important to show that our method can be applied to a wide variety of conformational search problems of relevance to industry. Second, the model systems are representative of a diversity of molecular graphs: systems A-C have star-like graphs, whereas D-I have linear graphs, albeit with different structures. Third, the model systems represent a significant variety of active torsions, ranging from very modest (e.g., A) to very substantial (e.g., F and I). Another motivation for choosing these systems was the existence of experimental data pertaining to their three-dimensional structure (see the Supplementary Information for details). energy model. In our experiments, the Lennard-Jones 6-12 potential is used to model the interaction energy of two atoms α and β as 12 6 where ε αβ is the depth of the potential well, σ αβ is the van der Waals bond length, and αβ r is the distance between the two atoms. An expression for Eq. (2) is then obtained by summing over all pairs of atoms. Values for the parameters ε αβ and σ αβ are taken from the UFF 36 without modification. One should note that while we use the Lennard-Jones potential in what follows to provide a proof of concept for the proposed VND method, other force fields, that allow for representing the molecular energy as Eq. (2), can also be used.

Results.
We present the results of solving the conformational search problem for the selected molecules using the proposed VND method. Since the D-Wave 2000Q quantum annealer is not guaranteed to find an optimal solution, we first present results for VND, where an exact QUBO problem solver is used. These results focus solely on the performance of the VND method as a conformational search approach by preventing any deterioration of the results attributable to the use of a device that has imperfections. For comparison purposes, we also present the results achieved by performing a local search heuristic (LS), and a hybrid of the two methods (LS-VND) in which a random conformation is first optimized using LS and then passed to VND for further optimization. The comparison with LS is helpful in understanding the improvement that can be achieved through the use of a more complex neighbourhood than what is used for LS. We also present the VND results with the quantum annealer used as the underlying QUBO problem solver. Finally, we compare the results with those found by our implementation of the PTMC algorithm for the conformational search problem, a state-of-the-art metaheuristic for conformational sampling. Details of the parameters used for VND, the quantum annealer, and PTMC are presented in the Supplementary Information.
Reference Conformations. We compare our results against reference conformations found using a PTMC conformational search method. More specifically, we use an initial conformation for each of the selected molecules as input to a PTMC algorithm and let it run with a sufficiently large number of sweeps, such that the resulting www.nature.com/scientificreports www.nature.com/scientificreports/ reference geometry can be assumed to represent the absolute minimum conformation with a high degree of confidence. Details of this procedure are given in the Supplementary Information. To ensure the fairness and accuracy of this comparison, the reference conformations are generated using the same potential energy model as the one we used for VND.
Performance Metrics. The following metrics have been used for performance evaluation.
• Success rate: the fraction of runs that found a conformation with an energy within 1 kcal/mol (roughly chemical accuracy) of the reference conformation's energy. • Number of energy evaluations: the number of molecular energy evaluations needed to arrive at the best found conformation, averaged over all runs. • Residual: the energy difference between the best conformation found in a run and the reference conformation. The normalized residual, when reported, refers to the ratio of the residual to the number of atoms in the system. • Time to solution (TTS): the time it took to find the best solution (conformation) in a single run.
Each run was terminated once it found a conformation within 0.1 kcal/mol of the reference conformation, or if it reached some other stopping criterion.
VND vs. LS vs. LS-VND. As a baseline for comparison, we used an exact QUBO problem solver to optimize the selected set of torsions at each VND iteration in runs of both VND and LS-VND. The success rate, residual, and TTS for each of the three methods, for all nine model systems, are presented in Table 4. Here, the TTS for VND includes the time needed for selecting the neighbourhoods, formulating their associated QUBO problems, solving the QUBO problems with an exact solver, and translating the solutions of the QUBO problems back to torsion vectors. The number of energy evaluations reported for the VND method is based on the number of energy evaluations required to find the θ θ U ( , ) ij k k i j and θ U( ) i k i coefficients of the QUBO problem in (9). As shown, whereas LS is faster than VND and LS-VND, its success rate and residual results are generally inferior to the other two methods. Our implementation of LS terminates when there is no neighbouring solution which has a lower energy value, ensuring that the algorithm terminates at a locally optimal solution with respect to the single-torsion neighbourhood. Further, since LS-VND combines both the LS and VND methods, it is not surprising that its success probability and residual are generally at least as good as those of LS and VND.
VND Using a Quantum Annealer. The results for VND using a quantum annealer as the underlying solver are presented in Table 5. Here, the TTS includes the time spent on all pre-and postprocessing steps needed to solve the QUBO problems on the quantum annealer. Note that the actual time spent using the quantum www.nature.com/scientificreports www.nature.com/scientificreports/ annealer, referred to as the annealing time, is much smaller than the reported TTS. In addition, to aid in visually assessing the found lowest-energy conformations, their graphical representations are provided in the Supplementary Information. A comparison of these results with the results of VND when using an exact QUBO problem solver (see Table 4) shows that those of the former are of lower quality. It is expected, however, that spending more effort on tuning the quantum annealer's parameters would improve these results (See Supplementary Information for more details). Another observation from Table 5 is that VND employing a quantum annealer to solve QUBO problems has a larger TTS than VND using an exact QUBO problem solver. The quantum annealer solves a QUBO problem much faster than the exact solver, so one might think that its TTS should also be lower.
To explain this observation, one should note that it is not sufficient to solve QUBO problems merely fast, because if they are not solved optimally by the quantum annealer, the VND method may take longer to converge. Further, a significant portion of time is spent on transforming the QUBO problems into Ising problems (see the Supplementary Information), communicating with the quantum annealer, and mapping the results from the quantum annealer back to the logical bits. It is worth noting that elapsed real time is a fair measure of the time required to solve actual problems using a quantum annealer. This is in contrast to the customary approach of reporting only the annealing time, which is very small in comparison (on the order of a few microseconds per annealing cycle).
There are several known factors that make the QUBO problems generated at the VND iterations challenging for the current generation of quantum annealers. First, the problem graphs are fully connected, whereas the quantum annealer's connectivity graph is extremely sparse, resulting in each logical bit being embedded onto chains of 17 qubits on the quantum annealer. It is difficult to maintain identical states for those qubits, resulting in a higher error rate. Second, these problems have a very large range of coefficients, due to the r 1/ 12 and r 1/ 6 terms in the Lennard-Jones 6-12 energy model. On the other hand, the couplers of the existing quantum annealers have a limited bit precision and a fixed range. The large range of coefficients results in a loss of precision, which manifests itself in a lower success probability. Third, the current generation of quantum annealers has a high level of noise, referred to as intrinsic control error (ICE), leading to a significant loss in precision. Future quantum annealers are expected to mitigate these factors, with more-dense hardware graphs, higher bit precision, and lower ICE levels.   Table 6. PTMC had a high success rate for all model systems except for system F. It is worth noting that PTMC's TTS and number of energy evaluations are generally significantly higher than that of VND.
Effect of the Neighbourhood Size on VND's Performance. In the VND experiments discussed above, we restricted the neighbourhood size, s, in each VND iteration to 63, the largest size of QUBO problem, whose underlying graph is complete, that can be solved using an equal-length embedding on the quantum annealer. In what follows, we report on the effect of s on the performance of VND. This is useful in predicting the performance improvement achievable by increasing either the number of qubits or their connectivity. Table 7 presents the results for molecules B and C for different neighbourhood sizes when an exact solver is used to solve the QUBO problem at each iteration. The reason for choosing these two molecules is that they have a star-like structure and, as discussed in Section 3.6, the advantage of using a quantum annealer over an exact solver is expected to be more pronounced for these molecules.
As shown, increasing the neighbourhood size from 30 to 60 and then from 60 to 90 noticeably improves the results for both molecules. However, the improvement exhibits diminishing returns when the neighbourhood size is increased beyond 90. We expect similar behaviour to occur at different neighbourhood sizes for different families of molecules.  Table 7. Effect of the neighbourhood size on the performance of VND. For each model system, we report the minimum, median, and 75th percentile of the residual, over 500 runs.  Table 5. Results for VND using a quantum annealer as the underlying QUBO solver. For each model system we report the minimum, median, and 75th percentile of the residual and TTS, over 25 runs.  Table 6. Results for PTMC. For each model system, we report the minimum, median, and 75th percentile of the residual and TTS, over 100 runs.

conclusion
In this paper, we have presented a variable neighbourhood descent (VND) method for conformational search. We introduced the concept of a rigid-body graph and used this simplified molecular structure to carefully define a neighbourhood structure to allow for efficient optimization using a binary quadratic optimizer. Based on current quantum annealing hardware, we selected a 2-torsion-dependent neighbourhood at each iteration such that finding the best solution in the selected neighbourhood could then be formulated as a QUBO problem. The size of the neighbourhood can be chosen such that the method can be adapted to the number of available qubits as well as to their connectivity on the quantum annealer. As a result, the proposed method is not only well-suited for current hardware, but can easily be adapted to take advantage of hardware improvements. Whereas the proposed method can be used as a standalone conformational search approach, it can be combined with existing conformational search methods for potentially improved performance. Beyond a simple presentation of the method, we also conducted a preliminary case study based on an implementation of the VND method using the D-Wave 2000Q quantum annealer for two families of molecules. In this exploration, we compared the results of our method with those of PTMC, a state-of-the-art solver for conformational search. To understand how much of the gap between the results of PTMC and those of VND used with a quantum annealer can be attributed to the imperfections of the quantum annealer, we replaced it with an exact QUBO problem solver. VND used along with the exact solver was able to find noticeably better conformations than those found using VND and the quantum annealer together. This observation points to the potential improvement achievable in the short term through more-advanced tuning of the existing quantum annealer, and in the long term using improved hardware.
This work suggests a number of possible future research directions. For example, investigating refined neighbourhood change functions rather than using a random ordering of the torsion angles to choose the neighbourhood could lead to significant improvements. This could involve further exploitation of the molecular graph or the solutions from previous VND iterations. In addition, an improved selection of the s discrete points at each VND iteration to account for hardware limitations could further improve the results. We leave these improvements for future work.
We believe that the proposed method, based on careful hardware-aware neighbourhood selection, holds the potential to provide promising solutions to important optimization problems. While PTMC shows better performance over the molecules studied, our method has opened a scalable path forward for leveraging emerging quantum technologies for conformational search, a critically important problem in the field of chemical and materials science.