Exploring the Rich Potential Energy Surface of (H2O)11 and Its

Subscriber access provided by READING UNIV

Article

Exploring the Rich Potential Energy Surface of (HO) and Its Physical Implications 2

11

Berhane Temelso, Katurah L. Klein, Joel W. Mabey, Cristóbal Pérez, Brooks H. Pate, Zbigniew Kisiel, and George C Shields J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00938 • Publication Date (Web): 12 Jan 2018 Downloaded from http://pubs.acs.org on January 12, 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 free 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 accessible to all readers and 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.

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

Exploring the Rich Potential Energy Surface of (H2O)11 and Its Physical Implications Berhane Temelso,*1,2 Katurah L. Klein,2 Joel W. Mabey, 2 Cristóbal Pérez,3,4 Brooks H. Pate,3 Zbigniew Kisiel,5 George C. Shields,*1,2

1 2

Department of Chemistry, Furman University, Greenville, SC 29613, USA

Dean’s Office, College of Arts and Sciences, and Department of Chemistry, Bucknell University, Lewisburg, PA 17837, USA 3

Department of Chemistry, University of Virginia, McCormick Rd., Charlottesville, VA 22904-4319, USA 4

Max Planck Institute for the Structure and Dynamics of Matter, Luruper Chausse 149, D-22761 Hamburg, Germany

5

Institute of Physics, Polish Academy of Sciences, Al. Lotników 32/46, 02-668 Warszawa, Poland

January 11, 2018

Corresponding authors: [email protected], [email protected]

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 41

ABSTRACT The rich potential energy surface of the water undecamer, (H2O)11, was explored with a basin hopping algorithm using a TIP4P potential and other methods followed by extensive ab initio MP2 minimizations and CCSD(T) corrections. This protocol yielded 17, 66, and 125 distinct isomers within 0.5, 1.0 and 2.0 kcal mol-1 of the complete basis set CCSD(T) global minimum, respectively. These isomers were categorized into fifteen different families based on their oxygen framework and hydrogen bonding topology. Determination of the global minimum proved challenging because of the presence of many nearly isoenergetic isomers. The predicted global minimum varied between ab initio methods, density functionals and model potentials, and it was sensitive to the choice of energy extrapolation schemes, higher-order CCSD(T) corrections and inclusion of zero-point vibrational energy. The presence of a large number of nearly degenerate structures and the isomerization between them has manifested itself in the anomalous broadening of the heat capacity curve of the undecamer in simulations around the melting region.

KEYWORDS Ab initio methods, MP2, electronic structure, water models, water potentials, global optimization, basin hopping, undecamer, broadband rotational spectroscopy, melting clusters


2

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


1. Introduction Finding the global minimum of large water clusters is very challenging because (a) the number of possible configurations (the configurational space) increases exponentially with the cluster size1 and (b) it is difficult to model the interactions between water monomers accurately.26

Many global optimization approaches including simulated annealing,7 evolutionary algorithm

(EA),8-10 genetic algorithm (GA),11,12 genetic evolutionary algorithm (GEA),13 basin hopping (BH)14,15, minima hopping16-18 and others19,20 have been developed to address the configurational search problem. The configurational search is usually performed on potential energy surface generated using empirical and ab initio-based model potentials. These models were developed to reproduce the properties of bulk water or high accuracy ab initio data for small water clusters.6,21 The most common models used in water cluster studies are SPC/E,22 TIP3P,23 TIP4P,23 TIP4P/2005,24 TIP4P-EW,25 TIP5P,26 TTM2,27 TTM2.1-F,28 MCY,29 AMOEBA30 and HBB2pol.31 Depending on which combination of global optimization algorithm and water model is used, one will often obtain different global minima, especially as the size of the water cluster increases. For a given water model, different global optimization algorithms applied to small water clusters, (H2O)n≤20, often yield the same global minima, but with varying computational efficiencies and scaling. For example, GA was found to be at least twice as fast as BH at finding global minima for (H2O)n, n=10-20.11 However, for a given global optimization algorithm, different water models often yield different global minima, because the potential energy surface of each water model is different, and this variation becomes more pronounced with increasing cluster sizes.11,15,18,32-35 Water models that produce global minima matching more accurate ab initio methods are considered more reliable.32,35-38 Therefore, choosing a good global


3


Page 4 of 41

optimization algorithm and a reliable water model results in a diverse set of low energy minima. This ensemble of low energy structures on the model potential energy surface can then be used as a starting point for optimizations on a quantum mechanical potential energy surface. The many model potentials and quantum mechanical methods used to study the undecamer predict different global minimum structures, as summarized in Figure 1 and Table 1.11,15,18,32-38 Holden and Freeman used Monte Carlo methods to investigate the classical thermodynamic properties of water clusters modeled by the TIP4P potential, and found that at 1 K all configurations quench to one distinct global minimum, and that at 9 K two minima are found that are nearly isoenergetic.39 Other works using model potentials include Guimarães et al. who used a modified genetic algorithm approach on TIP3P and TIP4P water models to determine a number of low energy structures.13 Do and Besley performed a basin hopping search combined with density functional theory to predict the minimum energy undecamer at the B3LYP-D/6-31+G* level38 while Sathyamurthy and co-workers used the Hartree-Fock (HF) method to locate lowlying undecamer structures.36 Belchior and co-workers used an approach based on a genetic algorithm search for conformers combined with a DFT calculation on certain minima.40 Sadlej used MP2 calculations with a double zeta basis set to study the sensitivity of the OH stretching frequency to hydrogen bonding.41 Kim and co-workers performed MP2 calculations with the TZ2P++ basis set to determine families of low energy conformers.32 Lu and co-workers employed MP2/6-31++G(d,p) calculations to locate two minima.42 One of the more comprehensive works was done by Bulusu et al. who employed a four-step searching/screening approach to determine the best candidates for the global minimum for the water undecamer.35 They used the minima hopping global optimization method16 together with four empirical water potentials, SPC/E,22 TIP3P,23

