AUTOMATON: A Program That Combines a Probabilistic Cellular

Dec 13, 2018 - Vicentina, Iztapalapa , C.P. 09340 Mexico City , Mexico. J. Chem. Theory Comput. , 2019, 15 (2), pp 1463–1475. DOI: 10.1021/acs.jctc...
0 downloads 0 Views 5MB Size
Subscriber access provided by AUSTRALIAN NATIONAL UNIV

Structure Prediction

AUTOMATON: a program that combines a probabilistic cellular automata and a genetic algorithm for global minimum search of clusters and molecules Osvaldo Yañez, Rodrigo Báez-Grez, Diego Inostroza, Walter Alfonso Rabanal-León, Ricardo Pino-Rios, Jorge Garza, and William Tiznado J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00772 • Publication Date (Web): 13 Dec 2018 Downloaded from http://pubs.acs.org on December 13, 2018

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 57 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

Journal of Chemical Theory and Computation

AUTOMATON: a program that combines a probabilistic cellular automata and a genetic algorithm for global minimum search of clusters and molecules Osvaldo Yañez1,2, Rodrigo Báez-Grez1,2, Diego Inostroza2, Walter A. Rabanal-León2,*, Ricardo Pino-Rios3, Jorge Garza4, and W. Tiznado2,* 1 Doctorado en Fisicoquímica Molecular, Facultad de Ciencias Exactas, Universidad Andres Bello,

República 275 (2do piso), Santiago 8370146, Chile. 2

Departamento de Ciencias Químicas, Facultad de Ciencias Exactas, Universidad Andres Bello,

República 275 (3er piso), Santiago 8370146, Chile. 3

Laboratorio de Química Teórica, Facultad de Química y Biología, Universidad de Santiago de

Chile (USACH), Av. Bernardo O´Higgins 3363, Santiago, 9170022, Chile. 4

Departamento de Química, División de Ciencias Básicas e Ingeniería, Universidad Autónoma

Metropolitana-Iztapalapa, San Rafael Atlixco 186, Col. Vicentina, Iztapalapa, C.P. 09340. Mexico-City, Mexico. KEYWORDS: Probabilistic cellular automaton, Genetic algorithms, Potential energy surface, Nanostructures, Cluster, Molecules.

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 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

Page 2 of 57

ABSTRACT

A novel program for the search of global minimum structures of atomic clusters and molecules in the gas phase, AUTOMATON, is introduced in this work. This program involves: first, the generation of an initial population, using a simplified probabilistic cellular automaton method, which allows to easily control the adequate distribution of atoms in space; second, the fittest individuals are selected to evolve, through genetic operations (mating and mutations), until the best candidate for a global minimum surfaces. In addition, we propose a simple way to build the descendant structures by establishing a ranking of genes to be inherited. Thus, by means of a chemical formula checker procedure, genes are transferred to the offspring, ensuring that they always have the appropriate type and number of atoms. It is worth nothing that a fraction of the fittest group is subject to mutation operations. This program also includes algorithms to identify duplicate structures: one based on geometric similarity and another on the similar distribution of atomic charges. The effectiveness of the program was evaluated in a group of 45 molecules, considering organic and organometallic compounds (benzene, cyclopentadienyl anion, ferrocene), Zintl ion clusters [Sn9-m-nGemBin](4-n)- (n = 1-4 and m = 0-(9-n)), star-shaped clusters (Li7E5+, E=BH, C, Si, Ge) and a variety of boron-based clusters. The global minimum and the lowestenergy isomers reported in the literature were found for all the cases considered in this article. These results successfully prove AUTOMATON´s effectiveness on the identification of energetically preferred structures of a wide variety of chemical species.

ACS Paragon Plus Environment

2

Page 3 of 57 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

Journal of Chemical Theory and Computation

Introduction Atomic clusters1–4, having dimensions between atoms and nanoparticles, are fascinating systems that have increasingly attracted the attention of the world of science and technology5–8. Such interest is due to their feasibility of tuning cluster properties for a variety of possible applications in various areas, such as the environment, energy, and biology7,9–14. A very active research field, with surprising results, is gas-phase cluster chemistry. Brand new gaseous species are always being reported in the literature contributing to enhance our understanding of the chemical structures and bonding of such exotic species. Their structural assignment is usually performed by comparing a measured gas-phase property, i.e. photoelectron spectra, with its computationally simulated counterpart for a set of plausible candidates15–19. Therefore, a realistic geometric and electronic structure characterization requires some reliable prediction of lower energy structures by using precise theoretical calculations. Chemical intuition, i.e. making use of Lewis electron dot diagram20 and valence-shell electron-pair repulsion (VSEPR) theory,21–23 is usually ineffective for predicting energetically preferred cluster structures. These predictions become more complicated as the number of atoms conforming the cluster increases, thus encouraging the development of modern techniques to sample these complex potential energy surfaces (PESs) as impartially, automatically and exhaustively as possible24,25. The quality and capability of the methods proposed to find the global minimum (GM) on the PES of clusters and molecules have increased significantly through the years, including those based on the classical Monte Carlo (MC) annealing26–30, genetic algorithm (GA)31–41, particle swarm42–45, basin hopping (BH)46–50, Kick method51–55, among others56–58. Particularly, genetic algorithm (GA) based

ACS Paragon Plus Environment

3

Journal of Chemical Theory and Computation 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

Page 4 of 57

methods have proven to be highly effective in identifying the GM (for both atomic and molecular clusters) on the PES37,40,59,60. These methods are characterized by transferring structural information (genetic code) from the parents of an initial population to successive descendent ones, until the best individual is found (candidate to GM). Some modern implementations include builtin procedures to guarantee that the final population is composed exclusively of true local minima. This can be achieved by combining geometric optimizations and vibrational frequency analyses. Thus, individuals with imaginary frequencies could be forced to a near local minimum by moving their atoms (as an independent procedure) in the direction of the imaginary frequency vibrational mode61. It is important to mention that when we refer to GM, we mean the electronic-ground states A familiar example of this approach is the ab initio gradient embedded genetic algorithm (GEGA), which was originally proposed in 200562. This method combines GA operations with ab initio geometry optimizations and vibrational frequency analyses of candidate structures. If a saddle point is identified, then the atoms are moved following the normal modes of the corresponding imaginary frequencies until the structure becomes a local minimum. In 2010, GEGA was updated for finding the global minimum on the PES of mixed clusters formed by atoms and molecules (i.e. hydrogen atom solvated by water molecules)63. In 2014, Kanters and Donald, introduced the CLUSTER program40, which includes methodologies analogous to GEGA. When compared to GEGA, this CLUSTER shows significant improvement by avoiding the creation of large generations and providing flexibility for the user to choose, among several alternatives, the ones that best suit the problem being addressed. More recently, in 2017, our group evaluated the effect of using an evolved starting population, built by chemical criteria (i.e. using Fukui function condensed values as local descriptor of reactivity) in GEGA performance64. Interestingly, this strategy nearly reduced in half the convergence time required by the canonical GEGA. These

ACS Paragon Plus Environment

4

Page 5 of 57 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

Journal of Chemical Theory and Computation

findings have motivated us to implement a new hybrid methodology, which consists of generating a population of diverse and evolved individuals, by initially controlling key structural parameters (bond distances, angles between consecutive atoms) to gain advantage in subsequent evolutionary processes. In view of all this, we introduce a structural search algorithm implemented in the new AUTOMATON program. The program consists of two main procedures: in the first one, a discrete population is generated combining rules of a simplified probabilistic cellular automaton method65 with geometrical optimizations (to the nearest stationary point) using an ab initio method. In the second one, this population is evolved through genetic operations followed by geometrical optimizations (to the nearest stationary point) using an ab initio method. Additionally, AUTOMATON includes a structure-recognition routine, which is used in different stages of the search process to identify and eliminate duplicates. Finally, AUTOMATON effectiveness has been evaluated by searching the lowest energy isomers in a diverse group of systems, including organic and organometallic molecules as well as experimentally identified inorganic three-dimensional and flat clusters. Finally, although AUTOMATON gives reliable structures, taking advantage of the use of quantum chemistry energy as cost function, this approximation limits the possibility of evaluating medium- and large-sized systems because of its highest computational cost.

Methodology Generation of the initial population a.

Cell selection procedure. For the construction of candidate structures, the AUTOMATON

code performs the procedure depicted in Scheme 1 (See also Scheme S1 in the SI), described as follows:

ACS Paragon Plus Environment

5

Journal of Chemical Theory and Computation 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



Page 6 of 57

A cubic space centered at the origin of coordinates is divided into regular cubes with sides of 0.3 Å. This value was chosen to be close to the smallest covalent radius (for hydrogen = 0.31 Å) found in the Periodic Table. AUTOMATON uses the new set of covalent atomic radii (deduced from crystallographic data)66. It is established that the side of the cube must be as close as possible (but never exceed) to twice the sum of the covalent radii of all the atoms that make up the system under study62. AUTOMATOM picks an "i" cell randomly, and its center defines the position of the first atom of the system.



