Computational Design of New Refrigerant Fluids Based on

Article pubs.acs.org/IECR

Computational Design of New Refrigerant Fluids Based on Environmental, Safety, and Thermodynamic Characteristics† Andrei Kazakov,* Mark O. McLinden, and Michael Frenkel Thermophysical Properties Division, National Institute of Standards and Technology, 325 Broadway, Boulder, Colorado 80305-3337, United States S Supporting Information *

ABSTRACT: We present a systematic search for new classes of refrigerants that would possess low values of global warming potential (GWP), along with low- to moderate flammability and suitable thermodynamic characteristics. We have developed new methods for estimating, solely from the molecular structure, the radiative efficiency (RE, a measure of radiative climate forcing) and atmospheric lifetime; the combination of RE and lifetime yields an estimate of the GWP. We also developed an estimate of the lower flammability limit (LFL) based on the enthalpy of formation. These estimation techniques, along with a previously developed technique for estimating critical temperature (Tc), are applied to a library of over 56 000 candidate molecules. We select fluids with GWP < 200; 300 K < Tc < 550 K; and LFL > 0.1 kg·m−3. Filters for toxicity and chemical stability based on functional groups are also applied to arrive at 1234 candidates for further study. The candidates that would be suitable for use in present types of refrigeration equipment (those having critical temperatures less than 400 K) are dominated by halogenated alkenes; additional chemical classes, including halogenated ethers and cyclic compounds, are identified among fluids with higher critical temperatures.

■

■

INTRODUCTION

ESTIMATION OF GWP Global warming potential was introduced in the first Intergovernmental Panel on Climate Change (IPCC) Assessment Report6 as a simple metric for comparison of potential greenhouse gases in terms of their impacts on the climate system. In the IPCC reports that followed,7−9 the GWP concept was further refined, and in spite of being criticized10,11 for having a number of significant drawbacks and limitations, it is widely used to guide and support decisions in industry and environmental policies. The definition of GWP is built upon the concept of radiative forcing due to a perturbation in or the introduction of a trace gas. Radiative forcing is defined as “the change in net irradiance ... at the tropopause after allowing stratospheric temperatures to readjust to radiative equilibrium, but with surface and tropospheric temperatures held fixed”.7 GWP, in turn, is defined as the ratio of radiative forcing, RF, due to a pulse release of 1 kg of trace compound integrated over a chosen time horizon, TH, to that of the reference gas, CO2:

Recently, the concerns of climate change have raised the issue of replacing widely used hydrofluorocarbon (HFC) refrigerants with low global warming potential (GWP) alternatives.1 Among the possibilities, hydrofluoroethers (HFEs) and, especially, unsaturated HFCs such as hydrofluoroolefins (HFOs) are being considered. Due to the presence of double bonds and, consequently, high reaction rates with atmospheric OH, HFOs have much lower atmospheric lifetime than the HFC-based refrigerants presently in use, which also results in much lower GWP. However, a more systematic search for new classes of refrigerants that, in addition to having low GWP, would also satisfy performance and safety requirements, has not been attempted. Initial proposals for design strategies seen in the literature are confined to qualitative considerations based only on radiative efficiency of potential candidates.2 On the other hand, large-scale virtual screening and computational design are actively used in drug discovery3 as well as explored4 and promoted at a policy level5 for design of novel materials. The search for efficient and environmentally friendly refrigerants is also within the scope of this methodology, although specific screening criteria and associated property estimation methods need to be developed. This work is an initial stage of a larger study that will involve not only property screening but also testing potential candidates in refrigeration cycle modeling. Here, we present our initial efforts on the development of computational tools for screening potential candidates for refrigerant fluids. Various screening criteria (GWP, flammability, critical temperature, toxicity, etc.) were applied based on candidates’ molecular structures. A screening example using a large diverse library of molecules is also presented. This article not subject to U.S. Copyright. Published 2012 by the American Chemical Society

TH

GWP =

∫0 RF dt TH

∫0 RFCO2 dt

(1)

If one adopts the global mean formulation for radiative forcing per unit increase in atmospheric concentration (radiative efficiency, RE) and a simple exponential decay function for Received: Revised: Accepted: Published: 12537

June 18, 2012 August 20, 2012 September 4, 2012 September 4, 2012 dx.doi.org/10.1021/ie3016126 | Ind. Eng. Chem. Res. 2012, 51, 12537−12548

Industrial & Engineering Chemistry Research

Article

beginning. Clearly, the indirect effects will be ignored during the initial screening. Second, the global-mean, well-mixed approximation is adopted throughout, even for short-lived compounds. This approach provides a consistent upper limit estimate, which should be sufficient for initial screening. By adopting these approximations, eq 2 can be used to estimate GWP given estimates for RE and τ (the first term the right hand side can be evaluated from documented model equations for CO217). Finally, for the remainder of this paper, we focus on the case of TH = 100 years, which is the value most commonly used in the policy framework. The GWP computed for a 100year time horizon will be referred to as GWP100. Estimation of Radiative Efficiency. The procedures for estimation of radiative efficiency available in the literature include the method proposed by Pinnock et al.12 and the group-contribution scheme by Young et al.25 The groupcontribution method25 is clearly a very computationally efficient option. However, it was developed mainly for fluorinated ethers and, as such, is of limited scope. Furthermore, Bravo et al.26 have noted degradation in its performance when applied even to perfluorocarbons that are structurally similar to the original target compounds. Therefore, the present efforts were focused on the method of Pinnock et al.12 This approach is also computationally efficient, exhibited very good performance in comparison with detailed models, and is widely used in the literature. In addition to the wellmixed approximation discussed above, the method also ignores the stratospheric temperature adjustment required by the RE definition. However, the error associated with this assumption is typically within 10−20%, which is comparable with experimental uncertainties associated with RE determination.9 In its original formulation, the method relies on a knowledge of the infrared (IR) spectrum obtained from experiment. Papasavva et al.27 have shown that Pinnock’s method also performs well when IR spectra predicted from quantumchemical calculations are used; this was further validated by subsequent theoretical studies.25,26,28 Good quantitative predictions were obtained with ab initio Møller−Plesset secondorder perturbation theory (MP2) and hybrid density functional theory (DFT) Becke, three-parameter, Lee−Yang−Parr (B3LYP) method using medium-to-large basis set sizes. These approaches are too computationally expensive for large-scale screening purposes; the present computational speed requirements do not permit levels of theory higher than those represented by the semiempirical family of methods. Semiempirical quantum chemical calculations for prediction of radiative efficiency were originally rejected for not being sufficiently accurate.29 This conclusion was made based on a comparison of computed frequencies and IR intensities with their experimental counterparts for a single compound, CF3CH2F, and no systematic assessment of RE predictions was made. On the other hand, semiempirical calculations can be performed very quickly even on modest modern computer hardware, and newer, improved methods have been developed over recent years.30,31 The data set for validation of estimation methods used in this study was assembled primarily from the most recent IPCC9,32 and WMO17 compilations that are normally used as the main reference sources for GWP-related data. The data from the reports were supplemented with data from the recent literature. The set contains slightly over 100 compounds, which are listed in the Supporting Information. While selecting the data and given an explicit choice, the preference was intentionally given

the time evolution of the trace gas pulse, the expression for GWP can be rewritten as 1

GWP =

∫0

TH RECO2cCO2(t ) WCO2

REτ[1 − exp(− TH/τ )] W dt (2)