TIP4P,23 and POL343 to create a database of 100 low-lying


4

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


isomers for each water model, 400 in total. They then applied the density functional tight-binding (DFTB)44 method to re-optimize the 400 lowest-lying minima and grouped them into eight distinct isomer families. A high-level work by Iwata et al.45 discovered other higher lying isomers of the undecamer using Monte Carlo Temperature Basin Paving method.46 They applied perturbation theory based on the ab initio locally-projected molecular orbitals (LPMO PT) to explain the relative electronic energy of 10 undecamers. In this work, we explore the rich potential energy surface of the water undecamer with particular emphasis on the challenge of picking a global minimum structure when many isomers are nearly degenerate.

2. Methodology 2.1.

Configurational Search

High temperature molecular dynamics (MD) simulations have been shown to successfully sample the configurational space of small water clusters (n = 7, 9, 10) to obtain a large number of low energy isomers.47-49 Here, we chose basin hopping (BH)14,15 as the global optimization method that explores the TIP4P model potential energy surface. BH or Monte Carlo50 minimization is equivalent to a canonical Monte Carlo simulation, except each step (translation and rotation of rigid water molecules in our case) is followed by a local optimization resulting in a minimum energy structure residing in a basin. This transforms a walk on the PES to barrierless hopping between basins which should sample all the local minima of a system over time. It has been successfully applied to atomic and molecular clusters, crystals and biomolecules.51 Wales and Hodges have used it to find the TIP4P global minima of (H2O)n, n < 21.15 TIP4P is a simple, rigid-monomer, non-polarizable empirical water model that can be used to explore the PES of


5


Page 6 of 41

water clusters inexpensively. It largely agrees with more elaborate water models and ab initio methods when it comes to the global minimum for small clusters with n ≤ 10.11,18,33,52 It is also found to be most compatible with higher-level QM methods for QM/MM simulations of liquid water.53 Beginning with random structures, 20 independent BH simulations were run at a constant temperature equivalent of 2 kJ mol-1 (240 K) for 20,000 steps each. Each Monte Carlo step corresponded to translation or rotation of the rigid TIP4P monomers. As suggested by Wales and Hodges,15 the steps were arranged to be a set of 100 translations, followed by 100 rotations, followed by 100 rotations or translations and a repeat of the same cycle. The maximum displacement and rotation were set to 0.30 Å and 30 degrees, respectively. Each step was followed by a minimization using the limited-memory Broyden-Fletcher-Golfarb-Shanno (LMBGFS) method.54 We used the above-mentioned basin hopping algorithm implemented by Stone et al. in ORIENT 4.7.55 For good measure, we confirmed that the TIP4P global minima we found for (H2O)n, n = 11 - 15 matched those reported by Wales and Hodges15 using the same BH algorithm implemented in an earlier version of ORIENT and a similar work by Kabrede and Hentschke.11 From the tens of thousands of unique isomers found by BH on the TIP4P PES of (H2O)11, we chose the 700 that were within 2 kcal mol-1 of the lowest energy isomer for further processing with the RI-MP2 (resolution-of-the-identity second-order Møller-Plesset perturbation theory) ab initio method. To ensure that the 2 kcal mol-1 energy cutoff on the TIP4P PES was safe, we investigated how well the final set of 102 ab initio isomers that originated from TIP4P basin hopping (BH) searches compare with the initial TIP4P isomers. Our analysis found that all 15 isomers within 0.5 kcal mol-1 of the ab initio global minimum isomer that originated from TIP4P


6

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


BH were within 1.1 kcal mol-1 of the TIP4P global minimum. The details of the analysis are summarized in the Supporting Information Section S1.

2.2.

Quantum Mechanical Calculations

RI-MP2 is an excellent method for determining the structures and energetics of water clusters and other systems bound by hydrogen bonds.47-49,56-59 The 700 isomers from TIP4P BH searches were minimized using the RI-MP2/6-31G* method and the resulting minima were subject to single-point energy calculations using RI-MP2/aug-cc-pVDZ (RI-MP2/aVDZ). This mapping from a TIP4P PES to the RI-MP2/aVDZ//RI-MP2/6-31G* PES resulted in 102 isomers within 2 kcal mol-1 of a new global minimum on the RI-MP2/aVDZ//RI-MP2/6-31G* surface. At this point, we combined these isomers with eighteen others found by adding a single donor-single acceptor (DA) water to pentagonal sides of a PP1 or PP2 water decamer47,48 and five unique structures found using an MD sampling approach described in our previous papers.47-49 The reasoning behind including structures beyond those arising from BH searches is to sample the configurational space more completely. On one extreme, the systematic inclusion of isomers derived by inserting water monomers to a PP decamer ensured that we did not miss highly symmetric structures that may be hard to capture from more stochastic methods like BH. At the other extreme, sampling structures randomly from a high-temperature MD simulation yields more disordered structures. These 125 isomers were then optimized and subject to harmonic vibrational frequency calculations using RI-MP2/aVDZ as implemented in ORCA 3.0.3.60 Successful geometry optimizations satisfied convergence criteria corresponding to energy, RMS (root mean squared) gradient, maximum gradient, RMS displacement, and maximum displacement thresholds of