The next cell “j” can be any of the cells intersected by the surface of a sphere of radius rij (covalent radius of atom in “i" + covalent radius of atom to be paced in “j”) centered in the “i” cell. In a two-dimensional representation, such as the one shown in Scheme 1, this set of possibilities corresponds to a Moore-like neighborhood67.



Once cell “j” is occupied, its Moore-like neighborhood is defined, but this time by a sphere of radius rjk (sum of covalent radii of atoms j and k) centered in cell “j”. In order to place the third and subsequent atoms, and to prevent using cells between two adjacent atoms, a restriction is established: the cells located in the overlap region between the two contiguous neighborhoods are not to be used. Once a third atom is placed in cell “k”, the process is repeated until the total number of atoms of the system is placed in viable cells.

ACS Paragon Plus Environment

6

Page 7 of 57 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

Journal of Chemical Theory and Computation

Scheme 1. Representation of the cell selection procedure for the construction of a four-atom cluster. For clarity, this scheme was reduced to a two-dimensional arranged cell, in AUTOMATON this process takes place in one, two and three dimensions.

b.

Automaton-population selection. A brief analysis of the influence of the initial

population size on the AUTOMATON effectiveness to identify the GM has been made and will be discussed at the beginning of the results section. Based on this analysis, the AUTOMATONbound suggests that the population size is equal or superior to 5N, where N is the number of atoms conforming the system. The 5N (or more) candidate structures, obtained in the step “Cell selection procedure”, are relaxed to the nearest local minimum using a gradient method. In this work, the local minimization of the cluster energy was performed using the "Berny algorithm" included in Gaussian 09 program64. The results discussed in this paper will not consider frequency analysis at this stage. Nevertheless, this option could be activated by the user, in such a way that each transition state structures are relaxed to their nearest local-minimum, similarly to GEGA approach.

ACS Paragon Plus Environment

7

Journal of Chemical Theory and Computation 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

c.

Page 8 of 57

Defining a fitness function. A fitness function is defined to select the best individuals in

a population. The fitness of the lowest and highest energy isomer in AUTOMATON is 1.00 and 0.05, respectively (this is because we used the exponential fitness function). Dynamic scaling is achieved by using a normalized value of the energy, , in fitness calculations: 𝐸𝑖 ― 𝐸𝑚𝑖𝑛

𝜌𝑖 = 𝐸𝑚𝑎𝑥 ― 𝐸𝑚𝑖𝑛

(1)

Where Ei, Emin and Emax correspond to the isomer being evaluated, to the lowest and to the highest energy isomers within the analyzed population, respectively. To obtain these results, the following exponential fitness function was used: 𝑓𝑖 = exp ( ― 𝑎𝜌𝑖) Considering a=361, all isomers with fi ≥ 0.5 (referred to as alpha individuals from now on) are selected to be subjected to the process of mating. In parallel, N of the alpha individuals (randomly selected) are selected to be mutated (vide infra). d.

Mating Procedure. AUTOMATON considers a variant of the “cut and splice” crossover

operator of Deaven and Ho68. To create a child, two parents are randomly chosen from the alpha individuals (vide infra).

ACS Paragon Plus Environment

8

Page 9 of 57 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

Journal of Chemical Theory and Computation

Scheme 2. Representation of the mating (crossover) operation.

ACS Paragon Plus Environment

9

Journal of Chemical Theory and Computation 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

Page 10 of 57

Each parent is centered (at the origin of coordinates) and subsequently randomly rotated40. Then, both parents are cut by a randomly oriented plane, , passing through the origin of coordinates (See Scheme 2). Next, the atomic distance from the cutting plane above (parent 1) and below (parent 2) are considered as a probability criterion to be inherited by the child (the greatest the atomic distance, the better chance to be transferred). Finally, the atoms are selected according to the formula verifier in Scheme 2. e.

Mutation Procedure. In order to introduce new genetic material and increase population

diversity, AUTOMATON considers three mutation schemes (See Scheme 3): i.

Atom displacement. N randomly selected isomers (from the alpha individuals) are mutated as follows: the atomic coordinates of 40% of their atoms are replaced by randomly generated values, varying each position coordinate to a value between -rα and +rα, where rα is the atomic radius of the atom to be displaced.

ii.

Atom permutation. Consisting of swapping positions between one or more pairs of atoms of different elements, without disturbing the structure of the system. This kind of mutation is used for hetero-elemental clusters, such as ionic and bimetallic clusters.

iii.

Cluster replacement. An entire cluster is replaced by a new one, which is generated in an identical way to the one used for the generation of the initial population. AUTOMATON generates N/2 linear, N/2 two-dimensional, and N/2 three-dimensional structures.

ACS Paragon Plus Environment

10

Page 11 of 57 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

Journal of Chemical Theory and Computation

Scheme 3. Representation of the mutation operation. All structures obtained by mutations are optimized to the nearest local minima. Then, these local minima, together with those obtained in other process steps (from its start) are ordered according to their energy and, the best 5N ones (lower energy) are chosen to begin a new exploration cycle that includes breeding and mutations. All low-lying isomers are detected and stored throughout the execution, and then reported to the user at the end of each cycle. AUTOMATON converges into a solution when the lowest energy isomer persists for 9 iterations (cycles). It is important to remark that there is no true convergence criterion for a non-deterministic global optimization procedure; consequently, the consideration of 9 cycles is merely an ad-hoc guess, the validity of which indeed also depends on cluster size and cluster type. This is evidenced by exploring the potential energy surface of Lennard-Jones clusters, which will be discussed in the final part of the results. We would also like to mention that these systems have not been considered in the initial version of our work. They were included as a suggestion by one of the reviewers and allowed us to detect that the effectiveness of our algorithm, with the characteristics presented here, is severely affected when the study system is large. f.

Similarity check algorithm. In AUTOMATON, after each process, where a set of

optimized structures is gathered, a similarity test is carried out to identify and eliminate structures corresponding to the same local minimum.

ACS Paragon Plus Environment

11

Journal of Chemical Theory and Computation 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

Page 12 of 57

After analyzing the effectiveness of some of the most popular similarity check algorithms (See Table S1 in the SI), the couple that proved to be the most successful one is implemented in our code: the first one was introduced by Rogan et al.69 as a modification to the proposal of Grigoryan and Springborg's (this algorithm is used by default)70–72 and it is based on the comparison of distances (Eq. 2): ℵ

d (α, β) =

[

N(N - 1)

(

dαn 2 ∑ 2 N(N - 1) n = 1 Dαave

dβn

- Dβ

2 1/2

)]

ave

(2)

Where 𝑑ℵ(𝛼, 𝛽) is a non-dimensional quantity, N is the number of the atoms in the system, 𝑑𝛼𝑛( 𝑑𝛽𝑛) and 𝐷𝛼𝑎𝑣𝑒(𝐷𝛽𝑎𝑣𝑒) are the ordered (from smallest to largest) interatomic distances and the average bond length between the atoms of system 𝛼(𝛽), respectively. The second one (Eq.3) is introduced in this work, replacing the distances, in Eq. 2, by atomic charges (e.g. Mulliken atomic charges):



q (α, β) =

[

N(N - 1)

(

qαn 2 ∑ 2 N(N - 1) n = 1 Qαave

qβn

- Qβ

2 1/2

)]

ave

(3)

Where 𝑞ℵ(𝛼, 𝛽) is a non-dimensional quantity, N is the number of the atoms in the system, 𝑞𝛼(𝛽) 𝑛 𝛼(𝛽) 𝛼 𝛽 = 𝑞𝛼(𝛽) 𝑖(𝑘) 𝑞𝑗(𝑙) (𝑖(𝑘) ≠ 𝑗(𝑙)) and 𝑄𝑎𝑣𝑒(𝑄𝑎𝑣𝑒), are the ordered (from smallest to largest) internuclear

charges product and the average of the charges product of system 𝛼(𝛽), respectively. Finally, the AUTOMATON program flowchart, depicted in Scheme 4, summarizes the sequence of the above described processes. It should be pointed out that frequency checking is included in the flow chart because it is available in the program. However, it is not used in the current work. It is also worth noting that local minimizations are performed in parallel; however, this collectand-compare stage acts as a serial bottleneck, severely degrading overall performance and scalability of the algorithm. A different way to introduce parallelization is changing the GA from

ACS Paragon Plus Environment

12

Page 13 of 57 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

Journal of Chemical Theory and Computation

the generational paradigm to the steady-state paradigm, which has shown to be clearly superior73,74. This will be implemented in future AUTOMATON updates

ACS Paragon Plus Environment

13

Journal of Chemical Theory and Computation 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

Page 14 of 57

Scheme 4. The AUTOMATON program flowchart.

Computational Details In the AUTOMATON processes, geometry optimizations and vibrational frequency calculations were carried out using the Gaussian 09 program75, employing the PBE0 hybrid functional76 and the Stuttgart/Dresden effective core potential (ECP) with its corresponding double- basis set77–81. Then, the lowest energy isomers delivered by AUTOMATON (in a range of 30 kcal.mol-1) were reoptimized at PBE0/Def2-TZVP82 level and, after verifying the absence of imaginary frequencies, the best minimum was identified as GM. Relativistic DFT calculations employing the zeroth-order regular approximation (ZORA)83,84 to the Dirac equation, including spin−orbit (SO) coupling and scalar relativistic effects, were also applied in specific cases. These calculations were performed using PBE (GGA) and PBE0 (hybrid) xc-functionals with slater type triple- plus two polarization functions (TZ2P) basis set85, as implemented in the ADF 2012.01 suit of programs86.

Results and Discussion To assess AUTOMATON’s performance and effectiveness, a set of 45 systems (including clusters, organic and organometallic molecules), already discussed in the literature, were considered. These molecules can be classified as follows: Zintl ion clusters [Sn9-m-nGemBin](4-n)- (n = 1-4 and m = 0-(9-n)); organic and organometallic molecules (benzene, cyclopentadienyl anion,

ACS Paragon Plus Environment

14

Page 15 of 57 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

Journal of Chemical Theory and Computation

and ferrocene); star-shaped clusters (Li7E5+, E=BH, C, Si, Ge); ionic bare boron clusters (B7-, B9-, B13+, B13- and B19-) and boron clusters doped with a single atom (CB62-, Co@B8- and Ru@B9-). With the purpose of defining an adequate population size for our explorations, we evaluated the of AUTOMATON’s effectiveness to identify the GM of the Ge5Bi4 system, using two initial populations. Both explorations, using 5N (45 individuals) and 3N (27 individuals), allowed us to find the GM. The one with the most significant population identified only one additional isomer in the range of 10 kcal.mol-1 above the lowest energy one (See Table S2 in SI). An initial population of 5N seems to be sufficient and adequate. Consequently, it is used in all searches reported in this work.

Zintl ion clusters. Fourteen of the thirty [Sn9-m-nGemBin](4-n)- (n = 1-4 and m = 0-(9-n)) possible combinations have been detected by electrospray mass spectrometry87, and candidates for the most stable structures of each stoichiometry have been previously reported in the literature88. They are, therefore, ideal candidates to be used as benchmark systems to test the performance of AUTOMATON. For each stoichiometry, a 5N population (45 individuals) was generated using the “cell selection procedure”. After optimizing these structures, duplicates were identified and eliminated, repeating the process until obtaining 45 different local minima. Then, the structures were re-optimized (see computational details) and sorted out at an increasing order of energy (See Figures S1-S22 in the SI). Table 1, summarizes the comparison between the most stable isomers found in this work with those reported by Muñoz-Castro and Sevov at PBE/TZ2P level (these calculations included scalar and spin-orbit relativistic effects using the zero-order regular approximation (ZORA) method)88.

ACS Paragon Plus Environment

15

Journal of Chemical Theory and Computation 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

Page 16 of 57

A check (), cross () or exclamation (!) mark means that the identified lowest energy isomer is the same, higher or lower in energy than the one of the reference, respectively. Interestingly, the same GM was identified in 20 cases and in two instances, a better candidate for GM was found ([Sn3Ge5Bi]3- and [Sn5GeBi3]-). However, at this stage, AUTOMATON could not find the GM for 8 reported cases. Table 1. Summary of the results obtained by the “cell selection procedure” (into the AUTOMATON code) used in the GM search of the Zintl ions of formula [Sn9-m-nGemBin](4-n)- for n = 1-4 and m = 0-(9-n). GMa

System

GMa

[Sn5Ge2Bi2] 2[Sn4Ge3Bi2] 2[Sn3Ge4Bi2] 2[Sn2Ge5Bi2] 2[SnGe6Bi2] 2[Ge7Bi2] 2[Sn8Bi] 3-

[Sn5Bi4] [Sn4GeBi4] [Sn3Ge2Bi4] [Sn2Ge3Bi4] [SnGe4Bi4] [Ge5Bi4] [Sn6Bi3][Sn5GeBi3]-

System

b

[Sn7GeBi] 3-

[Sn4Ge2Bi3][Sn3Ge3Bi3][Sn2Ge4Bi3]-

[Sn6Ge2Bi] 3[Sn5Ge3Bi] 3[Sn4Ge4Bi] 3-

[SnGe5Bi3]-

[Sn3Ge5Bi] 3-

[Ge6Bi3][Sn7Bi2]2[Sn6GeBi2]2-

[Sn2Ge6Bi] 3[SnGe7Bi] 3[Ge8Bi]3-

b

a

Comparison with the structures reported as global minima in reference 88.

b

A lower energy structure was found.

In order to compare the energies of the candidates with the literature, PBE- and PBE0-SOZORA/TZ2P level calculations were performed for four selected cases: [Sn3Ge2Bi4], [Sn5GeBi3]ACS Paragon Plus Environment

16

Page 17 of 57 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

Journal of Chemical Theory and Computation

, [Sn5Ge2Bi2]2- and [Sn3Ge5Bi]3- (one for each charge, including two cases where a best GM candidate was found, see Figure 1). These retained the lowest energy structures in a range of 10 kcal.mol-1, with small variations of the relative energy values (see Tables S3-S6, in SI). For instance, at PBE0/Def2-TZVP level, [Sn3Ge5Bi]3- and [Sn5GeBi3]- the lowest energy isomers are 0.37 and 0.21 kcal.mol-1 more stable than their second energy structures (reported as GMs by Muñoz-Castro et al.88). These small differences are shortened at PBE- and PBE0-SOZORA/TZ2P level. Then, it is more appropriate to consider these structures as energetically competitive.

Figure 1. Energy comparison between the GM found by “cell selection procedure” (into the AUTOMATON code) and those reported in reference 88.

ACS Paragon Plus Environment

17

Journal of Chemical Theory and Computation 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

Page 18 of 57

It is worth to note that AUTOMATON identified more isomers close in energy to the GM. Figure 2, shows a comparison of the lowest energy isomers identified, within a range of 10 kcal.mol-1, by bar chart. The red bars show isomers found solely in this work; the yellow bars only by reference 88; and, the orange bars by both, reference 88 and us. However, it is important to clarify that Muñoz-Castro et al. did not define 10 kcal.mol-1 as a sampling window, instead, their highest energy structures reported were 7.5 kcal.mol-1 above GM. The populations of the eight cases where the GM reported in the literature was not found (blue formulas in Figure 2) were submitted to the two additional procedures considered by AUTOMATON program, mating and mutation. Then, after further optimization (at PBE0/Def2TZVP level) of the delivered candidate structures, all GMs were found.

ACS Paragon Plus Environment

18

Page 19 of 57 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

Journal of Chemical Theory and Computation

Figure 2. Bar representation of the number of isomers under 10 kcal.mol-1 found just on this work by the “cell selection procedure” (red bar), just in reference 88 (yellow bar), and the isomers found by both works (orange bar). These results prove the advantage of involving additional evolutionary procedures on a partially evolved population. Besides, these refinements also enrich the variety of lower energy structures. For example, for the representative case [Sn3Ge3Bi3]-, fifteen local minima under 4.7 kcal.mol-1

ACS Paragon Plus Environment

19

Journal of Chemical Theory and Computation 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

Page 20 of 57

(depicted in Figure 3) were identified, out of which nine persist from the initial population, and only four have been previously reported in the literature (GM, 1, 2 and 3 structures). Further information can be found on the SI Figures S23-S30. Subsequently, the systematic study of Zintl ion clusters has allowed us to define the best strategy to use AUTOMATON in the search for the GM of the other molecular systems proposed in this work. For all the cases described below, the GM search was performed by sequentially executing all the AUTOMATON procedures, with an initial population of 5N.

Figure 3. The GM and the lowest energy isomers (within 5 kcal.mol-1) for the [Sn3Ge3Bi3]- anionic Zintl cluster identified after running all AUTOMATON processes. It is important to note that all considered Zintl clusters are isoelectronic with [Sn9]4- and [Pb9]4-. Consequently, these cage-like clusters should follow the Wade-Mingos rules89–91, adopting the

ACS Paragon Plus Environment

20

Page 21 of 57 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

Journal of Chemical Theory and Computation

nido-capped square antiprism (CSA) structure –which would have ideal C4v symmetry for E94-. However, it is known that the closo-deltahedral (triangulated) structure of the tricapped trigonal prism (TTP: the D3h symmetry ground state structure of [Sn9]2-) is also preferred by [Sn9]4- and [Pb9]4-anions, depending on the surrounding counterion92,93. TTP structure was also predicted (based on combined GA-DFT searching with electrostatic deflection measurements) as the GM for neutral Sn6Bi3 (isoelectronic with [Sn9]3-) in the work reported by Heiles et al.94, having full D3h symmetry, (Bi atoms preferentially cap the square faces of the trigonal prism). All the Zintl systems studied in this work [Sn9-m-nGemBin](4-n)- have GM geometries between TTP and CSA ideal structures, these have been carefully classified by Muñoz-Castro and Sevov88. Organic and organometallic molecules: C6H6, C5H5-, Fe(C5H5)2 The global minima of the compounds with the formula: C6H6 (benzene), C5H5- (cyclopentadienyl anion) and (C5H5)2Fe, (ferrocene), were chosen as search targets, due to their unique structural and electronic characteristics. In particular, they are highly symmetric structures with their atoms joined by covalent bonds (localized, delocalized, polar and nonpolar ones) and ionic bonds (in the case of ferrocene). Besides, they are essential compounds in organic and organometallic chemistry with several studies confirming their minimum energy structures unequivocally. In the three cases studied, their global minimum, as well as the second isomer at 34.1, 45.2 and 57.7 kcal.mol-1 for C6H6, C5H5- and Fe(C5H5)2, respectively (see SI Figures S31-S33) were identified (at the PBE0/def2-TZVP level). Figure 4, shows the three structures identified by the search with AUTOMATON, demonstrating its effectiveness. Also, it is important to mention that the second best isomer, for C6H6 stoichiometry, coincides with the one reported in reference 95.

ACS Paragon Plus Environment

21

Journal of Chemical Theory and Computation 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

Page 22 of 57

Figure 4. GM structures of organic and organometallic molecules found by AUTOMATON. Star-shaped clusters: Li7(BH)5+, Li7C5+, Li7Si5+, Li7Ge5+ Recently, some of us showed that Li7E5+ clusters (E = BH, C, Si, Ge) adopt star shapes in their minimum energy conformations (singlet states)96–100. These studies involved the use of different state-of-the-art methods for scanning the clusters’ PESs. Besides, all these species appear to be aromatic according to magnetic criteria. This valuable reference information makes these clusters suitable study cases for assessing AUTOMATON’s performance. In all the cases studied, AUTOMATON identified the star-shaped as the energetically preferred structures (Shown in Figure 5). The AUTOMATON search also allowed to identify some low energy local minima not yet reported in the literature. In the case of Li7(BH)5+, four additional local minima (∆E = 17.2, 21.0, 21.6 and 22.8 kcal.mol-1) were found between the GM and the isomer previously reported98 as the second one in energy (∆E = 25.4 kcal.mol-1), which was also found by AUTOMATON (See Figure S34 in the SI). For Li7C5+, the second best isomer corresponds to the one previously reported96 (See Figure S35 in the SI). However, the relative energies differ depending on the methods employed for

ACS Paragon Plus Environment

22

Page 23 of 57 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

Journal of Chemical Theory and Computation

calculations (∆E = 18.1 kcal.mol-1 at the PBE0/def2-TZVPP level, and ∆E = 10.5 kcal.mol-1 at CCSD(T)/def2-TZVPP//B3LYP/def2-TZVPP level). In the case of Li7Si5+, two additional isomers to those reported by Tiznado et al.97 were found considering the same energy range (above GM), these local minima were not found by Clemence et al. either59; but, they have been recently identified by Vásquez-Espinal et al.100 (See Figure S36 in the SI). The Li7Ge5+ lowest energy isomer was reported as a distorted Cs star-like structure at the PBE/TZ2P level99; the current calculations (at the PBE0/def2-TZVPP level) support this finding. However, Vásquez-Espinal et al. showed that D5h structure is energetically preferred at the CCSD(T)/def2-TZVP//PBE0-D3/def2-TZVP level 100. For this reason, we report the D5h structure in Figure 5. However, it is important to point out that this structure has two doubly degenerated imaginary frequencies (31i cm-1 and 19i cm-1), at the PBE0/def2-TZVP level. Additionally, all other low energy isomers (within an energy range of 10 kcal.mol-1 above GM) reported by Vásquez-Espinal et al. were also found by AUTOMATON (See Figure S37 in the SI).

ACS Paragon Plus Environment

23

Journal of Chemical Theory and Computation 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

Page 24 of 57

Figure 5. Star-shaped GMs of Li7(BH)5+, Li7C5+, Li7Si5+, and Li7Ge5+ clusters found by AUTOMATON.

ACS Paragon Plus Environment

24

Page 25 of 57 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

Journal of Chemical Theory and Computation

Boron clusters Boron clusters in gas phase have been extensively studied, by both theoretical and experimental methods in recent years2,101. This research has made it possible to identify different species, which promote the development of new theoretical tools102–104, increase knowledge of the chemistry of these compounds2, and establish several electronic and structural similarities between some clusters and the nanostructures embedded in solid phase materials105,106. Some recent examples of these developments are: the identification of the molecular Wankel motors (due to the high fluxionality of internal fragments that undergo internal rotation almost without energy barrier)107– 111

and the extension of concepts applied in organic chemistry, like aromaticity, to rationalize

boron clusters' structures and their stability112. For all these reasons, we have chosen a representative group of boron clusters to search their global minima, starting only from the data of their formula, employing the AUTOMATON program. The results are shown and discussed below.

Homo-atomic boron clusters In this series, B7-, B9-, B13+, B13- and B19- ionic species were evaluated. For anions, their GM geometries are well established in the literature by comparing their experimental photoelectron spectra with their theoretical counterpart. While for the case of B13+, there are numerous theoretical works, with different levels of theory and rigor, which provide certainty about the structure of the GM. After running AUTOMATON in the search for lower energy structures (at the PBE0/SDD level) and refining the geometries and energies (at the PBE0/def2-TZVP level), we found the previously reported global minimums for all the systems of this series (See Figure 6). It is important to

ACS Paragon Plus Environment

25

Journal of Chemical Theory and Computation 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

Page 26 of 57

highlight that the most significant lower energy isomers, due to their proximity to the global minimum, were also identified (See Figures S38-42 in the SI). In the case of B7-, the isomer C6v (triplet) was identified as GM, as well as two neighboring isomers at 5.4 (C2v, singlet) and 11.0 (C2v, triplet) kcal.mol-1 (at the PBE0/def2-TZVP level).

Figure 6. GM structures of small- and medium-size boron clusters successfully identified by AUTOMATON.

These results are in agreement with those reported in the literature. For B9-, the global minimum (D8h, singlet) and other three nearby energy isomers were identified at 3.7 (Cs, singlet), 6.0 (Cs,

ACS Paragon Plus Environment

26

Page 27 of 57 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

Journal of Chemical Theory and Computation

singlet) and 9.1 kcal.mol-1 (Cs, singlet), respectively (at the PBE0/def2-TZVP level). These energies change to 22.9, 22.8, 20.8 kcal.mol-1 at the B3LYP/6-311+G* level (the same used in reference

113),

thus demonstrating the dependence of the relative energies on the level of theory

used in the calculations101 (See Figure S39 in the SI). For B13+, we have identified the C2v structure, predicted by Ricca and Bauschlicher114, as a global minimum and, as the second best minimum, the structure Cs that Schleyer and co-workers115 also predicted (See Figures S40-S41 in the SI). For B19-, the global minimum, and a C2v isomer very close in energy, were identified (E= 1.8 kcal.mol-1 to CCSD(T)/def2-TZVP//PBE0/def2-TZVP), which is in total agreement with that reported in the literature116. For this case, it was not possible to find the third best reported minimum (∆E = 8.0 kcal.mol-1, at the B3LYP/6-311+G* level). However, we identified new local minima in the energy range considered in reference

116

(See Figure S42 in the SI).

CB62-, Co@B8- and Ru@B9- Boron based clusters The clusters CB62-, Co@B8- and Ru@B9- are systems whose interest focused in the types of bonds present in their minimum energy structures. In all cases, it has been proposed that the flat structures, with the heteroatom at the center and surrounded by a boron-ring, is the energetically favored one. However, the isomer containing the flat hexacoordinated carbon structure does not correspond to the global minimum. This fact is fascinating since this structure has been promoted as a "divining molecule" and stands out on the cover of Chemical & Engineering News117. Subsequent theoretical works showed that the global minimum, and other minimum energy structures, contain the carbon atom in the peripheral edge of planar structures118. Interestingly, the AUTOMATON results are in agreement with these

ACS Paragon Plus Environment

27

Journal of Chemical Theory and Computation 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

Page 28 of 57

predictions, both in the identification of the global minimum and of the other representative lowest energy isomers. It is also important to note, that AUTOMATON identified three new isomers between GM and the D6h isomer (See Figures. S43 in the SI). The preference of C to a lowercoordinate site, rather than sitting in the middle of a Bn template, is consistent with the known relative isomer stability of carbaboranes119.

Figure 7. GM and the seventh lowest energy isomer for the light heteroatomic CB62- cluster (top); GM structures of the Co@B8- and Ru@B9- clusters (bottom). All these structures were found by AUTOMATON.

ACS Paragon Plus Environment

28

Page 29 of 57 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

Journal of Chemical Theory and Computation

The clusters Co@B8- and Ru@B9- are examples of structures where the coordination number of an atom within a stable chemical system is taken to the limit, as highlighted by Heine and Merino some years ago120. The geometries of the global minimum of these species have been unequivocally assigned by comparing experimental photoelectronic spectra with those obtained for the best candidates employing theoretical calculations121. Finally, AUTOMATON also successfully found the isomers that contain the octa- and nona-coordinate plane metal as the global minima (reported in the literature121), as shown in Figure 7 (the most relevant isomers are shown in Figures S44-S45, in the SI).

Lennard–Jones clusters One of the reviewers asked us to explore the PES of Lennard-Jones clusters. Suggesting clusters with very challenging search spaces, for example N = 75, 76, 77, 98, 102, 103, 104, etc. To address the problem posed, we had to make some modifications to our methodology. Firstly, we changed the program for local optimizations by the LAMMPS software122, which includes the potentials for the Lennard-Jones (LJ) clusters. Besides, the generation of individuals in space 1D and 2D was suppressed, taking into account that the LJ clusters prefer 3D structures. With these modifications, some small clusters were explored (N = 10, 12, 13, 20, 30, 38). For all these cases, the minimum energy structures in agreement with those reported in the literature123–128 were identified. It is important to note that the cluster with N = 38, has been reported as a problematic case, with two competitive structures of minimum energy34. As can be seen in Table S7 (in the SI) AUTOMATON identifies both lowest-energy structures. However, despite these initial positive and encouraging results, the AUTOMATON lowest energy structure for N = 75 is very above in energy compared to the reported GM. With the aim

ACS Paragon Plus Environment

29

Journal of Chemical Theory and Computation 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

Page 30 of 57

to improve these results, we evaluated some additional modification in the AUTOMATON procedure: a) we changed the initial population to 10N, b) the percentage of atoms to be moved in the mutation operation was reduced to 10%, because the standard process (30%) provides very high energy individuals for this system, c) the process of identifying duplicates was suppressed (because given the size of the cluster it took a long time to compute it) and we finally increased the cycles from 9 to 50. Even with all these modifications, the minimum structure found is still above the GM reported in the literature (-391.09 vs. -397.49 in the GM128). Given these results, we can conclude that our proposal has severe limitations when dealing with problematic systems of more than 50 atoms (as the LJ clusters with complexes PESs), which we hope to overcome in the short and medium term.

