An Evolutionary Algorithm for the Global Optimization of Molecular

ARTICLE pubs.acs.org/JPCA

An Evolutionary Algorithm for the Global Optimization of Molecular Clusters: Application to Water, Benzene, and Benzene Cation J. L. Llanio-Trujillo,† J. M. C. Marques,*,† and F. B. Pereira‡,§ †

Departamento de Química, Universidade de Coimbra, 3004-535 Coimbra, Portugal Instituto Superior de Engenharia de Coimbra, Quinta da Nora, 3030-199 Coimbra, Portugal § Centro de Inform
atica e Sistemas da Universidade de Coimbra (CISUC), 3030-290 Coimbra, Portugal ‡

bS Supporting Information ABSTRACT: We have developed an evolutionary algorithm (EA) for the global minimum search of molecular clusters. The EA is able to discover all the putative global minima of water clusters up to (H2O)20 and benzene clusters up to (C6H6)30. Then, the EA was applied to search for the global minima structures of (C6H6)nþ with n = 2-20, some of which were theoretically studied for the first time. Our results for n = 2-6 are consistent with previous theoretical work that uses a similar interaction potential. Excluding the very symmetric global minimum structure for n = 9, the growth pattern of (C6H6)nþ with n g 7 involves the (C6H6)2þ dimer motif, which is placed off-center in the cluster. Such observation indicates that potentials commonly used in the literature for (C6H6)nþ cannot reproduce the icosahedral-type packing suggested by the available experimental data.

1. INTRODUCTION Clusters may be viewed as a bridge across many subjects in physical chemistry,1 ranging from gaseous matter (e.g., in atmospheric chemistry) to condensed phase. It is well-known that the structural optimization of atomic and molecular clusters is a major issue in the study of their properties (see, e.g., ref 2). For instance, global optimization algorithms are very useful to efficiently explore the potential energy surface and identify some important regions that can be subsequently searched at a higher theoretical level (e.g, by applying an accurate ab initio framework). Besides, the process of building up atomic or molecular cluster structures for being optimized with ab initio or DFT methodologies is a cumbersome task, because, in general, the best starting geometries cannot be easily obtained through chemical intuition. Thus, several methods have been proposed3-12 to help in the challenging task of finding out the structure of the cluster that corresponds to the global minimum energy. In recent years, we have developed completely unbiased evolutionary algorithms13-16 (EAs) for the global geometry optimization of atomic clusters. Such methods have been successfully applied to discover the global minima of various types of clusters: argon17 (up to Ar78), short-ranged Morse15 (up to 80 atoms), mixed argon-krypton,16 and binary Lennard-Jones16 (up to 50 atoms). For the latter, the EA was able to update the global minimum reported in the Cambridge Cluster Database18 for the cluster with 38 atoms whose size ratio between unlike atoms is 1.05.16 In comparison with atomic clusters, molecular aggregates present an additional complexity for the optimization procedure, because one needs to take into account the relative orientation of r 2011 American Chemical Society

molecules. In general, different orientations of a molecule may be associated with distinct local minima on the potential energy surface (PES) of the cluster. The modification of the orientational degrees of freedom has been discussed for benzene clusters by Takeuchi19 and, also, by other authors5,20,21 for various molecular systems. Additionally, for clusters such as water, whose structure is dominated by hydrogen bonds, it is known22 that different H-bond networks may be formed, which lead to distinct local minima (though, sometimes, with similar energies). Such knowledge about the system was, then, used by Takeuchi22 to propose specific operators that enhance the efficacy of his global search algorithm. In this work, we propose a completely unbiased hybrid EA to search for the lowest-energy structure of molecular aggregates. Indeed, the present algorithm does not encode any specific information about the system, and hence, we expect that it can be more generally applied to different molecular clusters. Here, the EA is tested for the global optimization of (H2O)n (n = 220) and (C6H6)n (n = 2-30), for which putative global minima have already been published in the literature. In fact, due to its ubiquitous presence in many physical chemistry phenomena, water aggregates are perhaps the most studied ones. Also, a considerable amount of work on neutral benzene clusters has appeared in the literature, since they constitute a prototype to study π-π interactions, which are particularly relevant for the structure of proteins.23,24 In contrast, there is only one theoretical Received: December 10, 2010 Revised: January 31, 2011 Published: March 03, 2011 2130

dx.doi.org/10.1021/jp1117695 | J. Phys. Chem. A 2011, 115, 2130–2138

The Journal of Physical Chemistry A

Figure 1. Flux diagram illustrating the EA for the molecular geometry optimization proposed in the present work. Comments included in dashed-line insets explain some steps of the algorithm.

study25 involving (C6H6)nþ clusters for n up to 6; in turn, we extend such study for the clusters (C6H6)nþ (n = 2-20). Accordingly, the plan of the paper is as follows. In section 2, we describe the molecular coordinates and the EA, while its application to benchmark systems (i.e., water and benzene clusters) is presented in section 3. In addition, results for cation-benzene clusters are given and discussed in section 4. Finally, some remarks concerning the present EA and future work are pointed out in section 5.

2. OPTIMIZATION ALGORITHM The goal of the hybrid EA is to discover the optimal configuration for a cluster with a given number (n) of molecules. It departs from an initial population of randomly generated solutions and iteratively seeks for the configuration with the lowest potential energy; see the flux diagram in Figure 1. A key issue in the design of an EA is the specification of an appropriate representation for the solutions being evolved. Thus, before describing the details of the EA, we present the coordinate frames used to represent the molecular cluster, as well as the gradient evaluation required for the local optimization. 2.1. Coordinates and Gradient Evaluation. The monomer (or molecule) is treated as a rigid-body with a fixed internal geometry, but it is allowed to translate and rotate (i.e., to orient) as a whole in order to lower the interaction energy of the molecular cluster during the optimization procedure. To