7


Page 8 of 41

1x10-7 atomic units (a.u.), 1x10-5 a.u. Å-1, 1x10-4 a.u. Å-1, 3x10-4 Å and 1x10-3 Å, respectively. Harmonic vibrational frequencies were calculated at the same level of theory from central finite differences of analytic gradients at 0.002 Å increments. The harmonic vibrational frequency calculations confirmed that they are local minima. Their uniqueness was determined at each level of theory by examining their rotational constants and relative electronic energy. In a few cases where the rotational constants and relative electronic energies were not enough to differentiate between isomers, structural alignments using ArbAlign61 and calculated principal dipole moment components were used to assess their uniqueness. ArbAlign61 uses the Kuhn-Munkres algorithm to optimally overlay any two arbitrarily ordered isomers and determines their similarity and uniqueness. When using a finite Gaussian basis set, basis set incompleteness (BSIE) and superposition (BSSE) errors need to be addressed properly to get accurate results. Using correlation consistent basis sets like the aug-cc-pVNZ

62,63

series allows the extrapolation of RI-MP2 energies to the

complete basis set limit (CBS) systematically. Of the various extrapolation schemes proposed in the literature, the 4-5 inverse polynomial extrapolation scheme64,65 has been used extensively for water clusters: RI−MP 2 ECBS = E NRI−MP 2 +

b c + 4 (N +1) (N +1)5

(1)

RI−MP 2 RI−MP 2 where E N is the RI-MP2/aVNZ//aVDZ energy, ECBS is the extrapolated RI-MP2/CBS

energy, N is the largest angular momentum number for the aVNZ basis set (N = 2, 3, 4 for N = D, T, Q, respectively), and b and c are fitting parameters. For the RI-MP2/aVNZ//aVDZ (N = D, T, Q) sequence, the CBS energy has the following form: ଵ

ோூିெ௉ଶ ‫ܧ‬஼஻ௌ = ଵଷଶ଴ ሾ243 ∗ ‫ܧ‬ଶோூିெ௉ଶ − 2048 ∗ ‫ܧ‬ଷோூିெ௉ଶ + 3125 ∗ ‫ܧ‬ସோூିெ௉ଶ ሿ

(2)


8

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


The RI-MP2/CBS binding energies are compared to benchmark explicitly correlated RIMP2-F1266 results. We chose to extrapolate the energies to their complete basis set limit and compare them with explicitly correlated RI-MP2-F12 results. The explicitly correlated RI-MP2F1266 method is implemented in Molpro 2010.1.67 We used the RI-MP2-F12 method employing the 3C(FIX) ansatz

68

with cc-pVNZ-F12 (VNZ-F12) orbital basis,69 the recommended VNZ-

F12/OPTRI auxiliary basis70 and VNZ/MP2FIT density fitting basis.71 To estimate the coupled-cluster with singles, doubles, and perturbative triple excitations ஼஼ௌ஽(்)

[CCSD(T)] CBS relative energy, we added CCSD(T) corrections (ߜெ௉ଶ

) to the RI-MP2/CBS

energy at the RI-MP2/aVDZ optimized geometry: ஼஼ௌ஽(்)

‫ܧ‬஼஻ௌ

஼஼ௌ஽(்)

ߜெ௉ଶ

஼஼ௌ஽(்)

ோூିெ௉ଶ ≈ ‫ܧ‬஼஻ௌ + ߜெ௉ଶ ஼஼ௌ஽(்)

= ‫ܧ‬௔௏஽௓

ெ௉ଶ − ‫ܧ‬௔௏஽௓

(3)

CCSD(T ) This δMP 2 correction is the energy difference between the conventional MP2 and

CCSD(T) energies, without any density fitting. Even though including small basis CCSD(T) correction to the binding energy (∆Ee) of hydrogen bonded systems is not always encouraged,7275

we investigated the effect of this correction on the relative energy (∆∆Ee) of the nearly

degenerate isomers and concluded that it is justified.76 Comparison of this effect using conventional

CCSD(T),77

and

explicitly

correlated

CCSD(T*)-F12a,

CCSD(T*)-F12,

CCSD(T**)-F12a and CCSD(T**)-F12b73,78-80 and corresponding MP2 and MP2-F12 methods yielded consistent corrections to the relative electronic energy across the different methods. Therefore, the inclusion of CCSD(T) correction is justified and essential.76 In the rest of this study, canonical MP2-CCSD(T) corrections, as opposed to explicitly correlated ones, will be used.


9


Page 10 of 41

RI-MP2/aVDZ zero-point vibrational energy (ZPVE) and finite temperature thermodynamic corrections to the enthalpy (H) and entropy (S) were added to the CCSD(T)/CBS electronic energy to obtain Gibbs free energies, G(T), at standard state conditions of 1 atmosphere pressure and given temperature, T. The thermodynamic corrections were derived from a canonical partition function built assuming ideal gas conditions with a rigid-rotor approximation for molecular rotations and a harmonic oscillator model for vibrations.81,82 The harmonic vibrational frequencies used to calculate ZPVE and vibrational thermodynamic corrections are not scaled.

2.3.

Classification and Nomenclature

