Tribal Competition Evolutionary Search as a Surrogate for Free

Subscriber access provided by UNIV OF LOUISIANA

A: New Tools and Methods in Experiment and Theory

"Teamwork Makes the Dream Work": Tribal Competition Evolutionary Search as a Surrogate for Free Energy Based Structural Predictions Troy David Loeffler, Henry Chan, Stephen K. Gray, and Subramanian K.R.S. Sankaranarayanan J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.9b00914 • Publication Date (Web): 02 Apr 2019 Downloaded from http://pubs.acs.org on April 3, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

”Teamwork Makes the Dream Work”: Tribal Competition Evolutionary Search as a Surrogate for Free Energy Based Structural Predictions Troy D. Loeffler, Henry Chan, Stephen Gray, and Subramanian K.R.S. Sankaranarayanan∗ Center for Nanoscale Materials, Argonne National Lab, Lemont, Illinois 60439, United States E-mail: [email protected],[email protected]

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Crystal structure prediction has been a grand challenge in material science owing to the large configurational space that one must explore. Evolutionary (genetic) algorithms coupled with first principles calculations are commonly used in crystal structure prediction to sample the ground state and metastable states of materials based on configurational energies. However, crystal structure predictions at finite temperature (T), pressure (P), and composition (X) require a free energy based search that is often computationally expensive and tedious. Here, we introduce a new machine learning workflow for structure prediction that is based on a concept inspired by evolution of human tribes in primitive society. Our tribal genetic algorithm (GA) combines configurational sampling with evolutionary optimization to accurately predict entropically stabilized phases at finite (T, P, X), at a computational cost that is an order of magnitude smaller than that required for a free-energy based search. In a departure from standard GA techniques, the populations of individuals are divided into multiple tribes based on a bond-order fingerprint and the genetic operations are modified to ensure that the cluster configurations are sampled adequately to capture entropic contributions. Team competition introduced into the evolutionary process allows winning teams (representing better set of individuals) to expand their sizes; this translates into a more expanded search of the phase space allowing us explore solutions near possible global minimum. Each team explores a specific section of the structural phase space and avoid bias on solutions arising from use of individual populations in a purely energy-based search. We demonstrate the efficacy of our approach by performing structural prediction of representative 2D two-body system as well as LJ clusters over a range of temperatures up to its melting. Our approach outperforms standard GA approaches and enables structural search under ”real non-ambient conditions” on both bulk systems as well as finite sized clusters.


Page 2 of 29

Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Structural searches 1–5 are a non-trivial problem due to the large configurational space that one must explore in order to find the stable structures. While many possible configurations exist, there are typically only a handful of stable structures for a given system. As such, a smart search algorithm is needed to sort through the massive amount of ”hay” to find the few ”needles” that are of interest to the researcher. In crystalline systems, the configurational space is usually defined by internal coordinates of the unit cell. Efficient procedures for exploring such space first involve generation of approximate structures followed by identification of regions that correspond to low-energy structures. In several instances, such procedures subsequently involve exploring the ’energy landscape’, although alternate approaches involve defining an objective based an energy function or other ’cost function’ or figure of merit. These can include descriptors such as bond lengths and coordination numbers or some physical constraints that are used to assess the viability of a structure. After defining the energy function or cost function, the crystal structure search usually involves searching the energy landscape using approaches such as simulated annealing, Monte Carlo basin hopping, genetic algorithms or more recently deep learning based search methods that are based on topological principles. 6,7 Among the various search methods, genetic Algorithms (GA) remain popular and have proven to be a useful optimization technique for high dimensional objective minimization. 8–12 GA works by creating pools of individual configurations/structures, also known as organisms, using operations such as the cross over and mutation. The cross over operation is where features of the two parent structures are taken and combined to form a new child structure. The mutation operation on the other hand takes a single parent structure and performs a slight perturbation. In the context of structural optimization, this can be moving an atom, changing unit cell dimensions, etc. In a typical GA run, these mutation and cross over



