Modeling the Binding Free Energy of Large Atmospheric Sulfuric Acid-Ammonia Clusters.

Sulfuric acid and ammonia are believed to account for a large fraction of new-particle formation in the atmosphere. However, it remains unclear how small clusters grow to larger sizes, eventually ending up as stable aerosol particles. Here we present the largest sulfuric acid-ammonia clusters studied to date using quantum chemical methods by calculating the binding free energies of (SA) n (A) n clusters, with n up to 20. Based on benchmark calculations, we apply the B97-3c//GFN1-xTB level of theory to calculate the cluster structures and thermochemical parameters. We find that the cluster structures drastically evolve at larger sizes. We identify that an ammonium ion is fully coordinated in the core of the cluster at n = 7, and at n = 13 we see the emergence of the first fully coordinated bisulfate ion. We identify multiple ammonium and bisulfate ions that are embedded in the core of the cluster structure at n = 19. The binding free energy per acid-base pair levels out around n = 8-10, indicating that at a certain point the thermochemistry of the clusters converges toward a constant value.

Introduction

Aerosols are airborne particles spanning a large range of sizes, from a few nanometers for freshly nucleated particles up to micrometer sizes for large cloud droplets. Aerosols have various compositions depending on both the chemical species present at formation and their growth via the condensation of low-volatile species onto the existing particles. Aerosols play an important role in relation to climate change both directly by scattering and absorbing light in the atmosphere and indirectly by acting as cloud condensation nuclei (CCN), which are necessary for initiating cloud formation.[1] These effects can have both negative and positive effects on the global energy budget depending on the exact composition of the aerosols. As designated by the recent sixth assessment report by the IPCC, aerosol–cloud interactions still contribute the largest uncertainty to climate estimation.[2]

Atmospheric new-particle formation (NPF) is initiated by the formation of stable atmospheric molecular clusters. NPF from gas-phase molecules is estimated to be the source of roughly half of all CCN,[3] with the other half being particulate matter emitted directly into the atmosphere, such as dust and sea spray. The majority of the clusters formed from gas-phase molecules are too small to initially act as CCN, and they are often unlikely to grow further before they are scavenged by larger particles. Therefore, elucidating the mechanisms leading to the successful growth of aerosol particles to CCN sizes is crucial to better understand the climate impact of aerosol particles. Several possible mechanisms have been proposed to explain particle formation, but it has been shown that nearly all NPF in the present-day atmosphere involves ammonia or biogenic organic compounds in addition to sulfuric acid.[4] Hence, as a starting point we herein focus on sulfuric acid–ammonia clusters.

Quantum chemical calculations have been useful for studying sulfuric acid (SA)–ammonia (A) cluster structures and thermochemistry. Early computational work was performed by Ianni et al. in 1999,[5] where they calculated the thermochemistry of small (SA)1–2(A) clusters with varying degrees of hydration at the B3LYP/6-311++G(2d,2p) level of theory. Through extensive contributions from several groups,[5−12] the size and studied composition of sulfuric acid–ammonia clusters have been expanded over the years. In 2012, Ortega et al.[13] reported the structures and thermochemistry of (SA)1–4(A)1–4 clusters. This was the first complete cluster set with all combinations of up to four acids and four bases. Such a data set allowed for subsequent simulations of the cluster kinetics using the atmospheric cluster dynamics code (ACDC).[14,15] The cluster size was substantially increased by DePalma et al.[16] in 2014, where they calculated the thermochemistry of large (SA)(A) clusters, with n up to 8, at the PW91/6-31++G(d,p) level of theory. It was shown that the free energy of the clusters decreased almost linearly with the system size. To the best of our knowledge, these are the largest sulfuric acid–ammonia clusters studied to date, and no attempts have been made to expand to larger sizes.

Atmospheric cluster dynamics simulations have given insight into the formation mechanism of sulfuric acid–ammonia clusters. It has been shown that the clusters with a 1:1 ratio of acids and bases are the most stable.[15,17] Furthermore, it was found that sulfuric acid–ammonia clusters, specifically (SA)3–4(A)3–4 cluster sizes, were already quite stable against evaporation.[17] Recently, Besel et al.[18] studied clusters with up to six sulfuric acid molecules and six ammonia molecules that were calculated at the DLPNO-CCSD(T0)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory. They found that if the simulated clusters were sufficiently large, the boundary conditions for outgrowing clusters only had a small influence on the simulated new-particle formation rates. This indicates that (SA)6(A)6 or larger clusters can be considered quite stable against evaporation. Here we attempt to push the limit of the cluster sizes modeled using quantum chemical methods by studying the cluster structures and thermochemistry of large (SA)(A) clusters, with n up to 20.