Of the 125 minima found to be within 2 kcal mol-1 of the CCSD(T)/CBS electronic energy global minimum, 17 were within 0.5 kcal mol-1 and 66 were within 1 kcal mol-1. Such a high density of low energy structures is unprecedented for any small water cluster or even larger (12 ≤ n ≤ 15) ones. These minima can be categorized into 15 different families based on their oxygen skeleton, rotational constants, and the number and type of hydrogen bonds. Since the heavy oxygen skeleton largely determines the rotational constants of a given isomer, we used the rotational constants to cluster the minima into families. Members of the same isomer family have rotational constants A, B, and C that fall within 0.03 GHz of the lowest energy member of an isomer family. The rotational constants for the lowest energy isomer of each family are given in Table 2. However, the differences in rotational constants were small and it was necessary to use the hydrogen bonding topology of each water molecule in the cluster to further group the isomers to families. Each family has a distinct hydrogen bond identification (HB-id) based on the number of singly coordinated (either D-donor, or A-acceptor), doubly coordinated (DD, AA, and DA), triply coordinated (DAA and DDA) and quadruply coordinated (DDAA, DDDA, DAAA) water molecules in the cluster. These HB-id values are reported in Table 2 as N1-N2-N3-N4, or more


10

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


specifically [D-A]-[DD-AA-DA]-[DAA-DDA]-[DDAA-DDDA-DAAA], so that N1 represents the number of singly coordinated waters (D and A), N2 represents the number of doubly coordinated waters (DD, AA, DA), N3 represents the number of triply coordinated waters (DAA, DDA), and N4 represents the number of quadruply coordinated waters (DDAA, DDDA, DAAA). While quadruply coordinated waters are generally of type DDAA, isomers in the 43’4A family have one water monomer of the rare DAAA type. This monomer is a strong HB donor, with DH…A distance of 1.60 - 1.65 Å, but a weak acceptor, with D-H…A distances of 1.95 - 2.06 Å. All four hydrogen bonds that the DAAA monomer in the 43’4A family members are involved in are well within the range of the geometric definition of a typical hydrogen bond and those criteria, as implemented in UCSF Chimera,83 are defined for different systems by Mills and Dean.84 The naming of these families was designed to be consistent with that of Bulusu et al.35 as much as possible since it is, along with the current work, the most complete published ab initio study of (H2O)11. Bulusu et al.35 reported nine isomer families named 515, 43’4, 55’1, 44’3’, 44’12, bow-tie, 41114, 41141. In this nomenclature, the individual digits indicate the number of water monomers on a given horizontal plane, and the presence of a prime or apostrophe (‘) signifies that those monomers do not form a closed hydrogen bonds ring. Such a naming and classification scheme makes studying and comparing such a complex system manageable. Since the energy differences between the clusters are very small, the binding and relative energies are rounded to the nearest 0.01 kcal mol-1 despite the fact that the error bars in our approach are expected to be much larger. The RI-MP2 geometry optimizations and frequency calculations were performed with the appropriate auxiliary basis85 using ORCA 3.0.3.60 RIMP2/aV[D,T,Q]Z, RI-MP2-F12/V[T,Q]Z-F12, CCSD(T)/aVDZ and CCSD(T)-F12/VDZ-F12


11


Page 12 of 41

single point energy calculations were all performed using Molpro 2010.1. Molecular graphics are generated with UCSF Chimera83 versions 1.8 and 1.9 using the default hydrogen bond distance and angle criteria84 relaxed by 20 degrees.

3. Results Many model potentials and quantum mechanical methods predict different global minimum structures as listed in Table 1 and displayed in Figure 1. The 125 isomers found in the current study are classified into 15 families and the lowest energy member of each family is shown in Figure 2 in order of increasing CCSD(T)/CBS electronic energy and the important signatures of each family are detailed in Table 2. The energy differences between the isomers are so small that their relative electronic energy (∆∆Ee) is extremely sensitive to the choice of basis sets and extrapolation scheme, as illustrated in Figure 3 and Figure 4. Moreover, Figure 5 demonstrates ஼஼ௌ஽(்)

that the higher-order electron correlation correction (ߜெ௉ଶ

) is significant enough to alter the

relative electronic energy of the most stable isomers. Table 3 highlights that the addition of ZPVE corrections and thermal corrections at to the enthalpies and Gibbs free energies at 100 K, 200 K, and 300 K leads to further re-ordering of their relative energies (∆∆Ee, ∆∆H, and ∆∆G). The seventeen lowest energy isomers within 0.50 kcal mol-1 of the CCSD(T)/CBS electronic energy global minimum isomer are displayed in Figure 6. At any temperature between 0 and 300 K, there are at least seventeen isomers that are within 0.5 kcal mol-1 of the Gibbs free energy global minimum isomer, but the precise isomer family they belong to changes substantially, as shown in Figure 7. This makes designating one isomer as the global minimum very challenging. This work concludes that the 43’4A-1 isomer is the most stable on the basis of its CCSD(T)/CBS electronic energy. Figure 7 shows that the inclusion of ZPVE and thermal corrections makes


12

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


515A-1 and 55’1A-1 isomers the most stable between 0 K and 150 K. The 55’1A-5 isomer is predicted to be the most stable up to 300 K even though other members 55’1A remain close. The Boltzmann populations of all 125 isomers at temperatues of 0 K, 50 K, 100 K, 150 K and 300 K are shown in Figure 8. A plot of the density of isomers at different temperatures in Figure 9 further illustrates the presence of many isomers with comparable relative Gibbs free energies. Even though many of these isomers have very distinct oxygen skeletons and hydrogen bonding networks, the fact that they remain energetically competitive at all temperatures is one of the more remarkable features of the water undecamer. Some of these structures are related by the nearly barrierless flipping of free hydrogens and other tunneling pathways that make isomerizations very likely.

