key: cord-0045662-0iy0zxjs authors: George, Jithin; Di, Zichao (Wendy) title: Trilateration-Based Multilevel Method for Minimizing the Lennard-Jones Potential date: 2020-05-25 journal: Computational Science - ICCS 2020 DOI: 10.1007/978-3-030-50426-7_13 sha: c32311b90d76c718c4b49e364fdfc482be63a2cf doc_id: 45662 cord_uid: 0iy0zxjs Simulating atomic evolution for the mechanics and structure of materials presents an ever-growing challenge due to the huge number of degrees of freedom borne from the high-dimensional spaces in which increasingly high-fidelity material models are defined. To efficiently exploit the domain-, data-, and approximation-based hierarchies hidden in many such problems, we propose a trilateration-based multilevel method to initialize the underlying optimization and benchmark its application on the simple yet practical Lennard-Jones potential. We show that by taking advantage of a known hierarchy present in this problem, not only a faster convergence, but also a better local minimum can be achieved comparing to random initial guess. Simulating interactions between atoms/molecules is particularly essential for understanding materials properties in materials science, chemistry, and biology. However, such problem presents an ever-growing challenge related to the huge number of degrees of freedom borne from the high-dimensional spaces in which increasingly high-fidelity material models are defined. One way is to simulate the arrangement of atoms by minimizing the potential energy. Among various mathematical models describing the interaction, the Lennard-Jones (LJ) potential [9] has attracted a lot of theoretical and computational attention due to its mathematical simplicity yet practical importance of discovering low energy configurations of clusters of atoms. There are two main difficulties in solving this problem. First, the non-convexity and highly non-linearity of many variables lead to a large number of local minima [3, 13] , and second, the intensive computational burden of many variables leads to slow convergence [3, 6] . Therefore, an acceleration of LJ potential with potentially better local minimum is desired. Given the domain-, data-, and approximation-based hierarchies present in this problem, it is natural to exploit them to speedup the convergence of LJ potential. One known property is that the Mackay icosahedron [4, 10, 19, 20] dominates the structures of LJ clusters at least in the size range of 10-150 atoms. In this work, we observe similar property as icosahedronal packing by examining the local minimum obtained by truncated-Newton method [12] using random initial guess. This observation inspires us to adapt the successive refinement [2, 15] developed in traditional multigrid method to accelerate the LJ convergence. Therefore, we propose a novel interpolation method as part of a successive refinement framework so that a better initial guess can be obtained for a larger system from the solution of a smaller system. We benchmark our method on various sizes of systems and show that the proposed method can lead to a much faster convergence. Furthermore, in our examples, we show that the proposed hierarchical approach can potentially lead to an overall better local minimum with a lower energy than using random initial guess. Atoms arrange themselves spatially to minimize a potential that is an accumulation of all the different kinds of attraction and repulsion at that scale. The Lennard-Jones potential is a simplistic attempt to model such interaction. Given a configuration of N atoms in a cluster V = {v 1 , . . . , v N }, its potential is given by where v i and v j are the locations of the ith and jth atoms, and d(v i , v j ) is the Euclidean distance between atoms v i , v j ∈ R 3 . The physical constants ε and σ are the depth of the potential well and inter-atom reaction limit, respectively, both depending on the type of atom. Because of its importance in computational chemistry, finding optimal configurations that locally minimize the potential energy (1) remains an active research area. For example, Maranas et al. [11] proposed an exotic optimization algorithm to find many stationary points, and Asenjo et al. [1] studied the mapping of the basins of attraction for various optimization algorithms. In this work, we follow the longstanding convention in [8] and consider the reduced-unit optimization problem The corresponding gradient is given by To avoid unnecessary degrees of freedom, we fix atom 1 at the origin, atom 2 on the x-axis and atom 3 on the x-y plane. Figure 1 illustrates the classical Leonard-Jones potential E 2 . We can see that even minimizing the smallest system can be challenging due to either the singularity as v i − v j → 0 or stationary but not local minimum as v i − v j → ∞ . We can further extend this visualization to E 3 , where we fix v 1 at the origin and v 2 on the y-axis at the optimal distance of r * = 2 1 6 (discussed more in Sect. 3) from the origin, then we examine the optimization performance by freeing v 3 on x-y plane. Again, through the paper, we use truncated-Newton method [12] to perform all optimization procedure due to its robust performance. Figure 2a shows the number of iterations (indicated by color intensity) needed for convergence at different initial locations of v 3 on x-y plane. Figure 2b shows the corresponding energy (indicated by color intensity) obtained while optimization converges. Again, even for the three atoms case, the basins of attractions are so rich in its complexity that a good initial guess is desired for a fast optimization process. One standard practice for initialization of the minimization of E N is by using random initial guess. The random locations of the atoms are obtained by uniform sampling in a box of size 2 √ 3r * [1] . As a side note, we compared it to another practice where the box sizes vary as N 0.8442 1/3 [21] (0.8442 is a specific density value). We observe that the fixed box has comparable performance for the range of number of atoms we study, therefore we use it through the paper. There is often visual modularity associated with the LJ minima. This fact inspired many researchers to leverage the geometric principles of LJ minima to generate good initial guesses (e.g., [17, 22] ). For example, for minimizing E N , instead of random initial guess, we initialize the iterate by simply adding an atom randomly to a local minimum of E N −1 . As we can see in Fig. 3 , compared to a random initial guess, using the initial guess built upon the local minimum of E N −1 requires much fewer number of iterations to converge, without sacrificing the quality of the local minimum. This observation prompts us to search for a multilevel approach to the LJ optimization, where one could generate good initial guesses by exploiting the hierarchical structure embedded in different systems, i.e., by interpolating solutions from small systems to serve as initial guess for large system, we could expect a more robust convergence. Let us first consider the pairwise interaction between atoms v i and v j , which is expressed by the following: It is easy to derive that the minimum of Eq. 4 is obtained when ||v i −v j || = r * = 2 1/6 , where r * is the ideal distance of separation between two atoms. Therefore, the global minima of E 3 and E 4 in 3D are equilateral triangle and tetrahedron, respectively, with pairwise distance r * . For larger systems, we first define the neighboring atoms of the atom v i as the ones connected with v i by edges of Delaunay triangulation [7] . Notice that the optimal pairwise distance of r * is not maintained between all neighboring atoms for N > 3 in 2D and N > 4 in 3D. Then an intuitive expectation is that an arrangement of the atoms, where all the pairwise neighboring distances are approximating r * , can be a good initial guess. This motivates the main idea we propose in this paper, which is to add new atoms to the boundary of a coarser atomic configuration (i.e., a smaller system with fewer atoms) at its local The number of iterations needed to find a local minimum for EN in 3D greatly reduces when a local minimum of En−1 in addition to a random point is used as the initial guess minimum so that the distance between the new atom and its nearest neighbours from amongst the existing atoms is r * . We employ the method generally known as trilateration to perform this action. Simply speaking, in 2D, we choose two neighboring atoms (i.e., connected by an edge in the Delaunay triangulation) and draw circles around them with radius r * , then the two intersection points of the circles provide two candidates as the new atoms to be added. For example, given two atoms at (0, 0) and (d, 0) where d ≤ 2r * , respectively, then the intersection points of trilateration . This idea is illustrated in Fig. 5 and can be extended to 3D. This trilateration technique of circle-circle intersection in 2D or spheresphere-sphere intersection in 3D has been used in various applications such as surveying and GPS systems [5, 14, 23] . After performing trilateration on every edge on the boundary of the atomic configuration, a list of candidate locations for the new atoms are obtained as shown in Fig. 6 . Each location is evaluated by adding an atom there and seeing if the resulted system has smaller energy, since the energies of stable systems should decrease as the number of atoms increase. If adding an atom at a location increases the energy of the system, this location is discarded. For simplicity, we summarize the detailed procedure (denoted as trilateration-based multilevel method) in 2D as follows (Fig. 4) , • Select edges that only have two atoms lying on it. • If an edge on the convex hull has more than 2 atoms on it, divide that edge into multiple edges comprising only of 2 atoms and select all the edges thus formed. • If the length of an edge is longer than 2r * , identify the triangular simplex that the edge belongs to. Select the other two edges in the triangular simplex instead of the original edge. If any of the edges formed this way have more than 2 atoms, divide it up into multiple edges with 2 atoms and select them all. -On each edge among the selected edges, use the trilateration method to find candidate locations to place the new atoms. -Take all the candidate locations found this way and add them one at a time to the original system of atoms. If the corresponding potential energy decreases, the new atom becomes part of the system. Otherwise, it is rejected. -Once the new configuration is obtained with m atoms where m > N, it is used as an initial guess for minimizing E m . One could also perform Delaunay triangulation and the trilateration procedure again on the new system to get a larger system of atoms. -This process can be repeated until we reach a desired system size. -Perform optimization on the final iterate giving us the generated initial guess to find a minimum. If the user desires a minimum for a particular number of atoms, this method can be applied to find a minimum for a system of size close to the desired number. Then, the remaining few atoms can be added either randomly or using trilateration as a heuristic. The resulted configuration can still serve as a good initial guess as suggested by Fig. 3 . The selection of the edge for trilateration is a critical step. For example, trilateration does not provide any points of intersection if the edge on the boundary (obtained from the convex hull) is longer than 2r * . If this fact is not taken into consideration, eventually the interpolation will stagnate as shown in Fig. 8 . Therefore, in the proposed algorithm, we provide a threshold to carefully choose the right boundary for trilateration, i.e., if the boundary edge is too long (e.g., >2r * ), the algorithm opts into the concave boundary as suggested in Fig. 7 . This guarantees our method to be scalable and avoids problems as shown in Fig. 8 where the boundaries are too long for the algorithm to perform trilateration anywhere. In terms of the computational complexity added by the proposed multilevel method, since the pairwise distance between all atoms is already available by calculating the potential, it is trivial to perform the threshold criteria for edge selection. The trilateration for each selected edge has a relatively simple analytical formula, therefore, this part of calculation is negligible as well. The main cost then is on the generation of the convex hull and the Delaunay triangulation which are O (N log N ) . Since these operations are one-time upfront cost, the overall added complexity by the proposed trilateration-based multilevel method for initialization is negligible comparing to one iteration cost of the underlying optimization. In this section, we demonstrate the performance of our proposed method for various sizes of systems primarily in 2D, with a preliminary exploration on 3D. Notice that in all experiments, we perform the optimization at the final step when the number of atoms reaches a desirable number. For the optimization step, we emphasize that any general-purpose optimization methods can be used in conjunction with the multilevel approach proposed in this paper. For our illustration purposes, we use the Truncated-Newton algorithm in the Scipy Python package [18] . Again, as for comparison, we evaluate the performance of our proposed method against randomly initialized guesses which is a standard practice for initialization. Fig. 9 . The multilevel process to extend a small system to a large system based on the proposed trilateration approach. Figure 9 demonstrates the iterates obtained at every step of the trilateration process with the initial system as the global minimum of E 3 . In Fig. 10 , we compare the number of iterations needed for the convergence of the optimization process for E N when the initial guess is chosen randomly versus if the initial guess is obtained using the trilateration procedure. We see dramatic savings in the number of iterations required. 11 . The energy of the local minimum obtained by using the initial guess from the trilateration process as opposed to taking a random initial guess. The question that remains is if the minimum found using the initial guess from trilateration is better than a random initial guess. Figure 11 shows that in average, the minimum obtained by the proposed trilateration has a lower energy than the one from averaged trials of random initial guesses. We also demonstrate a preliminary extension of the proposed method to 3D. As briefly shown in Fig. 12 , the reduced number of iterations needed for convergence by trilateration is maintained in the 3D case. In Fig. 13 , we use the energy of the global minimum extracted from the online database [20] as a reference, and compare it against the energies of the minima found by random initial guesses and those found through the trilateration method, respectively. Due to the increased geometric complexity in terms of locating the concave regions arising in 3D, the edge selection criteria in 2D can not be trivially extended to Fig. 13 . The energy of the local minimum obtained by using the initial guess from the trilateration process as opposed to taking a random initial guess in 3D. 3D, therefore, we only show the results up to 35 atoms. We leave the development of more optimal way for edge selection to the future. Although computational resources have grown significantly over the last few years, their growth alone will not suffice to address the volume and complexity encountered in applications such as molecular dynamic simulations. Multilevel methods presents one avenue for taking advantage of a known hierarchy to overcome the intensive computational burden bearing in such problems. In this work, we propose a trilateration-based multilevel method to speed up the convergence by providing a better initial guess, and prototype its benefit on the minimization of the Lennard-Jones potential. We exploit the hierarchies embedded in the atomic configuration with different sizes so that a good initial guess can be obtained by interpolating through trilateration from a local minimum of a smaller system. We observe that for 2D, the interpolated initial guess can save the number of iteration greatly comparing to random initial guess, without sacrificing the quality of the minimum. Although only preliminary results for 3D are shown for limited sizes of system, it is promising that the proposed method can be scaled to practical use. However, a more careful design of the framework is needed, for example, how to more optimally choose the edge to perform trilateration is critical for the scalability of the proposed method. Although the current work can not immediately address the challenge of finding global minimum, one simple way to take advantage of our trilaterationbased multilevel method is to plug in to a more complicated framework for global optimization such as the one proposed in [3] . For example, the random initial guess can be simply replaced by the proposed method in this work. Alternatively, one can utilize multistart framework for global optimization (e.g, [16] ) with the proposed trilateration-based multilevel method for each subproblem. One could also investigate better interpolation strategy to combine trilateration with physics to guarantee a better local minimum. We leave these questions to be the focus for the future. Visualizing basins of attraction for different minimization algorithms The cascadic multigrid method for elliptic problems Optimization of Lennard-Jones atomic clusters Structural optimization of Lennard-Jones clusters by a genetic algorithm Trilateration and extension to global positioning system navigation Size-temperature phase diagram for small Lennard-Jones clusters Bulletin de l'Académie des sciences de l'URSS. classe des sciences mathématiques et naturelles Physical cluster mechanics: statics and energy surfaces for monatomic systems On the determination of molecular fields, II: from the equation of state of a gas A dense non-crystallographic packing of equal spheres A global optimization approach for Lennard-Jones microclusters A survey of truncated-Newton methods Data compression for optimization of a molecular dynamics system: preserving basins of attraction Revisiting trilateration for robot localization Scatter search and local NLP solvers: a multistart framework for global optimization Packing schemes for Lennard-Jones clusters of 13 to 150 atoms: minima, transition states and rearrangement mechanisms SciPy 1.0-Fundamental Algorithms for Scientific Computing in Python Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms Global optimization of clusters, crystals, and biomolecules ESPResSo 4.0-an extensible software package for simulating soft matter systems Structural distribution of Lennard-Jones clusters containing 562 to 1000 atoms Beyond trilateration: on the localizability of wireless ad-hoc networks Acknowledgments. This work was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of two U.S. Department of Energy organizations (Office of Science and the National Nuclear Security Administration) responsible for the planning and preparation of a capable exascale ecosystem, including software, applications, hardware, advanced system engineering, and early testbed platforms, in support of the nation's exascale computing imperative. The material was also based in part on work supported by the U.S. Department of Energy, Office of Science, under contract DE-AC02-06CH11357.