University of Birmingham Structures and energy landscapes of hydrated sulfate clusters

: The sulfate ion is the most kosmotropic member of the Hofmeister series, but the chemical origins of this e ﬀ ect are unclear. We present a global optimization and energy landscape mapping study of microhydrated sulfate ions, SO 42 − (H 2 O) n , in the size range 3 ≤ n ≤ 50. The clusters are modeled using a rigid-body empirical potential and optimized using basin-hopping Monte Carlo in conjunction with a move set including cycle inversions to explore hydrogen bond topologies. For clusters containing a few water molecules ( n ≤ 6) we are able to reproduce ab initio global minima, either as global minima of the empirical potential, or as low-energy isomers. This result justi ﬁ es applications to larger systems. Experimental studies have shown that dangling hydroxyl groups are present on the surfaces of pure water clusters, but absent in hydrated sulfate clusters up to n ≈ 43. Our global optimization results agree with this observation, with dangling hydroxyl groups absent from the low-lying minima of small clusters, but competitive in larger clusters.


INTRODUCTION
The Hofmeister series ranks anions and cations according to their relative effectiveness at influencing a number of interesting phenomena in physical chemistry. 1,2 Within the series, ions can be characterized as being either kosmotropic (structure making) or chaotropic (structure breaking). Kosmotropes increase the surface tension of liquids 3 and the stability of proteins, 4 and decrease the solubility of hydrophobic particles 5 and the denaturation of proteins 6 (to name but a few properties), while chaotropes behave in the opposite manner. 7 However, despite the ubiquity of Hofmeister effects, their underlying chemical origins remain unclear, with theory and experiment seemingly providing evidence both for and against 3,8−10 the long-range structuring of solvent molecules around Hofmeister ions. The sulfate dianion sits at the far kosmotropic end of the Hofmeister series, and its solvation has been studied extensively by experiment 11−13 and theory, 14−18 due to its relevance in understanding the ordering behavior, as well as atmospheric chemistry, 19 and a number of industrial processes.
Infrared photodissociation (IRPD) spectroscopy of gas-phase hydrated sulfate clusters, SO 4 2− (H 2 O) n , in the size range 37 ≤ n ≤ 80 has revealed a size dependence of the appearance of OH bonds protruding from the surface of the clusters. The experimental signature, a small peak in the IRPD spectra around 3700 cm −1 , is only observed for clusters with n > 43, indicating that, for smaller clusters, all OH bonds participate in an OH···O hydrogen bond. 11 This result contrasts with the surface of bulk water 20 and small water clusters, 21 where many water molecules are able to orient themselves to expose a dangling OH bond, and suggests that the sulfate ion has a significant effect on the surrounding water structure, even beyond the second solvation shell. Here we present a systematic computational study of SO 4 2− (H 2 O) n , in the range 3 ≤ n ≤ 50, providing microscopic insight into the hydrogen bonding network of water molecules around the sulfate dianion and into the size-dependent emergence of dangling OH bonds.