Conclusions This paper describes the implementation of AUTOMATON, a novel program to explore potential energy surface. The effectiveness of the program is demonstrated by successfully identifying the global minima of various atomic clusters and molecules of different nature in the gas phase, assigned according to the data available in the literature. One major novelty of this program is that it uses an easy way to create individuals with appropriate geometries from the beginning of the search process (i.e., favoring the formation of bonds and avoiding both very short or long distances between nuclei), based on a probabilistic cellular automaton method. We assessed the effectiveness of this process in the Zintl type benchmark clusters, where a population of 5N individuals was generated for each stoichiometry (where N is the number of atoms that make up the cluster). Surprisingly, after optimizing the candidate structures, the global minima were identified for 73% of the combinations. We used the

ACS Paragon Plus Environment

30

Page 31 of 57 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

Journal of Chemical Theory and Computation

clusters, for which global minimum was not identified (27%), as testing systems for the subsequent set of the program’s processes. This consists in applying genetic operations to the individuals to make them evolve until the global minimum is reached. At the end of these processes, the program successfully found all the global minima. The geometries, of all the other cases studied, were explored using the complete sequence of processes, thus allowing us to find the global minimum and all the lowest energy isomers. Besides, AUTOMATON includes two algorithms to identify duplicate structures throughout all search processes: one based on comparing atomic distances and another on comparing atomic charge distribution. The use of these algorithms prevents computational expenses in repeated structures and allows to maintain the diversity among the individuals of each generation. The variety of molecular systems successfully analyzed, as well as the possibilities of updating the program continuously improving it, depicts AUTOMATON as a valuable handy computational tool for molecular design.

