Chemical Energetics and Atomic Charges Distribution of Variably Sized Hydrated Sulfate Clusters in the light of Density Functional Theory

Among the ions classified in the Hofmeister series, the firstly ranked divalent sulfate anion has the strongest hydrating and water-structure making propensity. This unique characteristic actually makes it kosmotropic which causes water molecules to interact each other and contributes to gain structural stability of its hydrated clusters [SO4(H2O)n] n = 1−40. In this study, few variably sized microhydrated sulfate clusters [SO4(H2O)n] n = 1−4, 16 are considered separately, and inquired their chemical energetics and atomic charge distributions through ab initio based theoretical model. The main objective of this insight is to specify and interpret their thermodynamic stabilities, binding energies, and specific bonding and electronic interactions quantum mechanically. An in-depth analysis of their change in relative ground state electronic energy with respect to hydration number indicates stronger affinity of the sulfate ion towards water molecules while attaining structural stability in any aqueous type solutions. The mathematically determined values of their binding energy (E) almost holds up the same with this structural stability order: [SO4(H2O)16] > [SO4(H2O)4] > [SO4(H2O)3] > [SO4(H2O)2] > [SO4(H2O)], as reliable as experimentally and molecular dynamics simulation predicted trend. Moreover, the Mulliken derived partial atomic charges feature qualitative charge distribution in them which not only depicts electronic interactions between the specific atoms but also exemplifies the involvement of central sulfate units in hydrogen bond formation with surrounding water molecules.