ARTICLE

establish the position and orientation of each molecule in the cluster, we define three different right-handed Cartesian coordinate frames: laboratory-fixed (LAB) frame, space-fixed (SF) frame, and body-fixed (BF) frame. The LAB frame is common to the whole cluster, while both SF and BF frames are attached to each individual molecule. The origin of the SF frame may be selected to coincide with the center-of-mass of each molecule, while their axes remain parallel to those of the LAB frame; hence, the SF frame is adequate to describe the translation of the corresponding molecule. In turn, the origin of the BF frame coincides with the SF frame, but it follows the molecular orientation. More specifically, the orientation of each molecule is defined by three Euler angles (φ, θ, ψ) associated with the rotation of the BF axes in the corresponding SF frame; see ref 26 for an alternative rigid-body coordinate system applied in the global optimization of clusters. In this work, the position of each molecule is defined by the center-of-mass Cartesian coordinates in the LAB frame, while Euler angles specify its orientation in the SF frame. Since the search of the cluster’s global minimum by the EA must be carried out in the space of the chosen coordinates, each solution encodes n six-tuples (xJ, yJ, zJ, φJ, θJ, ψJ): in each tuple, the first three variables define the Cartesian coordinates of the Jth molecule center of mass and the last three specify the corresponding Euler angles. In turn, the intermolecular potential (either isotropic or anisotropic) and its derivatives are usually given in terms of relative site-site coordinates. Therefore, the development of a general code for molecular optimization requires the transformation from the aforementioned coordinates used by the EA to those in which the potential is written (e.g., site-site coordinates). As described in the following, this transformation is done, first, to Cartesian coordinates from which the calculation of any kind of coordinates becomes easier; then, the potential is calculated and the derivatives needed for the local optimization are obtained by applying the chain rule. According to the present choice of coordinates, the LAB frame site-position vectors (rLAB i(J) ) of the i-site within the J-monomer can be expressed in terms of the corresponding ) and reference BF frame site-position center-of-mass (T LAB J vectors (rBF i(J)) as rLAB iðJÞ ¼ T

LAB J

þ rSF iðJÞ ¼ T

LAB J

þ R -1 ðφJ , θJ , ψJ ÞrBF iðJÞ

ð1Þ

where R-1 is the inverse of the rotation matrix; R-1 = RT, since R is an orthogonal matrix. Here, we assumed the XYZ-convention for R, i.e., the first rotation is about Z-axis by an angle φ (φ ∈ [0, 2π]), the second one is performed around the intermediary (new) Y-axis by an angle θ (θ ∈ [0, π]) and, finally, the third one corresponds to a change of the angle ψ in the interval [0, 2π] about the newest X-axis. The relation given in eq 1 allows for the transformation from the solution-string coordinates used in the optimization procedure (i.e., the n six-tuples mentioned above) to the LAB Cartesian coordinates of the molecular sites (e.g., the atoms of each molecule). Specifically, to compute the Cartesian components of the LAB frame site-position vectors for a solution resulting from the EA search, eq 1 uses the first three coordinates vectors of all molecules) and the of the n tuples (i.e., the T LAB J Euler angles given by the last three variables of the corresponding tuples. Then, the LAB Cartesian positions of the relevant molecular sites can be employed, in a straightforward way, to calculate the relative coordinates in which the cluster’s potential energy is represented. 2131

dx.doi.org/10.1021/jp1117695 |J. Phys. Chem. A 2011, 115, 2130–2138

The Journal of Physical Chemistry A

ARTICLE

Since the present hybrid EA incorporates the gradient-based L-BFGS local optimization method,27,28 one has to provide also the potential first derivatives with respect to the above-mentioned n six-tuples coordinates (i.e., center-of-mass positions of the molecules and corresponding Euler angles). This requires the calculation of the LAB Cartesian site-gradients from the intermolecular potential first derivatives with respect to the site-site coordinates, as well as the partial first derivatives of RT with respect to the corresponding Euler angles (φJ, θJ, ψJ). 2.2. EA Details. As the EA adopts a fully generational model, in each iteration the whole population is replaced by its offspring. An elitist operator ensures that the quality of the best solution does not deteriorate from one generation to the other. The iterative process ends when a predetermined number of generations is reached. In concrete, the generation of a new set of solutions is completed after performing the following steps (see also Figure 1): (i) Selection of the current best solutions from the population. Tournament selection, a probabilistic operator widely adopted by EAs, is used to choose the parents of the next population. (ii) Application of the genetic operators to the selected individuals. First, crossover is applied, with a given probability, to all pairs of parents. Afterward, mutation is applied, with a given probability, to each gene of the resulting individuals. Below we provide a detailed description of the selected operators. (iii) Local optimization and evaluation of the descendants. Each individual created by the genetic operators must be evaluated to access its relevance to the problem being solved. Before evaluation, the new solution is relaxed until convergence to the nearest local optimum. L-BFGS,27,28 a quasi-Newton local search procedure, is adopted in this step. In the present EA, we rely on state-of-the-art operators for numeric optimization: Simulated binary crossover (SBX)29 and sigma mutation.14 SBX has achieved good results in a number of real-valued problems of varying difficulty and dimensionality29-31 and belongs to the group of arithmetic crossover operators that have the important feature of allowing the modification of gene values. The new values that will compose the children solutions {C1, C2} are obtained from the parents {P1, P2} in the following way:30 (1) Select a random value μ ∈ [0, 1] (2) Calculate β ¼ 2μ1=ðη þ 1Þ ,

if

μ e 0:5

β ¼ ½1=2ð1 - μÞ1=ðη þ 1Þ ,

if

μ > 0:5

(3) Obtain children C1 ¼ 0:5½ð1 þ βÞP1 - ð1 - βÞP2  C2 ¼ 0:5½ð1 - βÞP1 þ ð1 þ βÞP2  where η should be a nonnegative real number. Large values of η increase the probability of creating descendants that are close to their parents, while smaller values allow distant solutions to be reached. As in our previous work,31 we adopt η = 3.0, which provides a good compromise between these two extremes.

Sigma mutation is applied to individual genes and performs a limited modification of the value encoded at that position. If applied to Cartesian coordinates, it slightly modifies the center of mass of a molecule, whereas the application to an Euler angle changes its orientation. The updated value G for a given gene is obtained according to the following expression: G r G þ KσNð0, 1Þ

ð2Þ

where N(0,1) represents a random value sampled from a standard normal distribution, K is the domain size of the variable undergoing mutation, and σ is a parameter from the algorithm that controls the magnitude of the change. This operator ensures that the new value of a mutated variable remains inside domain boundaries. The parameters of the EA have been empirically determined by performing some tests in the optimization of a limited number of water clusters instances. Results obtained were compared to determine a robust EA setting. The adopted values are the following: population size, 100; tournament size, 5; crossover rate, 0.75; mutation rate, 0.1; σ, 0.1.

3. APPLICATION TO BENCHMARK SYSTEMS As a validation test, we have applied the EA to the molecular clusters (H2O)n (with n = 2-20) and (C6H6)n (with n = 2-30), which can be considered benchmark systems for global optimization methods. The difference between these two types of systems assures a suitable evaluation of the robustness of the present EA, in what concerns the global optimization of molecular clusters. While water clusters exhibit dipole-dipole interactions with important hydrogen bonds, the structure of benzene aggregates is dominated by van der Waals interactions. The potentials employed in the global minimum search are the TIP4P32 for water and the one proposed by Williams and Starr33 for benzene clusters. Indeed, global minimum structures for such potentials have been already reported in the literature (see refs 18, 19, 22, 34, 35, and references therein), and hence, they can be directly compared with the results obtained by the EA previously described. While being a benchmark system for our purposes, it is worth noting, however, that the potential of Williams and Starr33 does not reproduce the geometry of the lowest-energy minimum of the benzene dimer obtained by high accurate electronic structure calculations available in the literature.36,37 The global minimum search was carried out by running the EA between 10 and 30 times for each cluster, while performing a total of 250 000 (50 000) evaluations per run for water (benzene). Despite the high number of evaluations granted to the EA, in most cases, the best solutions were discovered in the early stages of the optimization (particularly for smaller clusters). The energies and Cartesian coordinates of all the global minima found out with the new EA are available in the Supporting Information. The global minima were always discovered for both sets of clusters, and we report in Table 1 the difference between the corresponding energies obtained in this work and those from ref 18 for water or from refs 19 and 35 for benzene. We observe that the energy differences, ΔE, for water and benzene clusters are always close to zero: |ΔE| e 1.1  10-4 kJ mol-1 for water clusters and |ΔE| e 2.4  10-2 kJ mol-1 for the aggregates of benzene. Although not shown for the latter in Table 1, the comparison between our global minima energies and those obtained with other genetic algorithms35,38 reveals discrepancies for some cluster’s sizes. This happens for benzene clusters with n = 8, 9, 11, 13-15 in the results of Pullan35 and for n = 11, 13 minima 2132

dx.doi.org/10.1021/jp1117695 |J. Phys. Chem. A 2011, 115, 2130–2138

The Journal of Physical Chemistry A

ARTICLE

Table 1. Energy Difference (ΔE = En - Eref n ) and Rmsd between the Global Minima Obtained in the Present Work (En) and Those (Eref n ) from Ref 18 (Refs 19, 35) for Water (Benzene) Clusters; the Success Rate of the EA Is Also Displayed (see the text) (H2O)n -1

n

ΔE/kJ mol

2

5.6 (-6)

3

-2.3 (-5)

4 5

2.6 (-5) 3.5 (-5)

(C6H6)n -1

% success rate

ΔE/kJ mol

100

2.4 (-2)b

c

0.25 [1.5 (-3)]

100

2.0 (-3)

b

c

100

8.2 (-7) 3.0 (-6)

100 100

-6.0 (-4) 2.1 (-4)

8.8 (-5) 6.3 (-5)

100 100

a

rmsd /Å 4.5 (-7)

rmsda/Å

% success rate 100

6

-4.8 (-6)

1.3 (-6)

100

-5.5 (-4)

0.54 [7.5 (-5)]

100

7

-6.6 (-6)

0.47 [7.6 (-7)]

100

-3.3 (-4)

1.9 (-4)

100

8

-3.7 (-6)

1.3 (-6)

100

-4.8 (-4)

7.1 (-5)

100

9

2.3 (-6)

5.2 (-6)

100

-2.0 (-4)

1.5 (-4)

100

10

-1.9 (-6)

8.3 (-6)

100

-6.8 (-4)

1.50 [2.3 (-4)]

100

11

-3.1 (-5)

0.91 [2.0 (-5)]

100

-9.4 (-4)

9.5 (-5)

100

12 13

-4.7 (-5) -4.0 (-5)

9.4 (-6) 3.0 (-5)

100 100

-1.2 (-3) -1.6 (-3)

1.03 [6.5 (-4)] 2.5 (-4)

100 100

100

-9.6 (-4)

1.04 [2.0 (-4)]

100

77

-1.1 (-3)

1.10 [1.2 (-4)]

100

14

-4.2 (-5)

0.77 [4.3 (-5)]

15

-5.8 (-5)

3.0 (-5)

16

-8.9 (-5)

0.49 [6.0 (-5)]

40

-1.2 (-3)

2.7 (-4)

100

17

-8.6 (-5)

1.11 [4.2 (-5)]

27

-1.6 (-3)

2.7 (-4)

100

18

-8.0 (-5)

0.89 [1.6 (-4)]

3

-1.3 (-3)

1.6 (-4)

19

-1.1 (-4)

4.6 (-5)

27

-2.3 (-3)

1.0 [1.2 (-4)]

100

20 21

-7.9 (-5)

0.96 [7.3 (-5)]

20

-2.4 (-3) -2.3 (-3)

1.9 (-4) 1.17 [2.1 (-4)]

100 90

40

22

-1.9 (-3)

1.33 [1.6 (-4)]

70

23

-3.5 (-3)

1.4 (-4)

90

24

-2.3 (-3)

1.17 [1.3 (-4)]

25

-4.3 (-3)

3.5 (-4)

30

26

-2.4 (-3)

2.7 (-4)

50

100

27

-2.2 (-3)

1.40 [2.3 (-4)]

20

28 29

-3.8 (-3) -2.4 (-3)

1.32 [1.3 (-4)] 1.6 (-4)

80 8

30

-2.8 (-3)

3.0 (-4)

40

a

The rmsd values in square brackets are obtained from the best superposition between the mirror-image of the structure discovered in this work and the global minimum structure published in ref 18 (ref 19) for water (benzene). Since the rmsd values outside square brackets are non-negligible, all such structures discovered in this work are enantiomers (i.e., nonsuperimposing) with the corresponding ones previously reported.18,19 b Energies for (C6H6)2 and (C6H6)3 were taken from ref 35. c Since the coordinates of the original minima of (C6H6)2 and (C6H6)3 were not published,35 it was not possible to calculate the rmsd for these clusters.

obtained in the work of White et al.38 by employing their own genetic algorithm. In such cases, the energy difference is sufficiently large and negative (in relation with the above-mentioned ΔE values for the comparison with the global minima of Takeuchi) to indicate that, probably, the algorithms of Pullan35 and White et al.38 were not able to discover the corresponding optima. Moreover, we ensure that the structures achieved for both systems with the present EA are similar to those reported before18,19 by employing a superimposing method39 and calculating the root-mean-square distance (rmsd) between the corresponding pairs of clusters. It is worth noting that most of the rmsd values shown in Table 1 are less than 0.00035 Å. However, in the case of n = 3, 7, 11, 14, 16-18, and 20 for water and n = 6, 10, 12, 14, 15, 19, 21, 22, 24, 27, and 28 for benzene clusters, the rmsd values obtained after the best superposition are not negligible. This is illustrated in Figure 2 for the particular case of n = 14 for both water and benzene clusters. Since those clusters have essentially the same energy of the published putative global minima (cf. Table 1), it is important to verify whether the corresponding pairs of

structures are enantiomers or, simply, there is no geometric relation between them. To investigate this, we have employed a recently developed scheme,40 which basically tests whether the mirrorimage of a structure is able to overlap the one used as reference (i.e., the corresponding putative global minimum taken from the literature). As we can observe from the negligible rmsd values shown in square brackets in Table 1, the overlap between the pairs of structures is always achieved, and hence, one assigns the abovementioned clusters as enantiomers. This specular relationship between each pair of structures is also apparent from Figure 2a for (H2O)14 and Figure 2b for (C6H6)14, where the clusters were oriented in such a way that allows one to perceive immediately the corresponding mirror plane. Table 1 also reports the success rates (i.e., percentage of runs where the putative optimum was found) obtained by the hybrid EA. It is clear from this table that the success of the EA is notably high for up to n = 15 for water clusters and n = 2-17, 19-24, and 28 in the case of benzene. Due to the increasing ruggedness of the potential energy surface, the success rate tends to decrease with 2133

dx.doi.org/10.1021/jp1117695 |J. Phys. Chem. A 2011, 115, 2130–2138

The Journal of Physical Chemistry A

ARTICLE

Table 2. Energies (En) and Success Rates Achieved for the Putative Global Minima of (C6H6)nþ Clusters by Applying the New EA En/kJ mol-1

% success rate

2

-76.934 091

100

3

-153.983 212

100

4

-214.956 305

100

5

-278.731 570

100

6 7

-342.220 791 -406.704 348

100 100

8

-472.070 586

100

9

-530.401 178

100

Figure 2. Global minimum structures of (H2O)14 (panel a) and (C6H6)14 (panel b) clusters discovered in this work (black) overlapped with the corresponding structures (gray) reported in refs 18 and 19, respectively.

10

-576.141 453

100

11

-617.265 784

90

12

-651.636 729

60

13

-685.285 904

70

the dimension of the system. For atomic clusters, Wille and Vennik41 proved that the global minimization of the PES is NPhard, since the number of local minima grows exponentially with the number of particles. It is widely accepted20,21,42,43 that this result is also valid for molecular clusters, even though a formal proof of NP-hardness does not exist. It is also apparent from Table 1 that some of the clusters studied are particularly difficult instances for the global minimum search, e.g., (H2O)18 and (C6H6)29. Also, large water clusters are, in general, harder to optimize than benzene clusters with the same number of molecules. This has been noticed by applying other methods, and it appears to be associated with the different possibilities for the hydrogen-bond network in the case of water clusters.22 Just to give an example, Niesse and Mayne5 have used a genetic algorithm to optimize water cluster modeled with a different potential and concluded that the performance of the method in obtaining the global minimum begins to degrade for n = 7 and falls to 10% for the simple case of n = 8 (i.e., the largest size for which the success rate was reported).

14 15

-720.394 011 -755.742 431

16 3

16

-794.108 686

23

17

-833.497 096

10

18

-868.046 702

3

19

-903.452 718

3

20

-941.908 205

3

n

4. STUDY OF (C6H6)nþ CLUSTERS Clusters formed by a (C6H6)þ species surrounded by benzene molecules provide a good model to study organic ions solvation. We emphasize that, besides the van der Waals forces of the benzene clusters, induced-dipole interactions due to the presence of a cation are, now, expected to display an important role for the structure and energetics of (C6H6)nþ. Elucidating about the growing pattern of the (C6H6)nþ clusters may be an important additional test for the EA proposed in this work. Thus, we expect that the new results for the microsolvation of the benzene cation obtained with the EA might clarify some aspects concerning the knowledge about this system. In contrast to the large number of works on water and benzene aggregates, there is only one paper25 reporting global minima of the (C6H6)nþ clusters for n = 2-6. Then, we follow closely that previous work25 to establish the potential energy surface in the present study. They have applied25 the OPLS-AA force field for neutral benzene-benzene interactions,44 while the (C6H6)þ-(C6H6) potential energy was modeled as ! 12 12 qi qj e2 Aij Bij Vcation - neutral ¼ ∑ ∑ þ 12 - 6 ð3Þ rij rij rij i j

where the double summation runs over all the 12 sites (corresponding to the atomic positions) of both benzene and (C6H6)þ. In the Coulomb term, e represents the magnitude of the electron charge, while qi and qj are the partial charges assigned to each site (i and j) of neutral benzene and benzene cation molecules, respectively. As in ref 25, the OPLS-AA charges44 were used for the neutral benzene. However, the partial charges calculated for (C6H6)þ are not given in ref 25, and hence, we have used here our own values. These were calculated at the HF/ 6-31G* level of theory using the GAMESS package,45 according to a usual procedure for the OPLS-AA force-field, which involves the application of the RED III.4 program46,47 with an ESP-C1 charge model; the values of the charges for (C6H6)þ are given in Table S1 of the Supporting Information. In turn, values of Aij and Bij are obtained from the usual combining rules, e.g., Aij = (AiiAjj)1/2; Aii (Bii) are related to the corresponding LennardJones (LJ) parameters σ and ε through the expression Aii = 4εiσi12 (Bii = 4εiσi6). Here, we have used for all LJ parameters the values proposed in ref 25. Thus, the present (C6H6)þ-(C6H6) potential energy becomes -76.9 kJ mol-1 (i.e., slightly lower than the value -70.7 kJ mol-1 of ref 25), which is within the recommended range for this interaction (see ref 48 and references therein). In Table 2, we report the energies of the (C6H6)nþ (n = 220) putative global minima obtained with the new EA; see also the Supporting Information for the corresponding Cartesian coordinates. Also shown in Table 2 is the success rate, which accounts for the percentage of runs where the minimum-energy cluster was obtained; in general, the EA was run between 10 and 30 times per instance, but for the important n = 14 case (see below), we performed 110 runs to be more confident about the global minimum obtained. It is apparent from Table 2 that the probability of finding the putative global minimum decreases as the cluster size increases, which is similar to the behavior we have discussed above for water and benzene. Moreover, a mandatory 2134

dx.doi.org/10.1021/jp1117695 |J. Phys. Chem. A 2011, 115, 2130–2138

The Journal of Physical Chemistry A

ARTICLE

Figure 3. Global minimum structures of (C6H6)nþ clusters discovered with the present EA for n = 2-6. Carbon atoms of neutral benzene molecules are represented by gray spheres, while the corresponding ones for (C6H6)þ are in black. See the text.

test for the EA concerning the cation-benzene clusters is the search of the global minima for sizes up to n = 6, which can be compared with the results of Rusyniak et al.25 It is worth noting from Figure 3 that the global minimum structures up to n = 6 discovered by the EA resemble those presented in ref 25; for instance, in agreement with that work,25 we have obtained for the global minimum of (C6H6)2þ a sandwich-type structure skewed by 30
, where the centers of the two rings are separated by 3.63 Å (similar to the value of 3.65 Å reported in ref 25). An important issue for the study of clusters is the identification of a growing pattern. Monte Carlo simulations carried out by Rusyniak et al.25 with a potential similar to the present one allowed them to discover various minima structures of (C6H6)nþ for n = 2-6. All lowest-energy minima obtained for n > 2 show25 as fundamental motif a pancake trimer stack core, which is compatible with gas-phase ion mobility experiments performed in the same work. As displayed in Figure 3, the same pattern is achieved by our global optimization algorithm, which is an expected result, since the two intermolecular potentials might be similar. Alternatively, another growth sequence based on a sandwich dimer core (similar to the above-mentioned most stable stacked motif for n = 2) was also observed for the higher-energy minima structures obtained by Rusyniak et al.25 This growth pattern based on a dimer-core is consistent with the apparent theoretical evidence that the positive charge of (C6H6)3þ is localized around two stacked benzene molecules.49 Moreover, from the experimental side, such localized-charge model has been invoked to suggest an icosahedral packing around a dimer cation, which

Figure 4. As in Figure 3, but for (C6H6)nþ clusters with n = 7-20.

explains the magic numbers observed by Schriver et al.50 for (C6H6)14þ, (C6H6)20þ, (C6H6)24þ, (C6H6)27þ [and, more recently,25 also for (C6H6)30þ]; in particular, (C6H6)14þ has shown to be highly stabilized.51 As a continuation of Figure 3, the putative (C6H6)nþ global minimum structures discovered by our EA for n g 7 are shown in Figure 4. In addition, to help in the analysis of the growing pattern obtained for the (C6H6)nþ global minimum structures, we represent in Figure 5 all the distances between the center of the ring of each benzene molecule and the center of mass of (C6H6)nþ (hereafter designated as dCM) for each cluster size; points corresponding to the smallest distance of all sizes are linked by the solid line, which helps in the inspection of cation dimers that eventually arise in each global minimum structure. As already mentioned above, the growth pattern for clusters with n = 3-6 involves a trimer pancake (see Figure 3), which is also apparent from the two short distances shown by the points that 2135

dx.doi.org/10.1021/jp1117695 |J. Phys. Chem. A 2011, 115, 2130–2138

The Journal of Physical Chemistry A

ARTICLE

Figure 5. Scatter plot of the distances between the center of the benzene rings and the center of mass of the corresponding (C6H6)þ species as a function of the cluster size. A solid line links the smallest-distances points as a guide for the eye.

Figure 6. Contributions for the global minima energy of (C6H6)nþ clusters. The energy contribution from the ion-benzene interactions is represented by the dashed line, while energy from the benzenebenzene pair-potentials is shown as a dotted line. Also shown by the solid line is the energy of the global minima. See the text.

essentially coincide with the solid line in Figure 5. However, as the size of the cluster increases (i.e., for n g 7), only a dimer motif formed by the cation and a neutral benzene can be identified in the global minimum structures (see Figures 4 and 5). An exception to this pattern arises for n = 9 (note that the smallest-distance line in Figure 5 presents a maximum for this size), where the central cation is surrounded by eight neutral benzene molecules (four above and four below the cation’s plane), leading to a quite symmetric structure (cf. Figure 4). It is particularly interesting to notice in Figure 5 that n = 9 is the highest cluster size where the benzene molecules are accommodated around the cation with values of dCM less than ∼5 Å. For n > 9, larger values of dCM begin to appear; for example, the largest values of dCM for n = 12-13, 14-16, and 19 are ∼7 Å, ∼8 Å, and ∼9 Å, respectively. Such structural transition occurring for n > 9 is also apparent in Figure 6, where the contributions of the (C6H6)-(C6H6)þ and benzene-benzene

interactions for the global minima of the clusters are represented as a function of n. It is clear that the major contribution for the energy of the cluster comes from the (C6H6)-(C6H6)þ interactions for n e 9. Conversely, for n > 9, the prevailing contribution for the global minimum potential is due to the interactions between neutral benzene molecules, while the energy from all (C6H6)(C6H6)þ pairs keeps essentially the same value as for n = 9. This observation may be an indication that larger clusters tend to keep a similar distribution of C6H6 molecules around the cation from one size to the next, while the additional benzene molecules of the growing structure come to positions that are not close to (C6H6)þ, and hence, the corresponding interactions become negligible. It is worth noting that the ion-dimer motif shown for larger clusters (mainly for n > 10) is not placed in the center of the structure, and hence, the aforementioned icosahedral packing is not achieved for this potential. Indeed, most of the neutral benzene molecules are preferentially placed on the cation side of the dimer motif (cf. Figure 4), since the interaction involving (C6H6)þ is stronger than the neutral benzene-benzene one. As shown in Figure 6, the contribution of the (C6H6)-(C6H6)þ pairs to the total energy of the global minimum is always higher than the corresponding one from benzene-benzene interactions. In turn, the structural effect associated with the asymmetric position of the dimer motif in the cluster may be due to the fact that the present potential does not explicitly account for the polarization induced by the cation on the neutral benzene forming the stacked dimer. Indeed, the inclusion of many-body polarization terms in the potential has been shown to have important effects on the thermodynamics of ion-pair microsolvated clusters.52 Recently, Sherrill et al.53 have attributed to the missing polarization terms in popular force fields (e.g., OPLSAA) the disagreement observed in the energy curves of π-π prototype systems when compared with highly accurate quantum mechanical results. The main energetic features for (C6H6)nþ as a function of the cluster size are shown in Figure 7; also represented in this figure for comparison are the corresponding curves for benzene clusters. It is apparent from Figure 7a that the average binding 2136

dx.doi.org/10.1021/jp1117695 |J. Phys. Chem. A 2011, 115, 2130–2138

The Journal of Physical Chemistry A

Figure 7. Energetics of (C6H6)nþ (solid line) and benzene (dashed line) clusters: (a) average binding energy and (b) second energy difference. Also shown by the open squares in panel a are the values of the average binding energy for (C6H6)nþ (n = 2-6) reported in ref 25. See the text.

energy, i.e., ÆEnæ = -En/n, presents a maximum value around n = 8-9 for (C6H6)nþ, while it always increases for benzene clusters. Such behavior may be due to the strong interactions involved in the solvation of the cation whose first shell is completed for n = 8-9. Also shown in Figure 7a are the values of ÆEnæ obtained by Rusyniak et al.25 for (C6H6)nþ (n = 2-6). As expected, their values25 essentially resemble those of the present work. In turn, curves of the second difference in energy, i.e., Δ2E = Enþ1 þ En-1 - 2En, are shown in Figure 7b. This quantity is an approximation of the curvature of En (hence, related to the relative stability of each cluster) and it is usually compared with experimental mass spectral intensities. Then, the maxima arising in Figure 7b correspond to particularly stable structures, and the corresponding values of n are the so-called “magic numbers”. It is clear from Figure 7b that benzene clusters present magic numbers for n = 13, 19, 23, 26, and 28. This is compatible with an icosahedral packing already observed by Takeuchi19 for benzene clusters modeled with the potential of Williams and Starr.33 In contrast, magic numbers appear for n = 9, 11, and 17 in the case of (C6H6)nþ clusters; the corresponding structures, shown in Figure 4 are quite symmetric. As expected from experiment,25,50 we emphasize that a magic number would arise for (C6H6)14þ in case of icosahedral packing (i.e., for n þ 1 in relation to the magic numbers in benzene clusters), which would correspond to a central dimer (involving the cation) solvated by 12 equivalent nearest-neighbor molecules of neutral benzene. Although a small maximum appears in Figure 7b for (C6H6)14þ, this essentially corresponds to Δ2E = 0. In addition, the structure shown in Figure 4 for n = 14 is clearly nonicosahedral (as already mentioned above); also, the structure of (C6H6)20þ is not a double icosahedral arrangement (see Figure 4) as would be expected from experiment.25,50 Again, these findings corroborate the nonicosahedral growing pattern for (C6H6)nþ clusters, which is in disagreement with the experimental results.25,50,51

5. FINAL REMARKS We have proposed a new EA for the geometry optimization of molecular clusters that is able to reproduce the global minima for aggregates of both water [up to (H2O)20] and benzene [up to

ARTICLE

(C6H6)30]. It is worth noting that the EA does not incorporate any information about the chemical or physical properties of the cluster system and, hence, it is completely unbiased. A close inspection of the EA effectiveness reveals that the success rate tends to decrease as the clusters grow in size. In accordance, our future research efforts will aim toward enhancements that increase robustness and scalability, as this will allow the effective optimization of large molecular clusters. In addition, we have applied the EA to obtain the global minima for (C6H6)nþ (n = 2-20) clusters. The growth pattern for n = 3-6 involves a trimer pancake structure, where the distances between the center of the two benzene rings and the center of mass of the (C6H6)þ are similar to the corresponding one in the dimer cluster. Such results are in agreement with previous theoretical work25 that uses a similar interaction potential. The growth pattern of the (C6H6)nþ clusters for n g 7, in contrast, involves only a (C6H6)2þ dimer core; the very symmetric structure for the global minimum of n = 9 is, however, an exception to this behavior. Moreover, an important outcome from this work is concerned with the off-center position of the ion-dimer motif within the cluster structure, especially for n > 10. Such asymmetric position of the dimer core in the cluster explains the absence of icosahedral-type magic numbers for n = 14 and n = 20 that are suggested by the available experimental data. Such disagreement with experiment may be attributed to the potential model, which does not explicitly account for the polarization induced by the cation on the neutral benzene forming the stacked dimer. Thus, additional theoretical work on the PES of these clusters is welcome. Finally, we emphasize the importance of developing this kind of algorithms to discover low-energy minima structures that can be, then, reoptimized at a higher level of theory (e.g., by using ab initio or DFT methods). According to this line of research, the application to other systems (including ionic and binary molecular clusters) is underway in our group. Also, an easy-to-use version of the code is being prepared and it will be available from the authors upon request.

’ ASSOCIATED CONTENT

bS

Supporting Information. Supplement structures.zip contains files with the Cartesian coordinates and energies of all structures of water, benzene, and cation-benzene clusters obtained in this work by applying the new EA and a table with the values of the charges used for (C6H6)þ. This material is available free of charge via the Internet at http://pubs.acs.org/.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This work was supported by Fundac-~ao para a Ci^encia e Tecnologia (FCT), Portugal, under project PTDC/QUI/ 69422/2006, which is financed by Programa Operacional Factores de Competitividade (COMPETE) of QREN and FEDER programs. We are grateful to the John von Neumann Institut f€ur Computing, J€ulich, for the provision of supercomputer time on the Juropa-JSC (Project EPG01). We also acknowledge Dr. P. E. Abreu for helping in the calculation of the partial charges of the benzene cation that were used in this work. 2137

dx.doi.org/10.1021/jp1117695 |J. Phys. Chem. A 2011, 115, 2130–2138

The Journal of Physical Chemistry A

’ REFERENCES (1) Jena, P.; A. W. Castleman, J. Proc. Nat. Acad. Sci. U.S.A. 2006, 103, 10560–10569. (2) Wales, D. J. Energy Landscapes: With Applications to Clusters, Biomolecules and Glasses; Cambridge University Press: Cambridge, UK, 2003. (3) Deaven, D. M.; Ho, K. M. Phys. Rev. Lett. 1995, 75, 288–291. (4) Gregurick, S. K.; Alexander, M. H.; Hartke, B. J. Chem. Phys. 1996, 104, 2684. (5) Niesse, J. A.; Mayne, H. R. J. Comput. Chem. 1997, 18, 1233–1244. (6) Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A 1997, 101, 5111. (7) Roberts, C.; Johnston, R. L.; Wilson, N. T. Theor. Chem. Acc. 2000, 104, 123. (8) Leary, R. H. J. Global Optim. 2000, 18, 367. (9) Locatelli, M.; Schoen, F. Comput. Optim. Appl. 2003, 26, 173. (10) Shao, X.; Cheng, L.; Cai, W. J. Comput. Chem. 2004, 25, 1693. (11) Grosso, A.; Locatelli, M.; Schoen, F. Math. Program. Ser. A 2007, 110, 373. (12) Rossi, G.; Ferrando, R. J. Phys.: Condens. Matter 2009, 21, 084208. (13) Pereira, F. B.; Marques, J. M. C.; Leit~ao, T.; Tavares, J. Proceedings of the 2006 IEEE Congress on Evolutionary Computation; Vancouver, 2006; Vols. 1-6, pp 2270-2277. (14) Pereira, F. B.; Marques, J. M. C.; Leit~ao, T.; Tavares, J. Advances in Metaheuristics for Hard Optimization, Springer Natural Computing Series; Springer: Berlin, 2008; pp 223-250. (15) Pereira, F. B.; Marques, J. M. C. Evol. Intel. 2009, 2, 121–140. (16) Marques, J. M. C.; Pereira, F. B. Chem. Phys. Lett. 2010, 485, 211–216. (17) Marques, J. M. C.; Pereira, F. B.; Leit~ao, T. J. Phys. Chem. A 2008, 112, 6079–6089. (18) The Cambridge Cluster Database, http://www-wales.ch.cam. ac.uk/CCD.html, accessed in December, 2010. (19) Takeuchi, H. J. Chem. Inf. Model. 2007, 47, 104. (20) Wales, D. J.; Hodges, M. P. Chem. Phys. Lett. 1998, 286, 65–72. (21) Bandow, B.; Hartke, B. J. Phys. Chem. A 2006, 110, 5809–5822. (22) Takeuchi, H. J. Chem. Inf. Model. 2008, 48, 2226–2233. (23) Burley, S. K.; Petsko, G. A. Science 1985, 229, 23–28. (24) Hunter, C. A.; Sanders, J. K. M. J. Am. Chem. Soc. 1990, 112, 5525–5534. (25) Rusyniak, M. J.; Ibrahim, Y. M.; Wright, D. L.; Khanna, S. N.; El-Shall, M. S. J. Am. Chem. Soc. 2003, 125, 12001–12013. (26) Chakrabarti, D.; Wales, D. J. Phys. Chem. Chem. Phys. 2009, 11, 1970–1976. (27) Nocedal, J. Math. Comp. 1980, 35, 773–782. (28) Liu, D.; Nocedal, J. Math. Program. B 1989, 45, 503–528. (29) Deb, K.; Agrawal, R. B. Complex Systems 1995, 9, 115–148. (30) Deb, K.; Beyer, H.-G. Evol. Comput. 2001, 9, 197–221. (31) Marques, J. M. C.; Prudente, F. V.; Pereira, F. B.; Almeida, M. M.; Maniero, A. M.; Fellows, C. E. J. Phys. B: Atomic Mol. Opt. Phys. 2008, 41, 085103. (32) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (33) Williams, D. E.; Starr, T. L. Comput. Chem. 1977, 1, 173–177. (34) Kazachenko, S.; Thakkar, A. J. Chem. Phys. Lett. 2009, 476, 120–124. (35) Pullan, W. J. J. Chem. Inf. Comput. Sci. 1997, 37, 1189–1193. (36) Podeszwa, R.; Bukowski, R.; Szalewicz, K. J. Phys. Chem. A 2006, 110, 10345–10354. (37) Janowski, T.; Pulay, P. Chem. Phys. Lett. 2007, 447, 27–32. (38) White, R. P.; Niesse, J. A.; Mayne, H. R. J. Chem. Phys. 1998, 108, 2208–2218. (39) V
asquez-P
erez, J. M.; Martínez, G. U. G.; K€oster, A. M.; Calaminici, P. J. Chem. Phys. 2009, 131, 124126. (40) Marques, J. M. C.; Llanio-Trujillo, J. L.; Abreu, P. E.; Pereira, F. B. J. Chem. Inf. Model. 2010, 50, 2129–2140.

ARTICLE

(41) Wille, L. T.; Vennik, J. J. Phys. A: Math. Gen. 1985, 18, L419. (42) Stillinger, F. H. Phys. Rev. E 1999, 59, 48–51. (43) Doye, J. P. K.; Wales, D. J. J. Chem. Phys. 2002, 116, 3777–3788. (44) Jorgensen, W. L.; Severance, D. L. J. Am. Chem. Soc. 1990, 112, 4768–4774. (45) Schmidt, M. W.; Baldridge, K. K.; Boats, J. A.; Elbert, S. T.; Gorgon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyenn, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J., Jr. J. Comput. Chem. 1993, 14, 1347–1363. (46) Pigache, A.; Cieplak, P.; Dupradeau, F.-Y. Automatic and highly reproducible RESP and ESP charge derivation: Application to the development of programs RED and X RED. 227th ACS National Meeting, Anaheim, CA, March 28-April 1, 2008. (47) Dupradeau, F.-Y.; Pigache, A.; Zaffran, T.; Savineau, C.; R. Lelong, N. G.; Lelong, D.; Rosanski, W.; Cieplak, P. Phys. Chem. Chem. Phys. 2010, 12, 7821–7839. (48) Pieniazek, P. A.; Krylov, A. I.; Bradforth, S. E. J. Chem. Phys. 2007, 127, 044317. (49) Mine, M.; Mori, H.; Ogata, M.; Tanaka, S.; Tsutsui, T.; Miyoshi, E. Chem. Phys. Lett. 2007, 438, 157–161. (50) Schriver, K. E.; Paguia, A. J.; Hahn, M. Y.; Honea, E. C.; Camarena, A. M.; Whetten, R. L. J. Phys. Chem. 1987, 91, 3131–3133. (51) Beck, S. M.; Hecht, J. H. J. Chem. Phys. 1992, 96, 1975–1981. (52) Peslherbe, G. H.; Ladanyi, B. M.; Hynes, J. T. J. Phys. Chem. A 2000, 104, 4533–4548. (53) Sherrill, C. D.; Sumpter, B. G.; Sinnokrot, M. O.; Marshall, M. S.; Hohenstein, E. G.; Walker, R. C.; Gould, I. R. J. Comput. Chem. 2009, 30, 2187–2193.

2138

dx.doi.org/10.1021/jp1117695 |J. Phys. Chem. A 2011, 115, 2130–2138