ASSOCIATED CONTENT Supporting Information. The source code of the AUTOMATON program, as well as the xyz coordinate files for all the studied systems in this paper, can be found and downloaded from

the

following

url:

https://github.com/HumanOsv/Automaton.

Additionally,

supplementary material containing: the flowchart for the molecular candidate construction using the “cell selection procedure”; the GM and lowest in energy structures for the Zintl ion clusters (Figure S1-S30), the organic and organometallic molecules (Figure S31-S33), the star-shaped clusters (Figure S34-S37), the boron clusters (Figure S38-S45)

ACS Paragon Plus Environment

31

Journal of Chemical Theory and Computation 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

Page 32 of 57

are also reported. Furthermore, a similarity algorithm test (Table S1), the evaluation of the effectiveness of the AUTOMATON code varying the initial population (Tables S2), and the benchmark results for Zintl ions clusters (Table S3-S6) are also reported.

AUTHOR INFORMATION Corresponding Author * William Tiznado E-mail: [email protected] * Walter A. Rabanal-León E-mail: [email protected]

Funding Sources FONDECYT – Regular grant Nº 1181165 FONDECYT – Postdoctoral fellowship Nº 3160388 CONICYT – Doctoral fellowship Nº 21140667

ACKNOWLEDGMENT

ACS Paragon Plus Environment