where RE and RECO2 are radiative efficiencies, and W and WCO2 molecular weights of the trace and reference gases, respectively, τ the atmospheric lifetime of the trace gas, and cCO2(t) is the normalized concentration of the reference gas after time t from the initial pulse (described with a model that is more complex than a simple exponential decay9). In most GWP definitions found in the literature, RE is defined on a mass basis, while the reported RE values are always given in volumetric units. For consistency with reported values, we adopt volumetric units for RE here; consequently, molecular weights appear in eq 2 for the conversion back to mass units required by the GWP definition. In general, the evaluation of GWP requires detailed computer simulations describing atmospheric transport, radiative heat transfer, and multiple trace gas removal mechanisms.12−14 In practice, this has been done only for a small number of compounds. For the majority of the reported values, GWP is estimated based on eq 2, along with separately estimated values of RE and τ. The situation with GWP estimation for compounds with short atmospheric lifetimes (typically, less than a year13,15) is particularly problematic. When the rate of compound removal becomes comparable with the rates of atmospheric transport processes and there is not enough time for homogeneous mixing throughout the troposphere, the global mean assumption breaks down, and the effective value of GWP becomes dependent on the emission scenario (geographic location, time of the year, etc.). The GWP estimates for such systems based on detailed simulations are lower than those based on the global-mean, well-mixed approximations; differences in excess of 2 orders of magnitude were reported in extreme cases.14 Consequently, the IPCC8,9 and the World Meteorological Organization (WMO)16,17 reports refrain from providing GWP estimates for short-lived compounds. Approximate GWP values for short-lived compounds often appear in research literature, usually with the note that they represent an upper limit for actual values.18−24 Another issue of concern is the fact that the approach based on eq 2 accounts only for direct effects of trace gas emission, that is, the direct participation of the compound itself in atmospheric radiative heat exchange. Indirect effects include the impacts of the breakdown products produced during compound removal via atmospheric chemical reactions. Reaction products have radiative efficiencies and atmospheric lifetimes of their own that may exceed those of the initial compound, causing larger effective GWP values. Another indirect effect is a compound’s possible participation in the destruction of stratospheric ozone, which offsets radiative equilibrium toward cooling. Because the primary objective of the present study is to develop an efficient GWP screening procedure capable of processing a large number of compounds, the main requirements are computational speed, broad coverage, and consistency, followed by the accuracy of the predicted values. With detailed simulations being out of the question, several necessary compromises have to be accepted from the 12538

dx.doi.org/10.1021/ie3016126 | Ind. Eng. Chem. Res. 2012, 51, 12537−12548


Article

TRAN2008),41 Northwest Infrared Library,42 and the NIST Quantitative Infrared Database.43 These provided extra validation and allowed the comparison of the experimental and predicted values on the same basis. The original sources were investigated, and a search for recent updates was also conducted. As a result of these efforts, modifications to the data set were made, as described next. For cases with IR spectra available from the reference databases, the recomputed RE values resulted in much better agreement with the present theoretical estimates. The discrepancies were caused by various corrections introduced into the originally reported values and roundoff truncations for very small RE values. When available, the uncorrected values from the same sources as the reported values were taken. In some cases, however, inaccuracies in the IR spectra used to derive the reported values were the main cause for disagreement. For example, the value for CClF2CF3 (CFC-115) can be traced to the work of Myhre et al.44 who used an IR spectrum from an earlier version of the HITRAN database. The spectrum from the Northwest Infrared Library has a broader coverage and includes two additional absorption bands that are not present in HITRAN, thus resulting in a higher RE value that also agrees well with the present theoretical estimate. Another, less severe example is tetrachloromethane. The reported value comes from Jain et al.13 who used an unpublished IR spectrum quoted in ref 45. While the full spectrum was not provided, the given45 integrated absorption cross section deviates from those computed using the spectra from the reference databases by 22−27%. In two cases, CF3CF2CH2OCH3 (HFE-365mcf3) and CF3CH2OCH3 (HFE-263fb2), new experimental results for RE values emerged recently22,24 that agree with the present estimates. Three compounds, CHF2CF2OCH2CF3 (HFE347pcf2), c-C3F6, and CF3CH2CH2CH2OH, were removed from the data set. The origins of the value for HFE-347pcf2 appear untraceable. The RE value for perfluorocyclopropane, cC3F6, is erroneous: it originates from the early WMO compilation16 and was taken from the theoretical estimate27 for a different compound. As expected, the present estimate for c-C3F6 also agrees with the result of Bravo et al.26 based on B3LYP/6-31g(d,p) calculations. The value for CF3CH2CH2CH2OH comes from ref 20 and is inconsistent with the data for CF3CH2CH2OH from the same source: according to this report, RE for CF3CH2CH2CH2OH is almost a factor of 2 lower than that for CF3CH2CH2OH. The present theoretical estimate agrees well with the data point for CF3CH2CH2OH and suggests a very similar value for CF3CH2CH2CH2OH, which is reasonable considering their structural similarity. The comparison of predicted and experimental RE, after modifications, is shown in Figure 1. The computed estimates generally agree with the experimental data within estimated uncertainties. Next, the predictive capabilities of semiempirical methods were investigated. Four methods, AM1,46 PM3,47 RM1,30 and PM6,31 were tested following the same procedures outlined for B3LYP/Sadlej pVTZ. As expected, all methods performed significantly worse than the DFT. If based on the value of objective function alone, surprisingly, the oldest method, AM1, performed the best, while the latest one, PM6, exhibited the worst performance. However, with the exception of PM6, all methods yielded very low frequency scaling factors, below 0.7. Inspection of several representative cases by comparison with

