Monte Carlo Temperature Basin Paving with Effective Fragment

ARTICLE pubs.acs.org/JPCA

Monte Carlo Temperature Basin Paving with Effective Fragment Potential: An Efficient and Fast Method for Finding Low-Energy Structures of Water Clusters (H2O)20 and (H2O)25 Sudhanshu Shanker and Pradipta Bandyopadhyay* School of Computational and Integrative Sciences, Jawaharlal Nehru University, New Delhi, India 110067 ABSTRACT: Determining low-energy structures of large water clusters is a challenge for any optimization algorithm. In this work, we have developed a new Monte Carlo (MC)-based method, temperature basin paving (TBP), which is related to the well-known basin hopping method. In the TBP method, the Boltzmann weight factor used in MC methods is dynamically modified based on the history of the simulation. The states that are visited more are given a lower probability by increasing their temperatures and vice versa. This allows faster escapes from the states frequently visited in the simulation. We have used the TBP method to find a large number of low-energy minima of water clusters of size 20 and 25. We have found structures energetically same to the global minimum structures known for these two clusters. We have compared the efficiency of this method to the basinhopping method and found that it can locate the minima faster. Statistical efficiency of the new method has been investigated by running a large number of trajectories. The new method can locate low-energy structures of both the clusters faster than some of the reported algorithms for water clusters and can switch between high energy and low-energy structures multiple times in a simulation illustrating its efficiency. The large number of minima obtained from the simulations is used to get both general and specific features of the minima. The distribution of minima for these two clusters based on the similarity of their oxygen frames shows that the (H2O)20 can have different variety of structures, but for (H2O)25, low-energy structures are mostly cagelike. Several (H2O)25 structures are found with similar energy but with different cage architectures. Noncage structures of (H2O)25 are also found but they are 67 kcal/mol higher in energy from the global minimum. The TBP method is likely to play an important role for exploring the complex energy landscape of large molecules.

’ INTRODUCTION Finding low-energy structures of water clusters has generated considerable interest in last two decades or so. There are two primary reasons for this interest. One is the inherent importance of water cluster structures in understanding the transition of structure and bonding from the gas phase to liquid water. Both theoretical and experimental works are being pursued to know the details of water cluster structure and spectroscopy and their connection to liquid water.15 The other reason is the challenge it poses to the optimization algorithms to find a few low-energy structures from a rugged energy landscape comprising of astronomical numbers of minima. This is similar to finding the native structure of a protein from its extremely complex conformational space. It has been shown that finding the native state of a protein from its conformational space is a NP hard problem in the language of computational complexity theory by Unger and Moult.6 However, for water cluster optimization no such proof exists. Finding low-energy structures for floppy molecules like large water clusters is quite different from the protocols used for small rigid molecules. For small molecules typically local optimization methods, such as NewtonRaphson, are used to get stationary points starting from a good initial structure built from chemical intuition. This protocol can also be used for small water r 2011 American Chemical Society

clusters. However, as the size of the cluster increases, finding good starting structures using chemical intuition becomes difficult. Generally speaking, the total number of minima increases at least exponentially with the increase of water cluster size. Moreover, the importance of the lowest energy structure also diminishes, as there can be a large number of structures close in energy and a proper understanding of structures and spectroscopy would require knowledge of several low-energy structures. Another issue in optimization of water clusters is the use of an accurate potential energy function. Although a quantum mechanical description is preferable, due to its large computational cost various classical potential functions are typically used for water cluster optimization. Both fixed charge and polarizable potential functions have been used for water clusters. There has been a large number of investigations of water cluster structures using different optimization techniques and potential functions. In one of the pioneering works, Wales et al. determined the water cluster structures using the basin hopping (BH) algorithm with TIP4P potential up to size (H2O)21.7 Received: August 2, 2011 Revised: September 13, 2011 Published: September 19, 2011 11866

dx.doi.org/10.1021/jp2073864 | J. Phys. Chem. A 2011, 115, 11866–11875