4. Discussion 4.1.

Basis Set and Electron Correlation Effects

One of the fascinating aspects of the water undecamer is the unprecedented sensitivity of its isomers’ relative energy to basis set and electron correlation effects. Figure 3 demonstrates that RI-MP2/aVDZ, RI-MP2/aVTZ and RI-MP2/aVQZ energies yield differing relative electronic energies whereas extrapolating to the RI-MP2/CBS limit or employing explicitly correlated RIMP2-F12/VTZ-F12 and RI-MP2-F12/VQZ-F12 methods give similar and more reliable ordering of the isomers. Figure 4 illustrates the same point for the lowest energy members of the four most important isomer families. Table 3 contains the CCSD(T)/CBS relative energies for the 15 isomer families, including the electronic energies, the internal energies with ZPVE at 0 K, and the enthalpies and free energies at 100 K, 200 K, and 300 K. The significant effect of the CCSD(T) correction to MP2 energies on the relative electronic energy of different isomers is


13


Page 14 of 41

displayed in Figure 5. The inclusion of higher-order electron correlation correction changes the relative energy of the 515A-1, 43′4A- 1, 55′1A-1, and 44′3′A-1 isomers from 0.00, 0.03, 0.04, and 0.28 kcal mol−1 at RI-MP2/CBS level to 0.27, 0.00, 0.31, and 0.19 kcal mol−1 at CCSD(T)/CBS level, respectively. Just the addition of the CCSD(T) correction changed the predicted global minimum isomer from 515A-1 to 43′4A-1. A closer examination of Figure 5 reveals that the CCSD(T) correction has little to no effect on the relative electronic energy of 43’4A and 44’3’A families, but it destabilizes 515A and 55’1A isomers by 0.15−0.30 kcal mol−1 depending on the flavor of CCSD(T) correction chosen.

4.2.

Comparison with Previous Studies

Table 1 and Figure 1 summarize the literature on predictions of the global minimum for the undecamer from different model potentials and quantum mechanical methods. Different methods vastly disagree on the global minimum. The TIP4P minimum reported by Wales and Hodges and confirmed by by other research groups is the 44’3 isomer,13,15 which is composed of two tetramers and a trimer stacked on top of each other, with the sandwiched tetramer (designated as 4’) missing one hydrogen bond in the ring. The three-site model potentials (SPC/E, TIP3P and POL3) and five-site TIP5P potentials all predict the 55ʹ1-like global minimum while the four-site TIP4P and the flexible and polarizable ab initio-based potentials, TTM2-F and TTM2.1-F, predict that the 4 family is more stable. The 44’3’ motif consists of a bottom tetramer, and then a middle tetramer and a top trimer that are both missing one of the usual hydrogen bonds in each ring. Clearly, both ab initio and density functional theory (DFT) methods predict very different global minima. It is worth remembering that even in the present study, RI-MP2/CBS predicts 515A-1, 43’4A-1 and 55’1A-1 to be degenerate based on their electronic energy while inclusion


14

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


of a CCSD(T) correction stabilizes the 43’4A-1 isomer. Another ab initio study by Iwata et al.45 has compiled five low energy isomers from Bulusu et al.35 as well as five additional isomers they generated. A re-optimization of these ten isomers and their comparison with our set of isomers found that their eight lowest energy isomers are indeed in our set while the two highest ones are not. Those two isomers are 1.7 kcal mol-1 higher than our global minimum isomer in terms of CCSD(T)/CBS electronic energy. These comparisons are reported in Supporting Information Section S6. Kim et al.32 refer to the 55’1 family of isomers as Pr56 because insertion of a water monomer in between two cyclic pentamers can also be described as a cyclic hexamer stacked on a pentamer.

4.3.

Classification of Isomers

The 125 isomers found to be within 2 kcal mol-1 of the CCSD(T)/CBS electronic energy global minimum were categorized into 15 isomer families based on their a) oxygen skeleton, b) rotational constants, and c) hydrogen bonding characteristics of the individual water molecules in the clusters. Table 2 summarizes the important features used for this classification and we will use the the 44’3A family to demonstrate how the classification and naming scheme works. The family 44’3’A consists of 18 isomers with an electronic energy within 2 kcal mol-1 of the global minimum, and each member of this family has 17 hydrogen bonds. The rotational constants of the lowest energy isomer in the 44’3’A family are 0.606, 0.412, and 0.371 GHz, and the other 17 members have rotational constants within 0.03 GHz of those for the lowest energy isomer of that family. The 44’3’A family differs from the 44’3’B family in terms of its shape, which manifests itself in the rotational constants. The 44’3’ notation refers to a trimer with a non-cyclic hydrogen bond (3’) stacked on top of a tetramer with a non-cyclic hydrogen bond (4’) which is stacked on top of a tetramer with a cyclic hydrogen bond ring (4). The hydrogen bonding pattern of 00-001-


15


Page 16 of 41