operations are performed until the structure pool reaches a sufficiently large size. Once this size is reached a set of pruning operations is applied to eliminate structures from the pool. This elimination is primarily done according to the score of the structure or in other words how well it minimizes the objective function. This pruning is done to eliminate suboptimal structures which ultimately leaves only the structures which are close to the desired outcome. 13 The primary advantage GA has is that it is a highly flexible algorithm which allows the user to define a wide variety of mutation and cross over operations in order to efficiently sample the large number of configurations that are possible. This includes methods for getting over local energy barriers which commonly hinder local optimization techniques. As such it is an excellent choice of algorithm for structural optimization. In modern times, GA is commonly applied to predict the lowest energy structure since energy is easily computed via quantum mechanical or classical forcefields. 14,15 As such, GA algorithms perform exceptionally well when the most energetically stable structure equates the free energy wise most stable structure. However, it is not the case when the stability of a structure is due to entropic reasons rather than energetic reasons. The statistical nature of entropic stabilization makes predicting entropically stable structures a non-trivial problem in structural optimization. 16,17 This is because unlike the energy, free energy calculations are much more compute intensive and therefore are a much more challenging quantity to optimize toward. For systems where free energy calculations are readily available and cheap enough, it is possible to simply replace the energy with the free energy as the GA objective function. This can be done in cases such as bulk crystal structures where methods such as the Einstein crystal method exists. 18–20 However, it is not always the case that cheap and easy to implement free energy calculations are readily available for whatever system a researcher might be interested in, and even if they do exist, they might be too expensive or too unruly to use for optimization purposes. As such coming up with methods that allows one to correctly 4 ACS Paragon Plus Environment

Page 4 of 29

Page 5 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


predict states which are low in free energy without detailed free energy calculations would be ideal for a great number of applications. The primary challenge of this problem is dealing with situations where entropic stabilization is dominant. Entropic stabilization occurs when a single minimal energy state is overcome by a sufficiently large collection of thermally accessible states. While these states may be higher in energy, the shear number of them is enough to overpower the small number of lower energy states. One can effectively view this as the weaker states ”teaming up” to form a collective state that is much stronger than any of the states are by themselves. It is still possible at low temperatures for a single strong structure to win out over a group of average structures, but at higher temperatures the larger groups tends to win out. From this observation it is postulated that the problem of free energy is not about optimization of a single objective, but rather the optimization of a group of objective functions. This idea can be directly compared to observations in team sports. While having a single elite player can help a team win games, it is generally the teams who have multiple talented players and good team cohesion typically win championships. While the idea behind using a ”tribal unit” for optimization has been attempted before, 21 the previous attempts were used to speed up convergence of normal GA algorithms or were used for feature detection in pattern classification. 22 Here, we use the idea of group evolutionary optimization to predict structural configurations which are the most stable configuration by using configuration energies of teams as a proxy for the system free energy. The grouping is performed by using bond orientational order as a fingerprint that reflects the local symmetries of various close-packed particle clusters. Our approach allows us to predict optimal configurations at finite temperature, pressure, and composition (P, T, X) while being computationally much cheaper than structural search based on free energies. To demonstrate the suitability of our approach, we test the ability of our workflow to represent free energies for a two-body system. For a more rigourous test, we also perform structural prediction of representative LJ clusters over a range of temperatures 5 ACS Paragon Plus Environment


Page 6 of 29

up to its melting point. Our global optimization approach is quite general and can enable structural search of both bulk systems and finite sized clusters under ”real conditions”.