The Journal of Physical Chemistry A Tsai et al. used the Monte Carlo (MC) simulated annealing method to determine water cluster structures up to (H2O)20 with TIP4P potential.8 Day et al. used a different potential, effective fragment potential (EFP), to determine water cluster structures of size 620 using the BH algorithm.9 EFP method was also used by Kemp et al. to determine low-energy structures up to size 50.10 For even larger clusters, Kazaimarski and Buch used a combination of methods comprising of molecular dynamics simulation (MD), hydrogen bond network optimization and diffusion MC to determine water cluster structures of size 48, 123, and 293.11 Genetic algorithms and related methods were used by Bandow et al. and Kabrede et al.1214 Among the ab initio investigations, Maheshwari et al. determined water clusters up to size 20 using the HartreeFock and density functional theory.15 Goddard’s group determined the structures up to size 19 using DFT.16 In a recent work, Liu et al. investigated the energetic and fragmentation stability of water cluster of size 2 to 30 at the MP2 level of theory.17 In the past few years, several new optimization algorithms have been developed and used for larger water clusters. The minima hopping method, originally developed by Godecker has been extended by the group of Thakkar to determine water cluster structures up to size 37.18,19 In another development, Takeuchi used topology based investigation of water clusters and determined the low-energy structures up to (H2O)30.20 Bandyopadhyay used the basin paving (BP) algorithm to determine lowenergy structure of (H2O)20 and (H2O)50.21 Among the important techniques used for water cluster optimization, the BH method stands out because of its several successful applications. In this method, Markov Chain MC method is combined with local optimization method. This method effectively eliminates the barrier in the energy landscape since the MC acceptance/ rejection is done between the minima only. However, despite the success of the BH method for water clusters, it may not work well always. For instance, if the MC move gives a structure in the same potential well as the starting structure, then moving out of the well could be difficult with the BH method. An extension of the BH method is the BP method, where the acceptance/rejection of states are done based on the history of the simulation. In the current work, we have developed a new optimization algorithm, which is a variant of the BP algorithm but much more efficient. In this method, one energy range is defined in which simulation will be run. The total energy range is divided into different bins. Each bin is given a temperature and the temperature of a bin increases with each visit to the bin. Thus, states visited more frequently have higher temperatures than the states visited less frequently. If the system is stuck in a local minimum, the temperature of the state increases monotonically and this eventually takes it out of the minimum. It is found that the system can come out of any minimum very quickly in this method. By use of this method we consistently got low-energy structures of (H2O)20 and (H2O)25 in much smaller number of steps compared to some other algorithms used by previous workers. The new method will be denoted by temperature basin paving (TBP) in the rest of this paper. By use of this newly developed algorithm, we investigated the low-energy structures of two large water clusters, (H2O)20 and (H2O)25. The water molecules were represented by the effective fragment potential (EFP), a polarizable water potential.22 For (H2O)20, several low-energy structures have been reported by previous workers. Hence, (H2O)20 was used for benchmarking our algorithm. For (H2O)25, although the lowest-energy

ARTICLE

structure has been reported by previous workers, a thorough characterization of low-energy structures has not been done. We found more than a thousand low-energy distinct minima for (H2O)25. We investigated the distribution of these minima for these two clusters and found that for (H2O)20, low-energy structures could have very different oxygen frames with respect to the lowest-energy structure; however, for (H2O)25, the lowenergy structures were mostly cage types, similar to the lowestenergy structure for this cluster. The putative global minimum structure of (H2O)20 obtained by previous workers was found in this work. For (H2O)25, we found structures with essentially the same energy as the global minimum structure found by previous workers. We also found noncage structures for (H2O)25, which are of about 67 kcal/mol higher in energy than the global minimum structure. The efficiency of TBP method in locating low-energy structures was investigated by running a large number of trajectories. The manuscript is arranged in the following manner. The next section explains the TBP method. This is followed by the details of simulation. The analysis of structures has been discussed in Results and Discussion section. The manuscript ends with a conclusion.

’ METHODOLOGY The MC technique is one of the well-known classes of optimization methods for molecular systems. One of the most widely used MC techniques used for geometry optimization is the BH technique.23,24 In this method, at each MC step, the structure is modified in two steps. At first, the structure is changed randomly, then it is minimized to its nearest minimum. The MC acceptance/rejection is done between the minima only. In effect it eliminates the barriers in the energy landscape. However, in practice the BH method can still have problems in moving out of potential well and can stuck in a narrow region of energy surface. One improvement over the BH method is the BP method.25,26 In this technique, the search for the low-energy structures is done by modifying the Boltzmann factor based on the history of the simulation. One histogram of energy is generated during the simulation by defining an energy range. The histogram of energy is put at the exponent of the exponential of the Boltzmann factor making the probability of states more visited less and vice versa. The new weight factor W(E) is modified in the following way WðE, tÞ ¼ eβðE þ CHðE, tÞÞ