32

Page 33 of 57 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

Journal of Chemical Theory and Computation

The Automaton code was initially developed at Universidad Andres Bello (Chile) as part of the PhD. dissertation (2018) of O.Y. under the supervision of W.T and J.G. It was financially supported by the FONDECYT grant Nº 1181165 and the CONICYT (PCHA/Doctorado

Nacional/

2014



21140667)

doctoral

fellowship.

W.A.R.L.

acknowledges the financial support from CONICYT for his postdoctoral Project FONDECYT/Postdoctorado - 2016 Nº 3160388.

REFERENCES (1)

Jortner, J. Clusters as a Key to the Understanding of Properties as a Function of Size and Dimensionality. In Physics and Chemistry of Finite Systems: From Clusters to Crystals; Jena, P., Khanna, S. N., Rao, B. K., Eds.; Springer Netherlands: Dordrecht, 1992; pp 1–17.

(2)

Alexandrova, A. N.; Boldyrev, A. I.; Zhai, H. J.; Wang, L. S. All-Boron Aromatic Clusters as Potential New Inorganic Ligands and Building Blocks in Chemistry. Coord. Chem. Rev. 2006, 250, 2811–2866.

(3)

Malinowski, N.; Schaber, H.; Bergmann, T.; Martin, T. P. Electronic Shell Structure in NaO Clusters. Solid State Commun. 1989, 69, 733–735.

(4)

Wade, K. Structural and Bonding Patterns in Cluster Chemistry; Radiochemistry, H. J. E. and A. G. S. B. T.-A. in I. C. and, Ed.; Academic Press, 1976; Vol. Volume 18, pp 1–66.

(5)

Castleman, A. W.; Khanna, S. N. N.; Castleman Jr., A. W. W.; Khanna, S. N. N. Clusters, Superatoms, and Building Blocks of New Materials. J. Phys. Chem. C 2009, 113, 2664– 2675.

ACS Paragon Plus Environment

33

Journal of Chemical Theory and Computation 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

(6)

Page 34 of 57

Claridge, S. A.; Castleman Jr, A. W.; Khanna, S. N.; Murray, C. B.; Sen, A.; Weiss, P. S. Cluster-Assembled Materials. ACS Nano 2009, 3, 244–255.

(7)

Jena, P.; Castleman, A. W. Nanoclusters: A Bridge Across Disciplines; Elsevier, 2010; Vol. 1.

(8)

Castleman, A. W.; Jena, P. Clusters: A Bridge between Disciplines. Proc. Natl. Acad. Sci. 2006, 103, 10552 LP-10553.

(9)

Hosmane, N. S. Boron Science: New Technologies and Applications; CRC press, 2011.

(10)

Ammu, M.; Thalappil, P. Noble Metal Clusters: Applications in Energy, Environment, and Biology. Part. Part. Syst. Charact. 2014, 31, 1017–1053.

(11)

Xin Ting Zheng , Arundithi Ananthanarayanan , Kathy Qian Luo, and P. C. Glowing Graphene Quantum Dots and Carbon Dots : Properties , Syntheses , and Biological Applications Glowing Graphene Quantum Dots and Carbon Dots : Properties , Syntheses , and Biological Applications. Small 2014, 11, 1620–1636.

(12)

Shang, L.; Dong, S.; Nienhaus, G. U. Ultra-Small Fluorescent Metal Nanoclusters: Synthesis and Biological Applications. Nano Today 2011, 6, 401–418.

(13)

Daniel, M. C.; Astruc, D. Gold Nanoparticles: Assemble, Supramolecular Chemistry, Quantum-Size-Related Properties, and Applications toward Biology, Catalysis, and Nanotechnology. Chem. Rev. 2004, 104, 293–346.

(14)

Choi, S.; Dickson, R. M.; Yu, J. Developing Luminescent Silver Nanodots for Biological Applications. Chem. Soc. Rev. 2012, 41, 1867–1891.

(15)

Fan, J. Photoelectron Spectroscopy of Size-Selected Transition Metal Clusters : Fen- (n = 324). J. Chem. Phys. 1995, 102, 9480–9493.

(16)

León, I.; Yang, Z.; Liu, H. T.; Wang, L. S. The Design and Construction of a HighResolution Velocity-Map Imaging Apparatus for Photoelectron Spectroscopy Studies of Size-Selected Clusters. Rev. Sci. Instrum. 2014, 85, 83106.

ACS Paragon Plus Environment

34

Page 35 of 57 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

Journal of Chemical Theory and Computation