Method Tribal Genetic Algorithm (TGA) The underlying idea behind GA is still used in the TGA algorithm. However, unlike GA where each structure is weighted by its score, TGA gathers similar structures into a group that must compete against other groups to stay in the pool. Basic genetic operations of mating and mutation algorithms are largely the same as a standard GA algorithm. However, unlike traditional GA, the method by which organisms are removed from the pool is significantly different. Organisms can exit the pool by one of two events that are called at random. The first is what is called a famine event. In this, a single group is chosen at random and members are removed from the pool until the famine event ended. The purpose of the famine event is to prune groups so that structures with low scores are removed on a regular basis. If this is not done, a scenario where a single good structure protected by a giant army of bad structures can occur, which inhibits sampling. An apt analogy would be like playing chess against an opponent if the opponent was able to start the game by filling every empty spot on the chess board with a pawn. The second event is what is known as a war event where two groups are chosen from the pool and are pitted against each other. For the war event, the probability of a group’s victory is determined by the sum of their score functions.

wgroup =

X

e−Ei/kB T

i


(1)

Page 7 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


This weight is the same derived from statistical mechanics to estimate a partition function over a large number of samples. The probability for a given group i to win a competition against group j is given by

Pi =

wi wi + wj

(2)

This functional form is known as the symmetric acceptance rule which can be used to satistfy detailed balance as an alternative to the Metropolis acceptance criteria for Monte Carlo simulations. When a group either loses in a war or has been chosen for a famine event, one of the structures of the group is chosen for removal. This is done according to the probability:

ui = Stop − Si + u0

(3)

ui Pi = P j uj

(4)

where Si is the score of a given structure i, Stop is the score of the best structure in the group, Pi is the probability of selecting an object for removal, and u0 is the weight when Stop = Si . This is chosen such that the best structures in the pool with have a low yet non-zero removal probability. Only one structure at a time is purged from the group and the group itself is only removed when it has no remaining members in the pool. Since traditional GA is largely used for energy optimization, it is often the case that the best structure is retained in the pool. However, when it comes to free energy optimization it is often the case that the 7 ACS Paragon Plus Environment


Page 8 of 29

best structure is not a member of the most stable group at a given temperature. As such the TGA algorithm has been designed to allow for the possibility that the most energetically stable structure will be removed from the pool. Therefore the retention of a structure in the pool is more dependent on the probability of generating a structure and the probability of a group’s repeated victory. This is to mimic the same stocatistic effects found in Metropolis Monte Carlo sampling. In addition, the probability of selecting a group for either a famine or war event is weighted such that groups that are just starting up are less likely to be selected than the groups which have already matured. This is to prevent a ”bully” effect where a freshly created group is immediately knocked out by a large group before the smaller group has a chance to grow. For the purposes of this study the groups were selected inversely proportional to the number of group members.

Simulation Details - Representative Test Systems The algorithm was applied to several test systems in order to measure it’s effectiveness at predicting optimal free energy configurations. These systems were chosen because it is possible to compute the free energy by other methods and thus the exact solution is known. 2D Two-body system: The first representative system was a 2D two body system where the interaction between the two particles is described by the energy function:

E(r) = 20 (r − 1.0)2 ∗ (r − 5.0)2 + 0.2(r − 1.0)2

(5)

With r being the distance between the two particles. For simplicity the energy is given in units of the Boltzmann constant. (I.e., If E = 1 then E/kB = 1K) This equation was


Page 9 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


designed to have two energy minima at r = 1˚ A and r = 5˚ A with the well at r = 1˚ A being lower in energy. This two well system can be seen in Figure 1. Since this is a simple two body system, the exact free energy for a given range can be readily obtained by numerically integrating the equation

e

−∆G/kB T

R r2

re−E(r)/kB T dr = R rr1 max re−E(r)/kB T dr 0

(6)

This was chosen as a toy model because the system’s temperature dependence will behave in a predictable fashion. At low temperatures the well are r = 1˚ A will be the most stable position since it is lower in energy. However, the well at r = 5˚ A has a larger number of configurations associated with it since the number of states scales at a rate of 2πr. As such as one increases the temperature to intermediate values eventually the most stable configuration should shift to the well at r = 5˚ A and eventually as the temperature increases to an extreme value the two particles should begin to completely dissociate. Cluster structure predictions: The second system is a practical GA application which is predicting optimal cluster configurations as a function of temperature. This is an application for which traditional GA has been used extensively 23,24 and is also generalizable into other structural prediction calculations for which GA has been used for. This system was chosen for it’s simplicity and the ability to directly compute free energy related properties. For spherically symmetric potentials, mono-atomic clusters of 13 and 55 atoms have an icosahedral geometry for their minimal energy structure. This is the kind of structure that a standard GA algorithm will predict as the most stable structure. However, if thermal effects are accounted for these clusters will only appear at very low temperatures. For temperatures as low as T ∗ = 0.4 a more disordered cluster configuration will be the most stable. Thus it provides a good test case for predicting thermal related effects.