44-200 signifies that for the 11 waters in the cluster, no waters are singly coordinated (N1=0), one water is doubly coordinated (N2=1), eight waters are triply coordinated (N3=8), and two waters are quadruply coordinated (N4=2). In addition, the one doubly coordinated water is a Donor/Acceptor (DA), four of the triply coordinated waters are Single Donor/Double Acceptors (DAA), and the other four triply coordinated waters are Double Donors/Single Acceptor (DDA), while the two quadruply coordinated waters are both Double Donors/Double Acceptors (DDAA). The number of isomers within each family ranges from a low of two for 44’12o and 515D to a high of 32 for the 55’1A family. The 515B family has the fewest hydrogen bonds, 15 for each of the five structures within 2 kcal mol-1 of the minimum while the three members of the 443 family have the most hydrogen bonds, 18. The most stable isomer families, 515A, 43’4A, and 55’1A, all have 10 triply coordinated waters. These three isomer families are related to the (H2O)10 and (H2O)12 global minima by the addition or removal of one water molecule. Starting with the PP1 or PP2 decamer47,48 and inserting a water molecule between the two pentamers yields a 515A undecamer whereas adding one on top of a pentamer results in a 55’1A isomer.35 In this isomer two pentamers are stacked, with the 11th water interacting with the top pentamer, and that pentamer is designated as 5’ because it is missing one of the usual hydrogen bonds in the pentamer ring. The presence of 32 isomers in the 55’1A family within 2 kcal mol-1 of the 55’1A-1 isomer highlights the vast number of ways these hydrogen-bonded clusters can form. All those 32 isomers share 16 hydrogen bonds, and a 00-001-55-000 hydrogen bonding network and similar shapes.

4.4.

Relative Stability of Isomers


16

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


The difficulty of assigning a definitive global minimum for (H2O)11 marks it as one of the most interesting water cluster systems along with (H2O)6,86,87 where two- and three-dimensional clusters have competitive, and (H2O)17 where the first internally-solvated water cluster is reported to be the global minimum.88 It is clear from Table 3 that the relative electronic energy differences between 515A, 43’4A, and 55’1A are so close together that we would expect them to be within the error bars of the RI-MP2/CBS calculations. Adding CCSD(T) higher order electron correlation correction to the RI-MP2/CBS electronic energies affects the 43’4A isomers minimally, but destabilizes the 515A and 55’1A isomers by as much as 0.3 kcal mol-1 as shown in Figure 5. The magnitude and sign of the CCSD(T) correction on the relative electronic energy of isomers was heavily dependent on the isomer family. For example, the 515A and 55’1A isomer families were destabilized by an average of 0.29 and 0.27 kcal mol-1, respectively while the 43’4A and 44’3’A isomers are stabilized by an average of 0.01 to 0.09 kcal mol-1, respectively. This correction effectively breaks the degeneracy between the 515A-1, 43’4A-1, and 55’1A-1 at the RI-MP2/CBS level and deems the 43’4A-1 isomer the global minimum at CCSD(T)/CBS level with the 515A-1, and 55’1A-1 isomers lying 0.27 and 0.30 kcal mol-1 higher, respectively. One possible explanation for the stabilizing effect of the CCSD(T) correction on 43’4A-1 over 515A-1 and 55’1A-1 could be dispersion interactions. More compact structures like the 43’4A and 44’3’A isomers have a larger dispersion interaction than MP2 is capturing compared to more loose structures like the 515A and 55’1A isomers.45 This pattern is analogous to the water hexamer where the more compact Prism isomer is stabilized by a CCSD(T) correction relative to the two-dimensional Book, Cyclic Boat and Cyclic S6 isomers because of the larger dispersion component in its binding energy.89


17


Page 18 of 41

The inclusion of ZPVE destabilizes the 43’4A-1 isomer by about 0.3 kcal mol-1 relative to the 515A-1 and 55’1A-1 isomers, thereby bringing the three isomers within 0.1 kcal mol-1 of each other. Thus, the three isomers are essentially degenerate based on the CCSD(T)/CBS E(0 K). Turning to enthalpies and free energies, Table 3 shows that at 100 K, 200 K, and 300K, the 43’4A-1 isomer remains the enthalpic minimum, with the 55’1A-1 and 515A-1 structures being virtually isoenthalpic at the three different temperatures. At 100 K, the 515A-1 and 55’1A-1 structures are the free energy global minima, while the 5th isomer of the 55’1A family, 55’1A-5, becomes the free energy minimum between 150 K and 300 K. It is only 0.11 kcal mol-1 above the 515A-1 isomer at 0 K and its slightly higher entropy accounts for the decrease in its relative Gibbs free energy with rising temperature. Figure 7 displays the relative standard free energies of the lowest electronic energy isomers of the 15 families as a function of temperature. In addition, the 55’1A-5 isomer has been included because it is the G0(T) global minimum for 150 K < T < 300 K. The plot shows that the 515A-1 and 55’1A-1 isomers have the lowest Gibbs free energy at tempreatures below 150 K while the 55’1A-5 isomer, followed by 44’12o-1 and 44’12a-1 will be the lowest Gibbs free energy at higher temperatures up to 300K.

4.5.

Implications of the Rich Potential Energy Surface

Figure 8 reveals the distribution of isomers with increasing temperature. At 0 K only one isomer should be present by definition and that isomer is 515A-1. At 50 K we have a Boltzmann distribution of isomers dominated by the two lowest energy isomers, namely 515A-1 and 55’1A1, which account for about 40% of the population. At any higher temperature, there are many competing isomers with comparable populations. Figure 9 further highights the competition between isomers by showing the energy density of undecamer isomers as a function of temperature. More importantly, it shows which isomer families are prominent at a given


18

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


temperature. For example, the 43’4A isomer family, which is most stable in terms of electronic energy, becomes unfavorable once ZPVE and finite temperature corrections are included. On the other hand, the 515A and 55’1A families gain in relative Gibbs free energy with increasing temperature. The presence of a large number of nearly degenerate structures has been used to explain the anomalous broadening of the constant volume heat capacity, Cv(T) curve of the undecamer in the melting region.39 A pronounced peak in the Cv(T) curve is a signature of phase change because it indicates that the density of states (DOS) above this temperature are considerably larger than at lower temperatures.90 The absence of such a clear differences in the DOS for the undecamer leads to a broad Cv(T) curve without any sharp peaks. Holden and Freeman39 ran Monte Carlo simulations using the TIP4P water potential to investigate the constant volume heat capacity of small water clusters as a function of temperature, Cv(T). While most clusters exhibited a sharp "melting" peak corresponding to significant structural isomerizations at a certain low temperature, the undecamer had a broad curve over a large temperature range. This anomalous behavior is attributed to the interconversion between many nearly degenerate isomers that have very different structures. Likewise, Kaneko et al.'s91 investigation of the melting behavior of water clusters using multicanonical-ensemble molecular dynamics highlighted the peculiar features for the water heptamer and undecamer. The small energy difference between many isomers in each case prevents structures from being trapped in a given configuration, thereby leading to the unusual melting behavior of the heptamer and undecamer. Some of these isomers are related by the nearly barrierless flipping of free hydrogens which result in torsionally averaged structures as observed in the case of the water nonamer.92 This is illustrated in Figure S6 where the conversion of 55’1A-1 to 55’1A-4 isomer is displayed.


19


Page 20 of 41

Because of the low barrier to torsional motion around the nearest hydrogen bond, it is conceivable that both isomers are equivalently populated even at low temperatures. Similar motions allow interconversions between isomers such as 55’1A-2 and 55’1A-5. While we did not rigorously attempt to account for nuclear

quantum effects like zero-point motion and

tunneling, they should further facilitate isomerizations between different structures.49,93-95 Zeropoint motion and tunneling will affect the relative population and interconversion between different isomers at any temperature,96-98 but taking them into account is still conceptually and computationally difficult despite significant progress being made using ring-polymer instanton theory.99-101

5. Conclusions While water clusters up to the decamer, (H2O)10, and larger ones like the dodecamer, (H2O)12, have been known to have a small number of very stable isomers, the undecamer, (H2O)11, has a rich potential energy surface with a much larger number of stable isomers. Our exploration of the water undecamer PES with a basin hopping algorithm using a TIP4P potential and other configurational sampling approaches uncovered a vast number of distinct isomers within a small energy range. Ab initio MP2 minimizations and CCSD(T) energy corrections that followed the TIP4P basin hopping search protocol yielded 17, 66, and 125 distinct isomers within 0.5, 1.0 and 2.0 kcal mol-1 of the complete basis set CCSD(T) global minimum. Using a classification scheme incorporating their rotational constants, oxygen framework and hydrogen bonding topology, we categorized the 125 isomers into fifteen different families. Because of the presence of many nearly isoenergetic isomers, determination of the global minimum proved difficult for ab initio methods as well as model potentials. The predicted global minimum not only varied between ab


20

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


initio and model potentials, but it was highly sensitive to the choice of energy extrapolation schemes, inclusion of zero-point vibrational energy, and nature of the model potential. We have predicted that many isomers belonging to different isomer families will coexist even at the lowest temperatures. Also, isomerizations should be facilitated by torsional motions of free hydrogens and tunneling pathways relating different isomers. These predictions are in line with previous works reporting an anomalous broadening of the heat capacity curve of the undecamer in the melting region, which has been attributed to the presence of a large number of nearly degenerate structures. Attempts to compare these isomers with experimental observations are underway.

Acknowledgement.

The authors thank Prof. Anthony Stone for providing a

