A Comment on Nadytko et al., “Amines in the Earth’s Atmosphere: A Density Functional Theory Study of the Thermochemistry of Pre-Nucleation Clusters”. Entropy 2011, 13, 554–569

Nadykto, Yu, Jakovleva, Herb and Xu have recently reported a DFT study on the structure and formation thermodynamics of sulfuric acid-base-water clusters, with ammonia and a handful of amines as bases [1]. This study partially overlaps with our previous work [2], and a significant part of the discussion in their manuscript concerns differences between their results and ours. This comment is intended to address some issues related to that discussion. Specifically, it is shown that the errors related to basis-set effects in our calculations are very likely much smaller than claimed by Nadykto et al. [1]. Composite calculations including e.g., higher-level electron correlation also agree better


Introduction
Clusters of sulfuric acid with water and various base molecules are generally believed to be important for the first steps of atmospheric nucleation, and have therefore been extensively studied by both experimental and computational methods.Both Nadykto et al. [1] and our study [2] agree that the binding of most amines to a single sulfuric acid molecule is considerably stronger than that of ammonia.Primary and secondary amines also promote the addition of a second sulfuric acid to the cluster more effectively than ammonia, though the two studies disagree in the magnitude of the difference.Our study only computed (H 2 SO 4 )(X) + H 2 SO 4 addition energies and free energies for OPEN ACCESS X = ammonia and X = dimethylamine, while Nadykto et al. also considered methylamine and trimethylamine.They found that while tertiary amines such as trimethylamine bind very strongly to one acid, they are not able to stabilize a cluster containing two sulfuric acids.This is not surprising: since tertiary amines are only able to form one hydrogen bond, they are therefore less effective at binding two acid molecules together.This is also known from experiments [3][4][5]: while primary and secondary amines promote aerosol formation mainly through salt formation (clustering with acid molecules), tertiary amines promote aerosol formation via their oxidation products.Thus, the fairly low enhancement factors found in the experimental study by Erupe et al. [6] (reference [37] in Nadykto et al. [1]) using a tertiary amine are likely due to the choice of amine.Studies using primary or secondary amines typically find much larger enhancement factors, as illustrated e.g., by the recent results of Berndt et al. [7].As shown also by the thermodynamics computed by both of our groups [1,2], the enhancement effect depends strongly on the specific structure of the amine-a single study on a single amine will inevitably give only a partial picture.
The main difference between Nadykto et al. and our study is the difference in the computed free energies of the (H 2 SO 4 )(X) + H 2 SO 4 clustering reactions for X = ammonia and X = dimethylamine.Whereas our study predicted a reaction free energy difference of slightly less than 5 kcal/mol, Nadykto et al. report values of about 2.7 kcal/mol.This difference is significant, as the ratio of atmospheric concentrations of ammonia and amines is believed to lie in the 100-1000:1 region (though local variations in both directions may be much larger [5]).A straightforward application of the law of mass action indicates that a 1 kcal/mol difference in free energies has the same effect (at 298 K) on cluster concentrations as a 5.4-fold difference in reactant concentrations.Thus, 5 kcal/mol would correspond to a concentration difference of around 4500, while 2.7 would only correspond to a factor of 90 or so.However, it should be noted that the concentration of free X does not enter the mass-balance equation for the (H 2 SO 4 )(X) + H 2 SO 4 ⇔ (H 2 SO 4 ) 2 (X) 1 reaction.Instead, the concentration of the (H 2 SO 4 ) 2 (X) 1 product depends directly only on the concentrations of H 2 SO 4 and the (H 2 SO 4 )(X) clusters.As both our studies qualitatively agree, the difference in free energies for the H 2 SO 4 + X ⇔ (H 2 SO 4 )(X) reactions for X = ammonia and X = dimethylamine is fairly large (e.g., 3.6 kcal/mol according to Nadykto et al. and 7.0 kcal/mol according to our study), indicating that many (H 2 SO 4 )(amine) clusters may have equal or larger concentrations than the (H 2 SO 4 )(NH 3 ) clusters despite the difference in vapor concentrations.
Nonetheless, I do not disagree with Nadykto et al. that formation of (H 2 SO 4 ) 2 (amine) 1 clusters (or possibly even (H 2 SO 4 ) 3 (amine) 1 clusters) are likely the bottleneck reaction in amine-enhanced nucleation, at least if the absence of significant activation energies in the collision reactions is assumed.Correspondingly, the bottleneck for ammonia-enhanced nucleation is likely to be the formation of (H 2 SO 4 )(NH 3 ) clusters, which is fairly unfavorable at typical monomer concentrations (ppt-level for H 2 SO 4 and ppb-level for NH 3 ).Thus, in order to compare the effect of ammonia and amines, both the H 2 SO 4 + X and (H 2 SO 4 )(X) + H 2 SO 4 clustering reactions should be accounted for, in addition to all other possible relevant reactions.Larger clusters may also need to be considered, especially as secondary amines such as dimethylamine can only form two hydrogen bonds, while primary amines can form three and ammonia up to four.This indicates that ammonia and primary amines very likely have the advantage over secondary amines in promoting the growth of two-acid clusters to three-acid clusters.Whether or not this is significant in turn depends on the acid concentration, and the size of the critical cluster.In any case, full description of the cluster-forming process will require a dynamic cluster model using evaporation rates computed from quantum chemistry-a simple comparison of two (or even four) reactions will not conclusively answer any question about the process.However, such comparisons can give preliminary indications of the relative importance of different mechanisms, and the statements made by Nadykto et al. thus deserve some further comment.Nadykto et al. [1] claim that most of the difference in the computed free energies between their study and ours [2] is likely to be caused by basis-set issues in the MP2 energies.I disagree with this claim, and offer herein computational evidence to show otherwise.