ð1Þ

where β is 1/kT, k is the Boltzmann constant, and T is the temperature. H(E) is the histogram of energy, and t is the MC step number. C is a system specific parameter. This method has been used successfully for peptides, water clusters, RNA.21,25,27 In the current work, we have introduced a modified version of BP, TBP, which is more efficient than the original implementation. In the new method, an energy dependent temperature has been used. At the beginning of the simulation, all energy bins have same temperature. However, during the simulation, temperature is modified as TðE, tÞ ¼ Tinitial þ C0 Tinitial HðE, tÞ

ð2Þ

where Tinitial is the initial temperature, T(E,t) is the temperature for energy bin E at the MC step t. H(E,t) is the histogram of energy at MC step t. C0 is a normalization constant. If any energy 11867

dx.doi.org/10.1021/jp2073864 |J. Phys. Chem. A 2011, 115, 11866–11875

The Journal of Physical Chemistry A

ARTICLE

bin is frequently visited, the temperature for this energy bin will increase monotonically and the chance of moving out of the highly visited states increases as its temperature increases. The MC acceptance/rejection is done in the following manner. If energy of the new state is less than that of the old state, then the new state is always accepted. On the other hand, if the energy of the new state is more than that of the old state, then the new state is accepted with the probability minð1, expð βold ðEnew Eold ÞÞ

ð3Þ

Where βold is 1/kTold, Told is the temperature of the old state. Enew and Eold are the energies of new and old states, respectively. It is to be noted that the temperature of the new state is not used during the acceptance process. This is because the idea is that the trial MC move should be able to move out from the current state. If the temperature of the new state is taken in the acceptance/ rejection, the acceptance condition becomes expð βnew Enew min 1, expð βold Eold Þ So, if the temperature of the new state is low, it will make βnew high, and the ratio of probabilities from the exponentials will be low, and it will be difficult to move out from the current structure. For this reason, temperature of only old state is used for MC acceptance/rejction. Of course, this choice of acceptance condition does not satisfy the detailed balance condition and hence cannot be used for calculating thermodynamic averages. If the temperature of some highly visited energy bins continuously increases in the simulation, there could be two types of problem. The first is that acceptance of new structures with high energy also increases. The second one is that it essentially precludes the visit of the highly visited states in rest of the simulation. Both of these may reduce the efficiency of the simulation. To overcome this, a temperature cooling scheme was used, where effective temperature was slowly reduced using the following expression HðE, tÞ ¼ ðIntðHðE, tÞ CFHðE, tÞÞÞ

ð4Þ

Where Int converts the floating point value to the lowest integer and CF is a cooling factor constant.

’ EFP In the current work, water clusters are described by the EFP method. In this method, the weak intermolecular interactions are described by different contributions such as electrostatics, polarization, exchange-repulsion/charge-transfer. In the EFP1 method used in the current work, electrostatics and polarization was determined from first principle by performing ab initio calculations on a water molecule. The short-range repulsion has been constructed by a fitting precedure. EFP1 method has been used for a large number of applications and it has been found that for water clusters results are very close to the restricted Hartree Fock (RHF) calculation with a double-ζ basis set. ’ DETAILS OF SIMULATION The TBP method has been implemented in a local version of GAMESS program package.28 One important issue in MC technique is the kind of random moves to be given to a molecule. In literature different kind of moves were used for water clusters. In the investigation by the group of Wales,7 the moves were translation and rotation. However, for larger water clusters,