Fingerprinting:To ensure proper clustering during the simulation the Stillinger cluster criteria 25 is enforced during the TGA run. This is the same connectivity criteria used in nucleation studies. 26,27 To be consistent with previous studies a neighboring distance of 1.5σ was chosen. Grouping Criteria: In order to create teams a grouping criteria is needed. For the purposes of this study, a reaction coordinate based grouping method was used. This was done to facilitate the comparison of the TGA’s results to exact free energy calculations. While this is convenient in this case, the grouping is not limited to only reaction coordinates. For the two body system the inter-atomic distance was used since this is a natural choice for the reaction coordinate. For the cluster system, the Steinhardt q6 order parameter was chosen as a measure of quantifying the local orienational order of the close packed clusters. 28 In this context the parameter space is divided into discrete bins of 0.2 ˚ A in length for the two body system and blocks of 0.02 for the cluster system. Each of these bins make up a given group and thus any structure whose reaction coordinate falls within the range of a given bin is assigned to that group. Free energy calculation: The exact free energy for the two body system can be computed directly using Eq. 6 while the free energy The free energy of a given cluster with respect to the q6 parameter can be gathered using a modified version of the AVUS-HR method. 29

Results and Discussion The general workflow of our ML framework (tribal GA) for structural prediction is illustrated in Figure 2. The idea is to take advantage of the multi-population technique and grouping to naturally capture entropic effects arising from a large collection of thermally accessible states. Natural degrees of freedom such as vibrational modes can effectively be viewed as a collection


Page 10 of 29

Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of states which are a small perturbation from a given equilibrium position. From this one can see that free energy prediction is not about indivdual stability, but rather the stability of a group of states. In a departure from the use of individual populations in traditional GA, we introduce the concept of an inter-tribe competition to the evolutionary process to predict the elite tribe; the global optimum solution exists in the sub-spaces searched by these elite tribes. The size of the tribe reflects its search power and the grouping is representative of the number of closely accessible configurational states. During the evolutionary process, the size of the elite tribes are enlarged via mating to provide more search power whereas the size of the worst performing tribe is reduced via a famine event, thereby maintaining a constant population size. It is worth noting that the penalty and award strategy allows the tribal GA to not only search the structural phase space locally, but also converge quickly towards the global optimal. These genetic operations and competitive process is propagated through the generations until we converge to desired fitness constant or the required number of generations have reached. The tribal GA algorithm for structural prediction is benchmarked using two representative examples. In the first, we apply our workflow to predict the minimum free energy configurations at various temperatures (10K, 100K, and 3000K) in a 2-D two-body system. At low temperatures, it should be noted that the free energy is dominated by the energetic term and as such the free energy curve for this system will look very similar to the energy curve given in Figure 1. This is seen in panel 1 of Figure 3 where the global free energy minimal is located at r=1˚ A as expected. The TGA was able to correctly predict that the groups at r=0.8˚ A and r=1.0˚ A were the highest scoring groups while the groups the next tier of high ranking groups were the distributed about the two minima. This provides a rudimentary sanity check for the algorithm since this is a result that could have been predicted by traditional GA albeit requires more number of evaluations. The results for the intermediate temperature simulation are shown in panel 2 of Figure 3. In this situation the global free 11 ACS Paragon Plus Environment