Computational Details

Density functional theory (PW91,[19] M06-2X,[20] and ωB97X-D[21]) and semiempirical (PM7[22]) calculations were performed using the Gaussian 16 program.[23] We applied the Gaussian 09 default convergence criteria to allow for direct comparisons with values in the atmospheric cluster database (ACDB)[24] and the study by DePalma et al.[16] The domain-based local pair natural orbital DLPNO-CCSD(T0)[25,26] and DLPNO-MP2[27] calculations, as well as the empirically corrected B97-3c[28] and PBEh-3c[29] calculations, were performed using the ORCA program (ver. 4.2.1).[30] We used the semicanonical (T0) approximation, which neglected the off-diagonal Fock-matrix elements, to calculate the perturbative triple correction, as it has been shown that it leads to a performance similar to that of the recently introduced improved iterative (T) approximation[31] for atmospheric molecular clusters.[32] The GFN1-xTB[33] and GFN2-xTB[34] calculations were performed using the xTB program.[35]

Binding Free Energies

We calculate the binding free energy (ΔGbind) of the (SA)(A) clusters asAnalogously, we can define the cluster electronic binding energy (ΔEbind), the binding enthalpy (ΔHbind), and binding entropy (ΔSbind). The binding free energy can conveniently be written as the binding electronic energy plus a ”thermal” correction termIn this manner, the cluster structures and therefore the ΔGbind,thermal term can be calculated using a lower level of theory, and the binding energies can be calculated on top of the cluster structures using a higher level of theory. In all cases we report the calculated binding free energies at 298.15 K and 1 atm using rigid rotor and harmonic oscillator approximations.

Determination of Unique Structures

The initial cluster structures were generated using the ABCluster program[36,37] with the CHARMM force field.[38] As force field methods cannot account for bond breaking, we calculated the lowest 1000 minima,using ionic HSO4– and NH4+ monomers. This enforces a single proton transfer from sulfuric acid to ammonia, essentially leading to ammonium–bisulfate clusters. As proton transfer is always found in the lowest free energy (SA)(A) clusters with n > 1,[7,11,15,17,18] this approach should provide an adequate description of the cluster structures. The 1000 minima generated by ABCluster were initially optimized at the PM7 level of theory.

Some of the structures found in the force field calculations done by ABCluster will end up converging to the same structure when treated at the PM7 level. To eliminate duplicates, the cluster structures were aligned using the ArbAlign program,[39] and the root mean square deviations (RMSDs) between the clusters were pairwise calculated using the Kabsch algorithm. The minimum RMSD that allowed two structures to be significantly different was set to 0.38 Å based on previous experience with sulfuric acid–water clusters.[40,41] This reduced the number of structures approximately by a factor of 2–3 (see the SI). Even after removing duplicates we end up with 157–489 unique structures for each value of n. This implies that, despite the clusters being quite large, many of the ABCluster-generated local minima actually converge to identical structures. For the five clusters with the lowest free nergies at the ωB97X-D/6-31++G(d,p)//PM7 level of theory, we did a full geometry optimization and vibrational frequency calculation using PW91/6-31++G(d,p) for values of n up to 9. Figure presents the binding free energies of the calculated clusters compared to the ones obtained by DePalma et al. (also at the PW91/6-31++G(d,p) level of theory).

Figure 1

Calculated binding free energies (at 298.15 K and 1 atm) of the (SA)(A) clusters with n = 6–9 compared to the work of DePalma et al.[16] with n up to 8.

Using the outlined cluster configurational sampling technique, we obtained cluster structures that were lower in free energy compared to the ones previously reported. While this is encouraging, it should be mentioned that our approach is by no means an exhaustive search for the lowest free energy structures, implying that lower minima may exist. However, as the clusters increase in size, the number of low-lying configurations will increase. Therefore, we can assume that the presented clusters are quite close in free energy to the global free energy minimum. From Figure we see that the near-linear trend found by DePalma et al.[16] continues with n = 9. However, continuing to calculate the cluster structures and vibrational frequencies using PW91/6-31++G(d,p) will be prohibitively computationally expensive for larger systems. Hence, to target larger clusters we need to identify a methodology that can be applied for values n up to 20.

Results and Discussion

Benchmarking