Results and Discussion
Nadykto et al. [1] frequently refer to basis-set superposition error (BSSE) as a major error source in MP2 calculations.For weakly bound clusters, an estimate of the BSSE can be obtained by performing a counterpoise (CP) calculation, where the cluster is partitioned into its constituent monomers.There are, however, significant controversies concerning the accuracy and applicability of the CP correction, with no clear consensus on the issue.See for example reference [25] for a discussion on the subject.In some cases, the CP correction clearly improves agreement with basis-set limit binding energies, while in other cases it just makes things worse.This is a subtle issue: while the CP correction almost always makes the convergence toward basis-set limit values smoother, there is no guarantee that the actual value of the correction for any single basis set will improve the agreement.As stated by Dunning [26]: "It is quite possible, and even probable, that binding energies computed without the counterpoise correction are closer to the complete basis set limit than the uncorrected values".As an illustration of this, a recent study by de Lange and Lane [27] on weakly bound complexes found that the CP correction sometimes improves, and sometimes deteriorates, agreement with basis-set limit binding energies.Whether the CP correction is beneficial or not depends on the relative magnitudes of basis set superposition compared to other basis-set-related errors, e.g., basis-set incompleteness error, BSIE (or "basis set convergence errors" in the terminology of Dunning [15]).If these other error sources dominate, application of the CP correction may well result in much larger errors than neglecting it.This is further supported by some recent computational results.For example, Alvarez-Idaboy and Galano [28] found that "CP corrections systematically leads to results that differ from the CBS-extrapolated ones to a larger extension than the uncorrected ones", and hence recommended not including CP corrections in interaction energy calculations.Furthermore, in a study on molecular complexes and clusters of aerosol nucleation precursors, Zhao et al. [29] found that composite methods (which account for basis-set effects using extrapolation rather than explicit CP corrections) "likely produce the most reliable binding energies".
In clusters of strong acids such as sulfuric acid and strong bases such as an amine, one or more protons have typically been transferred from the acid to the base, forming an ion pair ("salt monomer").This is the case for example for the sulfuric acid-dimethylamine dimer cluster according to both Nadykto et al. [1] and our study [2].For such a cluster, the counterpoise correction cannot be performed properly, as the original neutral monomer subunits no longer exist in the product cluster (this is especially true for larger clusters, where the transferred proton may end up rather far from the original donating functional group-and in which the question of which proton is the transferred one may not even have a unambiguous answer).Performing the CP calculation with the anion and cation of the ion pair as subunits, in the cases where this can unambiguously be done, does not provide information on the BSSE in the formation energy from neutral monomers.A further issue is that proton transfer may not be a completely binary issue: also clusters of H 2 SO 4 with somewhat weaker bases like ammonia may have a "partially transferred" proton, in which case application of the CP correction may also be problematic.Unfortunately, the issue of CP corrections for clusters with partially transferred protons, while relevant and interesting, would likely require a separate study to assess completely.
The only way to estimate the magnitude of BSSE (and other basis-set-related errors, which may be of equal or greater importance) for proton-transferred clusters of strong acids and bases is to systematically perform calculations with larger and larger basis sets, and then extrapolate the basis-set limit binding energies if necessary.Various well-established methods exist for such extrapolations, see e.g., Dunning [15] or Jensen [30] for examples and details.Since our study [2] contained RI-CC2 and RI-MP2 energies computed using both the aug-cc-pV(D + d)Z and aug-cc-pV(T + d)Z basis sets, such extrapolations could be performed already with the data presented there-with the result that the basis-set limit binding energies are within about 1 kcal/mol of the aug-cc-pV(T + d)Z values.However, since the quality of extrapolations using the aug-cc-pV(D+d)Z data may sometimes be questionable, I have recalculated the RI-MP2 binding energies for the (H 2 SO 4 )((CH 3 ) 2 NH) cluster using basis sets up to aug-cc-pV(6 + d)Z.To keep the task feasible, I used the PW91/6-311++G(3df,3dp) geometries; the same level of theory used by Nadykto et al. [1].Results given below indicate that the difference between single-point energies computed using PW91/6-311++G(3df,3dp) and RI-MP2/aug-cc-pV(D + d)Z geometries are less than 0.1 kcal/mol.I also computed explicitly correlated DF-MP2-F12 binding energies up to the VQZ-F12 basis set.Explicitly correlated methods converge much faster toward the basis-set limit, and also provide an additional check of the reliability of the values.The results are given in Table 1.It can be seen from Table 1 that the aug-cc-pV(6 + d)Z and VQZ-F12 binding energies are both within 0.6 kcal/mol of the aug-cc-pV(T+d)Z values.Also, the difference between the aug-cc-pV(6+d)Z and aug-cc-pV(5 + d)Z binding energies is only 0.03 kcal/mol, demonstrating that the binding energies are essentially converged to the basis-set limit at the aug-cc-pV(6 + d)Z level, and further extrapolation is thus unnecessary if an accuracy of around ± 0.1 kcal/mol is sufficient.The absolute RI-MP2 energies, as well as the absolute RHF energies, and the Cartesian co-ordinates of the (H 2 SO 4 )((CH 3 ) 2 NH) cluster and free monomers used in the calculation, are given in a supplementary information file, and can be used to test the various available extrapolation schemes to obtain even more exact estimates of the RI-MP2 basis-set limit energies.In any case, it is clear that for the (H 2 SO 4 )((CH 3 ) 2 NH) cluster, RI-MP2 basis-set effects beyond aug-cc-pV(T+d)Z are below 1 kcal/mol, not several kcal/mol as claimed by the Nadykto et al. [1].Table 1.Formation energies (electronic energies, no zero-point corrections) for the (H 2 SO 4 )((CH 3 ) 2 NH) cluster, computed with various methods.All values in kcal/mol.(a) single-point energies at the RI-MP2/aug-cc-pV(D+d)Z geometries, [2]; (b) single-point energies at the PW91/6-311++G(3df,3pd) geometries, see supplementary information for co-ordinates.
This result does not mean that RI-MP2 energies or free energies are necessarily within 1 kcal/mol of the real values-both e.g., higher-order correlation and vibrational anharmonicity may easily change the free energies by 1 kcal/mol or more.In order to obtain more reliable values for the free energies including these effects, I further performed G2 and G3 calculations for the (H 2 SO 4 )((CH 3 ) 2 NH) and (H 2 SO 4 )(NH 3 ) clusters.The G2 and G3 composite methods have been found to yield reliable energies for water clusters by Dunn [31].I also performed G2MP2 and G3MP2 calculations for the larger (H 2 SO 4 ) 2 ((CH 3 ) 2 NH) 1 and (H 2 SO 4 ) 2 (NH 3 ) 1 clusters, which were too large for the full G2 and G3 calculations.G2MP2 and G3MP2 are approximate versions of the full G2 and G3 methods, which usually yield very similar results at much reduced cost.The PW91 geometries were used as input guesses for these computations; each composite method calculation then performs geometry optimizations and frequency calculations at levels specified in the corresponding method references.
The results are given in Table 2.It can be seen that the "Gn"-family of methods predict free energies that are in all cases fairly close to the values given by our study [2] and quite far off from the PW91 values [1].