(17)

Li, X.; Kuznetsov, A. E.; Zhang, H. F.; Boldyrev, A. I.; Wang, L. S. Observation of AllMetal Aromatic Molecules. Science (80-. ). 2001, 291, 859–861.

(18)

Li, J.; Li, X.; Zhai, H. J.; Wang, L. S. Au20: A Tetrahedral Cluster. Science (80-. ). 2003, 299, 864–867.

(19)

Ji, M.; Gu, X.; Li, X.; Gong, X.; Li, J.; Wang, L.-S. Experimental and Theoretical Investigation of the Electronic and Geometrical Structures of the Au32 Cluster. Angew. Chemie 2005, 117, 7281–7285.

(20)

Lewis, G. N. THE ATOM AND THE MOLECULE. J. Am. Chem. Soc. 1916, 38, 762–785.

(21)

TSUCHIDA, R. A New Simple Theory of Valency. Nippon KAGAKU KAISHI 1939, 60, 245–256.

(22)

Sidgwick, N. V.; Powell, H. M. Bakerian Lecture. Stereochemical Types and Valency Groups. Proc. R. Soc. A Math. Phys. Eng. Sci. 1940, 176, 153–180.

(23)

Gillespie, R. J. The Electron-Pair Repulsion Model for Molecular Geometry. J. Chem. Educ. 1970, 47, 18.

(24)

Baletto, F.; Ferrando, R. Structural Properties of Nanoclusters: Energetic, Thermodynamic, and Kinetic Effects. Rev. Mod. Phys. 2005, 77, 371–423.

(25)

Ferrando, R.; Jellinek, J.; Johnston, R. L. Nanoalloys:  From Theory to Applications of Alloy Clusters and Nanoparticles. Chem. Rev. 2008, 108, 845–910.

(26)

Metropolis, N.; Ulam, S. The Monte Carlo Method. J. Am. Stat. Assoc. 1949, 44, 335–341.

(27)

Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087–1092.

(28)

Vanderbilt, D.; Louie, S. G. A Monte Carlo Simulated Annealing Approach to Optimization over Continuous Variables. J. Comput. Phys. 1984, 56, 259–271.

(29)

Schmidt, N. Simulated Annealing. In Simulated annealing: Theory and applications;

ACS Paragon Plus Environment

35

Journal of Chemical Theory and Computation 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

Page 36 of 57

Springer, 1983; pp 1–14. (30)

Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by Simmulated Annealing. Science (80-. ). 1983, 220, 671–680.

(31)

Hartke, B. Global Geometry Optimization of Clusters Using Genetic Algorithms. J. Phys. Chem. 1993, 97, 9973–9976.

(32)

Hartke, B. Global Geometry Optimization of Clusters Using a Growth Strategy Optimized by a Genetic Algorithm. Chem. Phys. Lett. 1995, 240, 560–565.

(33)

Deaven, D. M.; Ho, K.-M. Molecular Geometry Optimization with a Genetic Algorithm. Phys. Rev. Lett. 1995, 75, 288.

(34)

Deaven, D. M.; Tit, N.; Morris, J. R.; Ho, K. M. Structural Optimization of Lennard-Jones Clusters by a Genetic Algorithm. Chem. Phys. Lett. 1996, 256, 195–200.

(35)

Alexandrova, A. N.; Boldyrev, A. I.; Fu, Y.-J.; Yang, X.; Wang, X.-B.; Wang, L.-S. Structure of the NaxClx+ 1−(X= 1–4) Clusters via Ab Initio Genetic Algorithm and Photoelectron Spectroscopy. J. Chem. Phys. 2004, 121, 5709–5719.

(36)

Davis, J. B. A.; Shayeghi, A.; Horswell, S. L.; Johnston, R. L. The Birmingham Parallel Genetic Algorithm and Its Application to the Direct DFT Global Optimisation of Ir N (N= 10–20) Clusters. Nanoscale 2015, 7, 14032–14038.

(37)

Johnston, R. L.; Mortimer-Jones, T. V; Roberts, C.; Darby, S.; Manby, F. R. Application of Genetic Algorithms in Nanoscience: Cluster Geometry Optimization. In Applications of Evolutionary Computing, Proceedings; Springer, 2002; Vol. 2279, pp 92–101.

(38)

Shayeghi, A.; Götz, D.; Davis, J. B. A.; Schäfer, R.; Johnston, R. L. Pool-BCGA: A Parallelised Generation-Free Genetic Algorithm for the Ab Initio Global Optimisation of Nanoalloy Clusters. Phys. Chem. Chem. Phys. 2015, 17, 2104–2112.

(39)

Vargas, J. A.; Buendía, F.; Beltrán, M. R. New AuN (N= 27–30) Lowest Energy Clusters Obtained by Means of an Improved DFT–Genetic Algorithm Methodology. J. Phys. Chem.

ACS Paragon Plus Environment

36

Page 37 of 57 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

Journal of Chemical Theory and Computation

C 2017, 121, 10982–10991. (40)

Kanters, R. P. F.; Donald, K. J. CLUSTER: Searching for Unique Low Energy Minima of Structures Using a Novel Implementation of a Genetic Algorithm. J. Chem. Theory Comput. 2014, 10, 5729–5737.

(41)

Rabanal-León, W. A.; Tiznado, W.; Osorio, E.; Ferraro, F. Exploring the Potential Energy Surface of Small Lead Clusters Using the Gradient Embedded Genetic Algorithm and an Adequate Treatment of Relativistic Effects. RSC Adv. 2018, 8, 145–152.

(42)

Eberhart, R.; Kennedy, J. A New Optimizer Using Particle Swarm Theory. In Micro Machine and Human Science, 1995. MHS’95., Proceedings of the Sixth International Symposium on; IEEE, 1995; pp 39–43.

(43)

Poli, R.; Kennedy, J.; Blackwell, T. Particle Swarm Optimization. Swarm Intell. 2007, 1, 33–57.

(44)

Zhan, Z.-H.; Zhang, J.; Li, Y.; Chung, H. S.-H. Adaptive Particle Swarm Optimization. IEEE Trans. Syst. Man, Cybern. Part B 2009, 39, 1362–1381.

(45)

Call, S. T.; Zubarev, D. Y.; Boldyrev, A. I. Global Minimum Structure Searches via Particle Swarm Optimization. J. Comput. Chem. 2007, 28, 1177–1186.

(46)

Li, Z.; Scheraga, H. A. Monte Carlo-Minimization Approach to the Multiple-Minima Problem in Protein Folding. Proc. Natl. Acad. Sci. 1987, 84, 6611–6615.

(47)

Wales, D. J.; Doye, J. P. K. Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms. J. Phys. Chem. A 1997, 101, 5111–5116.

(48)

White, R. P.; Mayne, H. R. An Investigation of Two Approaches to Basin Hopping Minimization for Atomic and Molecular Clusters. Chem. Phys. Lett. 1998, 289, 463–468.

(49)

Liberti, L.; Maculan, N. Global Optimization: From Theory to Implementation; Springer Science & Business Media, 2006; Vol. 84.

ACS Paragon Plus Environment

37

Journal of Chemical Theory and Computation 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

(50)

Page 38 of 57

Zhao, Y.; Chen, X.; Li, J. TGMin: A Global-Minimum Structure Search Program Based on a Constrained Basin-Hopping Algorithm. Nano Res. 2017, 10, 3407–3420.

(51)

Saunders, M. Stochastic Exploration of Molecular Mechanics Energy Surfaces. Hunting for the Global Minimum. J. Am. Chem. Soc. 1987, 109, 3150–3152.

(52)

Bera, P. P.; Schleyer, P. v R.; Schaefer III, H. F. Periodane: A Wealth of Structural Possibilities Revealed by the Kick Procedure. Int. J. Quantum Chem. 2007, 107, 2220– 2223.

(53)

Averkiev, B. Geometry and Electronic Structure of Doped Clusters via the Coalescence Kick Method, Utah State University, 2009.

(54)

Metha, G. F.; Addicoat, M. A. Kick: Constraining a Stochastic Search Procedure with Molecullar Fragments. J. Comput. Chem. 2009, 30, 57–64.

(55)

Zhai, H.; Ha, M. A.; Alexandrova, A. N. AFFCK: Adaptive Force-Field-Assisted Ab Initio Coalescence Kick Method for Global Minimum Search. J. Chem. Theory Comput. 2015, 11, 2385–2393.

(56)

Heiles, S.; Johnston, R. L. Global Optimization of Clusters Using Electronic Structure Methods. Int. J. Quantum Chem. 2013, 113, 2091–2109.

(57)

Zhang, J.; Dolg, M. Global Optimization of Clusters of Rigid Molecules Using the Artificial Bee Colony Algorithm. Phys. Chem. Chem. Phys. 2016, 18, 3003–3010.

(58)

Jackson, K. A.; Horoi, M.; Chaudhuri, I.; Frauenheim, T.; Shvartsburg, A. A. Unraveling the Shape Transformation in Silicon Clusters. Phys. Rev. Lett. 2004, 93, 13401.

(59)