energy minima is located not at the global energy minima at r=1˚ A, but rather at the second minima at r=5˚ A due to the higher entropic contribution. In a similar manner the TGA was able to correctly predict the groups at r=5˚ A and r=4.8˚ A to be the most stable. While the correct ordering could not be directly predicted by traditional GA one could reasonably argued that in this case that traditional GA could identify both minima and once discovered separate free energy calculations could then be performed to determine stability. Thus, one additional case was performed which was the high temperature simulation. In this case, the free energy minima is not located at either of the energetic minima, but rather at a value between 5.5˚ A and 6˚ A. This is a clear case where energy optimization alone is not able to correctly predict the free energy minima since this is a minima created by entropic stabilization. The TGA was correctly able to identify the entropically stable minima in our representative two-body system. As shown in Figure 4, for the two body system the top groups generally distinguished themselves within a very short time frame even as the scores of each group continue to rise. Eventually the scores for each group levels out once the groups reach an equilibrium of sorts. For this system the correct ordering can be obtained very rapidly. We next apply our TGA algorithm to a cluster structure prediction problem since it reflects a more practical free energy optimization problem. Specifically, we use the TGA workflow in Figure 2 to perform structural predictions of a LJ cluster from room temperature up to its melting. The results from each of these different TGA simulations are reported in Figure 5. Since the search space is significantly larger compared to the toy problem case, it took several generations to finally sort out the correct ordering. This is visualized in Figure 6. While it was found that even though the correct ordering may take some time, the most stable groups generally separated themselves relatively quickly from the unstable groups. The groups which usually ”battle it out” for the top spot are generally groups that are close in free energy to the global minima. Snapshots of the top group at each different temperature are shown in Figure 7. At lower temperatures, the TGA predicts structures which are close 12 ACS Paragon Plus Environment

Page 12 of 29

Page 13 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in shape to the icosahedral configuration which is lowest energy structure for this system. While they have many of the same characteristics of a perfect icosahedral, there are still small defects present at T*=0.2. At T*=0.2 the free energy profile is very similar to the energy profile shown in Figure 8 and as such one might expect that traditional GA would give a similar result to the TGA algorithm. However as the temperature is increased to T*=0.4 it is abundantly clear from the AVUS-HR calculations that the free energy minimal begins to appear at much higher q6 values compared to the low temperature case. On Figure 8 this region is a fairly high energy region and as such this minima is an effect of entropic stability instead of energetic stability. Traditional GA using a simple energy minimization, as shown in the supporting information, naturally struggle to predict these minima given the entropic stabilization. Even if one attempts to use a Boltzmann weight, this is effectively a futile exercise since the lowest energy structure will always have the best weight regardless of the choice of T. As such without either using free energy as the objective or performing a different style of sampling, traditional GA tends to miss out on the entropically stabilized structures (see supporting information - Figure S1 and Figure S2 for details). The two representative examples highlight the ability of the TGA to adequately predict structures and configurations by implicitly accounting for entropic contributions. Given the much lower compute cost of energy evaluations over explicit free energy calculations, the TGA offers the advantage of being computationally efficient for search at finite temperatures and pressures. This provides distinct advantage over the traditional GA which relies on individual populations with search criteria commonly based on evaluating ground state energies of the system. The algorithm is capable of handling multiple free energy minima simultaneously since even in the simple case there can be up to 15+ groups in the pool at a time. However, for larger problem spaces (e.g. proteins), the total number of structures retained in the pool will need to be increased in order to properly sort the various minima of interest. Note that having too few structures can result in a good minimum being pre13 ACS Paragon Plus Environment


maturely removed from the pool and hence one requires a large number of populations in such cases to ensure adequate sampling. The result of this is that it may take some time for that group to eventually return back into the pool. If there are a sufficient number of structures, then the various free energy minima have a very high retention time in TGA and it usually takes several crushing defeats for it to be removed completely. This will not normally happen if the group has a favorable score or a high rate of replacement. This procedure ensures that we adequately sample the various local minima. We note that while relaxation of generated structures using Molecular Dynamics or Monte Carlo is completely feasible, it was avoided in the context of this paper in favor of uniform generation schemes in order to show that the predictions were a result of the way the TGA sampling was done and not because one ”cheated” by introducing a generation bias using Boltzmann sampling. In future applications however the use of short Molecular Dynamics or Monte Carlo relaxations should assist the convergence rate by reducing the number of junk structures that are generated. It should be noted excessive use of local energy minimization schemes such as a Conjugated-Gradient approach can inadvertently bias the system toward energetically stable configurations and exclude entropically favored states. Therefore any relaxation schemes used in conjunction with TGA should have temperature related sampling incorporated into them. Since the TGA algorithm is still a GA based algorithm, many of the same sampling considerations need to be taken. For example, a poor choice of mating or mutation algorithm can fail to generate a diverse set of structures and as such may stifle the results of the algorithm. Finally, it is worth noting that our current study performs grouping using a single reaction coordinate that was representative of our chosen test systems. In theory any number of reaction coordinates can be used to perform the grouping. An ideal goal, however, would ultimately be to have the grouping be done on the fly without apriori knowledge of a reaction coordinate as this would improve the generality of the algorithm.


Page 14 of 29

Page 15 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Conclusion We have shown in this paper that the TGA algorithm is able to include some of the entropic contributions to structural optimization that are typically missing from traditional energy based GA. This allows one to optimize toward the free energy minima instead of simply the energetic minima. In addition we have presented a discussion which talks about extending the algorithm to more complicated problems of interest.

Supporting Information Available The results for the traditional GA simulations can be found within the supporting information. This includes information such as the workflow diagram, the resulting stable structures, and the methodology used.

Acknowledgement Funding: This work was supported by Argonne LDRD-2017-012-N0. Use of the Center for Nanoscale Materials and the resources of the Argonne Leadership Computing Facility was also supported by the U. S. Department of Energy (DOE), Office of Science, Office of Basic Science, under the contract no. DE-AC02-06CH11357. This research used resources of the National Energy Research Scientific Computing Center, a DOE office of science user facility supported by the Office of Science of the US Department of Energy under contract no. DE-AC02-05CH11231. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computational Facility, which is a DOE Office of Science user facility supported under contract no. DE-AC02-06CH11357.



References (1) Liu, T.-D.; Xu, L.-Y.; Shao, G.-F.; Tu, N.-N.; Tao, J.-P.; Wen, Y.-H. Structural Optimization of PtPdRh Trimetallic Nanoparticles Using Improved Genetic Algorithm. Journal of Alloys and Compounds 2016, 663, 466 – 473. (2) Zhao, X.; Nguyen, M. C.; Zhang, W. Y.; Wang, C. Z.; Kramer, M. J.; Sellmyer, D. J.; Li, X. Z.; Zhang, F.; Ke, L. Q.; Antropov, V. P. et al. Exploring the Structural Complexity of Intermetallic Compounds by an Adaptive Genetic Algorithm. Phys. Rev. Lett. 2014, 112, 045502. (3) Kaczmarowski, A.; Yang, S.; Szlufarska, I.; Morgan, D. Genetic Algorithm Optimization of Defect Clusters in Crystalline Materials. Computational Materials Science 2015, 98, 234 – 244. (4) Zhang, Y.-Y.; Gao, W.; Chen, S.; Xiang, H.; Gong, X.-G. Inverse design of materials by multi-objective differential evolution. Computational Materials Science 2015, 98, 51 – 55. (5) Jrgensen, M. S.; Groves, M. N.; Hammer, B. Combining Evolutionary Algorithms with Clustering toward Rational Global Structure Optimization at the Atomic Scale. Journal of Chemical Theory and Computation 2017, 13, 1486–1493, PMID: 28186745. (6) Vilhelmsen, L. B.; Hammer, B. A Genetic Algorithm for First Principles Global Structure Optimization of Supported Nano Structures. The Journal of Chemical Physics 2014, 141, 044711. (7) Wang, Z.; Zhou, X.-F.; Zhang, X.; Zhu, Q.; Dong, H.; Zhao, M.; Oganov, A. R. Phagraphene: A Low-Energy Graphene Allotrope Composed of 567 Carbon Rings with Distorted Dirac Cones. Nano Letters 2015, 15, 6182–6186, PMID: 26262429. 16 ACS Paragon Plus Environment