To identify an adequate level of theory for modeling the large (SA)(A) clusters structures with n = 6–20, we have to benchmark the cluster structures, the thermal contribution to the free energy (ΔGbind,thermal), and the electronic binding energies (ΔEbind). The aim is to identify a level of theory that can be applied across all the studied cluster systems. Besel et al.[18] reported the (SA)(A) clusters, where n = 1–6, at the DLPNO-CCSD(T0)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory. To the best of our knowledge, this set of clusters constitutes the highest level of theory applied for values n up to 6, and we will use it as a benchmark set for the structures as well as the ΔGbind,thermal and ΔEbind contributions.

Cluster Structures

We tested the semiempirical PM7, GFN1-xTB, and GFN2-xTB methods to obtain the molecular geometries of the (SA)(A) clusters with n = 1–6. We used the geometries obtained by Besel et al.[18] at the ωB97X-D/6-31++G(d,p) level of theory as a reference, and these structures were also used as the input for the semiempirical calculations. As a comparison, we also tested the GGA functional PW91 with the small 6-31+G(d) and 6-31++G(d,p) basis sets. Table presents the calculated RMSD values (using the ArbAlign Pprogram[39]) of the (SA)(A) clusters with n = 1–6 compared to those calculated at the ωB97X-D/6-31++G(d,p) level of theory.

Table 1

Root Mean Square Deviation (RMSD, Å) in the Cluster Geometries Compared to the ωB97X-D/6-31++G(d,p) Level of Theory

	PM7	GFN1-xTB	GFN2-xTB	PW91/6-31+G(d)	PW91/6-31++G(d,p)
(SA)1(A)1	0.28	0.18	0.23	0.04	0.06
(SA)2(A)2	0.40	0.32	0.22	0.10	0.11
(SA)3(A)3	0.57	0.23	0.52	0.23	0.23
(SA)4(A)4	0.69	0.55	0.57	0.25	0.25
(SA)5(A)5	0.85	0.33	0.55	0.18	0.18
(SA)6(A)6	0.59	0.27	0.57	0.24	0.25
mean	0.56	0.31	0.44	0.17	0.18

Unsurprisingly, we see that the PW91/6-31++G(d,p) and PW91/6-31+G(d) levels of theory yield structures quite similar to the ωB97X-D/6-31++G(d,p)-calculated structures, with RMSD values of 0.25 Å or below. In our similarity tests (see section 2.2), we treat structures with RMSD values of 0.38 Å or below as duplicates, indicating that the PW91 and ωB97X-D optimized (SA)(A) clusters with n = 1–6 clusters are in fact identical. Out of the semiempirical methods, PM7 performed the worst with RMSD values as high as 0.85 Å for the (SA)5(A)5 cluster. Surprisingly, GFN2-xTB performed significantly worse than GFN1-xTB, with RMSD values of 0.52 Å and above for n ≤ 3. From the data in Table , it is clear that GFN1-xTB might be an attractive alternative for obtaining the (SA)(A) cluster structures in cases where DFT is not applicable.

Thermal Contribution to the Free Energy

Figure presents the calculated deviations in the value of ΔGbind,thermal for the tested methods compared to those from the ωB97X-D/6-31++G(d,p) calculations. The numerical values for each cluster are shown in the SI.

Figure 2

Calculated thermal correction to the binding free energies (at 298.15 K and 1 atm) of the (SA)(A) clusters with n = 1–6. The calculations are compared to the ωB97X-D/6-31++G(d,p) values (taken from ref (18), harmonic approximation).

The error in the value of ΔGbind,thermal for PM7 increases linearly with the system size, with a maximum error of 28.7 kcal mol–1 for the (SA)6(A)6 cluster (see the SI). Such a catastrophic error implies that PM7 is unsuited for calculating the ΔGbind,thermal values for these systems. It should be noted that PM7 is in fact parametrized toward the heat of formation and not electronic energies, which might contribute to the calculated large error. Out of the tested methods, the PW91/6-31++G(d,p) level of theory shows the best agreement with the ωB97X-D/6-31++G(d,p) calculations, with a mean absolute error (MAE) of 0.8 kcal mol–1 and maximum error (MaxE) of 1.7 kcal mol–1. Lowering the basis set to 6-31+G(d) leads to a substantial increase in the errors (MAE = 2.2 and MaxE = 4.2 kcal mol–1). Interestingly, the GFN1-xTB and GFN2-xTB methods show error more or less similar to that of the PW91/6-31+G(d) level of theory. This implies that for (SA)(A) clusters where the PW91/6-31++G(d,p) level becomes too computationally expensive applying, either the semiempirical GFN1-xTB or GFN2-xTB method might be a valid choice.

Binding Energies