Fabrice, A.; Clemence, C. Identifying Clusters as Low-Lying Mimina—Efficiency of Stochastic and Genetic Algorithms Using Inexpensive Electronic Structure Levels. J. Comput. Chem. 2011, 33, 502–508.

(60)

Zhao, J.; Shi, R.; Sai, L.; Huang, X.; Su, Y. Comprehensive Genetic Algorithm for Ab Initio Global Optimisation of Clusters. Mol. Simul. 2016, 42, 809–819.

ACS Paragon Plus Environment

38

Page 39 of 57 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

Journal of Chemical Theory and Computation

(61)

Johnston, R. L. Evolving Better Nanoparticles: Genetic Algorithms for Optimising Cluster Geometries. Dalt. Trans. 2003, No. 22, 4193–4207.

(62)

Alexandrova, A. N.; Boldyrev, A. I. Search of the Lin0/+1/-1 (n = 5-7) Lowest-Energy Structures Using the Ab-Initio Gradient Embedded Genetic Algorithm (GEGA). Elucidation of the Chemical Bonding in the Lithium Clusters. J. Chem. Theory Comput. 2005, 1, 566.

(63)

Alexandrova, A. N. H·(H2O) n Clusters: Microsolvation of the Hydrogen Atom via Molecular Ab Initio Gradient Embedded Genetic Algorithm (GEGA). J. Phys. Chem. A 2010, 114, 12591–12599.

(64)

Alejandro, V.; Pino-rios, R. A Fukui Function-Guided Genetic Algorithm . Assessment on Structural Prediction of Sin (n = 12-20) Clusters. J. Comput. Chem. 2017.

(65)

Fernández, R.; Louis, P.-Y.; Nardi, F. R. Overview: PCA Models and Issues. In Probabilistic Cellular Automata: Theory, Applications and Future Perspectives; Louis, P.Y., Nardi, F. R., Eds.; Springer International Publishing: Cham, 2018; pp 1–30.

(66)

Cordero, B.; Gómez, V.; Platero-Prats, A. E.; Revés, M.; Echeverría, J.; Cremades, E.; Barragán, F.; Alvarez, S. Covalent Radii Revisited. Dalt. Trans. 2008, No. 21, 2832–2838.

(67)

Raabe, D. Cellular Automata in Materials Science with Particular Reference to Recrystallization Simulation. Annu. Rev. Mater. Res. 2002, 32, 53–76.

(68)

Deaven, D. M.; Ho, K. M. Molecular Geometry Optimization with a Genetic Algorithm. Phys. Rev. Lett. 1995, 75, 288.

(69)

Rogan, J.; Ramírez, M.; Varas, A.; Kiwi, M. How Relevant Is the Choice of Classical Potentials in Finding Minimal Energy Cluster Conformations? Comput. Theor. Chem. 2013, 1021, 155–163.

(70)

Grigoryan, V. G.; Alamanova, D.; Springborg, M. Structure and Energetics of Nickel, Copper, and Gold Clusters. Eur. Phys. J. D-Atomic, Mol. Opt. Plasma Phys. 2005, 34, 187– 190.

ACS Paragon Plus Environment

39

Journal of Chemical Theory and Computation 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

(71)

Page 40 of 57

Grigoryan, V. G.; Springborg, M. Structure and Energetics of Ni Clusters with up to 150 Atoms. Chem. Phys. Lett. 2003, 375, 219–226.

(72)

Grigoryan, V. ~G.; Springborg, M. Structural and Energetic Properties of Nickel Clusters: 2 ≤ N ≤ 150. Phys. Rev. B 2004, 70, 205415.

(73)

Bäck, T.; Fogel, D. B.; Michalewicz, Z. Evolutionary Computation 1: Basic Algorithms and Operators; CRC press, 2000; Vol. 1.

(74)

Oiso, M.; Yasuda, T.; Ohkura, K.; Matumura, Y. Accelerating Steady-State Genetic Algorithms Based on CUDA Architecture. In Evolutionary Computation (CEC), 2011 IEEE Congress on; IEEE, 2011; pp 687–692.

(75)

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision C.01. Gaussian 09, Revision C.01, Gaussian, Inc., Wallingford CT. Wallingford CT 2010.

(76)

Adamo, C.; Barone, V. Toward Realiable Density Functional Methods Without Ajustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158.

(77)

Fuentealba, P.; Preuss, H.; Stoll, H.; Von Szentpály, L. A Proper Account of CorePolarization with Pseudopotentials: Single Valence-Electron Alkali Compounds. Chem. Phys. Lett. 1982, 89, 418–422.

(78)

Bergner, A.; Dolg, M.; Kuchle, W.; Stoll, H.; Preuss, H. Ab-Initio Energy-Adjusted Pseudopotentials for Elements of Groups 13-17. Mol. Phys. 1993, 80, 1431–1441.

(79)

Küchle, W.; Dolg, M.; Stoll, H.; Preuss, H. Energy‐adjusted Pseudopotentials for the Actinides. Parameter Sets and Test Calculations for Thorium and Thorium Monoxide. J. Chem. Phys. 1994, 100, 7535–7542.

(80)

Fuentealba, P.; Von Szentpaly, L.; Preuss, H.; Stoll, H. Pseudopotential Calculations for Alkaline-Earth Atoms. J. Phys. B At. Mol. Phys. 1985, 18, 1287.

(81)

Küchle, W.; Dolg, M.; Stoll, H.; Preuss, H. Ab Initio Pseudopotentials for Hg through Rn:

ACS Paragon Plus Environment

40

Page 41 of 57 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

Journal of Chemical Theory and Computation

I. Parameter Sets and Atomic Calculations. Mol. Phys. 1991, 74, 1245–1263. (82)

Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305.

(83)

van Lenthe, E.; Ehlers, A.; Baerends, E.-J. Geometry Optimizations in the Zero Order Regular Approximation for Relativistic Effects. J. Chem. Phys. 1999, 110, 8943–8953.

(84)

Van Lenthe, E.; Van Leeuwen, R.; Baerends, E. J.; Snijders, J. G. Relativistic Regular TwoComponent Hamiltonians. Int. J. Quantum Chem. 1996, 57, 281–293.

(85)

Van Lenthe, E.; Baerends, E. J. Optimized Slater Type Basis Sets for the Elements 1-118. J. Comput. Chem. 2003, 24, 1142–1156.

(86)

Baerends, E. J.; Ziegler, T.; Atkins, A. J.; Autschbach, J.; Bashford, D.; Baseggio, O.; Bérces, A.; Bickelhaupt, F. M.; Bo, C.; Boerritger, P. M.; et al. ADF2012.01, SCM, Theoretical Chemistry, Vrije Universiteit. The Netherlands, https://www.scm.com.

(87)

Gillett-Kunnath, M. M.; Muñoz-Castro, A.; Sevov, S. C. Tri-Metallic Deltahedral Zintl Ions: Experimental and Theoretical Studies of the Novel Dimer [(Sn6Ge2Bi)2]4-. Chem. Commun. 2012, 48, 3524–3526.

(88)

Muñoz-Castro, A.; Sevov, S. C. Trimetallic Deltahedral Zintl Ions [Sn9-m-nGemBin](4-n)- for n = 1-4 and m = 0-(9 - n): A Theoretical Survey with Prediction and Rationalization of the Possible Structures. Phys. Chem. Chem. Phys. 2013, 15, 986–991.

(89)

Wade, K. The Structural Significance of the Number of Skeletal Bonding Electron-Pairs in Carboranes, the Higher Boranes and Borane Anions, and Various Transition-Metal Carbonyl Cluster Compounds. J. Chem. Soc. D Chem. Commun. 1971, No. 15, 792–793.

(90)

Mingos, D. M. P. A General Theory for Cluster and Ring Compounds of the Main Group and Transition Elements. Nature 1972, 236, 99–102.

(91)

Welch, A. J. The Significance and Impact of Wade’s Rules. Chem. Commun. 2013, 49,

ACS Paragon Plus Environment

41

Journal of Chemical Theory and Computation 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

Page 42 of 57

3615–3616. (92)

Scharfe, S.; Kraus, F.; Stegmaier, S.; Schier, A.; Faessler, T. F. Zintl Ions, Cage Compounds, and Intermetalloid Clusters of Group 14 and Group 15 Elements. Angew. Chemie Int. Ed. 2011, 50, 3630–3670.

(93)

Michael, D.; Mingos, P.; Johnston, R. L. Theoretical Models of Cluster Bonding. In Theoretical Approaches; Springer, 1987; pp 29–87.

(94)

Heiles, S.; Hofmann, K.; Johnston, R. L.; Schäfer, R. Nine‐Atom Tin‐Bismuth Clusters: Mimicking Excess Electrons by Element Substitution. Chempluschem 2012, 77, 532–535.

(95)

Dinadayalane, T. C.; Priyakumar, U. D.; Sastry, G. N. Exploration of C6H6 Potential Energy Surface:  A Computational Effort to Unravel the Relative Stabilities and Synthetic Feasibility of New Benzene Isomers. J. Phys. Chem. A 2004, 108, 11433–11448.

(96)