METHODOLOGY
2.1. Modeling the System. Due to the relatively large system sizes under consideration, we have modeled the hydrated sulfate clusters using an empirical potential with bond lengths and bond angles held rigid. The water molecules are represented by the four-site rigid-body TIP4P water potential, 22 which includes the two hydrogen atoms, the oxygen atom, and a lone pair site.
where i and j are the atom indices, q i is the partial charge on atom i, r ij is the distance between nonbonded atoms i and j, and σ and ε are Lennard-Jones parameters. Lennard-Jones parameters for interactions between unlike atom types were calculated using the Lorentz−Berthelot mixing rules. All energies in the paper are given in units of kcal mol −1 .
The sulfate potential that we have used was derived from Møller−Plesset MP4SDTQ level calculations by McCammon et al. and has been shown to reproduce experimental solution data. 14 Previous calculations of the microsolvation of the sulfate ion at the MP2 level also agree with the MP4 results. 17 McCammon et al. suggested two sulfate potentials, with the same Lennard-Jones parameters but different partial charges on the sulfur and oxygen atoms. In order to choose between these parameter sets, we performed short basin-hopping runs for SO 4 2− (H 2 O) n in the size range 3 ≤ n ≤ 7 for both potentials. The 40 lowest-energy isomers for each n were then reoptimized using density functional theory (DFT) level with the B3LYP exchange-correlation functional and a 6-311++G** basis set, as implemented within the NWChem package. 35 Little difference in isomer ordering was found when comparing the two potentials, so we chose to use the first of the two parameter sets suggested by McCammon et al. (Table 2) because it had a slightly better Spearman rank correlation with the DFT calculations.
As both the water molecules and the sulfate anion are treated as rigid bodies, we have used angle-axis variables to describe the rigid-body rotational coordinates. This representation is preferred to Euler angles, since it avoids the problems in geometry optimization that arise from singularities when these coordinates are used. 36 (H 2 O) n were located using the basin-hopping algorithm 39−41 implemented in the GMIN 42 and pele 43 software packages. For each cluster in the range 3 ≤ n ≤ 50, eight searches were conducted starting from different random geometries. The basin-hopping algorithm works as follows: (1) With an initial potential energy minimum as a starting point, the geometry is perturbed. (2) A local minimization of the geometry is performed with respect to the potential energy. (3) The new minimum is accepted or rejected according to a Metropolis criterion. 44 The geometry was perturbed according to three distinct move classes: rigid-body translations of the center-of-mass of the water molecules and sulfate ion, rigid-body rotations about the oxygen atom (for the water molecules) or sulfur atom (for the sulfate ion), and cycle inversion moves (equivalent to Takeuchi's closed chain perturbations 24 ). A cycle in a directed network is a closed loop of edges with the direction of each edge pointing the same way around the loop. 45 Such a loop can be inverted by pointing edges in the opposite direction. The cycle inversion move class searches for simple cycles (those in which no node appears twice) in the directed network of water−water hydrogen bonds using the simple_cycles method in networkx, 46 a Python library for the creation, manipulation, and study of networks. In a cycle inversion move, a cycle is chosen at random from the hydrogen bond network and inverted. If no cycles are present, a rigid-body translation or rotation is performed instead. A water molecule i is considered to donate a hydrogen bond to water molecule j if their oxygen− oxygen distance is less than 3.50 Å and the OH bond axis of molecule i is within 30°of the displacement vector between the oxygen atoms in i and j. Some structures with positively charged atoms close to the TIP4P lone pair underwent cold fusion. These structures were discarded from the basin-hopping search.
Global optimization studies of pure water clusters have suggested that the use of block moves, where one type of move class is used exclusively for a set number of steps before switching to another, is an effective strategy for finding lowenergy structures. 30 In this study, we used blocks of 100 steps for each type of move. To ensure reasonable parameters for searching across a variety of n, we tuned the size of the translation step and the temperature used in the Metropolis acceptance criterion. This tuning was achieved by running 160 000 basin-hopping steps for each parameter combination for every tenth cluster size starting from n = 5 (i.e., n = 5, 15, 25, ...). The combination of temperature and translational step size that consistently found the lowest-energy minima was considered to be optimal, and was assumed to be transferable to neighboring sizes.
2.3. Energy Landscape Mapping. Transition states connecting minima on the potential energy surface of small hydrated sulfate clusters were located using the pele software package. 43 A doubly nudged elastic band method 47 was used to find transition state candidates, with interpolation between end points using a linear interpolation of the translational coordinates of the rigid-body molecules and spherical linear quaternion interpolation of the rotational coordinates. 48 Transition state candidates were further optimized with hybrid eigenvector-following. 49,50 Initially, connection attempts were made between pairs of minima chosen at random, or between a randomly chosen minimum and the global minimum. Later, we adopted a strategy to remove artificial frustration in the network by attempting to connect minima to the global minimum, according to a weighting proportional to the barrier separating the two minima divided by their difference in energy. 51 The energy landscapes were visualized as disconnectivity graphs 52,53 using the PyConnect software package. 54 Analysis of clusters with n = 9 and n = 12 is presented here, and disconnectivity graphs for other clusters in the range 3 ≤ n ≤ 12 are available in the Supporting Information.

RESULTS AND DISCUSSION
For clusters with n ≤ 18, all eight basin-hopping runs located the same lowest common energy minimum, and for n = 19, seven of the eight basin-hopping runs located the same lowest common energy minimum. For n = 19, approximately 1.2 million energy minimizations were conducted, requiring ≈3.4 × 10 9 energy evaluations. These structures are probably good candidates for global minima (Figure 1). For n ≥ 20, none of the eight basin-hopping runs locate a common lowest-energy structure, and so the putative "global minima" presented at this size have been observed by only one of the independent basinhopping searches. Hence, they may be representative of lowlying minima, rather than the true global minimum. This is a sharp decrease in search success. To investigate further, at each cluster size in the range 16 ≤ n ≤ 21 (i.e., 3 water molecules either side of the transition) we performed 100 basin-hopping searches of 300 000 steps each and found that the effect persists over the larger number of searches. Studying the geometry of low-energy minima on either side of the drop-off does not reveal any obvious feature that might inhibit search success, and so the sharp fall could be a result of the increasing dimensionality of the search space (for a rigid-body molecular system such as this, the number of degrees of freedom increases by 6, not 3, for each water molecule added). The statistics for the successful searches are reported in the Supporting Information. It should be noted that the calculations reported here do not include zero-point vibrational energies, and it has been shown that some reordering of the isomer energies can occur when they are considered. 15,55 Global minima of SO 4 2− (H 2 O) n for n ≤ 7 have been presented in a previous study by Head-Gordon et al. 15 using a combined empirical potential/DFT search method. Our global minima for n = 4 and 5 agree with the DFT results (minima 4.5.3 and 5.7.3 in ref 15, respectively), and the DFT global minima for n = 3 and 6 (minima 3.6.0 and 6.7.5 in ref 15, respectively) are low-energy isomers (but not global minima) for the present potential. We were unable to reproduce the structure previously proposed as the global minimum for n = 7 (minimum 7.9.5 in ref 15). However, when treated at the DFT level using the B3LYP exchange-correlation functional and 6-311++G** basis set, this structure rearranged to a different geometry, previously reported as the third minimum, which we did find in our global optimization study. Our empirical global minimum for n = 7 is also the lowest energy structure when relaxed and evaluated at the DFT level, and corresponds to the second lowest DFT isomer previously reported (minimum 7.8.6 in ref 15). To the best of our knowledge, global minima for n ≥ 8 have not been reported in the literature, so for clusters above this size we have no previous results for a comparison. The evolution of the energy per water molecule, U(n)/n (Figure 2), suggests that n = 50 is not sufficient to have reached the asymptotic limit of the bulk, infinitely dilute, solution.
We can identify particularly stable cluster sizes using the central difference approximation to the second derivative of U(n), Figure 3: The more positive Δ 2 U(n) is, the more stable that structure is relative to its next nearest neighboring sizes. In the size range 3 ≤ n ≤ 15, the more stable structures are those that contain a high proportion of trimeric water rings (n = 7, 9, 12, 15). For 24 ≤ n ≤ 50 we observe alternating even−odd behavior, with structures containing even numbers of water molecules being more stable. This result contrasts with the behavior of homogeneous TIP4P water, 21 where even−odd behavior is observed for smaller clusters, before it disappears at larger sizes.
In Figure 4 we plot the mean hydrogen bond length between pairs of water molecules, r OH ww , and between a water molecule and one oxygen of the sulfate ion, r OH ws , of the global minima as a function of n. r OH ww is consistent with the mean hydrogen bond length in pure TIP4P water clusters, and r OH ws is ≈4% shorter than r OH ww for all n. An exception to this trend occurs at n = 5, where the mean is skewed by the long hydrogen bonds between

Journal of Chemical Theory and Computation
Article the two water molecules bonded along the edge of the sulfate. The comparative shortening of hydrogen bonds between water molecules and the sulfate is presumably due to the larger partial charge on the oxygen atom in the sulfate ion than on the lone pair in the TIP4P water. The shortening of r OH ws in the range 13 ≤ n ≤ 25 is probably due to the compression of the inner water shell in order to allow for tangential hydrogen bonding between water molecules that might otherwise be too far apart. Fluctuations between 12 and 14 coordination of the sulfate ion are responsible for the oscillatory behavior in r OH ws above n = 25. Here, we define the length of a hydrogen bond to be the OH···O distance.
In order to study the size-dependent appearance of dangling OH bonds, we examined the hydrogen bonding network of water molecules and the sulfate ion. We deemed structures with a water molecule donating only a single hydrogen bond, either to another water molecule or to the sulfate ion, to contain a dangling OH bond. Dangling OH bonds protruding radially from the cluster were observed in the global minimum structures for the larger hydrated sulfate clusters at n = 43, 45, and 47. In the smaller clusters n = 14 and 17, we observe semidangling OH bonds, which do not engage in a hydrogen bond, but do not protrude radially from the cluster surface. Subsequent local minimizations of these structures at the DFT level (B3LYP, 6-311++G**) show that these semidangling OH groups remain and are not an artifact of the potential. We define f ̅ as the Boltzmann-weighted mean number of dangling OH bonds where f i and ΔU i are the number of dangling OH bonds and energy above the global minimum for minimum i, respectively, and β ≡ 1/k B T. The sum is taken over all unique structures found during any one of the independent basin-hopping global optimization searches or (for 3 ≤ n ≤ 12) in transition state searches. In Figure 5 we plot f ̅ as a function of n, weighted by temperature, T = 130 K, to be consistent with experiment, 11 and assume that at this temperature the free energy can be estimated from the potential energy alone (i.e., that the entropic contribution to the free energy from each isomer is approximately the same). For the size range 3 ≤ n ≤ 19 there are three peaks in f ̅ at n = 8, 14, and 17 water molecules, due to the global minima of n = 14 and 17 exhibiting a dangling OH bond each. For SO 4 2− (H 2 O) 8 , there exist a number of lowenergy isomers with one or two dangling OH bonds, which are stabilized by the other water molecules participating in two trimeric water rings. For 20 ≤ n ≤ 39 we observe no energetically relevant structures exhibiting a dangling OH bond. For 40 ≤ n ≤ 50 we begin to observe system sizes with prevalent numbers of dangling OH bonds more frequently, with f ̅ ≈ 0.9 at n = 43, 45, and 48. At this stage, it is not obvious whether the pseudo-odd−even behavior of f ̅ in this size range is physical, or an artifact due to the difficulty of sampling structures at this size. Provisionally, the data is consistent with the IRPD spectra of size-selected hydrated sulfate clusters, 11 which indicate that the sulfate ion suppresses the appearance of dangling OH bonds until n ≈ 43.
This suppression is not currently well-understood, but we can begin to rationalize it as follows: a water molecule is capable of donating and accepting up to two hydrogen bonds (i.e., four hydrogen bonds in total). In order for a watercontaining system to totally inhibit the appearance of dangling OH bonds, each water must donate both OH bonds into a hydrogen bond. For systems solely composed of water molecules, this can only be achieved in highly coordinated environments such as bulk ice. At surfaces or in clusters, such coordination numbers are unachievable, and thus dangling OH bonds appear. In hydrated sulfate clusters, the sulfate anion occupies the center of the cluster, favoring highly coordinated

Journal of Chemical Theory and Computation
Article sites. When the first solvation shell is filled, it is able to accept 12−14 hydrogen bonds without donating any back into the system. In this manner, the sulfate anion acts as a net sink for hydrogen bonds, relaxing the requirement that the mean number of donated and accepted hydrogen bonds must be the same. We define r ̅ ad to be the ratio of the Boltzmann-weighted mean number of hydrogen bonds accepted by and donated to a water molecule where a i and d i are the number of hydrogen bonds accepted by and donated to a water molecule in minima i respectively, and ΔU i is the energy of minimum i above the global minimum. In Figure 6 we plot r ̅ ad as a function of the number of water molecules n, and observe that for small n, r ̅ ad is approximately 0.5 (i.e., a water molecule on average donates twice as many hydrogen bonds as it accepts). As the system size increases, the significance of the sulfate ion as a net hydrogen bond acceptor decreases, with r ̅ ad rising to ≈0.9 at n = 50, implying that 90% of hydrogen bonds are being accepted by water molecules. Note that, for systems composed entirely of water molecules (or in the limit of an infinitely dilute solution), r ̅ ad = 1. Though this analysis provides some insight into the size-dependent inhibition of dangling OH bonds, it is still not clear why they should begin to appear around 40 water molecules and above, and this will be the focus of future work.
3.1. Structural Analysis of Global Minima and Energy Landscapes for 3 ≤ n ≤ 12. Common structural motifs for the structures of the global minima of SO 4 2− (H 2 O) n for 3 ≤ n ≤ 12 are trimeric water rings about the open faces of the sulfate ion, with the water molecules in n = 3, 6, 9, 12 all engaged exclusively in such rings. There is some dispute regarding the size-dependent appearance of the second solvation shell, with estimates placed between 8 and 12 water molecules. 56,57 If we define the second solvation shell as beginning when a water molecule hydrogen bonds only with other water molecules (and not the sulfate ion), we find that the onset is ambiguous. In the lowest energy structures with n = 10 and 11, we observe a water molecule which hydrogen bonds only with other water molecules, and participates in a tetrameric ring instead of bonding onto the exposed face of the sulfate ion. At n = 12 all water molecules return to sharing a hydrogen bond with the sulfate ion, but for sizes larger than this, at least one water molecule is hydrogen bonded to other waters only.
We built databases of minima and transition states for clusters in the size range 3 ≤ n ≤ 12, and present results for n = 9 and n = 12, which contain 12 094 minima and 144 677 transition states and 14 419 minima and 228 415 transition states, respectively, and display interesting kinetic features. Data for other sizes can be found in the Supporting Information. Figures 7 and 8 show the disconnectivity graph and two lowestenergy isomers of SO 4 2− (H 2 O) 9 , respectively. The two isomers both have three trimeric hydrogen bonded water rings (i.e., the same oxygen skeleton), but differ in the relative directionality of the hydrogen bonded rings, resulting in an energy difference of ΔU = 0.22 kcal mol −1 . In spite of the energetic and structural similarity, the two isomers are separated by a barrier of ≈6 kcal mol −1 , leading to a frustrated landscape. The relationship between low-energy isomers that share an oxygen skeleton but differ in hydrogen bond directionality is observed for clusters of other sizes. The disconnectivity graph for SO 4 2− (H 2 O) 12 is shown in Figure 9, with the five lowest-energy isomers labeled and their structures shown in Figure 10. The isomers are organized into two oxygen skeletons, differing in the relative directionality of the hydrogen bonding in the rings alone. Isomers a, d, and e have four trimeric hydrogen bonded water rings, and isomers b and c have structures with one trimeric ring and a nine-membered water cycle. Both oxygen skeletons contain 12 water−sulfate and 12 water−water hydrogen bonds. As with the n = 9 system, large energetic barriers exist between the isomers, which will lead to frustrated kinetics. In both systems, the energetic gap (≈2.5 kcal mol −1 and ≈1.5 kcal    9 global minimum (a) and the next lowestenergy isomer (b), as labeled in the disconnectivity graph in Figure 7. ΔU is the energy of a given isomer above the global minimum. Energies are in kcal mol −1 .

Journal of Chemical Theory and Computation
Article mol −1 for n = 9 and n = 12, respectively) between the labeled minima and the next lowest-energy isomer suggests that the oxygen skeleton geometry is the key factor in determining the energy of hydrated clusters in this size range.
3.2. Structural Analysis of Global Minima for 13 ≤ n ≤ 19. For clusters with n ≥ 13, we begin to see 3-fold symmetric concave caps of six water molecules about the vertex of the sulfate (three water molecules hydrogen bond to an oxygen atom on the sulfate and the other three hydrogen bonds between them, see n = 18 for an illustration). One and two concave caps are observed for n = 15 and 18, respectively, and feature in the global minima and low-energy isomers for many of the larger cluster sizes. Search methods that bias toward such caps may prove to be a useful optimization strategy for future studies of hydrated sulfate clusters. Clusters containing n = 14 and n = 17 waters feature a single semidangling OH bond, and are the only two global minimum structures found in our study below n = 43 to do so. In both cases, the structures are similar to their n + 1 neighboring global minimum, but with a water molecule missing from the concave cap, causing one of the water molecules to have a dangling OH bond. It should be noted that the dangling OH bond belongs to a water molecule in the first solvation shell, and donates a hydrogen bond directly to the sulfate ion. Analysis of the normal-mode frequencies of these unusual structures confirms that they are minima on their respective potential energy surfaces.
The largest cluster size at which a trimeric water-ring appears around the open face of the sulfate ion is n = 15. For n ≥ 16 (with the exception of n = 17 which has a dangling water), the graph of hydrogen bonded water molecules is connected, and there are no self-contained subnetworks of water molecules that hydrogen bond to either themselves or the sulfate (for example, the four trimeric water rings in n = 12 are each disconnected from the other). Water "cubes" are observed, 31,58,59 which share an edge with the sulfate ion (see n = 16 and n = 19). The largest cluster for which the same putative global minimum is found consistently by the majority of independent basinhopping searches is n = 19.
3.3. Structural Analysis of Hydrated Sulfate Clusters for 20 ≤ n ≤ 50. For each cluster with n ≥ 20, the putative global minimum is found in only a single search, so we are not confident that they are the global minima, but we believe that they are structurally representative of low-energy isomers. Structures of putative global minima for 3 ≤ n ≤ 50 are given in the Supporting Information. Minima containing one, two, or three concave caps are a recurring feature, along with structures that appear for the homogeneous TIP4P water clusters, including fused cubes and edge-sharing pentagonal prisms. 30 The sulfate remains fairly central within the cluster throughout. The most important feature is that we do not see the emergence of dangling OH bonds in any of the global minima or low-energy isomers until n ≈ 40, with f¯growing from 0 to ≈0.9 ( Figure 5). The suppression of dangling OH bonds, and the cluster size at which they begin to appear, is roughly equivalent to what is found in the experimental IRPD spectra.
3.4. Conclusions. Microhydrated ions are interesting systems both experimentally and computationally for investigating the influence of ions on local water structure. We have shown that rigid-body modeling of hydrated sulfate clusters is capable of replicating the physical chemistry of these systems. For clusters containing a few water molecules (n ≤ 6), we are able to reproduce ab initio global minima, either as global minima of the empirical potential, or as low-energy isomers. In the size range 3 ≤ n ≤ 12, the structures of low-energy minima exhibit water molecules engaging in trimer water rings. The global minima for n = 3, 6, 9, and 12 water molecules are particularly stable, with the water molecules participating in such rings exclusively. For these systems, the lowest energy isomers typically have the same oxygen skeleton, but with different relative directionality of the hydrogen bonds, as explored with cycle inversion moves. Landscape analysis reveals that although these structures differ only slightly in energy, the barriers to interconversion can be very large, suggesting that such systems will display frustrated kinetics. At larger cluster sizes we note that dangling OH bonds are unfavorable up to n ≈ 43, consistent with IRPD spectra of size-selected clusters. We suggest that this is due to the sulfate ion acting as a net acceptor of hydrogen bonds, which allows water molecules in the system to accept only a fraction of the hydrogen bonds that they donate.
Future work will involve global optimization of hydrated sulfate clusters directly at the ab initio level using the BCGA-DFT code 60,61 and a combined DFT/basin-hopping 62 approach. We will also extend this work to study other microhydrated Hofmeister ions.    12 global minimum (a) and the next four lowest-energy isomers (b−e), as labeled on the disconnectivity graph in Figure 9. ΔU is the energy of a given isomer above the global minimum. Energies are in kcal mol −1 .

Journal of Chemical Theory and Computation
Article ■ AUTHOR INFORMATION Corresponding Author *E-mail: r.l.johnston@bham.ac.uk.