While the applicable methods for obtaining the large cluster structures and vibrational frequencies are limited, there are more possibilities for calculating the electronic binding energies. Here we tested a range of semiempirical methods (PM7, GFN1-xTB, and GFN2-xTB), empirically corrected DFT methods (B97-3c and PBEh-3c), density functionals (PW91, M06-2X, and ωB97X-D with the 631++G(d,p) basis set), and DLPNO methods (DLPNO-CCSD(T0) and DLPNO-MP2 with the aug-cc-pVDZ basis set). Figure presents the calculated errors in the binding energy (ΔEbind) of the tested methods compared to those from the DLPNO-CCSD(T0)/aug-cc-pVTZ calculations reported by Besel et al.[18]

Figure 3

Calculated binding energies of the (SA)(A) clusters with n = 1–6. The calculations are compared to the DLPNO-CCSD(T0)/aug-cc-pVTZ values (taken from ref (18)).

The lowest error was obtained at the DLPNO-CCSD(T0)/aug-cc-pVDZ level of theory. However, it is unlikely that the DLPNO-CCSD(T0)/aug-cc-pVDZ level of theory can routinely be applied to very high values of n due to the steep memory requirements of the triples corrections in the DLPNO-CCSD(T0) methods. The DLPNO-MP2/aug-cc-pVDZ level of theory yield a large MAE of 7.0 kcal mol–1, with maximum errors up to 14.8 kcal mol–1, and in general does not seem to outperform the DFT/6-31++G(d,p) calculations. Out of the tested DFT functionals, the PW91/6-31++G(d,p) level of theory exhibited the lowest errors, with a MAE of 3.1 kcal mol–1 and a MaxE of 6.8 kcal mol–1. Both M06-2X and ωB97X-D present significantly larger deviations. Interestingly, the ωB97X-D/6-31++G(d,p) level of theory presents linearly increasing errors with system size (R2 = 0.99), indicating that even though it presents the largest errors out of the tested DFT functionals the errors might be more systematic. The semiempirical methods also all present rather large MAE and MaxE values but are on par with (and not significantly worse than) the DFT/6-31++G(d,p) levels. Interestingly, the empirically corrected B97-3c method exhibits low errors of, with a MAE of 2.1 kcal mol–1 and a maximum error of 3.4 kcal mol–1.

Based on the preceding three sections, we can conclude that the B97-3c//GFN1-xTB level of theory appears as the most cost efficient approach to obtain relatively accurate binding free energies of (SA)(A), with n = 6–20 cluster systems. While this level of theory works well for these systems, it cannot necessarily be expected that it is transferable to other cluster compositions.

Cluster Thermochemistry

We calculated the binding free energies of the unique (SA)(A) clusters structures with n = 6–20 at the B97-3c//GFN1-xTB level of theory. Figure present the identified lowest free energies as a function of the number of acid–base dimers (n). The inclusion of one or two H2SO4 and SO42– pairs in the sampling of the largest cluster with n = 20 was also tested. However, including one or two H2SO4 and SO42– pairs in the sampling yielded structures at least 9.94 and 3.53 kcal mol–1 higher in free energy, respectively (calculated at the GFN1-xTB level of theory). Nevertheless, including such H2SO4 and SO42– pairs in the sampling might be important for clusters even larger than those studied here or clusters that are composed of more bases than acids.

Figure 4

Calculated lowest binding free energies of the (SA)(A) clusters. The n = 1–6 clusters were taken from Besel et al.,[18] with the red dot data points calculated at the DLPNO-CCSD(T0)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory and the green dot data points recalculated at the B97-3c//GFN1-xTB level of theory. The blue dot data points are the extension for n up to 20 calculated at the B97-3c//GFN1-xTB level of theory. The black line is a linear least-squares fit to the data points.

The data from Besel et al., which was calculated at the DLPNO-CCSD(T0)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory (red dots), matches quite well with the data recalculated at the B97-3c//GFN1-xTB level of theory (green dots). This further illustrates the applicability of the B97-3c//GFN1-xTB level of theory in obtaining the binding free energies of the clusters. The binding free energy is seen to almost linearly decrease as a function of the number of acid–base pairs (n). Hence, the almost linear trend originally observed by DePalma et al.[16] for the (SA)(A) clusters, with n up to 8, is here shown to continue up to n = 20. When reaching n = 15, 16, n = 17, 18, and n = 19, 20 there is a stepwise function in the free energy, with the jump from even to odd being significantly more favorable than the jump from odd to even. This could indicate that there is an even–odd preference for larger clusters. However, we were not able to deduce what might cause this effect based on the cluster structures.

Evaporation Free Energies