Perez-Peralta, N.; Contreras, M.; Tiznado, W.; Stewart, J.; Donald, K. J.; Merino, G. Stabilizing Carbon-Lithium Stars. Phys. Chem. Chem. Phys. 2011, 13, 12975–12980.

(97)

Tiznado, W.; Perez-Peralta, N.; Islas, R.; Toro-Labbe, A.; Ugalde, J. M.; Merino, G. Designing 3-D Molecular Stars. J. Am. Chem. Soc. 2009, 131, 9426–9431.

(98)

Torres-Vega, J. J.; Vásquez-Espinal, A.; Beltran, M. J.; Ruiz, L.; Islas, R.; Tiznado, W. Li7(BH)5+: A New Thermodynamically Favored Star-Shaped Molecule. Phys. Chem. Chem. Phys. 2015, 17, 19602–19606.

(99)

Contreras, M.; Osorio, E.; Ferraro, F.; Puga, G.; Donald, K. J.; Harrison, J. G.; Merino, G.; Tiznado, W. Isomerization Energy Decomposition Analysis for Highly Ionic Systems: Case Study of Starlike E5Li7+ Clusters. Chem. Eur. J. 2013, 19, 2305–2310.

(100) Vásquez-Espinal, A.; Palacio-Rodríguez, K.; Ravell, E.; Orozco-Ic, M.; Barroso, J.; Pan, S.; Tiznado, W.; Merino, G. E5M7+(E=C-Pb, M=Li-Cs): A Source of Viable Star-Shaped Clusters. Chem. - An Asian J. 2018, 13, 1751–1755. (101) Sergeeva, A. P.; Popov, I. A.; Piazza, Z. A.; Li, W.-L.; Romanescu, C.; Wang, L.-S.;

ACS Paragon Plus Environment

42

Page 43 of 57 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

Journal of Chemical Theory and Computation

Boldyrev, A. I. Understanding Boron through Size-Selected Clusters: Structure, Chemical Bonding, and Fluxionality. Acc. Chem. Res. 2014, 47, 1349–1358. (102) Oña, O. B.; Torres-Vega, J. J.; Torre, A.; Lain, L.; Alcoba, D. R.; Vásquez-Espinal, A.; Tiznado, W. Chemical Bonding Analysis in Boron Clusters by Means of Localized Orbitals According to the Electron Localization Function Topology. Theor. Chem. Acc. 2015, 134, 1–9. (103) Olson, J. K.; Boldyrev, A. I. Electronic Transmutation: Boron Acquiring an Extra Electron Becomes ‘Carbon.’ Chem. Phys. Lett. 2012, 523, 83–86. (104) Zubarev, D. Y.; Boldyrev, A. I. Comprehensive Analysis of Chemical Bonding in Boron Clusters. J. Comput. Chem. 2007, 28, 251–268. (105) Galeev, T. R.; Chen, Q.; Guo, J. C.; Bai, H.; Miao, C. Q.; Lu, H. G.; Sergeeva, A. P.; Li, S. D.; Boldyrev, A. I. Deciphering the Mystery of Hexagon Holes in an All-Boron Grapheme α-Sheet. Phys. Chem. Chem. Phys. 2011, 13, 11575–11578. (106) Carter, T. J.; Mohtadi, R.; Arthur, T. S.; Shirai, S.; Kampf, J. W. Boron Clusters as Highly Stable Magnesium Battery Electrolytes. Angew. Chemie Int. Ed. 2014, 53, 3173–3177. (107) Jiménez-Halla, J. O. C.; Islas, R.; Heine, T.; Merino, G. B19-: An Aromatic Wankel Motor. Angew. Chemie - Int. Ed. 2010, 49, 5668–5671. (108) Zhang, J.; Sergeeva, A. P.; Sparta, M.; Alexandrova, A. N. B13+: A Photodriven Molecular Wankel Engine. Angew. Chemie - Int. Ed. 2012, 51, 8512–8515. (109) Merino, G.; Heine, T. And Yet It Rotates: The Starter for a Molecular Wankel Motor. Angew. Chemie Int. Ed. 2012, 51, 10226–10227. (110) Liu, L.; Moreno, D.; Osorio, E.; Castro, A. C.; Pan, S.; Chattaraj, P. K.; Heine, T.; Merino, G. Structure and Bonding of IrB12- : Converting a Rigid Boron B12 Platelet to a Wankel Motor. RSC Adv. 2016, 6, 27177–27182. (111) Cervantes-Navarro, F.; Martínez-Guajardo, G.; Osorio, E.; Moreno, D.; Tiznado, W.; Islas,

ACS Paragon Plus Environment

43

Journal of Chemical Theory and Computation 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

Page 44 of 57

R.; Donald, K. J.; Merino, G. Stop Rotating! One Substitution Halts the B19-Motor. Chem. Commun. 2014, 50, 10680–10682. (112) Boldyrev, A. I.; Wang, L. S. Beyond Organic Chemistry: Aromaticity in Atomic Clusters. Phys. Chem. Chem. Phys. 2016, 18, 11589–11605. (113) Zhai, H. J.; Alexandrova, A. N.; Birch, K. A.; Boldyrev, A. I.; Wang, L. S. Hepta- and Octacoordinate Boron in Molecular Wheels of Eight- and Nine-Atom Boron Clusters: Observation and Confirmation. Angew. Chemie - Int. Ed. 2003, 42, 6004–6008. (114) Ricca, A.; Bauschlicher, C. W. The Structure and Stability of Bn+ Clusters. Chem. Phys. 1995, 208, 233–242. (115) Gu, F. L.; Yang, X.; Tang, A.-C.; Jiao, H.; von R. Schleyer, P. Structure and Stability of B13+ Clusters. J. Comput. Chem. 1998, 19, 203–214. (116) Huang, W.; Sergeeva, A. P.; Zhai, H.; Averkiev, B. B.; Wang, L.; Boldyrev, A. I. A Concentric Planar Doubly π-Aromatic B19- Cluster. Nat. Chem. 2010, 2, 202–206. (117) Kemsley, J. Molecules That Could Be. C&E News 2007, Auguts 13, 17. (118) Averkiev, B. B.; Zubarev, D. Y.; Wang, L.-M.; Huang, W.; Wang, L.-S.; Boldyrev, A. I. Carbon Avoids Hypercoordination in CB6− , CB62− , and C2B5− Planar Carbon−Boron Clusters. J. Am. Chem. Soc. 2008, 130, 9248–9250. (119) Jemmis, E. D. Overlap Control and Stability of Polyhedral Molecules. Closo-Carboranes. J. Am. Chem. Soc. 1982, 104, 7017–7020. (120) Heine, T.; Merino, G. What Is the Maximum Coordination Number in a Planar Structure? Angew. Chemie Int. Ed. 2012, 51, 4275–4276. (121) Romanescu, C.; Galeev, T. R.; Li, W.; Boldyrev, A. I.; Wang, L. Aromatic Metal‐Centered Monocyclic Boron Rings: Co@B8- and Ru@B9-. Angew. Chemie Int. Ed. 2011, 50, 9334– 9337. (122) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput.

ACS Paragon Plus Environment

44

Page 45 of 57 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

Journal of Chemical Theory and Computation

Phys. 1995, 117, 1–19. (123) Hoare, M. R.; Pal, P. Physical Cluster Mechanics: Statics and Energy Surfaces for Monatomic Systems. Adv. Phys. 1971, 20, 161–196. (124) Hoare, M. R.; Pal, P. Statics and Stability of Small Cluster Nuclei. Nat. Phys. Sci. 1971, 230, 5. (125) Hoare, M. R.; Pal, P. Geometry and Stability of “Spherical” Fcc Microcrystallites. Nat. Phys. Sci. 1972, 236, 35–37. (126) Northby, J. A. Structure and Binding of Lennard‐Jones Clusters: 13≤ N≤ 147. J. Chem. Phys. 1987, 87, 6166–6177. (127) Pillardy, J.; Piela, L. Molecular Dynamics on Deformed Potential Energy Hypersurfaces. J. Phys. Chem. 1995, 99, 11805–11812. (128) Doye, J. P. K.; Wales, D. J.; Berry, R. S. The Effect of the Range of the Potential on the Structures of Clusters. J. Chem. Phys. 1995, 103, 4234–4249.

ACS Paragon Plus Environment

45

Journal of Chemical Theory and Computation 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

Page 46 of 57

GRAPHIC FOR MANUSCRIPT

It is introduced AUTOMATON, a novel program for the search of global minimum structures, specially proposed to cover sub-nano materials, where realism on the computations strongly depend of the level of theory.

ACS Paragon Plus Environment

46

Page 47 of 57 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

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

ACS Paragon Plus Environment

Page 48 of 57

Page 49 of 57 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

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

ACS Paragon Plus Environment

Page 50 of 57

Page 51 of 57 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

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

ACS Paragon Plus Environment

Page 52 of 57

Page 53 of 57 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

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

ACS Paragon Plus Environment

Page 54 of 57

Page 55 of 57 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

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

ACS Paragon Plus Environment

Page 56 of 57

Page 57 of 57 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

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment