Structures, energetics and vibrational spectra of (H2O)32 clusters: a journey from model potentials to correlated theory

Empirical model potentials are found to be very useful for generating most competitive minima of large water clusters, whereas correlated (e.g. second order-Møller–Plesset perturbation (MP2) theory or higher) calculations are necessary for predicting their accurate energetics and vibrational features. The present study reports the structures and energetics of (H2O)32 clusters at MP2 level using aug-cc-pvDZ basis set, starting with low-lying structures generated from model potentials. Such high-end and accurate calculations are made feasible by the cost-effective fragment-based molecular tailoring approach (MTA) in conjunction with the grafting procedure. The latter is found to yield electronic energies with a sub-millihartree accuracy with reference to their full calculation counterparts. The vibrational spectra of nine low-lying (H2O)32 isomers are obtained from the corresponding MTA-based Hessian matrix. All these low-lying isomers show almost similar spectral features, which are in fair agreement with the experiment. The experimental spectrum of (H2O)32 is thus better understood from the vibrational features of this set of very closely spaced isomers. The present case study of (H2O)32 clearly demonstrates the efficacy in obtaining accurate structures, energetics and spectra at correlated level of theory by combining model potential-based structures with fragmentation methods.


Introduction
Being the significant chemical constituent of our planet's surface, water and hence water clusters have received more scientific and technological attention than any other substance [1]. Vast amount of work has already been devoted towards the understanding and modelling the role of water clusters in several chemical and biological processes [2][3][4][5]. Numerous fundamental questions arise for the understanding of the unusual properties of the extended hydrogen bond networks in water clusters in both pure and solvent form. For example, to quote Buck et al. [6] verbatim: (1) At what cluster size, typical structures of the bulk-like completely four-coordinated molecules appear? (2) In which size region, the transition to completely crystalline behaviour does occur? Additionally, there are still unresolved mysteries of the extra stability of hydrogen bond networks, appearance of internally solvated water molecules etc. in aqueous chemistry. Thus, the structure and dynamics of water clusters are still the crux of many theoretical and experimental investigations.
Recent developments in the high-resolution and sizeselective spectroscopic techniques such as infra-red (IR), matrix isolation IR spectroscopy and broadband rotational spectroscopy have accelerated the understanding and unearthing of the structural information and size-dependent * Corresponding author. Email: gadre@iitk.ac.in properties of water clusters [6][7][8][9][10][11][12][13][14]. Especially, the sizeselective experiments enable accurate exploration of the microscopic behaviour of water. Several experimental studies [8][9][10][11][12][13][14] have been carried out for probing the molecular structures and properties of water clusters of increasing size. These experimental spectra point to the existence of several energetically close water cluster isomers and an explanation in terms of a unique isomer is often not possible. The proper separation and characterisation of individual structures within experimental conditions is still an arduous task, owing to the rapid increase in the possibility of the number of iso-energetic isomers with the size of cluster. For example [15], the potential energy surface of (H 2 O) 25 cluster has around 2.6 × 10 7 hydrogen bond topologies whereas for (H 2 O) 32 , this number is approximately equal to 1.8 × 10 9 . In view of this, many empirical model interaction potentials have been developed and applied for searching favourable hydrogen-bonded networks of water clusters.
In the past few years, several investigations have appeared on medium-sized water clusters (n > 15), using different optimisation techniques and potential functions [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. In some of the earlier works, Wales et al. [16] and Tsai et al. [17] employed the empirical rigid-body TIP4P potential for searching the global minimum structures C 2015 Taylor & Francis of water clusters (H 2 O) n for n ≤ 21. Some other works include the rigid body diffusion Monte Carlo calculation by Buch and co-workers [18], empirical point-charge models for water due to Kabrede and Hentschke [19,20], TIP5P and polarisable model potential by Wales and co-workers [21,22], path-integral molecular dynamics and centroid molecular dynamics methods by Paesani et al. [23], and full-dimensional ab initio potential energy and dipole moment surfaces based on many body potentials [24]. Recently, simulated annealing with effective fragment potential [26] and Monte Carlo algorithm-based basin paving/hopping and temperature basin paving method [27,28] along with minima hopping method [29][30][31] have been widely used for locating the lowest energy structures for the medium-sized water clusters. In an extensive series of works, Kazachenko and Thakkar [32] reported new putative global minimum geometries of water clusters containing up to 55 monomers employing five model potentials, viz. TIP4P, TIP4P-Ew, TIP4P/2005, TTM2.1_F and AMOEBA. In spite of the ability to identify and to introduce stochastics in the generation of water clusters of well-characterised families of favourable hydrogen bond networks, the potential-based techniques have their own limitations. It has been recognised that the ad hoc parameterisation of empirical potentials for water results into an incomplete treatment of many-body interactions. This may, in turn, lead to an inconsistent prediction of the energetics and the relative rank ordering of various cluster isomers. Therefore, obtaining the structures of low-lying candidates using reliable ab initio methods is necessary for assessing the accuracy of empirical force fields for water.
The Hartree-Fock (HF) and density functional theories have been widely employed for studying medium-as well as large-sized water clusters [15,[33][34][35][36]. However, for an accurate description of cooperativity and hydrogen bond energy, methods that include electron correlation, such as the second-order Møller-Plesset perturbation (MP2) theory or coupled cluster with singles, doubles and perturbative triple replacements (CCSD(T)) are recommended. Especially, due to effective parallelisation and low scaling, viz. O(N 5 ) (N, being the number of basis functions), MP2 method is still widely used for studying water clusters. Furthermore, this is the least expensive one among the well-known post-HF methods that includes electron correlation. In the recent years, the availability of highly parallel codes as well as high-end computational architectures have made such calculations feasible for medium-sized (15 < n < 25) water clusters. Xantheas and co-workers [37][38][39] reported a series of such high-level computational studies on water clusters in this size-range. For instance, in 2009, the CCSD(T)/aug-cc-pvTZ single point energy full calculation (FC) [40] for the (H 2 O) 24 cage was performed on all 223,200 processors of Oak Ridge National Laboratory's Jaguar Supercomputer with a usage of a total of 150 Terabytes of memory and a demonstrated sustained performance of 1.39 Petaflops/s in double precision. However, due to the requirement of such heavy computational resources as well as the existence of enormous number of nearly iso-energetic isomers, finding the global minimum structures of large water clusters employing correlated ab initio methods remains out of reach for most of the workers in this area.
In this context, the introduction of fragment-based approaches [41] has proven to be a boon for making such calculations feasible, especially for the geometry optimisation and vibrational spectra of large-sized water clusters. In the recent years, a significant amount of effort has been devoted towards the development of efficient algorithms enabling a reduction in the computational cost for largesized systems. Some notable studies include the works by Gadre and co-workers [42,43] employing molecular tailoring approach (MTA), Li and co-workers [44,45] through generalised energy-based fragmentation (GEBF) method, Federov and Kitaura [46] using fragment molecular orbital (FMO) methods and Gordon and co-workers [47] employing both FMO and systematic molecular fragmentation techniques. Other significant works in this domain include the electrostatically embedded many-body expansion by Truhlar and co-workers [48], multi-centred integrated QM:QM method due to Tschumper and co-workers [49] and recent ones by Raghavachari and co-workers employing dimer-based two-body interaction model [50]. Among the fragment-based methods, mainly MTA, FMO and GEBF have been widely applied to water clusters of various sizes for the estimation of the relative stability and structural parameters of the extensive hydrogen bond networks. Especially, MTA assisted by the newly proposed grafting procedure [42] (a two-level calculation technique) has been extensively used for geometry optimisation followed by accurate energetics and vibrational spectra of medium-sized water clusters during the past five years.
In summary, the empirical potentials in conjunction with stochastic methods such as basin paving, genetic algorithms, etc. are suitable for searching low-energy isomers. On the other hand, the utilisation of the fragmentationbased techniques such as MTA, is crucial for carrying out computationally expensive ab initio calculations with minimal off-the-shelf hardware. The judicious combination of these two can be profitably used for the high-end correlated calculations of large water clusters. In view of the current experimental interest in vibrational spectra of water clusters [6], the present study reports the MP2/aug-cc-pvDZ (MP2/aDZ) level geometry optimisation and vibrational spectra of (H 2 O) 32 clusters employing the most recent version of MTA.

Methodology and computational details
The present section gives a brief introduction of the MTA method and the grafting procedure for accurate energetics and vibrational calculations followed by the detailed computational procedure employed for the calculations.
As indicated in the Introduction, MTA [51][52][53], the in-house fragmentation-based method from the authors' group, is mainly meant for the high-level quantum chemical calculations of spatially extended molecular systems. Within this, a large molecular system under study is cut into smaller fragments on which ab initio calculations are carried out. The estimation of an appropriate electronic property, P, of the large system is achieved by patching the results from individual fragments by means of the set inclusion-exclusion (cf. Equation (1)) principle [54,55]. In other words, within MTA framework, the value of any property, P, of a large system can be approximated as the sum (and difference) of the corresponding values of the fragments. Thus far, P has been chosen to be the electronic energy, energy gradients, Hessian elements as well as tensors, such as the dipole derivatives and polarisability derivatives.
where P F i denotes the value of P in the ith fragment, P F i ∩F j , stands for the value of the binary overlap between the ith and jth fragments and so on. K represents the order of overlap between the fragments, e.g., K = 1 for binary overlaps and K = 2 for ternary overlaps of fragments. The details of the working principle and the controlling parameters, viz. R-goodness, average size of fragments, etc., needed for the MTA-based computation are described in detail in the previous publications from the authors' group, and further summarised in Reference 52. The present form of the MTA-Plugin enables the geometry optimisation and calculations of Hessian matrix elements, dipole derivatives, polarisability derivatives and molecular orbital energies of large molecular clusters/systems employing offthe-shelf hardware [56][57][58][59][60]. Due to its demonstrated ability to work with Gaussian, GAMESS and NWChem software packages at the back-end, MTA has emerged as a powerful tool for ab initio treatment of large molecular systems. With increasing size of the systems, and especially for MP2 level geometry optimisation, MTA shows a significant time advantage vis-à-vis the corresponding FCs even with limited hardware.
It has been demonstrated in our earlier works that the MTA-based geometry optimisation is able to provide good geometries for large molecular clusters. For instance, in a recent study on (H 2 O) 25 clusters from the authors' group [43], it has been found that the maximum, average and RMS differences between the full and MTA optimised geometries amount to 0.06, 0.02 and 0.02Å for the hydrogen bond H···O lengths and to 0.006, 0.001 and 0.002Å for the covalent O-H bond lengths, respectively. However, due to the approximations involved, the electronic properties derived from MTA always carry some error with respect to their FC counterparts. In the last few years, significant effort has been made for minimising the errors through a grafting procedure [42,43,[61][62][63]. This hinges on the observation of the basis set independent nature of the errors between the FC and MTA-based calculations [64]. The correct estimation of the molecular energies is effected by adding a correction (termed as 'Grafting') estimated with a smaller basis (sb) to that of the computed quantity at a larger basis (lb) [42,43,[61][62][63]. Such a correction for large molecular systems is also recently re-examined by Raghavachari and co-workers [65,66].
Simulation techniques based on molecular dynamics studies incorporating basin-hopping, genetic algorithm and simulated annealing have been crucial for introducing stochastics as well as identifying the energetically favourable isomers of water clusters. The current study employs starting geometries that are derived from such empirical potentials. The initial structures are taken from the existing literature, viz. from the works due to Hartke and co-workers [6], Kazachenko and Thakkar [32], and Lenz et al. [33], generated by using a variety of empirical model potentials and/or low-level ab initio calculations. Out of a total of 45 starting geometries used in the current study, 38 are derived from the global geometry optimisation using genetic algorithm by Hartke and co-workers [6] whereas other structures [32] are the putative global minima for various intermolecular potential energy models. These models include, three empirical, pairwise-additive potential energy surfaces, viz. TIP4P, TIP4P-Ew and TIP4P/2005 and other two, TTM2.1_F and AMOEBA, are polarisable models including non-additive inductive effects.
All these isomers are subjected to MTA-based MP2/aDZ level geometry optimisation followed by the MTA-based single point energy calculations at MP2/aTZ level of theory. This is followed by the accurate estimation of the electronic energies at both MP2/aDZ and MP2/aTZ levels employing the grafting procedure as shown in Equations (2) and (3), respectively. For the grafting method at MP2/aDZ level, 6-31 + G * is chosen as the sb whereas aDZ basis is taken as lb. Later, for the purpose of testing and benchmarking, the grafted values at this level are then compared with their FC counterparts. Similarly, for the estimation of accurate energetics at MP2/aTZ level of theory, aDZ and aTZ are chosen as the sb and lb, respectively. Equation (2) hinges on the fact that the major contribution to the error in energy comes from the HF level rather than from the second-order perturbative correction i.e. E2. The grafting correction at MP2/aTZ level is performed via Equation (3) assuming that the error within MTA energy is approximately constant for similar quality basis sets [43]. The details of the grafting procedure for the estimation of accurate energies are given elsewhere [43,62,63].
where E GMTA@MP2/aDZ and E GMTA@MP2/aTZ stand for the grafted energies at MP2/aDZ and MP2/aTZ levels of theory, respectively. E FC@HF/aDZ , E FC@MP2/aDZ and E MTA@MP2/aDZ denote the corresponding electronic energies at aDZ basis set. On the other hand, E2 FC@MP2/6-31 + G * , E2 MTA@MP2/6-31 + G * and E2 MTA@MP2/aDZ represent the second-order perturbative correction (correlated energy) part from the corresponding MP2 level calculations. It may be pointed out that Equation (2) requires a FC at HF level using lb, which could be cumbersome for a large cluster. Thus, the use of Equation (3), being devoid of a FC at HF/lb, is more useful for enhancing the balance between the computational cost and accurarcy [43]. The correct estimation of electronic energies even at MP2/aQZ, MP2/a5Z or higher levels of theory can be easily computed within Equation (3), employing aDZ basis as sb set. Later, the cluster binding energies at both MP2/aDZ and MP2/aTZ levels are calculated as shown in Equation (4).
The distinct hydrogen bond networks of the isomers are analysed through the connectivity, type and number of water monomers as well as through the number of hydrogen bonds present in the isomers. The geometry optimisation is followed by the estimation of the Hessian matrix elements and dipole derivative tensors, leading to the vibrational spectrum of an isomer. In a very recent study [62], the present authors have reported MP2 level harmonic vibrational spectra of neutral and charged (both cationic and anionic) water clusters, employing MTA. This procedure is applied to the nine low-lying isomers of (H 2 O) 32 (W32) cluster at MP2/aDZ level of theory. The distinct features of the calculated vibrational spectra of the isomers are compared to the recent experimental spectrum reported by Hartke and coworkers [6]. All the MTA-based calculations are performed using Gaussian 09 package [67] at the back-end.

Results and discussion
Due to the extensive compactness and randomness in the orientation of the hydrogen bond networks, the fragmentation of even a medium-sized water cluster is still a challenging task. Therefore, the automatic fragmentation procedure developed in our group [52] is employed for carrying out MTA-based geometry optimisation. For assuring uniformity of treatment, almost similar (both qualitatively and quantitatively) set of fragments is used for all of the isomers. These schemes yield close to 25 main fragments with typical size in between 10 and 12 water monomers for all of the isomers. As an example, a brief account of the fragments used for the W32-TIP4P-Ew isomer, in terms of the monomer indices (only for main fragments), is included in the supplementary material.
All of the 45 isomers that are derived from the global geometry optimisation using genetic algorithm (38 in number) by Hartke and co-workers [6] and the corresponding five intermolecular potential energy models (TIP4P, TIP4P-Ew, TIP4P/2005, TTM2.1_F and AMOEBA) by Kazachenko and Thakkar [32] are optimised at MP2/aDZ level of theory employing MTA. This is followed by a FC of the single point energy at MP2/aDZ level of theory at the corresponding MTA-optimised geometry. The detailed structures and energetics of nine low-lying MTA-optimised geometries are presented here. For the sake of uniformity, the nomenclature of the isomers used in Ref. [32] is followed in the present work. For instance, W32-AMOEBA and W32-TIP4P refer to the global minima geometries [32] obtained employing the AMOEBA and TIP4P water models. The corresponding MP2/aDZ-level MTA-optimised geometries for the above five isomers are displayed in Figure 1. In a similar manner, W32-G1 and W32-G2 are derived from the global geometry optimisation due to Buck et al. [6] whereas other two isomers W32-T1 and W32-T7 are derived from some other potential models (see Figure 1). The cartesian atomic coordinates (inÅ) of the MTA-optimised geometries of all these nine isomers are given in the supplementary information file. Table 1 lists out the interaction energies as reported by Thakkar et al. [32] using model potentials along with the corresponding MP2/aDZ-level energy from the FCs at the initial geometries derived directly from these potentials. How does the relative energy rank ordering of these isomers change upon MTA-based geometry optimisation? For answering this query, the single-point FC energies for the MTA-optimised geometries are also included in Table 1. A change in the relative ab initio energy rank ordering of the isomers is noticed upon the geometry optimisation although the 'best' two isomers (W32-AMOEBA and W32-TTM2.1_F) still remain at top. The W32-TIP4P/2005 global minimum isomer (numerically largest [32] interaction energy) being fourth before geometry optimisation becomes the third most stable after optimisation at MP2/aDZ level of theory. Similarly, the W32-TIP4P and W32-TIP4P-Ew isomers are seen to switch their order from fifth and third to fourth and fifth, respectively, after effecting the optimisation process. In view of this, it can be concluded that the relative energy rank ordering of the isomers at the potential model geometries, in general changes upon the geometry optimisation at ab initio level of theory. Thus, for the ac- curate energetics and benchmarking with the experimental observations, correlated level investigations are warranted for the structures derived from various model potentials.
Later, the grafting procedure developed in our group [42,43,61] has been implemented for obtaining more reliable energetics of the isomers within MTA. Table 2 reports the MTA-and grafted MTA-based (calculated via Equation (2)) energies for all the nine isomers at MP2/aDZ level of theory. It may be seen that the MTA energies at MP2/aDZ level of theory before the grafting procedure are in error between 7 and 35 mH from the respective FC values. Additionally, the MTA-based energy rank ordering among the isomers is quite different from their FC counterparts. However, upon grafting, the errors in MTA are reduced to a maximum of 1 mH (W32-G1) and the relative energy ordering of the isomers is restored with respect  to the FCs. The most stable isomer [W32-AMOEBA] at MP2/aDZ level of theory resembles the global minimum geometry obtained from the AMOEBA potential model by Kazachenko and Thakkar [32] whereas the isomers W32-TTM2.1_F, W32-TIP4P/2005, W32-TIP4P-Ew and W32-TIP4P correspond to the global minimum geometries from the respective models. On the other hand, isomers W32-G1 and W32-G2 correspond to two isomers from the global geometry optimisation [6]. Within the potentialbased global geometry optimisation, W32-G2 turns out to be the most stable isomer [6] whereas after MP2/aDZ level geometry optimisation, it turns out to be energetically unfavourable as compared to other isomers (cf. Table 2). The optimised isomers W32-T1 and W32-T7 are found to be more stable than W32-TIP4P/2005, W32-TIP4P-Ew and W32-TIP4P isomers. Unlike the MP2/aDZ-level grafted-MTA energy, the grafting correction at MP2/aTZ level is performed using Equation (3) assuming that the error within MTA energy is approximately constant for similar quality basis sets. Thus, the difference in the energy between the MTA and FC values at MP2/aDZ level is grafted [43] on to the gross MTA values at MP2/aTZ level in order to obtain the corrected energies at the MP2/aTZ level of theory. Due to the unavailability of adequate computational power, the FC at MP2/aTZ level of theory (total 2944 basis functions for W32) is not possible. However, in accordance with the earlier work [43] from our group, the FC counterparts and grafted-MTA values at MP2/aTZ level of theory are expected to be in good agreement with each other. At this level of theory, a switching of the relative rank order between isomers W32-TIP4P/2005, W32-G1 and W32-TIP4P is noticed. However, the top four isomers at MP2/aDZ level retain their place at MP2/aTZ level of theory. Isomers W32-T1 and W32-T7 turn out to be practically isoenergetic at MP2/aTZ level of theory (cf. Table 2). Overall, the energy rank ordering between the grafted values at both aDZ and aTZ basis sets is found to match well. The isomers W32-AMOEBA and W32-TIP4P-Ew retain their positions at the top and bottom of stability chart at both levels (see Table 2).
Energies of all the isomers fall within 2-3 kcal/mol of each other (see Table 3). This indeed brings out the ability of the potential models for the generation of distinct and yet energetically very close isomer structures. On the other hand, the energy window of the interaction energies of the isomers (cf. Table 1) corresponding to the global minima of the potential models [32] is seen to be much wider as compared to that of the binding energies at MP2/aDZ and MP2/aTZ levels of theory (see Table 3). In a similar fashion, the structures derived by Buck et al. [6] from the global optimisation employing the genetic algorithm followed by the MTA-based geometry optimisation at MP2/aDZ level of theory lie within 3-4 mH (∼2 kcal/mol) of each other although they are energetically less favourable as compared to W32-AMOEBA, W32-TTM2.1_F and W32-TIP4P/2005 isomers [32]. On the other hand, W32-G1 isomer is energetically more favourable than the TIP4P isomer at this level. Besides these, the MP2/aDZ level geometry optimisation of the geometries obtained from the global optimisation procedure produces around 10 new isomers (less favourable) that are placed within an energy window of 1.2 kcal/mol from the W32-TIP4P-Ew isomer. Furthermore, MTA-based geometry optimisation is able to generate three more isomers that are almost isoenegetic to W32-G2 isomer. However, due to limitation of space, all those isomers are not included in the present study. In summary, the potential model-based studies are crucial for generating most competitive minima with extensive hydrogen networks. The MP2 level computation gives the accurate energetics of the clusters. All the isomers have somewhat similar networks containing 56 hydrogen bonds each, with two internally solvated water molecules (see Figure 1). Table 1 reports the total number of hydrogen bonds along with the type and number of water monomers in all the isomers derived from the model potentials [32]. Isomers W32-AMOEBA, W32-T1, W32-TIP4P/2005, W32-G1 and W32-G2 contain 8 DAA-, 8 DDA-and 16 DDAA-type (D = proton Donor, A = proton Acceptor) water molecules, i.e. 16 monomers connected via four hydrogen bonds and other 16 monomers with three hydrogen bonds to nearest neighbours. On the other hand, isomers W32-TTM2.1_F, W32-T7, W32-TIP4P and W32-TIP4P-Ew have one bi-coordinated (DA-type) water molecule in the hydrogen bond networks. These three isomers are endowed with 17 DDAA-, 7 DDA-and 7 DAAtype of water molecules. Thus, all the isomers have eight dangling -OH (not participating in hydrogen bond formation) groups in their networks. With very little distinguishing feature of these networks, the prediction of the relative stability is quite difficult. Furthermore, the average similarity index [32] among the isomers is found to be very close to 1. It may hence be anticipated that the characterisation and separation of such isomers is prohibitively difficult even with the advanced modern spectroscopic tools.
The MTA-based vibrational spectra [56,62] of all the isomers (expect W32-T7) show no imaginary frequencies at MP2/aDZ level of theory, confirming their minimal nature. However, the only imaginary frequency in W32-T7 isomer corresponds to 10 cm -1 with a negligible intensity of 0.029 KM/mol and hence may be ignored. We would like to point out that this is a first-ever study of the MP2 level vibrational spectra for such large water clusters. Some earlier examples of such large calculations include the MP2/aDZ-level vibrational analysis of (H 2 O) 20 clusters due to Pradzynski et al. [13] and (H 2 O) 25 clusters by the present authors [43].
The vibrational spectrum of water clusters is mainly divided into three distinct regions: (1) libration modes (intermolecular region less than <1100 cm −1 ) and the other two being the intramolecular regions associated with the (2) HOH bending (1600-1750 cm −1 ) and (3) O-H stretching vibrations (>3000 cm −1 ). The latter region displays distinctively different features for the HB networks and is regarded as the 'fingerprint' region for water clusters [13]. A scaling factor of 0.959 is used for scaling of the calculated frequencies at MP2/aDZ level of theory. Figure 2 displays the MTA-based convoluted spectra (Lorentzian function with full-width at half-maximum (FWHM) of 10 cm −1 ) of W32-AMOEBA, W32-TTM2.1_F and W32-T1 isomers at MP2/aDZ level of theory. The vibrational spectra of other six isomers are included in the supplementary information file (see Figures S1-S6). In the O-H stretching region, the peaks beyond 3700 cm −1 appear due to free O-H stretching vibrations of 3-coordinated (DAA-) and 2-coordinated (DA-) type water molecules, the latter being more blue shifted as compared to the earlier type. The peaks below 3600 cm −1 are attributed to the hydrogen bond-bonded O-H stretching vibrations. Out of these, O-H vibrations associated DDA-type water molecules show peaks in 3450-3600 cm −1 regions whereas further red-shift (generally more intense) is observed in those associated with DAA-and DDAA-type of water molecules. Due to the presence of eight dangling O-H groups, the spectra of all the isomers show a peak in the region from 3700 to 3725 cm −1 (cf. Figure 2). The most intense peaks are attributed to those of the hydrogen-bonded O-H groups associated mainly with DAA-and DDAA-type water molecules.
In a very recent study, Buck et al. [6] have simulated the vibrational spectra of W32 clusters employing the global geometry optimisation and photoionising sodium doped clusters close to the ionisation threshold. From this study [6], no spectrum corresponding to a single-structure from the top 38 unique structures (generated from global optimisation) was able to provide a good approximation to the experimental spectrum. All the structures were found to have two interior water molecules within a potato-shaped cage; a cuboid protrusion is located edge-down above the midpoint of an imaginary line joining the two interior water molecules and are found to share same general structural features, but differ in finer details. Very similar to the observation to Hartke and co-workers [6], all the isomers examined in the present study are structurally similar to one another with two internally solvated water molecules and with 56 hydrogen bonds in each.
A qualitative comparison of the calculated spectra with the experimental observed spectrum follows. The experimental spectrum provided in Reference 6 is indeed very broad and does not show very well-resolved peaks. It is very difficult to see if any of the spectra in Figure 2 has/have a similar intensity distribution to the experimental data. Simulation of the theoretical spectra with a relatively larger line width is necessary to have a better representation of these calculated spectra. Figure 3 displays the MP2/aDZlevel vibrational spectrum for the W32-AMOEBA, W32-TTM2.1_F and W32-T1 isomers in the region from 2800 to 3800 cm −1 with a trial FWHM value of 40 cm −1 chosen for this purpose. The experimental spectrum due to Buck et al. [6] shows some intense peaks in the region from 3300 to 3400 cm -1 and from 3420 to 3550 cm -1 region, some other peaks around 3150 cm -1 with moderate intensity along with a less intense region around 3700-3720 cm -1 . Our calculated spectra of all the isomers of course are endowed with features similar to those displayed by the experimental one. For example, the unconvoluted spectrum of W32-AMOEBA isomer has its most intense peaks at 3319, 3241, 3374, 3178, 3425 and 3465 cm −1 . Some other peaks with relatively less absolute intensities near 3260, 3118, 3120 and 3172 cm −1 are also noticed for this isomer. Similarly, the unconvoluted spectrum of W32-TTM2.1_F isomer exhibits the most intense peaks around 3500, 3420, 3400, 3465, 3425 and 3135 cm -1 whereas the most intense peaks for the W32-T1 isomer correspond to 3478, 3541, 3308, 3581 and 3209 cm -1 , respectively. Furthermore, all the isomers show distinct peaks between 3700 and 3720 cm -1 which are attributed to the dangling O-H groups. Similar features are also noticed after the convolution of the calculated spectra (cf. Figures 2 and 3). The convoluted spectrum of W32-AMOEBA isomer shows five intense regions around 3175, 3240, 3317, 3420 and 3550 cm −1 , respectively. Similarly, the W32-TTM2.1_F and W32-T1 isomers have these regions, respectively, around 3192, 3283, 3361 and 3515 cm −1 and 3167, 3331, 3397 and 3483 cm −1 . The most intense regions in the experimental spectrum reported by Buck et al. [6] happen to be approximately around 3290, 3370, 3390, 3450 and 3540 cm −1 . Furthermore, the O-H stretching region of the calculated spectra spread over a region from 3000 to 3730 cm -1 as compared to 3000-3750 cm -1 due to the observed spectrum. All these features are found to be in a qualitative agreement with the experimental observed spectrum, although a single isomer is inadequate for interpreting the complete experimental spectrum. Overall, after the convolution, the vibrational spectra of the isomers are found to be qualitatively similar to each other and are found to share common features (see the supplementary material for vibrational spectra of other isomers). Thus, it is apparent that the possibility to assign only one or two structures for explaining the features of the experimental spectrum would diminish for larger clusters [6].

Concluding remarks
The number of local minima in water clusters, with very diverse hydrogen bond networks, increases very sharply with the size of the cluster. The combination of stochastic optimisation methods with different model potentials is hence crucial for the exhaustive search of the minimum-energy structures of these clusters. However, the ad hoc parameterisation of the empirical interaction potentials may lead to an inconsistent prediction of the energetics and hence the resulting relative rank ordering of various cluster isomers. That is where high-level ab initio methods are indispensable for obtaining the structures and accurate energetics of low-lying isomers. From the current study, it is noticed that the energy rank ordering of the potential-derived isomers of (H 2 O) 32 cluster indeed changes upon the MP2-level geometry optimisation. In view of the requirement of large computational resources, finding the global minimum-energy conformers of such a large-sized cluster is out of reach for even the least expensive correlated level of theory, viz. MP2. In this context, the fragmentation-based methods such as MTA are immensely beneficial. The present study has reported the accurate structures, energetics and vibrational spectra of various isomers of (H 2 O) 32 cluster (derived from a variety of potential models) at MP2/aug-cc-pvDZ and MP2/aug-cc-pvTZ levels of theory using MTA.
The diversity in the hydrogen bond networks of (H 2 O) 32 cluster is brought in by introducing various potential models. The starting structures are taken from the global minimum geometries employing AMOEBA, TTM2.1_F, TIP4P/2005, TIP4P and TIP4P-Ew potential models as well as from the global geometry optimisation using genetic algorithm. A two-level calculation technique (termed as the grafting procedure) is found to yield reliable electronic energies, with sub-millihartree accuracy, at MP2/aug-cc-pvDZ level of theory while retaining the relative energy rank ordering of the isomers vis-à-vis the corresponding FC results. A total of nine low-energy conformations (each containing 56 hydrogen bonds with two internally solvated water molecules) within an energy window of 3 kcal/mol at MP2/aug-cc-pvDZ level are subjected to MTA-based vibrational analysis, the results of which are compared to the recent experimental data. Indeed the present work is a firstever study of the MP2 level vibrational spectra for such large water clusters. For all of the nine isomers, the MTAbased vibrational spectra in the O-H stretching are found to be in a fairly good agreement with the observed spectrum. The experimental spectrum can be better understood from the vibrational features of a set of very closely placed isomers rather than a single isomer.
In summary, due to the appreciable savings in computational time and resources by MTA, the present methodology is expected to serve as a useful tool for such art-ofthe-possible calculations, which are otherwise prohibitively difficult even on huge computational hardware. It is hoped that such investigations on the stability and energetics of large-sized clusters would ultimately provide a link to the corresponding bulk behaviour.