The evaporation free energy of an acid–base pair can be calculated asThe calculated evaporation free energies as a function of n are plotted in Figure . The n ≤ 6 data were calculated at the DLPNO-CCSD(T0)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory (data were taken from Besel et al., ref (18)), and the n ≤ 7 data were calculated at the B97-3c//GFN1-xTB level of theory.

Figure 5

Calculated dimer (n) evaporation free energy of the (SA)(A) clusters.

An oscillatory behavior is visible in the evaporation free energies of the clusters with values between −11.3 and −39.5 kcal mol–1. Such low free energies should translate into quite low evaporation rates, implying that the evaporation of an acid–base dimer is highly unlikely. This is consistent with previous studies[15,17,42] on the evaporation kinetics of acid–base clusters, where it was found that evaporation predominantly occurred via the monomers. Hence, the erratic oscillatory behavior could indicate that it is important to further investigate the off-diagonal clusters such as (SA)(A) and (SA)(A) to identify whether the evaporation of single SA or A components is prominent. It should be noted that the calculations were performed at 298.15 K and 1 atm, corresponding to boundary layer conditions. At higher altitudes (i.e., lower temperatures) the free energies will be even lower and thus evaporation will be further suppressed.

Evolution of the Cluster Structures

The cluster structures drastically evolve as a function of the number of acid–base pairs. For n ≤ 6, all the molecules in the clusters are exposed to the exterior. At n = 7, the lowest free energy cluster structure has a single ammonium ion in the core of the cluster that is fully coordinated to all the surrounding molecules. A similar effect was observed by DePalma et al.[43] in positively charged (SA)(A)+ clusters for n = 7–10. This “encapsulation” effect was argued to make the core ions inaccessible to substitution with stronger bases such as alkylamines.[43] At n = 13, we observed multiple ions that were fully coordinated (one bisulfate ion and two ammonium ions). Hence, the structure almost resembles a fully coordinated cluster in the “particle” environment. At n = 19, we observed the first emergence of two fully coordinated bisulfate ions together with four fully coordinated ammonium ions.

The linear trend observed in Figure does not yield much information about the free energy change of the clusters as the value of n increase. As additional acid–base dimers are added to the cluster, the free energy gain per added n will decrease. This can be illustrated by plotting the free energies per acid–base dimer n (see Figure ). The n ≤ 6 data were taken from Besel et al.,[18] where they were calculated at the DLPNO-CCSD(T0)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory. The n ≤ 7 data were calculated at the B97-3c//GFN1-xTB level of theory. The free energies were calculated at 298.15 K and 1 atm.

Figure 6

(Left) Calculated average thermal contribution to the free energy of the (SA)(A) clusters per n. (Right) Calculated binding free energy of the (SA)(A) clusters per acid–base pair. The red dot data points were calculated at the DLPNO-CCSD(T0)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory (taken from ref (18)), and the blue dot data points were calculated at the B97-3c//GFN1-xTB level of theory.

By inspecting the thermal contribution to the free energy, we find that it does not vary much. We find that the average thermal contribution to the free energy rapidly reaches a value of 24 kcal mol–1 per n at n = 6 and more or less converges at ∼26 kcal mol–1 per n at n = 20.

The cluster stabilization (i.e., the the free energy per n) rapidly diminishes and levels out with a value around −25 kcal mol–1 around n = 8–10. This clearly suggests that there is a large stabilization for small clusters, and the free energy gain per increasing n reaches a constant value as the clusters grow to larger sizes.

Conclusions

We have presented the largest (SA)(A) cluster structures studied to date using quantum chemical methods. We tested the performance of several semiempircal methods (PM7, GFN1-xTB, and GFN2-xTB) to obtain the cluster structures and the thermal contribution to the free energy and compared the data to literature ωB97X-D/6-31++G(d,p) calculations. We further tested the performance of several methodologies to calculate the binding energies of the clusters compared to high- level DLPNO-CCSD(T0)/aug-cc-pVTZ calculations. We identified the B97-3c//GFB1-xTB level of theory as an efficient low-cost methodology that could be applied to very large clusters. It should be further tested whether this methodology also yields satisfactory results for other cluster systems.

Applying the identified methodology, we studied the binding free energies of (SA)(A) clusters, with n = 6–20. The free energy of the cluster structures was found to decrease almost linearly as a function of n. Considering the free energy gain per n, we see that the cluster stabilization rapidly decreases as a function of n and levels out around n = 8–10. This work is the first to study the free energy surface of massive atmospheric molecular clusters and will in the future be extended to include the off-diagonal (SA)(A) and (SA)(A) clusters to further explore the free energy surface.