Page 16 of 29

Page 17 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(8) Graser, J.; Kauwe, S. K.; Sparks, T. D. Machine Learning and Energy Minimization Approaches for Crystal Structure Predictions: A Review and New Horizons. Chemistry of Materials 2018, 30, 3601–3612. (9) Ryan, K.; Lengyel, J.; Shatruk, M. Crystal Structure Prediction via Deep Learning. Journal of the American Chemical Society 2018, 140, 10158–10168, PMID: 29874459. (10) Ziletti, A.; Kumar, D.; Scheffler, M.; Ghiringhelli, L. M. Insightful Classification of Crystal Structures Using Deep Learning. Nature communications 2018, 9, 2775. (11) Deringer, V. L.; Proserpio, D. M.; Csnyi, G.; Pickard, C. J. Data-Driven Learning and Prediction of Inorganic Crystal Structures. Faraday Discuss. 2018, 211, 45–59. (12) Kinaci, A.; Narayanan, B.; Sen, F. G.; Davis, M. J.; Gray, S. K.; Sankaranarayanan, S. K.; Chan, M. K. Unraveling the Planar-Globular Transition in Gold Nanoclusters Through Evolutionary Search. Scientific reports 2016, 6, 34974. (13) Custdio, F. L.; Barbosa, H. J.; Dardenne, L. E. A multiple minima genetic algorithm for protein structure prediction. Applied Soft Computing 2014, 15, 88 – 99. (14) Chen, S.; Zheng, F.; Wu, S.; Zhu, Z. An Improved Genetic Algorithm for Crystal Structure Prediction. Current Applied Physics 2017, 17, 454–460. (15) Wu, S.; Ji, M.; Wang, C.-Z.; Nguyen, M. C.; Zhao, X.; Umemoto, K.; Wentzcovitch, R.; Ho, K.-M. An Adaptive Genetic Algorithm for Crystal Structure Prediction. Journal of Physics: Condensed Matter 2013, 26, 035402. (16) Oganov, A. R. Crystal Structure Prediction: Reflections on Present Status and Challenges. Faraday Discuss. 2018, 211, 643–660.



(17) Oganov, A. R.; Lyakhov, A. O.; Valle, M. How Evolutionary Crystal Structure Prediction Worksand Why. Accounts of Chemical Research 2011, 44, 227–237, PMID: 21361336. (18) Freitas, R.; Asta, M.; de Koning, M. Nonequilibrium Free-Energy Calculation of Solids using LAMMPS. Computational Materials Science 2016, 112, 333 – 341. (19) Frenkel, D.; Ladd, A. J. C. New Monte Carlo Method to Compute the Free Energy of Arbitrary Solids. Application to the FCC and HCP Phases of Hard Spheres. The Journal of Chemical Physics 1984, 81, 3188–3193. (20) Vega, C.; Sanz, E.; Abascal, J. L. F.; Noya, E. G. Determination of Phase Diagrams via Computer Simulation: Methodology and Applications to Water, Electrolytes and Proteins. Journal of Physics: Condensed Matter 2008, 20, 153101. (21) Lin, Y.; Li, J.-J.; Zhang, J.; Wan, M. A Tribal Ecosystem Inspired Algorithm (TEA) for Global Optimization. Proceedings of the 2014 Annual Conference on Genetic and Evolutionary Computation. 2014; pp 33–40. (22) Ma, B.; Xia, Y. A Tribe Competition-Based Genetic Algorithm for Feature Selection in Pattern Classification. Applied Soft Computing 2017, 58, 328 – 338. (23) Deaven, D. M.; Ho, K. M. Molecular Geometry Optimization with a Genetic Algorithm. Phys. Rev. Lett. 1995, 75, 288–291. (24) Daven, D.; Tit, N.; Morris, J.; Ho, K. Structural Optimization of Lennard-Jones Clusters by a Genetic Algorithm. Chemical Physics Letters 1996, 256, 195 – 200. (25) Stillinger, F. H. Rigorous Basis of the FrenkelBand Theory of Association Equilibrium. The Journal of Chemical Physics 1963, 38, 1486–1494.