topological moves such as changing the hydrogen bond network have also been used by different groups.11,19 We tried these two kinds of move both independently and together. We found that simple translational and rotational moves, when implemented in our algorithm gave structures close to global minimum consistently. Although topological moves also gave a few structures low in energy, generally speaking these moves were found to be less efficient than translation and rotation. We used a temperaturedependent translational move. When the energy of the structures came to a high temperature region, the translational move was increased and vice versa. The rotational moves were given in terms of Euler angles of the system. The cooling factor was used after every 20 MC steps. The energy range initially chosen for (H2O)20 was 192 to 165 kcal/mol. The energy range chosen for (H2O)25 was 246 to 176 kcal/mol. The bin size was taken as 0.5 kcal/mol. After some trials, C0 value in eq 2 was taken as 0.003. The cooling factor value used was 0.001. In some runs the number of water molecules to move were kept fixed at 2 and 3 for (H2O)20 and (H2O)25, respectively. For some other runs, the number of water molecules to move was taken randomly between 1 and 5. We did not see significant difference in results for these two choices. For classifying the large number of minima obtained from simulation, we did two types of analysis. In rootmean-square distance (rmsd) analysis, each pair of minima were superimposed using the algorithm of Kabasch and rmsd between the two molecules were calculated.29 In dissimilarity index (DSI) analysis, used by Buch and co-workers,10 all oxygenoxygen pair wise distances were calculated for a minimum. All oxygenoxygen distances were sorted from maximum to minimum. The set of distances with n(n1)/2 elements (where n is the number of oxygens) for one structure can be thought as a vector. The Euclidian distance between a pair of minima (i.e., between two vectors) gave the DSI between the two structures. No Hessian calculation was done for the minima found. All the calculations were done without any symmetry imposed on water clusters.

’ RESULTS AND DISCUSSION Figure 1 shows representative trajectories for BH, BP, and TBP methods for (H2O)20 starting from the same structure. All three runs were started at temperature 300 K. The initial energy was at a high value of 174 kcal/mol. It can be seen from the figure that both BP and TBP method consistently gave lower energies than the BH method. The BH method got trapped in a narrow region of energy several times. The lowest energy obtained with the BH method was about 187.0 kcal/mol. The BP method performed better than the BH method; however, it also got stuck in some energy region several times. It usually took between 500 and 1000 MC steps to come out from a confined region. The lowest energy obtained with the BP method was about 1 kcal/mol less than that found by the BH method. The TBP method was found to be most efficient among the three methods. It can come out of any region even it gets stuck there temporarily. Usually it took less than 100 MC steps to come out from a potential well. Although, it sampled high energy regions after coming out from the potential well, it could reach structures with energy lower than those found by the other two methods consistently. It reached energy below 188 kcal/mol four times. It is to be noted that it is possible to get low-energy structures with both BH and BP methods as shown for the BH method by the group of Wales and for the BP method by Bandyopadhyay. However, it is clear that the new method moved to the 11868



ARTICLE

Figure 1. Comparison of BH, BP, and TBP for a representative MC trajectory of (H2O)20.

low-energy region faster than these two methods. Both BP and TBP also sampled energies much higher than the BH method. This is due to the non-Boltzmann dynamic weight factor used in these methods, which increased the probability of states less visited. This is one of the reasons for better performance of both BP and TBP method over the BH method. The TBP method changes the weight factor more efficiently than the original BP method and hence can locate low-energy minima much faster. Similar findings were obtained for (H2O)25 also. Figure 2 shows the variation of energy and temperature for the TBP method for (H2O)25 for one representative trajectory. This run was started with temperature 500 K. It can be seen that at the beginning of the figure the energy was stuck around 234 kcal/mol. The temperature increased from 500 K to about 900 K monotonically in the beginning of the trajectory. At this point, the system came out of that minima and reached a state, where the temperature of the bin, where the new state lies, was still 500 K. There are other regions in the plot, where the temperature increased if the system could not escape from a minimum. The energy plot between about 300 and 500 steps shows that the system was getting stuck and coming out and again getting stuck. The temperature plot in the same region shows the same feature. The temperature increased and then decreased and then increased again. It can be seen that in most of the cases, the system could come out of the potential well within 100 MC steps. Table 1 shows the statistics coming from 17 independent production runs for (H2O)25 with the TBP method. In the cases where energies of the starting structures were the same, other simulation parameters such as random number seed and temperatures were varied to make the trajectories independent. In the table reaching energy lower than 240 kcal/mol has been