Table 2.
Free energies (at 298 K and 1 atm reference pressure) for the addition of sulfuric acid to various clusters.All values in kcal/mol.(a)RI-CC2/aug-cc-pV(T+d)Z//RI-MP2/augcc-pV(D+d)Z, [2]; (b)PW91/6-311++G(3df,3p), [1].The good agreement between G2, G3 and the more approximate G2MP2 and G3MP2 methods for the smaller (one-acid) clusters indicate that the results for the larger (two-acid) clusters are very likely also quite reliable.It should be noted that while the Gn family of methods account for higher-order correlation, their treatment of vibrational anharmonicity is based on the use of single scaling factors, which may be insufficient for treating clusters.Especially as the clusters grow to larger sizes, more accurate treatment of anharmonicity such as hindered rotor analysis and frequency-dependent scaling factors may be required.As clusters are more anharmonic than free molecules, these corrections would lead to even more negative cluster formation free energy values, and thus further increase the difference between the computed free energies and the PW91 results [1].

Conclusions
While I strongly agree with Nadykto et al. [1] that further research on ammonia-and amine-containing sulfuric acid and sulfuric acid-water clusters is required, especially concerning hydrated clusters and clusters with more than two acids, the calculations presented here show that (1) basis-set superposition error is extremely unlikely to be the main explanation for the differences between RI-CC2/RI-MP2 and PW91 results on sulfuric acid-amine clusters, and (2) higher-level calculations agree significantly better with RI-CC2/RI-MP2 results than with PW91 results.Supplementary Table 2. Cartesian co-ordinates (in Ångström) for sulfuric acid, dimethylamine and their dimer cluster, optimized at the PW91/6-311++G(3df, 3pd) level.Note: these PW91 geometries were optimized with Gaussian 09.Nadykto et al. [1] have not stated which program they have used for their calculations (presumably Gaussian 03 or Gaussian 09).Especially if a different version has been used, there may be minor differences between these geometries and those used in Nadykto et al.This will not affect any of the conclusions as the effect on energies is small (<0.1 kcal/mol).