> Li + > Na + = K + > NH4 + ) [3], [4] but also existing predominantly in different sized hydrated clusters [SO4 2− (H2O)n], n ≤ 50 [5] in (i) sea water (885 mg/L) [6]; (ii) flow battery systems [7], [8]; (iii) human body [9], [10]; (iv) beverages such as beer, juice, milk (coconut, human, animal, and infant formula), wine, mineral and distilled water etc.; (v) foods and vegetables such as almonds, broccoli, bread, cauliflower, dried fruits, sausage, peanuts etc.; (vi) soil (850 mg/kg); (vii) laxatives (Na2SO4 solution) used in medical treatments etc. In particular, a spontaneous hydration of the SO4 2− ions in the human body which serves as an essential step in facilitating their gastrointestinal absorption, biosynthesis of PAPS (3'-phosphoadenosine-5'-phosphosulfate: the most common coenzyme in sulfotransferase reactions), and various metabolic activities [10] has triggered theoretical, experimental and biochemists very much in order to investigate ground state molecular geometries, thermodynamic stabilities, quantum mechanical perspectives of the chemical bonding, step-by-step hydration process, chemical energetics, binding energies, atomic charges distribution, kinetic stabilities and molecular reactivity etc. of their variably sized hydrated clusters deeply. In the same context, present author has already reported an in-depth theoretical investigation of the sequential hydration followed by ground state molecular geometries of the microhydrated sulfate ions [SO4 2− (H2O)n] n = 0−4, 16 elsewhere [11]. In this contribution, the detailed chemical energetics including binding energy and atomic charges distribution of each hydrated cluster:  16] plus their role in stabilizing these distinctly hydrated clusters are described on the basis of quantum mechanical information derived from density functional theory (hereafter, DFT) model.
In physical chemistry, the hydrated molecular/ionic specimens are only said to be in thermodynamically stable states whenever they are in the lowest energy states or in the ground states. In other words, thermodynamic stability of any chemical systems can be interpreted on the basis of their chemical and structural stability which in turn relates directly to their energetic stability and ground state energy terms. Meanwhile, the energetic stability of any ionic or molecular specimen increases with decreasing potential energy itself is a function of these typical molecular aspects: bond strength (lengths and angles), chemical compositions, charge/electronic distributions, states of matter, intermolecular forces, ion-water interactions, hydration number etc. Therefore, the major impacts of these features must be incorporated sufficiently while computing ground state electronic energy/structure of the individual clusters [SO4 2− (H2O)n] n = 1−4, 16 quantum mechanically [12]. To address the same, a theory of electronic structure: the DFT model is quite popular among the researchers as it is based on electron density distribution function and Hohenberg-Kohn paradigm [13]. Moreover, determining binding energy (E) of each and every hydrated sulfate clusters from their DFT derived electronic energy is another quite essential task to understand their additional chemistry. In aqueous type solutions, the E of any hydrated ions is defined by an energy required to dismantle the hydrated molecular structures into their constituent molecules and ions [1] [2]. It is mathematically formulated with the total electronic energy of the whole hydrated sulfate structure ( ), and that of the unhydrated sulfate ( ) and free water molecule ( ) as: Where, n stands for the hydration number of each hydrated state of the sulfate ion.
Additionally, as the computational procedure expands molecular orbitals (hereafter, MOs) in terms of atomic orbital (hereafter, AO) basis functions and hence, are called electron rich mathematical functions, they can be used to depict electronic location and wave like behavior in the molecular/ionic specimens [14], [15], [16], [17]. But, while extracting electronic charges distributions associated with MOs, AOs, and atomic centers, population analysis methods [18] must be used. Among them, the most frequently used is Mulliken population analysis method [19] due to its highest efficiency and simplicity in estimating relative atomic charges of the molecular /ionic specimens theoretically through post-processing type quantum mechanical procedures. This method follows some arbitrary and orbital based schemes while partitioning the electron density (molecular wavefunction ( )) mathematically. Despite its recognition as a basis set dependent method, it is quite useful in examining electronic charge redistribution and bond populations in the bulk materials as well as in simple molecular/ionic specimens [19]. This is why, it is implemented here via the LCAO-Kohn-Sham ansatz based DFT model [20], [21], [22] for extracting Mulliken atomic charges related information from the DFT derived molecular wavefunctions of each cluster [SO4 2− (H2O)n] n = 1−4, 16. The numerical values of the Mulliken atomic charges assigned to each atom of these clusters would be very much needful for getting a rough idea of the charge distribution in them that in turn can quantify and compare their specific type bonding interactions more closely.
In short, this article emphasizes on quantification and interpretation of chemical energetics, binding energies, partial atomic charges distribution, and specific bonding plus electronic interactions in the hydrated sulfate clusters of the form [

II. COMPUTATIONAL DETAILS
Prior to determine the binding energy and partial atomic charges distribution of the variably sized hydrated sulfate clusters [SO4 2− (H2O)n] n = 1−4, 16, the ground state electronic energies and the concerned electronic structures were computed individually. For it, the computationally designed molecular models were optimized to their low energy states through the DFT model. While modelling the clusters [SO4 2− (H2O)n] n = 1−4, a divalent SO4 2− ion with a central S atom holding four O atoms in 3D axis was built as a structural module through GaussView visualization software [23], and then nH2O (n = 1, 2, 3, 4) molecules were added one by one to its immediate environs. But, in the case of [SO4 2− (H2O)16] cluster, the trial structure was displayed directly by using Monte Carlo simulation produced X, Y, Z coordinates available in the supplementary materials of the research article available elsewhere [24]. The symbol 'n' in the [SO4 2− (H2O)n] n = 1−4, 16 actually tells the hydration number (size of the hydrated cluster): monohydrated (n = 1, one molecule equivalent of water), dihydrated (n = 2, two molecules equivalent of water), trihydrated (n = 3, three molecules equivalent of water), tetrahydrated (n = 4, four molecules equivalent of water), and hexadecahydrated (n = 16, sixteen molecules equivalent of water). While adding one or more H2O molecules to the central SO4 2− unit, all the four O atoms of the latter ion were left unbound to neighboring H2O molecules themselves were completely free and disconnected to each other. This molecular model actually forces the DFT algorithms for realizing its H2O molecules as the solvated water molecules around the central SO4 2− computationally. All the Gaussian keywords and methodologies related to compute ground state electronic structure of each hydrated sulfate were selected according to the instructions given in Gaussian 09 manuals [25] [26]. The B3LYP hybrid functional with 6−31G (d, p) basis set was used i.e. the methodology used was DFT: B3LYP / 6−31G (d, p). Under it, a selfconsistent field (SCF) with default SCF procedure (SCF=Tight) and Berny algorithm were implemented for solving electronic Schrodinger equation iteratively while optimizing the given trial structure to the global minimum respectively. The Gaussian output files such as .log, checkpoint (.chk), and formatted checkpoint (.fchk) files were read via GaussView, and all the required three dimensionally displayed chemical data were extracted. Since the current interest of this author is limited to explore chemical energetics, binding energy, and partial atomic charge distribution in each [SO4 2− (H2O)n] n = 1−4, 16 cluster separately, the ground state electronic energy term of each of them were extracted carefully followed by the mathematical calculations of the respective binding energy with the help of the equation formulated in Eq. 1, and all the Mulliken derived partial atomic charges of each cluster were displayed on the respective DFT derived ground state electronic structure on the basis of their charge density (color indices and ranges).

Chemical Energetics and Binding Energy
As reported earlier by the same author elsewhere [11], the DFT derived 3D geometry of each hydrated and unhydrated sulfate clusters [SO4 2− (H2O)n], n = 0−4, 16 followed the criteria of "no imaginary frequencies" while computing ground state frequencies at their DFT produced atomic Cartesian coordinates. It reconfirmed us that the DFT produced geometry of each of them represent a true minimum of their own potential energy surface or these structures serve as their most stable ground state electronic structures. This conclusion was made on the basis of this pre-existing principle: absence of negative frequencies in the Gaussian output file means the atomic Cartesian coordinates of the theoretically converged geometry do not dislocate at all from their initial positions while the frequency calculations are on the fly. In terms of quantum mechanical sense in physical chemistry, such global minimum electronic structure always means its thermodynamically stable state which further represents ground state electronic geometry or energy, i.e. a thermodynamic stability only occurs whenever a chemical system is in its lowest energy state. Hence, determining low energy system (more stable structure) theoretically from the concerned high energy (less stable) system is one of the most significant thermodynamically favorable steps. Here, the DFT model is employed to determine ground state electronic energies from the related high energy state (trial structure) of each and every hydrated sulfate ion clusters. All these DFT derived energies in atomic unit are presented in Table 1. For easy comparison, the relative electronic energies of them are also determined and tabulated in the fourth column of the same table. A graphical plot between the relative electronic energies of them and their respective hydration number (n) (Figure 1) clearly shows that greater the hydration number, lesser is the electronic energy (more stable) possessed by the sulfate ion. This conclusion agrees well with the trends derived from the advanced experimental tools and  [27]: smaller the sulfate-water clusters, more is the electronic energy possessed by them; and the clusters having an average n of ~12 (experiment) and ~11 (MD simulation) would have the least electronic energy than any other smaller clusters. Unfortunately, the electronic energies for the sulfate-water clusters with n = 11, and 12 were not computed here due to the shortage of their X, Y, Z Cartesian coordinates. As a whole, the decreasing trend of the ground state electronic energies with increasing hydration numbers ( Figure 1) clarified that the sulfate ions always want to make clusters with variable number of H2O molecules during the course of attaining stabilization in aqueous type solutions without which their predominant existence in several aqueous type media could not be recognized. The same is the reason why sulfate ions prefer to attract large number of oppositely charged poles of water molecules in their vicinity just like other Hofmeister ions. Moreover, while plotting the positions of these variably sized hydrated sulfate clusters on an energy  Figure 2, where energy gap between each cluster with respect to monohydrated (n = 1) sulfate ion is clearly mentioned and logically arranged.   Similarly, the calculated values of the binding energy (E) for each sulfate ion cluster are listed in Table 2. As a rule, greater the binding energy, higher is the structural stability, and more negative are its values, more stable are the hydrated structures. It further stresses that the relative stability of the sulfate-water clusters increases with increasing hydration number n, but there is a specific value of n in the most stable cluster. As sketched in Figure 3, the binding energy for the [SO4 2− (H2O)16] cluster is very high but that for the lower hydrated clusters are almost equal. It may be due to the instability of the lower hydrated sulfate clusters in the aqueous type solutions, further stressing a strong hydrating tendency of the sulfate ion. While arranging all of these hydrated sulfate clusters on the basis of their ground state electronic energies and the calculated E values, the approximate stability order would be

Atomic Charge Distribution
In theoretical research, determining numerical values of the atomic charges and analyzing their distribution in the molecular/ionic specimens are quite significant tasks while exploring electronic interactions, and quantifying and comparing specific type atomic interactions. One of the widely used methods to accomplish them scientifically is population analysis method, a computationally cheap yet reliable technique to derive numerical values of the atomic charges [24]. This method actually divides up the density matrix and overlap matrix terms of the orbital occupation ( ) expression ( ( ) = ), and predicts  tentative positions of the electron density lying inside the orbital surfaces of the concerned molecular/ionic specimens. More specifically, the Mulliken population analysis is quite popular method among researchers for estimating partial atomic charges and approximating their distribution in the molecular/ionic systems qualitatively. Mathematically, it involves summation of all the populations for all the atomic orbitals on a single atom followed by the substraction of the nuclear charge to give the resultant partial atomic charge on each atom. The DFT-implemented Mulliken derived atomic charge distribution in variably sized hydrated sulfate clusters [SO4 2− (H2O)n], n = 1−4, 16 are summarized in Table 3 and 4, where all the numerically labeled atoms are listed in accordance with the atomic label of Figure 4(a-e). On the basis of the magnitude of these Mulliken charges, each and every atoms of each hydrated sulfate ion cluster are colored and displayed individually as DFT-generated ground state electronic structures (Figure 4), where color indices and ranges are shown in corresponding insets. From the color indices and color intensity of each of these optimized structures, the qualitative charge distribution in the entire ionic specimens can be understood: all the red colored (O atoms), the most green colored (S atoms), and the dark green colored (H atoms) spheroids represent the most negatively, most positively, and less positively charged atoms respectively. The same conclusion can be drawn by analyzing the respective numeral Mulliken atomic charges listed in Table 1. In particular, the central S atom of each SO4 2− unit of each cluster is comparatively the most positively charged due to its involvement in bonding with more electronegative O atoms themselves are involved in interacting with the surrounding H2O molecules through the hydrogen bond (hereafter, H−bond) as can be seen clearly in Figure 4(a-e). Beside this, all the H−bonded hydrogen atoms of the solvated H2O molecules are more positively charged than the central O atoms, since the latter atoms are more electronegative they would prefer attracting the electrons towards their vicinity. Again, while summing up all the atomic charges distributed in each hydrated sulfate cluster, the net Mulliken atomic charges are found as around −2.00 a.u, exactly equal to the charge unit carried by the central sulfate ion during hydration process with neutral H2O molecules. In conclusion, such a significant inhomogeneous charge distribution in each of these ionic specimens is mostly arose due to a strong ability of the O atom to attract a bond pair towards itself as its electronegativity (µO = 3.44 in Pauling scale) is significantly greater than that for the H (µH = 2.20 in Pauling scale) and S (µS = 2.58 in Pauling scale) atoms. This effect is more pronounced in the H−bonded oxygen and hydrogen atoms as observed in each and every hydrated sulfate clusters. An over interpretation of the Mulliken derived partial atomic charges is not recommended as the Mulliken method is highly sensitive to the types of basis sets used.

IV. CONCLUSION
In this contribution, the DFT computed ground state electronic energies of the variably sized hydrated clusters of the sulfate ion [SO4 2− (H2O)n], n = 1−4, 16 followed by the mathematical determination of the minimum amount of energy required to dismantle them into their constituent molecules and ions (binding energy) are presented. Apart from this, the Mulliken derived partial atomic charges distributions in them and their consequences in bonding/electronic interactions are also explained logically.
While analyzing the graphical plot between the DFT derived ground state electronic energies (relative) of each and every hydrated sulfate clusters and their respective hydration number (n), it was noticed that greater the hydration number, lesser (higher) is the electronic energy (structural stability). This trend ultimately clarified why sulfate ions tend to undergo sequential hydration in the course of attaining thermodynamic stability in any aqueous type solutions. This conclusion was found to be quite reasonable with the step-by-step hydrating nature of other Hofmeister ions (Hofmeister series of anions: SO4 2− > HPO4 2− > CH3COO − > Cl − > NO3 − and, cations: Mg 2+ > Li + > Na + = K + > NH4 + ). Accordingly, the stability order [SO4 2− (H2O values of the binding energy (E). Besides this, the Mulliken derived partial atomic charges were found to feature qualitative charge distribution in each and every clusters: the O atoms are most negatively charged, the S atoms are most positively charged, and the H atoms are least positively charged. More particularly, this atomic charge distribution not only depicted the electronic and bonding interactions that take place in between more electropositive S atom and more electronegative O atoms of each SO4 2− unit of each cluster but also described the involvement of these four O atoms in hydrogen bond formation with H atoms of the surrounding H2O molecules.