defined as reaching low-energy structure. The lowest-energy structure of (H2O)25 found by Thakkar’s group18 when the TIP4P potential was reoptimized with EFP1 potential, and its energy was found to be 243.79 kcal/mol. It can be seen from the table that in all 17 trajectories, the TBP method finds lowenergy minima. Ten runs found energies below 242 kcal/mol, and 2 runs found energy below 243 kcal/mol. In most of the cases low-energy structures (below 240 kcal/mol) were found within 22000 steps (with three exceptions). The lowest-energy structures found in the trajectories were generally obtained in less than 100000 steps except for two cases. This can be compared with more than one million steps used in ref 18. It is to be noted that the structures were taken as fully random structures as can be seen from their very high initial energies. Performance can become even better if structures are prepared using chemical intuition. Since, the purpose of the current work is to check the efficacy of the algorithm, no chemical intuition was used to prepare the starting structures. The last column of the table shows how many times the system moved to high energy structures and again returned to low-energy structures, which can be defined as tunneling events. The high energy is defined as 5 kcal/mol more than that of the low-energy structure, i.e., 235 kcal/mol. It can be seen that all the runs moved from high to low-energy structures, the maximum being 16 times. It is an indicator that, even if the system goes out of the low-energy minima during the MC search, it can return to lowenergy minima. This is expected from the TBP method, since it is designed to escape from the frequently visited states. We have checked these multiple low-energy structures; they were found to be mostly different indicating the algorithm is not finding the same structure again and again. 11869



ARTICLE

Figure 2. Temperatures and energies for a representative MC trajectory for (H2O)25.

Table 1. Statistics of TBP Runs for (H2O)25 serial

starting

total no. of

first time low-energy structure

lowest-energy structure

steps at which lowest-energy structure was

tunneling

no.

energy

steps

found

found

found

events

1 2

219.50 219.50

105231 114563

19643 21946

242.30 241.31

68229 21982

3

219.50

102834

14571

242.16

15505

7

4

219.50

77537

3133

241.86

24691

11

7 5

5

229.24

104911

79767

240.34

79767

1

6

229.24

115142

9912

243.31

17835

14 16

7

229.24

103222

1927

242.80

93378

8

223.03

26401

20885

242.03

23838

2

9 10

223.03 223.03

119546 119984

39895 19873

243.63 242.05

40896 111327

8 8

11

223.03

117198

6699

240.66

7523

8

12

225.48

108624

36884

241.82

45551

8

13

221.88

110174

11977

242.01

109554

9

14

226.25

110012

10854

241.00

14833

4

15

220.02

109879

22975

242.00

44115

12

16

229.48

114999

5376

241.95

18913

14

17

229.48

109050

3437

242.94

63726

13

Overall Features of the Minima. As mentioned in the Method section, both rmsd and DSI analysis were done for finding the unique structures from a large pool of minima. Figure 3a shows the distribution of structures for (H2O)20 in terms of the DSI with respect to the lowest-energy structure found by Thakkar’s group (shown in Figure 3b). The structure pool was reduced by eliminating similar structures using a DSI threshold of 1 Å before comparing with Thakkar’s structure. This reduced the total number of about 19000 minima found from the simulations to about 150 unique minima. It is to be noted that the DSI deals only with the oxygen frame of the structures. It can be seen from the

figure that the DSI values range from close to zero to close to 12 Å indicating large variation of the oxygen frames of (H2O)20. As can be seen from Figure 3b, the lowest-energy structure has two planes on top of each other. The oxygen frames of the low DSI structures are like two planes superimposed with each other. However, there are other varieties of structure, which are of cage type, and oxygen frames are significantly different from the plate like structures. We have found the global minimum structure found by Thakkar;18 however, this structure is not included in Figure 3a. Figure 4a shows the DSI distribution for (H2O)25 structures with respect to the lowest-energy structure found by Thakkar’s 11870



ARTICLE

Figure 3. (a) Plot of DSI (see text) vs energy for unique minima of (H2O)20. (b) The global minimum structure found in ref 18 for (H2O)20. The DSI values shown in part a are calculated with respect to this structure.

group (the structure is shown in Figure 4b). The number of structures was reduced before comparing with Thakkar’s structure using the same procedure used for (H2O)20. It reduced about 300000 minima found from the simulations to about 3000 minima. The reason for getting many more minima for (H2O)25 compared to (H2O)20 is that 25 sized cluster has been investigated more thoroughly in this work. The distribution of minima for (H2O)25 is found to be very different from that of (H2O)20. It is skewed to the lower value of the DSI (less than 5 Å). This suggests that mostly the (H2O)25 structures are cagelike, similar to the structure shown in Figure 4b. Figures 3 and 4 clearly show the