modified version of ORIENT 4.7., and Volodymyr Babin and Francesco Paesani for sharing their HBB2-pol code. C. P. acknowledges a Research Fellowship from the Alexander von Humboldt Foundation. Acknowledgment is made to the NSF, Bucknell University and Furman University for their support of this work. This project is supported in part by NSF grants CHE-1213521 and CHE-1508556. Computational resources were provided by NSF grant CHE-1229354 as part of the MERCURY high-performance computer consortium (http://www.mercuryconsortium.org). This research used the NSF XSEDE resources provided by the Texas Advanced Computing Center (TACC) under grant TG-CHE120025 and resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under contract DE-AC02-05CH11231.


21


Page 22 of 41

Supporting Information Available. An expanded methodology, justification of choices in configurational sampling, the optimized Cartesian coordinates, the rotational constants, harmonic vibrational frequencies, absolute, relative and binding energies of all 125 isomers. This information is available free of charge via the Internet at http://pubs.acs.org/. The Cartesian coordinates

and energies are also presented

in a convenient form at

http://www.molecularclusters.xyz


22

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


Table 1. Global Minimum of (H2O)11 Predicted by Different Model Potentials and Quantum Mechanical (QM) Methods. Model Potentialb

Global Minimumb

Reference

QM Method

Global Minimum

Reference

MCY

55'1(Pr56)

Lee, 200132

HF/6-31G(d,p)

55'1

Maheshwary,200136

SPC/E

55'1

Kabrede, 200311

B3LYP/aVDZ

515

Lenz, 200937

TIP3P

55'1

Wales, 199815

B3LYP+D/6-31+G*

515

Do, 201238

TIP4P

44'3

Wales, 199815

B3LYP/6-311++G**

55'1(Pr56)

Lee, 200132

TIP5P

55'1

James, 200534

B3LYP/6-311+G(d,p)

55'1

Bulusu, 200635

HBB2-Pol

44'3

Current

MP2/6-311++G(d,p)

515

Liu, 201142

TIP4P-Ew

44'3'

Kazachenko, 201318

MP2/6-311++G(d,p)// HF/6-311G(d,p)

551

Liu, 201142

TIP4P-2005

44'3'

Kazachenko, 201318

MP2/TZ2P++/TZ2P

55'1 (Pr56)

Lee, 200132

TTM2-F

44'3

Hartke, 200333

MP2/aVDZ

43'4

Bulusu, 200635

TTM2.1-F

44'3'

Hartke, 200333

MP2/aVTZ

43'4

Bulusu, 200635

AMOEBA

515

Kazachenko, 201318

MP2/aVQZ//aVTZ

43'4

Bulusu, 200635

POL3

55'1

Bulusu, 200635

MP2/aV5Z//aVTZ

43'4

Bulusu, 200635

MP4(SDQ)/6311++G(2d,2p)// B3LYP/6-311+G(d,p)

515

Bulusu, 200635

a

- This notation is intended to be consistent with Bulusu et al.35 Some refer to the 55’1 isomer as Pr56.

b

- References to the original MCY,29 SPC/E,22 TIP3P,23 TIP4P,23 TIP5P,26 HBB2-Pol,31 TIP4P-Ew,102 TIP4P2005,24 TTM2-F,103 TTM2.1-F,28 AMOEBA,30 and POL343 definitions.


23


Page 24 of 41

Table 2. Definition of Isomer Families Based on Oxygen Skeleton and Nature of Hydrogen Bonding. Isomer Familya

Rotational Constants [GHz] b A B C

N1-N2-N3-N4c [D-A]-[DD-AA-DA]-[DAA-DDA]-[DDAA-DDDA-DAAA]

N[HBs] d

Counte

43'4A

0.534 0.468 0.406

00-000-46-001

17

6

44'3'A

0.606 0.412 0.371

00-001-44-200

17

18

44'3'B

0.648 0.386 0.353

00-001-44-200

17

12

515A

0.538 0.418 0.399

00-001-55-000

16

17

55'1A

0.495 0.446 0.386

00-001-55-000

16

32

515B

0.557 0.424 0.372

01-001-45-000

15

5

515C

0.574 0.430 0.373

00-002-44-100

16

5

43'4B

0.583 0.396 0.370

00-001-55-000

16

5

551B

0.553 0.407 0.380

00-001-44-200

17

5

443'

0.665 0.395 0.357

00-002-33-300

17

3

44'12a

0.581 0.377 0.329

00-003-33-200

16

7

515D

0.584 0.404 0.360

10-000-54-100

16

2

44'12o

0.544 0.371 0.350

00-003-33-200

16

2

4X

0.601 0.438 0.349

00-002-44-100

16

3

443

0.718 0.370 0.360

00-000-44-300

18

3

a-

- This notation is consistent with that of Bulusu et al.35

b

- The rotational constants of the lowest energy isomer of each family. All members of an isomer family have A, B, C values within 0.03 GHz of that of the lowest energy isomer. c

- Nn = number of n-coordinated waters in the cluster.

d

- Number of hydrogen bonds.

e

- The number of isomers with ∆∆Ee < 2 kcal mol-1 belonging to each isomer family.


24

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


Table 3. CCSD(T)/CBSa Relative Energiesb for the Lowest Energy Members of the 15 Isomer Familiesc CBS

0K

100 K

200 K

300 K

Isomer

∆∆Ee[CCS D(T)]

∆∆E

∆∆H

∆∆G

∆∆H

∆∆G

∆∆H

∆∆G

43'4A-1

0.00

0.09

0.00

0.28

0.00

0.68

0.00

1.19

44'3'A-1

0.19

0.23

0.17

0.35

0.23

0.62

0.28

0.97

44'3'B-1

0.26

0.27

0.22

0.35

0.29

0.58

0.35

0.87

515A-1

0.27

0.00

0.01

0.00

0.09

0.08

0.12

0.23

55'1A-1

0.31

0.01

0.03

0.00

0.13

0.05

0.20

0.16

515B-1

0.42

0.33

0.34

0.37

0.44

0.49

0.51

0.66

515C-1

0.59

0.29

0.30

0.29

0.41

0.37

0.50

0.49

551B-1

0.64

0.36

0.41

0.35

0.56

0.35

0.67

0.38

443'-1

0.79

0.88

0.82

1.00

0.91

1.27

0.99

1.59

43'4B-1

0.82

0.48

0.51

0.40

0.61

0.38

0.67

0.42

44'12a-1

0.84

0.53

0.60

0.37

0.76

0.21

0.85

0.07

515D-1

1.03

0.75

0.81

0.71

0.96

0.66

1.07

0.65

44'12o-1

1.09

0.69

0.79

0.48

0.95

0.24

1.05

0.02

4X-1

1.15

1.27

1.18

1.49

1.27

1.90

1.36

2.35

443-1

1.21

1.50

1.40

1.67

1.44

2.06

1.50

2.51

55'1A-5d

0.56

0.11

0.16

0.04

0.29

0.00

0.37

0.00

a

- RI-MP2/aVDZ//aVDZ *, RI-MP2/aVTZ// aVDZ*, RI-MP2/aVQZ// aVDZ * energies extrapolated using eq 1 and corrected for higher-order electron correlation according to eq. 3.

b

- See Table S3 for absolute (E, H, G) and binding energies (∆E, ∆H, ∆G).

c

- All energies are in kcal mol-1. Global minima shown in bold.

d

- This isomer is included because it is the G0(T) global minimum for 150 K < T

Exploring the Rich Potential Energy Surface of (H2O)11 and Its

Recommend Documents