to the values obtained in a manner consistent with the Pinnock’s method (i.e., without corrections for stratospheric temperature adjustment or lifetime corrections13,15). This deliberate decision was made to separate errors caused by quantum-chemical methods (the subject of this analysis) from the deviations caused by various, often empirical, corrections that make the reported value inconsistent with predictions by the Pinnock’s method. All quantum-chemical calculations in the current study were performed with the Gaussian 09 package.33 The initial three-dimensional molecular structures were produced following an improved and extended version of the procedure used previously34 with sampling, where applicable, of multiple conformations and stereoisomers. A final single structure was chosen based on the lowest free energy at standard conditions obtained at the level of theory used. Before proceeding to semiempirical calculations, to establish a baseline and to verify the consistency of the data, B3LYP calculations with the Sadlej pVTZ basis set35,36 were performed for all compounds in the set. The Sadlej pVTZ basis set was optimized for electrical properties, and the B3LYP/Sadlej pVTZ combination has shown superior performance in predicting vibrational frequencies and IR intensities as compared with more expensive, higher-level methods.37,38 The computation of RE using predicted IR spectra followed the procedures of Papasavva et al.;27 that is, individual band broadening was not considered. The use of band broadening was tested and did not show significant improvements, consistent with the findings of Bravo et al.26 When computing RE with quantum chemical methods, a single adjustable parameter, a uniform vibrational frequency scaling factor, was used. Unlike the previous studies that selected a frequency scaling factor based on a comparison between the computed and experimentally measured vibrational frequencies, the scaling factors here were determined by minimizing the deviations between the predicted and reported RE. Scaling factors are known to show mild frequency dependence,37 and minimization conducted in this manner provides proper weighting for the frequencies contributing to RE the most. The minimization procedure was carried out with an in-house program based on a differential evolution algorithm.39 The stochastic and derivative-free nature of this optimizer allows handling objective functions that exhibit discontinuities present in Pinnock’s method. An objective function based on the Cauchy distribution40 was chosen, as it offers some resistance to anticipated outliers. The uncertainties in the listed RE values were estimated based on the number of significant digits provided in compilations and a typical value of 15%.9 The optimized scaling factor for B3LYP/Sadlej pVTZ vibrational frequencies was found to be 1.0345, and the corresponding RE predictions are also given in the Supporting Information. The agreement between the literature results and the predictions based on B3LYP/Sadlej pVTZ IR spectra is excellent, although some outliers outside of estimated uncertainties were present. The observed level of agreement suggests the possibility that these theoretical results can also be used as a tool for the identification of erroneous and inconsistent RE data. To illustrate this point and to improve the data set consistency, a detailed outlier analysis was performed for the cases in which current calculations deviate from the reported data in excess of 30%. Several venues were pursued. When available, RE values were computed using the Pinnock’s model and experimental IR spectra from several reference databases: High-Resolution Transmission Molecular Absorption Database (HI12539



Article

Figure 1. Comparison of radiative efficiencies predicted using B3LYP/ Sadlej pVTZ IR spectra with the experimental values. Computed vibrational frequencies were scaled by 1.0345. Open symbols: shortlived compounds. Closed symbols: long-lived compounds. Crosses: compounds excluded from the set (see text).

Figure 2. Comparison of radiative efficiencies predicted using PM6 IR spectra with the experimental values. Computed vibrational frequencies were scaled by 0.9069. The shaded area represents the logarithmic RMSD.

photolysis can play an important role for long-lived compounds. A complete consideration of all these factors during virtual screening is clearly not feasible. For this initial effort, we focus only on one loss mechanism, removal via reactions with OH. This decision is based on the following considerations. From eq 2, GWP depends on lifetime via an integrated normalized concentration,

the baseline B3LYP/Sadlej pVTZ results revealed that abnormally low frequency scaling factors are not due to overestimated vibrational frequencies but rather due to high IR intensities being exhibited by incorrect vibrational modes at much lower wavenumbers. Although this correction appeared fairly systematic, it was believed to be due to insufficient molecular diversity of the present data set dominated by very similar fluorinated compounds. This was further tested by using subsets of the present data set with higher molecular diversity sampled with data balancing techniques48 using Soergel molecular similarity distance metrics as computed by the OpenBabel package.49 Again, with the exception of that for PM6, the scaling factors exhibited substantial instabilities as the diversity of the sample increased. Based on this analysis, PM6 appears to offer better consistency and extendability in predicting RE as compared to other semiempirical methods. These features are more critical for the present purposes, even at the expense of some accuracy loss for the particular data set used here. Therefore, PM6 was selected for RE screening calculations. The comparison of PM6-calculated RE using the final optimized frequency scaling factor of 0.9069 with the experimental data for the present data set is shown in Figure 2. As noted, the agreement is substantially worse than that for the DFT calculations, with several prominent outliers present. The outliers mainly include small highly symmetric compounds such as CH3Cl, CH3Br, CH3CCl3, and NF3. The logarithmic root-mean-squared deviation (RMSD) computed for the entire set corresponds to a factor of 1.84. This level of agreement, however, is very reasonable for virtual screening applications, especially considering the several orders of magnitude difference in computational effort between PM6 and B3LYP for typical cases. Estimation of Lifetime. Computationally efficient estimation of atmospheric lifetime is an even more challenging problem than that for RE. The first issue arises from the fact that there are multiple mechanisms of compound removal from the atmosphere.17 Trace compounds with short lifetimes are primarily removed in the troposphere via chemical reactions (with OH, ozone, NO3, and Cl), ultraviolet (UV) photolysis, and hydrolysis-related processes such as wet deposition, rainout, and washout. Stratospheric loss processes via reactions with OH, ozone, and singlet oxygen atoms and via UV

c ̃ = τ[1 − exp( −TH/τ )]

(3)

c̃ approaches τ for τ ≪ TH and becomes less sensitive to τ as it increases, approaching TH for τ ≫ TH, that is, becoming independent of τ. Therefore, the accuracy in estimating τ is important mainly for the cases when τ ≤ TH. For this range of lifetimes with the chosen TH (i.e., 100 years), reaction with OH appears to be the dominant loss mechanism for the majority of compounds for which atmospheric lifetime data are available (with a few notable exceptions to be discussed later). Estimation of τ based on the reaction with OH, in turn, requires an estimate for the corresponding rate coefficient, kOH, which poses a challenge of its own. The methods for estimation of rate coefficients for bimolecular reactions range from fundamental approaches based on transition state theory to empirical group-contribution techniques.50 Application of transition state theory requires the knowledge of the potential energy surface, with details and accuracy achievable only with very high-level ab initio methods, and, as such, is not feasible for the present framework. Computational time restrictions confine our choices to empirical estimation methods. Among those specifically developed for atmospheric oxidation by OH, there are the widely used group-contribution method of Atkinson51 and correlations based on descriptors computed from quantum chemical calculations proposed by Klamt52,53 and further extended by Böhnhardt et al.54,55 (MOOH method). The MOOH method was reported to have better performance than the Atkinson scheme and can be used with computationally inexpensive semiempirical methods.54 However, to-date, the Atkinson group-contribution scheme has wider coverage than MOOH, while the differences in predictive abilities are relatively modest for the present purposes. Additionally, the Atkinson method has a readily available reference implementation,56,57 and it was chosen for this study. An extended57 implementation of the Atkinson scheme was used to predict kOH with one additional modification. By design, the Atkinson scheme cannot reproduce the trends in 12540



Article

kOH for fluorinated and chlorinated ethylenes that exhibit nonmonotonic behavior with the degree of halogenation.58 The substituent factors for halogen atoms attached to doublebonded carbon were, therefore, taken at compromise values that do not work for the entire series of fluorinated and chlorinated ethylenes. These choices, on the other hand, affect the assignments for the substituent factors for other groups attached to double-bonded carbons, including halogenated methyl groups that appear in some modern refrigerants. In its original formulation,51 the scheme had a substituent factor only for −CH2Cl derived from the data for a single compound. The same value, 0.76, was later extended to all −CH2X groups,57 where X is F, Cl, Br, or I. However, no substituent factors were applied for −CHX2 or −CX3 groups, which makes the estimation procedure somewhat inconsistent. Our tests using the same base group assignments and recent experimental data58 indicate that the substituent factors for −CH2F and −CF3 have similar values between 0.5 and 0.6. As a compromise, during the present estimation procedure, the substituent factors for all −CH2X, −CHX2, and −CX3 groups were assigned the single value of 0.6. Having the estimated kOH values, the corresponding atmospheric lifetime was computed as

τ=

α k OH[OH]

reproduced well over several orders of magnitude. A cluster of outliers appears at the predicted c̃ of 100 years, an upper limit defined by eq 3. These are compounds for which the Atkinson scheme produced no rate estimate and an infinite lifetime was assumed. These are also compounds that are mainly destroyed by UV photolysis, and the central assumption of the present estimation procedure fails. UV photolysis is typically a slow process; thus, the assumption of an infinite lifetime results in a relatively modest error in c̃ for the majority of cases. A particularly drastic exception is demonstrated by an outlier in the right bottom corner of Figure 3: c̃ for this case is overpredicted by several orders of magnitude. This point corresponds to CF3I. The C−I bond is very weak, which makes UV photolysis extremely effective for this compound; the lifetime of CF3I is among the shortest in the data set (see the Supporting Information). This problem is expected to occur for other iodine-containing compounds. Excluding these outliers, the logarithmic RMSD for the data set corresponds to a factor of 2.7. This represents a much higher error than reported for kOH for the original method.51 However, it was also noted that halogenated compounds that dominate the present data set exhibit larger than average deviations. Overall, the observed performance is acceptable for virtual screening applications. Estimation of GWP100. Having established the framework for estimation of RE and τ as described in the previous sections, GWP100 can now be estimated using eq 2 and compared with the reported values listed in the Supporting Information. The results are presented in Figure 4. In spite of the significant

(4)

where [OH] is the annual global average OH concentration taken to be 1 × 106 molecules cm−3,59 and α is the empirical adjustable parameter introduced to provide an additional crude correction for various approximations introduced in the estimation procedure (for example, while the average temperature for estimation of τ is typically taken at around 270−275 K, the Atkinson scheme provides estimates for kOH at 298 K). The parameter α was determined by fitting the predicted integrated normalized concentrations defined by eq 3 to the corresponding values derived from the reported τ listed in Supporting Information. As with RE, this was done to impose a proper weighting for actual evaluation of GWP100 via eq 2. The resulting value of α was found to be 2.1 and did not show significant sensitivity to molecular diversity. The comparison of predicted and reported c̃ is shown in Figure 3. As seen, the results exhibit significant scatter, but the overall trend is

Figure 4. Comparison of estimated and reported GWP100. The shaded area represents the logarithmic RMSD.

scatter, the overall performance of the present estimation procedure is quite reasonable, covering the range of over 4 orders of magnitude. The uncertainty mainly comes from lifetime estimates (see Figure 3), and the worst outlier in the bottom right corner is CF3I, as discussed in the previous section. Excluding this outlier, the logarithmic RMSD for GWP100 estimate corresponds to a factor of 3.0.

■

SCREENING FOR LOW GWP100 The fast estimation procedures described in the previous sections allow use of large molecular libraries for screening. We start with the public domain PubChem database60 of the National Institutes of Health. The number of compounds in PubChem is reduced to a more tractable 56 000 candidates by considering only molecules with 15 or fewer atoms and comprised only of the elements C, H, F, Cl, Br, N, and S.

Figure 3. Comparison of integrated normalized concentrations predicted using eq 4 (α = 2.1) and derived from the reported values of τ for TH = 100 years. The shaded area represents the logarithmic RMSD. 12541



Article

compounds from the present collection have an estimated GWP100 below 200. This can be explained by considering another feature seen in Figure 5. With few exceptions, variations in RE/W are confined to a relatively narrow range, and significant variations in GWP100 are mainly controlled by atmospheric lifetime, consistent with recent suggestions.1 When large numbers of compounds are examined, constrained only by valency rules, the resistance toward reactions with OH is a relatively rare feature, and the majority of compounds exhibit short atmospheric lifetimes. On the other hand, high reactivity toward OH is expected to correlate with many undesirable properties (i.e., toxicity, flammability, tendency to polymerize, lack of thermal stability, incompatibility with materials, etc.). Also of note is the fact that in the data set used for validation (see the Supporting Information), which contains presently used as well as newly proposed refrigerants, only 29% of compounds have a GWP100 below 200. It appears that the selection of refrigerants based on avoiding various undesired properties (in addition to their performance in a refrigeration cycle) is more likely to result in long-lived, high-GWP compounds. The present result also suggests that screening efforts focused exclusively on radiative efficiency2 can produce only limited improvements, considering the statistical dominance of the effect associated with lifetime. An analysis of compounds with the lowest values of GWP100 shows two main groups. Several compounds have zero values of GWP100 (not shown in Figure 5 due to logarithmic scaling) caused by effectively zero values of RE. These are diatomic molecules assembled from the elements considered (i.e., N2, O2, HF, F2, etc.). Nongreenhouse air components such as N2 are clearly superior to any other options from the environmental point of view, but they make very poor refrigerants in the traditional vapor compression refrigeration cycle. The same applies to all other diatomics that also exhibit additional undesirable features (e.g., hydrohalogenic acids). The next group of compounds that follows the diatomics (the scattered points in the bottom left corner of Figure 5, outside of the main cluster) include cyclic molecules containing multiple sulfur atoms: tetrasulfur, hexathiane, 2,4,6,8,9,10-hexathiaadamantane, 1,3,5-trithiane-2,4,6-trithiol, 1,2,5,6-tetrathiocine, etc. From a formal point of view, these compounds have the most desirable combination of properties to achieve low values of GWP100: the presence of multiple connected sulfur atoms moves vibrational frequencies outside of the atmospheric window and increases molecular weight, both resulting in a lower RE/W. Additionally, the Atkinson scheme suggests very high reaction rates with OH for sulfur sites. As a result, these compounds exhibit the lowest values of GWP100 after the diatomics. However, all of them are clearly not practical as refrigerants. In addition to obvious environmental concerns, they are expected to have very high melting temperatures and to become unstable at elevated temperatures. Filtering Based on Additional Property Constraints. Based on the presented evidence, GWP100 alone does not provide a sufficient constraint for compound selection. Filtering based on other properties needs to be applied to obtain a pool of potential candidates of manageable size. In fact, virtual screening would be more efficient if more constraining and less computationally expensive filters are applied first; screening for low GWP should be one of the last filters in future studies of this kind to avoid molecular geometry optimization for all the candidates from the library. Several additional filters were

Radicals, ions, ionic compounds, and compounds containing specific atomic isotopes were excluded. The refrigerants in current commercial use are all small molecules. The largest molecule listed in the American Society of Heating, Refrigerating and Air-Conditioning Engineers (ASHRAE) classification standard for refrigerants,61 for example, is pentane, with 17 atoms, but this fluid is used primarily as a minor constituent in blends. A study by McLinden62 provides a thermodynamic basis for preferring small molecules: the optimum range of critical temperature was 340 K < Tc < 470 K, with values of the ideal-gas heat capacity in the range 35 J mol−1·K−1 < c0p < 210 J mol−1·K−1 (depending on which variation on the basic vapor−compression cycle was considered). Such values are observed for alkanes with three to five carbons or halogenated alkanes with three or fewer carbons. The choice of elements traces back to Midgley,63 who observed that only a small portion of the periodic table would form compounds sufficiently volatile to serve as a refrigerant. It is likewise confirmed by the ASHRAE standard,61 where the only exceptions among the listed refrigerants are helium, neon, and argon, which are used in cryogenic refrigeration systems but have boiling points and critical temperatures that are too low for normal refrigeration applications. For all compounds in the collection, RE and τ were estimated as described above. The results are presented in Figure 5. The estimated values are plotted in logarithmic

Figure 5. Estimated GWP100 for the library of over 56 000 compounds collected from PubChem.60

coordinates of RE/W and c̃ (see eqs 2 and 3). In this coordinate system, GWP100 isolines become simple straight lines, as shown in Figure 5. Following a GWP100 isoline from left to right leads to a decrease in the lifetime with simultaneous increase in RE/W, thus keeping GWP100 constant. This presentation clarifies the relative effects of these two factors defining GWP. An immediate and somewhat surprising observation can be made: an overwhelming majority of compounds have low GWP100 values (note that values of GWP100 below 1 have to be viewed with skepticism, because the indirect effects associated with oxidation products are neglected when making estimates). We selected a target value of GWP100 < 200 based on the current policy outlook,64 including the so-called Mobile directive regulation of the European Union,65 which mandates fluids with GWP100 < 150 in automotive air-conditioning and a joint United States/Canada/Mexico proposal to the Montreal Protocol66 that would phase down the use of HFCs to 15% of 2004−2006 consumption levels by 2033. Almost 94% of 12542



Article

applied sequentially, as described next, and the progress of compound elimination from the library is given in Table 1. Table 1. Filtering Criteria and the Corresponding Decrease in the Number of Candidate Compounds step

applied filter

constraint

resulting compound count

0 1 2 3 4 5

GWP100 toxicity flammability critical temperature stability

GWP100 < 200 see text LFL > 0.1 kg·m−3 300 K < Tc < 550 K see text

56203 52565 30135 20277 1728 1234

Figure 6. Correlation between LFL and heat of combustion computed at the PM6 level. Points: experimental data.72 Line: empirical fit, LFL[kg·m−3] = 2.38 × (ΔHc°[MJ·kg−1])−1.19. Shaded area represents the logarithmic RMSD.

Safety codes (e.g., ref 67) require refrigerants of low toxicity in residential and most commercial applications. A cursory toxicity filtering was performed based on a list of markers and associated rules for the elimination of compounds that contain them compiled by Lagorce et al.68 Although this test is by no means comprehensive, it allows elimination of some obviously toxic compounds and can be performed very rapidly with an inhouse utility built upon the OpenBabel library.49 A significant reduction in the number of compounds was observed, as indicated in Table 1. The next filtering was based on each candidates’ flammability, as expressed by the lower flammability limit (LFL). The ideal refrigerant would be nonflammable, but moderate flammability may be tolerated. We take as “moderate” a lower flammability limit greater than 0.1 kg m−3, which corresponds to the boundary of the “Class 2” flammability rating in the ASHRAE refrigerant classification standard.61 A number of methods for estimation of LFL of varying scope and accuracy are available in the literature.69 As mentioned previously, the computational speed and the scope of coverage take the precedence in choosing the estimation method for LFL. McLinden70 has shown that experimental LFL correlates well with experimental heat of combustion when both are expressed on a mass basis. While the fact that LFL and heat of combustion are correlated is well-known, both are usually given in volumetric units, and in this case, correlations often involve the molecular weight as an additional variable.71 An advantage of using only the heat of combustion for LFL estimation is that there are heats of formation computed at the PM6 level already available for all compounds in the library (obtained in the course of optimization for IR spectra calculations), and PM6computed heats of combustion can be obtained with essentially no additional computational cost. Although these values are not of chemical accuracy, they may, nevertheless, be sufficient for the use in a correlation for LFL. To verify this, a set of experimental LFLs was collected from the DIPPR 801 compilation72 (468 total, only values labeled as experimental were taken). For all selected compounds, heats of formation were obtained at the PM6 level. The corresponding heats of combustion were evaluated based on combustion products as specified in the ASHRAE standard61 and using heats of formations for O2 and reaction products also computed using PM6. The experimental lower flammability limits are plotted against computed heats of combustion in Figure 6. As seen, the correlation observed for experimental heats of combustion70 also holds for the PM6-predicted ones. The experimental points are reasonably described with an empirical power-law function, LFL[kg·m−3] = 2.38 × (ΔHc°[MJ·kg−1])−1.19, with a

logarithmic RMSD corresponding to a factor of 1.24. Several outliers are present; some of them, such as hydrogen sulfide and carbon disulfide, also exhibited deviations from a correlation when experimental heats of combustion were used.70 When the correlation described was applied to filter out flammable compounds with the cutoff value of 0.1 kg m−3, the size of the candidate pool was reduced by almost 10 000 (Table 1). The next step of compound elimination was based on a constraint on the critical temperature, Tc. We use the critical temperature as the primary thermodynamic criterion for fluid selection. The lower bound is 300 K; this is the approximate lowest operating temperature of the condenser in the vapor compression cycle; fluids with lower values of Tc would operate in a supercritical cycle. Supercritical cycles are under active investigation, typically using carbon dioxide as the refrigerant, but these systems operate at very high pressures and require sophisticated control strategies to avoid efficiency penalties. As the refrigerant Tc increases, the cycle operates at lower reduced temperatures and pressures; this generally results in higher efficiencies,62 but the volumetric capacity (refrigeration effect per unit of volume of refrigerant vapor flowing into the compressor) decreases, requiring physically larger equipment. Operation at subatmospheric pressures is problematic, because air leakage into the sealed system would reduce heat transfer coefficients and, for flammable refrigerants, create an explosion hazard. For these reasons, the most commonly used refrigerants are in the range 343 K < Tc < 385 K. Systems with centrifugaltype compressors (i.e., turbomachinery) often use fluids having higher critical temperatures, up to about 471 K. We adopt 550 K as a generous upper limit on critical temperature, so as not to prematurely exclude otherwise promising candidates. The quantitative structure−property relationships (QSPR) method developed previously34 utilizing the NIST ThermoData Engine73,74 technology was used to estimate Tc for the remaining compounds. As a result of this filtering, the number of remaining compounds was also reduced quite significantly, by more than an order of magnitude (Table 1). Although the total count is still fairly large, it is of a size amenable to more direct analysis. Visual inspection of the remaining compounds revealed that there was a significant portion of them that may have stability problems. A refrigerant is typically expected to operate for many years in a sealed system, so unstable compounds are not 12543



Article

suitable. To remedy this, additional filtering rules were applied. Compounds with any of the following functional groups were eliminated: (1) triple bond, (2) nitro (NO2) group, (3) peroxide (OO) group, (4) 3- and 4-member rings, (5) disulfide (SS) group, (6) CS group, (7) linear diene (CCCC), (8) NO and NN groups (any bond order), (9) NX and OX groups, where X = F, Cl, or Br, (10) C group, and (11) groups exhibiting keto−enol tautomerism. This final round of filtering resulted in a list of 1234 compounds.

high GWP100 and high Tc, occupying the upper right corner of the graph. This is the result of competing constraints applied during filtering. To qualitatively understand this phenomenon, a representative example of a homologous series of fluorinated propanes is used (top panel of Figure 8). Shown in the figure is

■

RESULTS AND DISCUSSION The complete list of 1234 compounds remaining after filtering is given in the Supporting Information. Here, to facilitate the discussion, the compounds are divided into chemical classes, and the results are presented in Figure 7 in Tc-GWP100

Figure 8. Estimated critical temperatures for a homologous series of fluorinated propanes (top) and propenes (bottom). ○ compounds with GWP100 < 200, × compounds with GWP100 > 200; shaded areas represent cutoffs due to flammability constraint, LFL > 0.1 kg·m−3.

the estimated critical temperature as a function of the number of hydrogen atoms in the molecule. As seen, the global trend exhibits an initial increase of Tc with the number of hydrogens going from C3F8, reaching a maximum, and, finally, a decrease to C3H8. The final part is truncated due to the LFL constraint, leaving only the initial, increasing part. To have a low GWP value, the compounds need to have as many C−H bonds as possible to increase the number of sites available to react with OH, thus reducing the lifetime. According to the present estimate, fluorinated propane should have at least 3−4 hydrogens to meet GWP100 < 200 criterion (Figure 8). This, in turn, causes an increase in the critical temperature (in excess of 400 K for this example). As a result of these two competing processes, the remaining halogenated alkanes have both high GWP100 and high Tc. This “hydrogen penalty”, that is, the increase in Tc caused by the increase in a number of hydrogen atoms, is likely to be one of the main factors controlling the chemical nature of halogenated low-GWP refrigerants. Cyclic alkanes and aromatics are represented by halogenated derivatives of toluene, benzene, and cyclopentane. They have, on average, slightly lower GWP100 and higher critical temperatures as compared to linear alkanes. Halogenated derivatives of linear olefins is the largest group among the potential candidates, over a third of the entire set. Because the dominant mechanism of their reaction with OH is via addition to a double bond rather than H-abstraction, their GWP100 values are not affected by the “hydrogen penalty” observed for alkanes. As shown in the bottom panel of Figure 8, all fluorinated propenes, including C3F6, meet the GWP100 < 200 criterion. Consequently, critical temperatures of linear olefins span over the entire range of interest. Simultaneously, they also show lower values of GWP100 as compared to alkanes and aromatics. This explains the present industry focus on

Figure 7. Compounds remaining after filtering. Top panel, halogenated hydrocarbons: red circle, linear alkanes; green triangle, cyclic alkanes and aromatics; blue square, linear olefins; black triangle, cyclic olefins. Middle panel, halogenated oxygenates: red circle, linear ethers; green triangle, cyclic ethers; blue square, alcohols, black triangle, contain both −O− and −OH groups. Bottom panel, halogenated compounds that contain N or S: red circle, contain N (linear); green triangle, contain N (cyclic or aromatic); blue square, contain S; black triangle, contain both N and S.

coordinates. The vast majority of the candidates are halogenated; all of the hydrocarbons, as well as the unsubstituted ethers and alcohols, were filtered out, mainly by the flammability criteria. Only six compounds out of the 1234 do not contain halogens; all of them contain N and/or S, and all have critical temperatures in excess of 470 K. Over 60% of the halogenated compounds contain only fluorine, because the addition of heavier Cl or Br atoms generally causes a large increase in critical temperature, often resulting in violation of the constraint on Tc. The compounds from the first and most abundant chemical class (over a half of the entire set), halogenated hydrocarbons, are presented in the top panel of Figure 7. The class is further subdivided into linear alkanes, cyclic alkanes and aromatics, linear olefins, and cyclic olefins. Linear alkanes generally exhibit 12544



Article

HFOs, since they are capable of achieving lower Tc compatible with current refrigeration equipment while maintaining acceptable GWP100. The last subclass, cyclic olefins, is mainly composed of halogenated derivatives of cyclopentene and generally shows lower GWP100, but higher Tc, as compared with their linear counterparts. The next chemical class of compounds present in the set is halogenated oxygenates, plotted in the middle panel of Figure 7. This class is further subdivided into linear ethers, cyclic ethers, alcohols, and compounds exhibiting both −O− and −OH functional groups. Linear ethers are most abundant within this class, covering a wide range of acceptable GWP100 and reaching Tc below 400 K. However, an examination of molecular structures of ethers with low Tc reveals that they are also olefins. In general, the low GWP100 values exhibited by ethers are expected to be due to higher H-abstraction rates from carbons adjacent to oxygen (as reflected by a substituent factor of 6.1 per C−H site in the Atkinson scheme used51). This is apparently not the case for the low-Tc ethers that are mostly fully fluorinated (the lowest estimated Tc = 352 K for ethers is shown by CF3OCFCF2 identified previously75 as a potential low-GWP refrigerant), and it is the presence of the double bond that causes their high reaction rate with OH and low GWP100. Therefore, these compounds are effectively olefins from the point of view of the present analysis. The lowest Tc for a linear ether without a double bond is 413 K, for CF3OCH2CH3. A similar situation occurs for cyclic ethers: many of them have double bonds, either in the ring or in the side chain. The presence of an −OH group in the two remaining subclasses results in both high Tc and GWP100, displaying the pattern similar to that observed for linear alkanes. These compounds are highly polar, which tends to produce elevated Tc. The lowest Tc among alcohols is exhibited by CF3OH, about 424 K. The last chemical class shown in the bottom panel of Figure 7 includes compounds that contain nitrogen or sulfur. Apart from ammonia (NH3), which is very significant industrial refrigerant, none of the common refrigerants contain N or S, yet these heteroatoms offer another alternative to increase rates of compound degradation, without the “hydrogen penalty” or having double bonds. Reaction rates of OH with nitrogen and, especially, sulfur sites are very high.51 On the other hand, the presence of these heteroatoms also tends to increase critical temperature. The subclasses in this case include compounds containing N (linear), containing N (cyclic or aromatic), containing S, and containing both N and S. Compounds containing N in a linear configuration are most prevalent in this class and include mostly amines. They show a wide range of acceptable GWP100 values and somewhat elevated Tc, as expected. Two of them, N(CF3)(CF3)(CFCF2) and N(CF3)(CF3)(CF2H) that have Tc < 400 K feature a tertiary amine group and are highly fluorinated (thus avoiding the “hydrogen penalty”); one of them also has double bond. Similar to other chemical classes, cyclic N-containing compounds generally have higher Tc and include five- or six-member rings, both aliphatic and aromatic. Sulfur-containing compounds are represented by thiols and thioethers. As noted previously, OH reaction rates with these S sites are the highest among all other functional groups encountered. Consequently, one can achieve very low GWP100 values just by adding one of these groups. As with nitrogen, three compounds with Tc < 400 K are fully fluorinated: one thiol, CF3SH, and two thioethers, CF3SCF3 and CF3SCF2CF3. In fact, this subclass contains a significant

number of perfluorocompounds that are modified with either a −SH or a −S− group, which results in low GWP100 values. Nine compounds have both N and S present; on average, they exhibit higher Tc values, as expected. Finally, it should be noted that compounds with N and S are generally viewed as problematic due to obvious environmental concerns, materials compatibility (amines and thiols are known to be corrosive), and other issues (e.g., the odor, in case of thiols). The present results do not consider the performance (i.e., efficiency and capacity) of a refrigerant in the vapor compression cycle, except for the very coarse screening parameter of critical temperature, and none of the candidates identified here should be rejected until these aspects are explored. The present work is the first phase of a continuing effort, and future work will include estimating additional thermodynamic parameters (including critical pressure and ideal-gas heat capacity), so that the performance of these fluids can be estimated. These results will be combined with those of our companion study,76 which will determine the optimal values of fundamental thermodynamic parameters, to identify a much smaller set (order of 20) fluids for detailed study and, perhaps, experimental verification.

■

SUMMARY

A method for estimation of global warming potential for a time horizon of 100 years was developed and successfully validated against available literature data. The method relies on separate estimates for radiative efficiency and atmospheric lifetime, both derived from molecular structures, and provides accuracy and computational speed, which are sufficient for virtual screening applications. When applied to a library of over 56000 compounds, it was found that a screening based on GWP alone does not allow reducing the candidate pool to a reasonable size, and additional filtering criteria were needed. When filtering based on toxicity, flammability, critical temperature, and stability were added, the number of remaining candidates was decreased to slightly over 1200. Investigation of the chemical nature of the remaining compounds revealed the following. The requirement of low flammability generally results in highly halogenated species. Their atmospheric destruction via reactions with the hydroxyl radical can proceed via hydrogen abstraction, addition to a double bond or aromatic ring, or attack on a heteroatom (N or S) site. A high reaction rate is required to achieve low GWP values. If H-abstraction is a dominant mechanism, a significant number of C−H sites must be present in the molecule. For highly halogenated compounds, the addition of hydrogen also results in a “hydrogen penalty” on the critical temperature; that is, an increase in Tc often makes the compound incompatible with current refrigeration equipment that can only accommodate fluids with Tc below about 400 K. Introduction of CC bonds and activation of a destruction mechanism via OH addition to a double bond eliminates the “hydrogen penalty” and makes halogenated olefins most promising candidates for the use with current refrigeration equipment. Introduction of heteroatoms, N and, especially, S, also results in potentially acceptable candidates. Future work will include estimating additional thermodynamic parameters so that the performance of these fluids can be estimated. 12545



■

Article

(11) Manning, M.; Reisinger, A. Broader perspectives for comparing different greenhouse gases. Phil. Trans. R. Soc. A 2011, 369, 1891− 1905. (12) Pinnock, S.; Hurley, M. D.; Shine, K. P.; Wallington, T. J.; Smyth, T. J. Radiative forcing of climate by hydrochlorofluorocarbons and hydrofluorocarbons. J. Geophys. Res. 1995, 100, 23227−23238. (13) Jain, A. K.; Briegleb, B. K.; Minschwaner, K.; Wuebbles, D. Radiative forcings and global warming potentials of 39 greenhouse gases. J. Geophys. Res. 2000, 105, 20773−20790. (14) Patten, K. O.; Khamaganov, V. G.; Orkin, V. L.; Baughcum, S. L.; Wuebbles, D. J. OH reaction rate constant, IR absorption spectrum, ozone depletion potentials and global warming potentials of 2-bromo3,3,3-trifluoropropene. J. Geophys. Res. 2011, 116, D24307. (15) Sihra, K.; Hurley, M. D.; Shine, K. P.; Wallington, T. J. Updated radiative forcing estimates of 65 halocarbons and nonmethane hydrocarbons. J. Geophys. Res. 2001, 106, 20493−20505. (16) Scientific Assessment of Ozone Depletion: 1998. World Meteorological Organization Global Ozone Research and Monitoring ProjectReport No. 44; World Meteorological Organization: Geneva, Switzerland, 1999. (17) Scientific Assessment of Ozone Depletion: 2010. Global Ozone Research and Monitoring ProjectReport No. 52; World Meteorological Organization: Geneva, Switzerland, 2011. (18) Nielsen, O. J.; Javadi, M. S.; Andersen, M. P. S.; Hurley, M. D.; Wallington, T. J.; Singh, R. Atmospheric chemistry of CF3CFCH2: kinetics and mechanisms of gas-phase reactions with Cl atoms, OH radicals, and O3. Chem. Phys. Lett. 2007, 439, 18−22. (19) Orkin, V. L.; Martynova, L. E.; Ilichev, A. N. High-Accuracy measurements of OH reaction rate constants and IR absorption spectra: CH2CFCF3 and trans-CHFCHCF3. J. Phys. Chem. A 2010, 114, 5967−5979. (20) Jiménez, E.; Antiõ lo, M.; Ballesteros, B.; Martínez, E.; Albaladejo, J. Atmospheric lifetimes and global warming potentials of CF3CH2CH2OH and CF3(CH2)2CH2OH. ChemPhysChem 2010, 11, 4079−4087. (21) Cometto, P. M.; Taccone, R. A.; Nieto, J. D.; Dalmasso, P. R.; Lane, S. I. Kinetic study of OH radical reactions with CF3CClCCl2, CF3CClCClCF3 and CF3CFCFCF3. ChemPhysChem 2010, 11, 4053−4059. (22) Thomsen, D. L.; Andersen, V. F.; Nielsen, O. J.; Wallington, T. J. Atmospheric chemistry of C2F5CH2OCH3 (HFE-365mcf). Phys. Chem. Chem. Phys. 2011, 13, 2758−2764. (23) Baasandorj, M.; Ravishankara, A. R.; Burkholder, J. B. Atmospheric chemistry of (Z)CF3CHCHCF3: OH radical reaction rate coefficient and global warming potential. J. Phys. Chem. A 2011, 115, 10539−10549. (24) Østerstrøm, F. F.; Nielsen, O. J.; Andersen, M. P. S.; Wallington, T. J. Atmospheric chemistry of CF3CH2OCH3: reaction with chlorine atoms and OH radicals, kinetics, degradation mechanism and global warming potential. Chem. Phys. Lett. 2012, 524, 32−37. (25) Young, C. J.; Hurley, M. D.; Wallington, T. J.; Mabury, S. A. Molecular structure and radiative efficiency of fluorinated ethers: A structure−activity relationship. J. Geophys. Res. 2008, 113, D24301. (26) Bravo, I.; Aranda, A.; Hurley, M. D.; Marston, G.; Nutt, D. R. Infrared absorption spectra, radiative efficiencies, and global warming potentials of perfluorocarbons: Comparison between experiment and theory. J. Geophys. Res. 2010, 115, D24317. (27) Papasavva, S.; Tai, S.; Illinger, K. H.; Kenny, J. E. Infrared radiative forcing of CFC substitutes and their atmospheric reaction products. J. Geophys. Res. 1997, 102, 13643−13650. (28) Bravo, I.; Díaz-de-Mera, Y.; Aranda, A.; Smith, K.; Shine, K. P.; Marston, G. Atmospheric chemistry of C4F9OC2H5 (HFE-7200), C4F9OCH3 (HFE-7100), C3F7OCH3 (HFE-7000) and C3F7CH2OH: Temperature dependence of the kinetics of their reactions with OH radicals, atmospheric lifetimes, and global warming potentials. Phys. Chem. Chem. Phys. 2010, 12, 5115−5125. (29) Papasavva, S.; Tai, S.; Esslinger, A.; Illinger, K. H.; Kenny, J. E. Ab initio calculations of vibrational frequencies and infrared intensities

ASSOCIATED CONTENT

S Supporting Information *

Data set used for validation of the methods for GWP estimation and final list of compounds with associated screening parameters. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes †

Contribution of the U.S. National Institute of Standards and Technology and not subject to copyright in the United States. Trade names are provided only to specify procedures adequately and do not imply endorsement by the National Institute of Standards and Technology. Similar products by other manufacturers may be found to work as well or better. The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Dr. L. S. Rothman of Harvard-Smithsonian Center for Astrophysics and Dr. S. W. Sharpe of Pacific Northwest National Laboratory for providing access to HITRAN2008 and Northwest Infrared Library, respectively. Helpful discussions with Dr. V. Diky and Dr. J. A. Widegren of NIST are also greatly appreciated. This work was supported by the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy, under Grant No. DE-EE0002057.

■

REFERENCES

(1) Velders, G. J. M.; Ravishankara, A. R.; Miller, M. K.; Molina, M. J.; Alcamo, J.; Daniel, J. S.; Fahey, D. W.; Montzka, S. A.; Reimann, S. Preserving Montreal Protocol climate benefits by limiting HFCs. Science 2012, 335, 922−923. (2) Bera, P. P.; Francisco, J. S.; Lee, T. J. Design strategies to minimize the radiative efficiency of global warming molecules. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 9049−9054. (3) Shoichet, B. K. Virtual screening of chemical libraries. Nature 2004, 432, 862−865. (4) O’Boyle, N. M.; Campbell, C. M.; Hutchison, G. R. Computational design and selection of optimal organic photovoltaic materials. J. Phys. Chem. C 2011, 115, 16200−16210. (5) Materials Genome Initiative for Global Competitiveness; National Science and Technology Council: Washington, DC, 2011. (6) Climate Change: The IPCC Scientific Assessment; Houghton, J., Jenkins, G., Ephraums, J., Eds.; Cambridge University Press: New York, 1990. (7) Climate Change 1995: The Science of Climate Change; Houghton, J., Meira Filho, L., Callander, B., Harris, N., Kattenberg, A., Maskell, K., Eds.; Cambridge University Press: New York, 1996. (8) Climate Change 2001: The Scientific Basis. Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change; Houghton, J. T., Ding, Y., Griggs, D. J., Noguer, M., van der Linden, P. J., Dai, X., Maskell, K., Johnson, C. A., Eds.; Cambridge University Press: New York, 2001. (9) Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change; Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., B., K., Tignor, M., Miller, H. L., Eds.; Cambridge University Press: New York, 2007. (10) Fuglestvedt, J. S.; Berntsen, T. K.; Godal, O.; Sausen, R.; Shine, K. P.; Skodvin, T. Metrics of climate change: Assessing radiative forcing and emission indices. Clim. Change 2003, 58, 267−331. 12546



Article

structure−reactivity relationship: An update. Atmos. Environ. 1995, 29, 1685−1695. (52) Klamt, A. Estimation of gas-phase hydroxyl radical rate constants of organic compounds from molecular orbital calculations. Chemosphere 1993, 26, 1273−1289. (53) Klamt, A. Estimation of gas-phase hydroxyl radical rate constants of oxygenated compounds based on molecular orbital calculations. Chemosphere 1996, 32, 717−726. (54) Böhnhardt, A.; Kühne, R.; Ebert, R.-U.; Schüürmann, G. Indirect photolysis of organic compounds: Prediction of OH reaction rate constants through molecular orbital calculations. J. Phys. Chem. A 2008, 112, 11391−11399. (55) Böhnhardt, A.; Kühne, R.; Ebert, R.-U.; Schüürmann, G. Predicting rate constants of OH radical reactions with organic substances: Advances for oxygenated organics through a molecular orbital HF/6-31G** approach. Theor. Chem. Acc. 2010, 127, 355−367. (56) Meylan, W. M.; Howard, P. H. Computer estimation of the atmospheric gas-phase reaction rate of organic compounds with hydroxyl radicals and ozone. Chemosphere 1993, 26, 2293−2299. (57) U.S. Environmental Protection Agency. The Atmospheric Oxidation Program for Microsoft Windows (AOPWIN) v1.92a; U.S. EPA: Washington, DC, 2010. (58) Sander, S. P.; Abbatt, J.; Barker, J. R.; Burkholder, J. B.; Friedl, R. R.; Golden, D. M.; Huie, R. E.; Kolb, C. E.; Kurylo, M. J.; Moortgat, G. K.; Orkin, V. L.; Wine, P. H. Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies, Evaluation No. 17, JPL Publication No. 10-6; NASA Jet Propulsion Laboratory: Pasadena, CA, 2011. (59) Montzka, S. A.; Spivakovsky, C. M.; Butler, J. H.; Elkins, J. W.; Lock, L. T.; Mondeel, D. J. New observational constraints for atmospheric hydroxyl on global and hemispheric scales. Science 2000, 288, 500−503. (60) Bolton, E.; Wang, Y.; Thiessen, P. A.; Bryant, S. H. PubChem: Integrated platform of small molecules and biological activities. In Annual Reports in Computational Chemistry; American Chemical Society: Washington, DC, 2008; Vol. 4, Chapter 12. (61) ANSI/ASHRAE Standard 34-2010: Designation and Safety Classification of Refrigerants; American Society of Heating, Refrigerating, and Air-Conditioning Engineers, Inc.: Atlanta, GA, 2010. (62) McLinden, M. O. Optimum refrigerants for non-ideal cycles: An analysis employing corresponding states. USNC/IIR-Purdue Refrigeration Conference and ASHRAE-Purdue CFC Conference, West Lafayette, IN, 1990; pp 69−79. (63) Midgley, T. From the periodic table to production. Ind. Eng. Chem. 1937, 29, 241−244. (64) Mascarelli, A. L. A bright future for the Montreal Protocol. Environ. Sci. Technol. 2010, 44, 1518−1520. (65) Directive 2006/40/EC of the European Parliament and of the Council of 17 May 2006 relating to emissions from air conditioning systems in motor vehicles and amending Council Directive 70/156/EEC; 2006. Available online: http://eur-lex.europa.eu/JOHtml.do? uri= OJ:L:2006:161:SOM:EN:HTML (accessed March 28, 2012). (66) Proposal by United States, Canada and Mexico to amend the Montreal Protocol; 2009. Available online: http://ozone.unep.org/ Meeting_Documents/mop/21mop/ MOP-21-3-Add-1E.pdf (accessed March 28, 2012). (67) ANSI/ASHRAE Standard 15-2010: Safety Standard for Refrigeration Systems; American Society of Heating, Refrigerating, and Air-Conditioning Engineers, Inc.: Atlanta, GA, 2010. (68) Lagorce, D.; Sperandio, O.; Galons, H.; Miteva, M. A.; Villoutreix, B. D. FAF-Drugs2: Free ADME/tox filtering tool to assist drug discovery and chemical biology projects. BMC Bioinf. 2008, 9, 396. (69) Vidal, M.; Rogers, W.; Holste, J.; Mannan, M. A review of estimation methods for flash points and flammability limits. Process Saf. Prog. 2004, 23, 47−55. (70) McLinden, M. O. unpublished results, 1987. (71) Suzuki, T.; Ishida, M. Neural network techniques applied to predict flammability limits of organic compounds. Fire Mater. 1995, 19, 179−189.

for global warming potential of CFC substitutes: CF3CH2F (HFC134a). J. Phys. Chem. 1995, 99, 3438−3443. (30) Rocha, G. B.; Freire, R. O.; Simas, A. M.; Stewart, J. J. P. RM1: A reparameterization of AM1 for H, C, N, O, P, S, F, Cl, Br, and I. J. Comput. Chem. 2006, 27, 1101−1111. (31) Stewart, J. J. P. Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements. J. Mol. Model. 2007, 13, 1173−1213. (32) IPCC Fourth Assessment Report. Climate Change 2007: Working Group I: The Physical Science Basis. Errata.; 2011. Available online: http://www.ipcc.ch/publications_and_ data/ar4/wg1/en/ errataserrata-errata.html (accessed March 28, 2012). (33) Frisch, M. J. et al. Gaussian 09, Revision B.01; Gaussian, Inc.: Wallingford, CT, 2010. (34) Kazakov, A.; Muzny, C. D.; Diky, V.; Chirico, R. D.; Frenkel, M. Predictive correlations based on large experimental datasets: Critical constants for pure compounds. Fluid Phase Equilib. 2010, 298, 131− 142. (35) Sadlej, A. J. Medium-size polarized basis sets for high-level correlated calculations of molecular electric properties. Collect. Czech. Chem. Commun. 1998, 53, 1995−2016. (36) Sadlej, A. J. Medium-size polarized basis sets for high-levelcorrelated calculations of molecular electric properties II. Second-row atoms: Si through Cl. Theor. Chem. Acc. 1991, 79, 123−140. (37) Halls, M.; Velkovski, J.; Schlegel, H. Harmonic frequency scaling factors for Hartree-Fock, S-VWN, B-LYP, B3-LYP, B3-PW91 and MP2 with the Sadlej pVTZ electric property basis set. Theor. Chem. Acc. 2001, 105, 413−421. (38) Zvereva, E. E.; Shagidullin, A. R.; Katsyuba, S. A. Ab initio and DFT predictions of infrared intensities and Raman activities. J. Phys. Chem. A 2011, 115, 63−69. (39) Price, K. V.; Storn, R. M.; Lampinen, J. A. Differential Evolution: A Practical Approach to Global Optimization; Springer-Verlag: Heidelberg, 2005. (40) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipies: The Art of Scientific Computing, 3rd ed.; Cambridge University Press: New York, 2007. (41) Rothman, L. S.; et al. The HITRAN 2008 molecular spectroscopic database. J. Quant. Spectrosc. Radiat. Transfer 2009, 110, 533−572. (42) Sharpe, S. W.; Johnson, T. J.; Sams, R. L.; Chu, P. M.; Rhoderick, G. C.; Johnson, P. A. Gas phase databases for quantitative infrared spectroscopy. Appl. Spectrosc. 2004, 58, 1452−1461. (43) Chu, P. M.; Guenther, R. R.; Rhoderick, G. C.; Lafferty, W. J. The NIST quantitative infrared database. J. Res. Natl. Inst. Stand. Technol. 1999, 104, 59−81. (44) Myhre, G.; Highwood, E. J.; Shine, K. P.; Strodal, F. New estimates of radiative forcing due to well mixed greenhouse gases. Geophys. Res. Lett. 1998, 25, 2715−2718. (45) Fisher, D. A.; Hales, C. H.; Wang, W.-C.; Ko, M. K. W.; Sze, N. D. Model calculations of the relative efects of CFCs and their replacements on global warming. Nature 1990, 344, 513−516. (46) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. Development and use of quantum mechanical molecular models. 76. AM1: A new general purpose quantum mechanical molecular model. J. Am. Chem. Soc. 1985, 107, 3902−3909. (47) Stewart, J. J. P. Optimization of parameters for semiempirical methods. 1. Method. J. Comput. Chem. 1989, 10, 209−220. (48) Vladislavleva, E.; Smits, G.; den Hertog, D. On the importance of data balancing for symbolic regression. IEEE Trans. Evol. Comput. 2010, 14, 252−277. (49) O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchinson, G. R. Open Babel: An open chemical toolbox. J. Cheminf. 2011, 3, 33. (50) Sumathi, R.; Green, W. H., Jr. A priori rate constants for kinetic modeling. Theor. Chem. Acc. 2002, 108, 187−213. (51) Kwok, E. S. C.; Atkinson, R. Estimation of hydroxyl radical reaction rate constants for gas phase organic compounds using a 12547



Article

(72) Evaluated Standard Thermophysical Property Values, DIPPR 801; Design Institute for Physical Properties/American Institute of Chemical Engineers: New York, 2011. (73) Frenkel, M.; Chirico, R. D.; Diky, V.; Muzny, C. D.; Kazakov, A. F.; Magee, J. W.; Abdulagatov, I. M.; Kroenlein, K.; Kang, J. W. NIST Standard Reference Database 103b. NIST ThermoData Engine Version 6.0: Pure Compounds, Binary Mixtures, Ternary Mixtures, and Chemical Reactions; Standard Reference Data Program; National Institute of Standards and Technology, Gaithersburg, MD, 2012. (74) Diky, V.; Chirico, R. D.; Kazakov, A.; Muzny, C. D.; Frenkel, M. ThermoData Engine (TDE): software implementation of the dynamic data evaluation concept. 4. Chemical reactions. J. Chem. Inf. Model. 2009, 49, 2883−2896. (75) Li, Z.; Tao, Z.; Naik, V.; Good, D. A.; Hansen, J. C.; Jeong, G.R.; Francisco, J. S.; Jain, A. K.; Wuebbles, D. J. Global warming potential assessment CF3OCFCF2. J. Geophys. Res. 2000, 105, 4019−4029. (76) McLinden, M. O.; Domanski, P. A.; Kazakov, A.; Heo, J.; Brown, J. S. Possibilities, limits, and tradeoffs for refrigerants in the vapor compression cycle. ASHRAE/NIST Refrigerants Conference, Gaithersburg, MD, 2012.

12548


Computational Design of New Refrigerant Fluids Based on

Recommend Documents