transition from different type of structures to cage like structures for these clusters. Specific structures are discussed later. Figure 5a shows the histogram of unique structures of (H2O)20 within about 4 kcal/mol in energy from the structure obtained by Thakkar’s group by the rmsd analysis. The energy range was divided into 10 equally spaced bins. The rmsd was calculated by taking all atoms of the molecule, including the hydrogens. The uniqueness of the structures is defined by the rmsd values of 0.4 Å or more. As expected the number of unique structures increases with increase of energy. Five unique structures have energy less than 189 kcal/mol including the lowest-energy 11871



ARTICLE

Figure 4. (a) Plot of DSI (see text) vs energy for unique minima of (H2O)25. (b) The global minimum structure found in ref 18 for (H2O)25. The DSI values shown in part a are calculated with respect to this structure.

structure obtained, which has energy 191.13 kcal/mol (which is same as the Thakkar’s structure18). Figure 5b shows the same plots for (H2O)25 using the same definition of bin as used for (H2O)20. Here rmsd threshold used was 0.2 Å to determine unique structures. A smaller threshold was used for (H2O)25 than (H2O)20, since (H2O)25 structures are very similar to each other. Here, 3 unique structures were found below 243 kcal/mol and 6 unique structures were found below 242 kcal/mol. Discussion on Specific Structures. Figure 6 shows the lowest energy structure of each of the 6 classes of structures obtained from the simulations for (H2O)20. These structures have at least 0.6 Å rmsd between any pair. It is to be noted that there are other classes of structures reported in the literature.7,30,31 For (H2O)20 there are 4 major classes of structures termed as edge-sharing pentagonal prism (structure e in the figure), face sharing

pentagonal prism, dodecahedron, and fused cubes. We did not attempt to find structures from all the classes for (H2O)20. It can be seen from the figure that there is remarkable variation in structures, but their energies are only within 2.5 kcal/mol from the lowest-energy structure reported by Thakkar’s group. Four out of six structures are of fused ring types. One is the cage type, and the other is the plate like. The six classes of structures reported here are named as (a) distorted fused rings, (b) cage, (c) fused 454 parallel, (d) fused 545 rings, (e) plate, and (f) fused 454 perpendicular. The meaning of 454 is three rings with 4, 5, and 4 members are fused. The difference between 454 parallel and 454 perpendicular is that the middle ring has different orientation in these two classes. The plate structure has also been described by edge sharing pentagonal prisms by previous workers. All the structures have 34 hydrogen bonds. 11872



ARTICLE

Figure 6. Six low-energy distinct structures of (H2O)20. The numbers in the figure represent their relative energy in kcal/mol from the global minimum structure.

Figure 5. (a) Energy histogram for low-energy unique structures for (H2O)20. (b) ) Energy histogram for low-energy unique structures for (H2O)25.

Figure 7. Six low-energy distinct structures of (H2O)25. The numbers in the figure represent their relative energy in kcal/mol from the global minimum structure.

Figure 7 shows the lowest energy structures of each of the six unique classes structures for (H2O)25. Any pair of structures among these six has at least 0.4 Å rmsd between them. All the structures shown in Figure 7 have 42 hydrogen bonds and are cage types. The structures are within 3.3 kcal/mol from Thakkar’s lowest energy structure. Even lower energy can be obtained by optimizing the hydrogen positions keeping the oxygen frame intact. However, this was not done, since the purpose of our work is to show unique classes of structures that are very close in energy. Also, using more sophisticated energy calculations such as quantum mechanical calculations may change the relative ordering of structures. The lowest-energy structure found by us essentially has same energy as Thakkar’s structure at the EFP level. However, these two structures differ in their cage architecture. Among the six structures, c and f look very similar; however structure f is more symmetric, while structure c is less symmetric. The general feature of structure f can be described as the following. It has two planes, each having three fused rings of size 6, 5, and 5. However, to accommodate three more waters, one 5-membered ring from each plane and the 6-membered ring from each plane move away from each other. In that place two water molecules come and form another five membered ring between the two rings. One more water molecule is also accommodated in the space created by the movement of two planes. As it is very difficult to understand the structural features of the complicated architecture of (H2O)25 structures, parts a and b of Figure 8 show the structure f schematically. It can be