Page 18 of 29

Page 19 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(26) Loeffler, T. D.; Henderson, D. E.; Chen, B. Vaporliquid Nucleation in Two Dimensions: On the Intriguing Sign Switch of the Errors of the Classical Nucleation Theory. The Journal of Chemical Physics 2012, 137, 194304. (27) Loeffler, T. D.; Chen, B. Surface Induced Nucleation of a Lennard-Jones System on an Implicit Surface at Sub-Freezing Temperatures: A Comparison with the Classical Nucleation Theory. The Journal of Chemical Physics 2013, 139, 234707. (28) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Bond-Orientational Order in Liquids and Glasses. Phys. Rev. B 1983, 28, 784–805. (29) Chen, B.; Siepmann, J. I.; Klein, M. L. Simulating VaporLiquid Nucleation of Water: A Combined Histogram-Reweighting and Aggregation-Volume-Bias Monte Carlo Investigation for Fixed-Charge and Polarizable Models. The Journal of Physical Chemistry A 2005, 109, 1137–1145, PMID: 16833423.



Figure 1: The two body energy given by Eq. 5 is plotted here as a function of radius. The system is effectively a two well harmonic potential with the global minima located a r=1˚ A ˚ and a second local minima at r=5A. The two wells differ in energy by approximately 5 meV with approximately a 30 meV barrier in the middle.


Page 20 of 29

Page 21 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: A flow chart which outlines the overall structure of the TGA algorithm is shown here. For the sake of comparison, we also provide the flowchart showing the evolutionary scheme in a standard GA algorithm based on evolution of individual populations.



Figure 3: Results from the Tribal GA simulations and the exact free energy calculations are presented here. From left to right the 10K, 100K, and 3000K two body systems are shown. The solid black curve represents the free energy obtained from Eq.6 The dashed red line shows the position of the group which had the best score at the end of the simulation while the dash-dot green lines illustrate the positions of the 2nd, 3rd, and 4th place groups. The dashed blue lines show the 5th-8th place positions.


Page 22 of 29

Page 23 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4: The group scores from the low temperature two body system are given as a function of the number of Generations are shown here. Each color on the line represents a different group contained within the structure pool. Group score is computed from Eq. 1



Figure 5: Results from the Tribal GA simulations and the exact free energy calculations are presented here. The solid black curve represents the free energy obtained from the Monte Carlo simulations. The dashed red line shows the position of the group which had the best score at the end of the simulation while the dash-dot green lines illustrate the positions of the 2nd, 3rd, and 4th place groups.


Page 24 of 29

Page 25 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6: Results from the T*=0.4 cluster simulation showing the total group score as a function of the number of generations are presented here. Each color on the line represents a different group contained within the structure pool. Due to the large number of groups the legend is omitted. The y-axis is plotted using a log-scale.



Figure 7: Representative clusters of the most stable groups at different temperature are presented here. At lower temperatures both TGA and standard GA can correctly predict the low energy icosahedral geometries. At high temperatures, tribal GA correctly predicts entropically stable clusters whereas standard GA still predicts energetically stable icosahedral cluster.


Page 26 of 29

Page 27 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8: The average energy for a cluster of 13 LJ particles is shown as a function of the q6 order parameter.




Page 28 of 29

Page 29 of 29

Tribal GA Best T* = 0.05

Entropic Effects

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


T* = 0.4

T* = 0.8

ACS Paragon Plus Environment

Traditional GA

Best

Tribal Competition Evolutionary Search as a Surrogate for Free

Recommend Documents