seen from Figure 8a that the structure is composed of three parts, two planes with six-, five-, and five-membered rings. In between two planes, there are three more waters, which also form other rings. Figure 8b shows them from a different angle. It can be seen from Figure 8b that the six- and five-membered rings move away from each other making room for the three other water molecules to come. The other five-membered ring does not move away from each other. Structure e is different from the other five structures. It has one four membered ring coming out of the main cage. In short, although all the six structures are of cage type, there are distinct differences in the formation of cages. It is clear there can be several other cage architectures with energies very close to each other. Figure 9a shows four other structures, which have different kind of cages and has quite different oxygen frames from the structures shown in Figure 7. These structures were detected from their high DSI values shown in Figure 4a. Structures a, b, and c have at least one (two for structure a) four membered ring coming out of the main cage. However, the cages excluding the four membered ring(s) have different architectures. Structure d has a different architecture than the other 3 structures. Figure 9b shows a noncage structure of (H2O)25, which is about 6.1 kcal/ mol higher in energy than in Thakkar’s structure. This has two planes, similar to the plate like structures for (H2O)20. As far as we know, not all the structures shown in Figures 7 and 9 have been reported before. The algorithm developed in this work has been quite successful in locating low-energy structures, including structures with 11873



ARTICLE

Figure 8. (a) Schematic illustration of structure f shown in Figure 7. It shows the three layers of the cage. The middle layer is shown with dotted lines. (b) The structure illustrated in part a is shown from a different angle. It shows how the two planes move away from the each other to accommodate 3 more water molecules.

same energy as the known global minimum, for very complicated water clusters of size (H2O)20 and (H2O)25. For (H2O)25 our investigation clearly shows that it is important to talk about structures close in energy rather than a single global minimum, since there are several structures within 12 kcal/mol from the lowest-energy structure. A proper understanding of bonding and spectroscopy of large water clusters would require consideration of low-energy structures. Our algorithm could locate both cagelike and noncagelike structures for (H2O)25. Although Figure 1 shows its superiority over the BH and BP method, it is not possible to say if this is better than other algorithms used for water clusters. Since, most of the studies of water clusters have been done with TIP4P potential, it is not possible to compare the efficiencies of those algorithms with ours. However, this algorithm should be one of the most efficient algorithms currently in use for large water cluster optimization. Our algorithm is also different from many of the previous algorithms, which used multiple methods to determine the structures. In our case, only MC-based simulations were used. Results can be improved even more by performing additional operations such as hydrogen positions optimization. We also comment that it is equally important to find optimized MC moves along with the MC algorithm for an efficient optimization technique. One issue is that whether there is any way to check if the structures described do contribute to spectroscopy. This, in principle, can be done by calculating vibrational frequencies and comparing with the free OH frequencies measured by the group of Fujii.3,4 Another issue is how the relative energies for different structures of (H2O)20 and (H2O)25 change if another potential function is used. We plan to investigate this thoroughly in a future publication. In the current work, we have investigated the relative energies of four major classes structures of (H2O)20 (as mentioned before and given in Figure 1 of ref 30) with EFP1 method and compared the results with Fanourgakis et al.’s MollerPlesset perturbation calculation at complete basis set limit (MP2/CBS method) as given previously.30 The relative

Figure 9. (a) Four distinct low-energy structures of (H2O)25 with high DSI values. The numbers show their relative energy in kcal/mol and DSI value with respect to the global minimum structure. (b) One noncage structure of (H2O)25. The numbers show its relative energy in kcal/mol and DSI value with respect to global minimum structure.

ordering of energies of these four structures is found to be same by both EFP1 and MP2/CBS method. The plate structure (edgeshearing pentagonal prism) is the lowest-energy structure in both the methods. The relative energies of face sharing pentagonal prism, fused cubes, and dodecahedron from the global minimum for EFP1 (MP2/CBS) method in kcal/mol are 2.3 (2.9), 2.7 (5.3), and 9.7 (17.8), respectively. EFP1 results were obtained by optimizing the structures at the EFP1 level. As can be anticipated, high level quantum chemical calculation may change the relative energies of other structures as well. That is why it is important to report several low-energy structures and not just one, which can be the global minimum in one potential function but not with another potential function.

’ CONCLUSION A robust and fast MC-based optimization algorithm, temperature basin paving has been developed to determine low-energy structures of large water clusters. It has been shown that this method finds low-energy structures of large water clusters of size 20 and 25 faster than the basin hopping and basin paving methods. This method can come out of any minima within 50100 MC steps. Out of 17 independent trajectories run for cluster size 25, this method can locate low-energy structures consistently. It finds the global minimum for both the clusters. The average number of steps it requires to find low-energy 11874


The Journal of Physical Chemistry A structures is much less than some of the previously used algorithms. A large number of minima obtained from the simulation are classified by rmsd calculation and DSI analysis. For (H2O)20, low-energy structures can have very different oxygen frames. However, for (H2O)25, all low-energy structures are of cage types. This efficient algorithm can be used for even larger water clusters.

ARTICLE

(31) Fanourgakis, G. S.; Apra, E.; De Jong, W. A.; Xantheas, S. S. J. Chem. Phys. 2005, 122, 134304.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This work is funded by a Department of Science and Technology grant awarded to P.B. (No. SR/S1/PC-45/2009). ’ REFERENCES (1) Clark, G. N. I.; Cappa, C. D.; Smith, J. D.; Saykally, R. J; HeadGordon, T. Mol. Phys. 2010, 108, 1415. (2) Ohno, K.; Okimura, M.; Akaib, N.; Katsumotoa, Y. Phys. Chem. Chem. Phys. 2005, 7, 3005. (3) Mizuse, K.; Mikami, N.; Fujii, A. Angew. Chem., Int. Ed. 2010, 49, 10119. (4) Hamashima, T.; Mizuse, K.; Fujii, A. J. Phys. Chem. A 2011, 115, 620. (5) Buch, V.; Bauerecker, S.; Devlin, J. P.; Buck, U.; Kazimirski, J. K. Int. Rev. Phys. Chem. 2004, 23, 375. (6) Unger, R.; Moult, J. Bull. Math. Biol. 1993, 55, 1183. (7) Wales, D. J.; Hodges, M. P. Chem. Phys. Lett. 1998, 286, 65. (8) Tsai, C. J.; Jordan, K. D. J. Phys. Chem. 1993, 97, 5208. (9) Day, P. N.; Pachter, R.; Gordon, M. S.; Merrill, G. N. J. Chem. Phys. 2000, 112, 2063. (10) Kemp, D. D.; Gordon, M. S. J. Phys. Chem. A 2008, 112, 4885. (11) Kazimirski, J. K.; Buch, V. J. Phys. Chem. A 2003, 107, 9762. (12) Hartke, B. Z. Phys. Chem. 2000, 214, 1251. (13) Kabrede, H. Chem. Phys. Lett. 2006, 430, 336. (14) Kabrede, H.; Hentschke, R. J. Phys. Chem. B 2003, 107, 3914. (15) Maheshwari, S.; Patel, N.; Sathyamurthy, N.; Kulkarni, A. D.; Gadre, S. R. J. Phys. Chem. A 2001, 105, 293. (16) Su, J. T.; Goddard, W. A., III. J. Phys. Chem. A 2004, 108, 10518. (17) Liu, X. J.; Lu, W. C.; Wang, C. Z.; Ho, K. M. Chem. Phys. Lett. 2011, 508, 270. (18) Goedecker, S. J. Chem. Phys. 2004, 120, 9911. (19) Kazachenko, S.; Thakkar, A. J. Chem. Phys. Lett. 2009, 476, 120. (20) Takeuchi, H. J. Chem. Inf. Model 2008, 48, 2226. (21) Bandyopadhyay, P. Chem. Phys. Lett. 2010, 487, 133. (22) Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. J. Phys. Chem. A 2001, 105, 293. (23) Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A. 1997, 101, 5111. (24) Li, Z.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1987, 84, 6611. (25) Zhan, L.; Chen, J. Z.Y.; Liu, W. K. Phys. Rev. E 2006, 73, 015701(R). (26) Hansmann, U. H. E.; Wille, L. T. Phys. Rev. Lett. 2002, 88, 068105–1. (27) Bandyopadhyay, P.; Kharerin, H. Chem. Phys. Lett. 2011, 502, 130. (28) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus, T. L. J. Comput. Chem. 1993, 14, 1347. (29) Kabsch, W. Acta Crystallogr. 1976, A32, 922. (30) Fanourgakis, G. S.; Apra, E.; Xantheas, S. S. J. Chem. Phys. 2004, 121, 2655. 11875


Monte Carlo Temperature Basin Paving with Effective Fragment

Recommend Documents