Application of the Generalized Connectivity-Based Hierarchy to

Article pubs.acs.org/JPCA

Application of the Generalized Connectivity-Based Hierarchy to Biomonomers: Enthalpies of Formation of Cysteine and Methionine Raghunath O. Ramabhadran, Arkajyoti Sengupta, and Krishnan Raghavachari* Department of Chemistry, Indiana University, Bloomington, Indiana 47405, United States S Supporting Information *

ABSTRACT: Computational challenges toward an accurate determination of the enthalpies of formation of amino acids are partly due to the nonavailability of systematic error-canceling thermochemical procedures for such biomonomers. Recently, we developed the connectivity-based hierarchy (CBH) to accurately compute the enthalpies of formations of organic molecules composed of main group elements. Advancing the applicability of CBH to biologically relevant molecules, we have computed the enthalpies of formation of the naturally occurring sulfur-containing amino acids cysteine and methionine which act as fertile testing grounds for the error-canceling ability of thermochemical schemes for biomolecules. We establish herein using the sophisticated error-canceling isoatomic scheme (CBH-2) that relatively inexpensive computational methods with modest basis sets can be used to accurately obtain the enthalpies of formations of the amino acids. Overall, we recommend the use of the isoatomic scheme over the currently popular isodesmic bond separation scheme in future applications in theoretical thermochemistry.

1. INTRODUCTION

To improve upon the isodesmic bond separation scheme and to construct a unique as well as generalized thermochemical hierarchy for organic compounds, we developed the connectivity based hierarchy (CBH) in 2011.12 The hierarchy solves the important problem of creating a reliable and automated computational procedure to accurately obtain the enthalpies of formation of all closed shell organic molecules. Furthermore, CBH has been shown to provide accurate enthalpies of formations with both density functional methods12 as well as wave function-based methods using modest basis sets.13 Further advancement in the applicability of a thermochemical hierarchy from simple organic molecules14,15 to larger biomolecules necessitates the demonstration of superior errorcancellation abilities even when a molecule has a number of different heavy (non-hydrogen) atoms. Biomonomers such as the sulfur containing naturally occurring amino acids cysteine and methionine serve as a useful testing ground in this regard. Besides their obvious biological significance,1 these biomonomers contain four different heavy atoms (C, N, O, and S), the maximum among all of the naturally occurring amino acids. All three heteroatoms (N, O, and S) are located within 2−3 carbon atoms resulting in significant interactions that have to be included. Roux et al.7,8 have discussed the enthalpy of formation of these molecules taking multiple low-energy conformations into account, thereby providing a suitable

The structures, conformations, and thermochemical properties of cysteine, methionine, and other amino acids in the solid phase, in solution, and in the gas phase have been of great scientific interest for decades.1 The solid and solution phase studies play a key role in understanding the crystal structures of various proteins and their physiological functions.1 In the gas phase, as Alonso and co-workers have rightly pointed out recently,2 studies on amino-acids provide valuable structural and thermochemical information on their neutral forms that are free from intermolecular interactions and zwitterionic forms commonly observed in the condensed phases. The derivation of accurate gas phase thermochemical values of the amino acids is challenging due to the lack of reliable experimental data.2−6 Additionally, the conformational flexibility of these molecules leads to a variety of local minima that have to be considered. The absence of gas phase enthalpies of formation for cysteine and many other amino acids in the hugely popular NIST chemistry WebBook underscores this problem. Addressing this issue, Roux et al. and Dixon and coworkers, have recently obtained accurate enthalpies of formations for the sulfur-containing amino acids in the gas phase. 7−9 In Roux et al.’s work, they compare their experimental results with the corresponding values from the popular G3, G3(MP2), and G4 methods10 via atomization and isodesmic bond separation schemes.11 Dixon has used the G3(MP2) method (via atomization and by constructing a clever thermochemical scheme) to obtain the enthalpies of formation for all of the amino acids in the gas phase. © XXXX American Chemical Society

Received: March 29, 2013 Revised: May 20, 2013

A

dx.doi.org/10.1021/jp403123c | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

study the effect of larger basis sets. We have intentionally used these relatively modest sized, practical basis sets to evaluate if they provide sufficient error cancellations. In addition, consistent with our approach of balancing the computational cost of CBH with the desired level of accuracy, we avoid directly performing the computationally expensive CCSD(T)/ 6-311++G(3df,2p) calculations. Instead, we use an additivity approximation to get the CCSD(T)/6-311++G(3df,2p) energies and subsequently use these energies in the CBH schemes to obtain the enthalpies of formations, i.e.

reference point. All of these points put together, these two amino acids provide key insights regarding the errorcancellation ability of CBH for future applications to more complicated biomolecules. This work features the first application of an automated thermochemical hierarchy to compute the enthalpies of formations of amino acids. On the basis of superior error cancellation obtained using the isoatomic scheme (CBH2),16,17 we demonstrate herein that accurate enthalpies of formation of formation of such biologically relevant mediumsized molecules can be obtained with relatively inexpensive computational methods using modest-sized basis sets. The effect of the level of theory employed for geometry optimizations as well as the role played by different low-lying conformations are also discussed. We finally recommend the use of the isoatomic scheme in place of the widely used isodesmic bond separation scheme in future applications in theoretical thermochemistry. For ease of comprehension, a succinct account of CBH is provided herein. More details have been provided elsewhere.12,13 CBH has several rungs, denoted CBH-1, CBH-2, CBH-3, etc. The hierarchy depends only on the connectivity of the atoms and can be derived by an automated computational procedure for any organic molecule. CBH-1 is formally the same as the isodesmic bond separation scheme. At the isoatomic CBH-2 rung, the immediate chemical bonding environment of all the heavy atoms is preserved. At CBH-3, the immediate chemical environment of all of the heavy atom bonds is preserved. Additional higher rungs can also be defined in a similar manner, if required. However for reasons mentioned in the ensuing paragraphs (vide infra), we restrict the discussions to the CBH-1 and CBH-2 rungs, with a specific emphasis on the CBH-2 (isoatomic) rung. It is evident that the basic strategy behind the construction of CBH is to increasingly preserve the chemical environment of all the heavy-atoms present in an organic molecule on ascending up the rungs of the hierarchy. When the bond-types and hybridizations are increasingly matched at the higher rungs, we expect superior error cancellation leading to accurate energies for the reaction schemes. Using these reaction energies and Hess’s law, accurate enthalpies of formation can be computed, provided accurate experimental data are available for all of the reference molecules.

E(CCSD(T)/6‐311++G(3df,2p)) = E(CCSD(T)/6‐31+G(d,p)) + E(MP2/6‐311++G(3df,2p)) − E(MP2/6‐31+G(d,p))

For the purpose of determining the enthalpies of formations, the single point energies obtained using these different methods were added to the thermally corrected enthalpies at 298.15 K. To take multiple conformations into account, the single point energies obtained using these different methods were added to the thermally corrected Gibbs free energies at 298.15 K.

3. RESULTS AND DISCUSSIONS We start the analysis with the lowest energy conformations of cysteine and methionine (Figure 1) optimized at the B3LYP/6-

2. COMPUTATIONAL DETAILS All of the calculations have been performed using the Gaussian 0918 suite of programs. The computational protocols involved (a) geometry optimizations at the B3LYP/6-31G(2df,p) level or MP2/6-31G(2df,p) level of theory. The associated harmonic frequencies were scaled by a factor of 0.9854 in the case of the B3LYP optimizations and 0.9445 in the case of the MP2 optimizations. (b) Geometry optimization was followed by single-point calculations with seven different density functional methods BPW91,19−21 BMK,22 B3LYP,23,24 M05-2X,25 M062X,26 TPSSh,27 and B2PLYP28 corresponding to different rungs of exchange and correlation energy functionals, as well as three wave function-based methods HF, MP2, and CCSD(T).29 Even though it is well-known that the HF method is severely deficient33 and is thus expected to yield poor results, we still include the HF computations for the sake of completeness. Two different Pople-style basis-sets (6-311++G(3df,2p) and 6-31+G(d,p)) and Dunning-style basis sets (aug-cc-pVTZ and aug-cc-pVDZ) were used in the single point computations to

Figure 1. Lowest energy conformers of cysteine (bottom) and methionine (top).

31G(2df,p) level (same level as in G4 theory) with subsequent single-point calculations at various levels of theory (see Computational Details). Our structures are in agreement with the lowest energy conformer of cysteine previously obtained by Allen, et al.,30 and the lowest energy conformer of methionine obtained by Notario et al.8 The CBH-1 schemes obtained by executing isodesmic bond separation reactions on cysteine and methionine, and the isoatomic CBH-2 schemes obtained by preserving the immediate chemical environment of all of the atoms in cysteine and methionine are given below (vide supra, see references 12 and 13 for more details on constructing the schemes). Cysteine: B



Article

Figure 2. Pictorial representation of how the reaction energies (Y axis, kcal/mol) vary as we go from the isodesmic bond separation baseed CBH-1 rung to the isoatomic CBH-2 rung using different wave function-based as well as DFT methods for both cysteine and methionine. The vertical order in which the legend is given is the same as the horizontal order in which the cones are given from left to right. The 6-311++G(3df,2p) basis set was used throughout (extrapolated in the case of CCSD(T), see Computational Details). Figure 1 contains the data corresponding to the plots. The geometries were optimized at the B3LYP/6-31G(2df,p) level of theory.

Table 1. Reaction Energy (kcal/mol) Data Corresponding to Figure 2a

a

molecule

CCSD(T)

MP2

BPW91

BMK

B3LYP

M05-2X

M06-2X

TPSSh

B2PLYP

HF

cysteine (CBH-1) cysteine (CBH-2) methionine (CBH-1) methionine (CBH-2)

58.1 1.9 63.0 1.8

62.8 2.1 68.8 1.5

55.7 1.4 57.2 −1.1

56.8 0.7 60.0 −0.4

53.7 0.1 56.3 −1.4

59.1 1.3 63.4 0.1

58.2 1.4 62.5 0.2

54.7 1.8 56.3 −0.4

57.0 0.7 60.9 −0.3

47.9 −2.9 52.2 −2.1

For details of the level of theory, see Figure 2.

Table 2. Reaction Energies (kcal/mol) for the Lowest Energy Conformers of Cysteine and Methionine with Different Basis Functions for the MP2 Methoda

a

molecule

6-311++G(3df,2p)

6-31+G(d,p)

aug-cc-pVTZ

aug-cc-pVDZ

cysteine (CBH-1) cysteine (CBH-2) methionine (CBH-1) methionine (CBH-2)

62.8 2.1 68.8 1.7

60.6 2.1 65.5 1.3

61.5 2.1 66.7 1.1

60.3 2.2 66.2 1.8

The geometries were obtained at the B3LYP/6-31G(2df,p) level of theory.

CBH-1:

C5H11NO2 S + 3C2H6 + H3CSH

C3H 7NO2 S + 5CH4

→ C2H5SH + C3H 9N(isopropyl amine) + H3CCOOH

→ 2C2H6 + H3CNH 2 + H3CSH + H 2CO + H3COH

+ C3H8 + H3CSCH3(dimethyl sulfide)

CBH-2:

3.1. Reaction Energies. At the CBH-1 level, with both wave function-based and density functional methods, the reaction energies for cysteine as well as methionine are quite large (of the order of 50−70 kcal/mol; Figure 2 and Table 1). In our previous works,12,13 the mean absolute reaction energies at CBH-1 for a test set comprising 20 organic molecules at various levels of theory was between 20 and 35 kcal/mol. However, no molecule in the test set contained as many as four different heavy atoms. Therefore the larger calculated reaction energies for cysteine and methionine at CBH-1 are not unexpected.

C3H 7NO2 S + 2C2H6 → C2H5SH + C3H 9N(isopropyl amine) + H3CCOOH

Methionine: CBH-1: C5H11NO2 S + 6CH4 + H 2S → 3C2H6 + H3CNH 2 + 2H3CSH + H 2CO + H3COH

CBH-2: C



Article

Table 3. Enthalpies of Formation (kcal/mol) of the Lowest Energy Conformer of Cysteine and Methionine at the CBH-2 Runga molecule

CCSD(T)

MP2

BPW91

BMK

B3LYP

M05-2X

M06-2X

TPSSh

B2PLYP

HF

cysteine methionine

−96.1 −104.3

−96.4 −104.2

−95.7 −101.6

−95.00 −102.2

−94.3 −101.2

−95.5 −102.7

−95.6 −102.9

−96.1 −102.2

−95.0 −102.3

−91.3 −100.6

a

The geometries were obtained at the B3LYP/6-31G(2df,p) level of theory. 6-311++G(3df,2p) basis set was used throughout. CCSD(T)/6-311+ +G(3df,2p) energies obtained using an additivity approximation.

Table 4. Enthalpies of Formation of (kcal/mol) for the Lowest Energy Conformers of Cysteine and Methionine with Different Basis Sets Using the MP2 Methoda

a

molecule

6-311++G(3df,2p)

6-31+G(d,p)

aug-cc-pVTZ

aug-cc-pVDZ

cysteine methionine

−96.4 −104.2

−96.4 −104.0

−96.3 −103.7

−96.4 −104.5

The geometries were obtained at the B3LYP/6-31G(2df,p) level of theory.

Table 5. Enthalpies of Formation (kcal/mol) of the Lowest Energy Conformer of Cysteine and Methionine at the CBH-1 Runga molecule

CCSD(T)

MP2

BPW91

BMK

B3LYP

M05-2X

M06-2X

TPSSh

B2PLYP

HF

cysteine methionine

−93.7 −101.4

−98.4 −107.1

−91.3 −95.5

−92.4 −98.3

−89.3 −94.6

−94.7 −101.7

−93.8 −100.8

−90.3 −94.6

−92.6 −99.2

−83.5 −90.5

a

The geometries were obtained at the B3LYP/6-31G(2df,p) level of theory. 6-311++G(3df,2p) basis set was used throughout. CCSD(T)/6-311+ +G(3df,2p) energies obtained using an additivity approximation.

simple Hartree−Fock (HF) method itself, all the density functionals as well as the MP2 method provide enthalpies of formations within 2 kcal/mol (for cysteine) and 3 kcal/mol (for methionine) in comparison to the accurate CCSD(T) method (Table 3). For cysteine, the enthalpy of formation varies between −94.3 and −96.4 kcal/mol. Similarly for methionine, the enthalpy of formation falls between −101.2 and −104.3 kcal/mol. At room temperature for cysteine, Roux et al. recently reported7 an experimental enthalpy of formation of −91.4 ± 0.44 kcal/mol. Weighted over ten different conformations, they obtained −94.3 kcal/mol using the G3(MP2) method, −94.0 kcal/mol using the G3 method, and −92.6 kcal/mol using an empirical correction32 on top of the G3(MP2) method. For the lowest energy conformation, Dixon and co-workers report9 values of −94.4 and −93.6 kcal/mol for the G3(MP2) method via atomization and by a designed thermochemical scheme, respectively. Finally, the popular NIST Web site does not provide any gas-phase thermochemical information on cysteine. For methionine, the NIST WebBook provides a value of −98.8 ± 0.98 kcal/mol. Notario, et al. take into account ten different conformations and based on the G3, G4 and experimental methods recommend8 a value of −102.8 ± 2.4 kcal/mol. For the lowest energy conformer, Dixon and coworkers report −102.5 and −101.6 kcal/mol using the G3(MP2) method via atomization and by a designed thermochemical scheme, respectively. The values obtained for the lowest energy conformations of cysteine and methionine using commonly used density functionals (Table 3) are thus very similar (within 1−2 kcal/ mol) to those obtained by Notario et al., using G3, G3(MP2), and G4 methods7,8 as well as the results obtained by Dixon and co-workers.9 Finally, as in the case of the reaction energies, the enthalpies of formations also do not vary much with changing basis sets (Table 4 and SI). To demonstrate that an additive approximation (see Computational Details) to the CCSD(T)/6-311++G(3df,2p) level of theory based on CCSD(T)/6-31+G(d,p) produces the same quality of results, we computed the enthalpies of

It was also expected that the CBH-2 reaction energies for cysteine and methionine will be smaller relative to the CBH-1 energies. But the magnitude of the fall in the reaction energies on going from CBH-1 to CBH-2 (Figure 2 and Table 1) is dramatic. The reaction energies for these two amino acids are very small (of the order of 1−3 kcal/mol) at CBH-2 and are, in fact, in the same range as those for the test set containing much simpler organic molecules in refs 12 and 13. This clearly indicates that the isoatomic CBH-2 rung provides significantly sophisticated bond matching for systems containing several different heavy atoms as well. The same trend, i.e., 50−70 kcal/ mol reaction energies at CBH-1 and a much smaller CBH-2 reaction energy, is noticed irrespective of the basis set (Poplestyle or Dunning) or the method (DFT or wave function-based method) used (Table 2 and SI), thus demonstrating the robustness of CBH. 3.2. Enthalpies of Formation of the Lowest Energy Conformations of Cysteine and Methionine. Having computed the reaction energies, we can now derive the enthalpies of formation of cysteine and methionine at room temperature (298.15 K) using the experimental enthalpies of formation for the components involved. The experimental enthalpies of formations used at the isoatomic rung (kcal/mol) are as follows: ethane, −20.08 ± 0.1; methanethiol, −5.46 ± 0.14; acetic acid, −103.40 ± 0.4; isopropylamine, −20.02 ± 0.19; propane, −25.00 ± 0.1; ethanethiol, −11.06 ± 0.14; and dimethylsulfide, −8.96 ± 0.14. These values are obtained from the NIST WebBook, ref 14b or Pedley et al.’s thermochemical tables.31 Given that CBH-2 provides error cancellation much more effectively than CBH-1 (refs 12 and 13), we focus on the enthalpies of formation obtained at CBH-2. It is important to note that we have not considered rungs higher than CBH-2 in this work as the experimental enthalpies of formation of several of the individual components involved in the higher rungs are not known. Table 3 provides the computed enthalpies of formation for the lowest energy conformations of cysteine and methionine at CBH-2 using various methods. For both molecules, barring the D



Article

formation of cysteine and methionine at CCSD(T)/6-311+ +G(3df,2p) as well. The enthalpies of formation of the lowest energy conformations of cysteine and methionine are −96.1 and −104.1 kcal/mol, respectively (compared with very similar numbers obtained via an additivity approximation, see Table 3). For the sake of completeness, the enthalpies of formation at the CBH-1 rung are presented in Table 5. On comparing the enthalpies of formations obtained using CBH-1 (isodesmic bond separation, Table 5) and those obtained using CBH-2 (isoatomic, Table 3), it is readily seen that MP2 and the DFT methods perform uniformly better in CBH-2 (within 1−3 kcal/ mol relative to the accurate CCSD(T) method) whereas larger deviations as high as 5−7 kcal/mol are seen with CBH-1. The relatively poor error canceling ability of the CBH-1 rung (isodesmic bond separation scheme) vs the advanced error canceling capabilities of CBH-2 (isoatomic scheme) are clearly shown by these results. 3.3. Role Played by the Level of Theory Involved in Optimization. Thus far, the geometries of cysteine, methionine and other species involved in the CBH-2 scheme were optimized using the B3LYP density functional. However, in the case of amino acids and other biomolecules, weak interactions such as hydrogen bonds or dispersion effects may play a significant role.33,34 To consider these effects, we optimized the geometries of the lowest energy conformers of cysteine, methionine and other constituents in the CBH-2 scheme at the MP2/6-31G(2df,p) level of theory which is known to provide a good treatment of such effects, and then carried out single point computations at different levels of theory. The enthalpies of formation obtained from the MP2 geometry are summed up in Table 6 for two different post-HF

A cautionary note however is that, for larger peptides and polypeptides containing secondary or tertiary interactions, dispersion effects may be more pronounced such that the level of theory employed for getting the geometries may be important. This aspect will be a subject of future investigation. For the amino acids however, B3LYP geometries suffice and we continue to use them in this study. 3.4. Enthalpies of Formations taking Multiple Conformations into Account. Given the conformational complexity of the amino acids, a legitimate comparison of the enthalpies of formation obtained using the isoatomic scheme with the experimental values is possible only when the issue of conformational flexibility of cysteine and methionine is aptly taken into account. To address the conformational issue in the present work, in addition to the lowest energy conformations of both the amino acids, we choose the next five most stable conformations of cysteine and the next four most stable conformers of methionine from the studies by Notario et al. (Figure 3).7,8

Table 6. Enthalpies of Formation (kcal/mol) of the Lowest Energy Conformer of Cysteine and Methionine at the CBH2 Rung Using MP2/6-31G(2df,p) Geometries for Two Different Density Functional Methods and Two Post-HF Methodsa molecule

CCSD(T)

MP2

M05-2X

B2PLYP

cysteine methionine

−96.1 −104.2

−96.4 −104.4

−95.5 −102.0

−95.7 −102.0

a

6-311++G(3df,2p) basis set was used throughout. CCSD(T)/6-311+ +G(3df,2p) energies obtained using an additivity approximation.

Figure 3. Other conformers of cysteine (bottom) and methionine (top) considered in this study.

methods and two density-functionals. A comparison of the results in Table 6 with the appropriate columns in Table 3 reveals that very similar enthalpies of formation are seen for both geometries. This similarity again illustrates that sufficient error-cancellation is achieved with the isoatomic CBH-2 scheme. Additionally, we have also added Grimme’s empirical D3 dispersion corrections to the B3LYP geometries. The corrections resulted in the following enthalpies of formation for the lowest energy conformation of cysteine: −95.6 kcal/mol (B3LYP), −95.7 kcal/mol (B2PLYP), and −95.7 kcal/mol (M06-2X), and the following enthalpies of formation for the lowest energy conformation of methionine: −103.9 kcal/mol (B3LYP), −103.9 kcal/mol (B2PLYP), and −103.1 kcal/mol (M06-2X). The general trend after the addition of dispersion correction is that, the enthalpies of formation still fall in same range as shown in Table 3 and the enthalpies of formations obtained with the density functionals fall within a kcal/mol of the corresponding values obtained with CCSD(T) (Table 3).35

We restrict ourselves to the top six and five conformations of cysteine and methionine, respectively, based on the fact that the rest of the conformers have a population of less than 5% and contribute very little to the overall enthalpy of formation. Throughout, it is always noticed that the most stable conformation in their work and ours remains the same though certain minor variations are noticed in the populations of the other conformers (SI). Coordinates of the optimized geometries of all of the conformers used are provided in the SI. Similar to Roux et al.’s approach, we first calculate the individual enthalpies of formations for each conformation (ΔHi) and then obtain the weighted enthalpy of formation (ΔH) as i=1

ΔH =

∑ XiΔHi n

where n is the number of conformations considered (6 for cysteine and 5 for methionine). Xi refers to the weight of the ith E



Article

conformer and is obtained by Boltzmann-weighting the relative free energies (referenced to the lowest energy conformer) as

Table 9. Enthalpies of Formation (kcal/mol) Taking into Account Multiple Conformers of Cysteine and Methionine at the CBH-2 Rung Using MP2/6-31G(2df,p) Geometries for Two Different Density Functional Methods and Two Post-HF Methodsa

i=1

Xi = {e−ΔGi / RT /∑ e−ΔGi / RT } n

Full details of the weights Xi and the various ΔHi values for both cysteine and methionine obtained using CBH-2 are provided in Table 7 and SI. The overall weighted ΔH values are

molecule

CCSD(T)

MP2

M05-2X

B2PLYP

cysteine methionine

−96.1 −104.2

−96.4 −104.4

−95.5 −102.7

−95.7 −102.0

a

6-311++G(3df,2p) basis set was used throughout. CCSD(T)/6-311+ +G(3df,2p) energies obtained using additivity approximation.

Table 7. Weights and Individual Enthalpies of Formation (kcal/mol) for All of the Conformers of Cysteine and Methionine Considered in This Study with the CCSD(T) Method at CBH-2a molecule

conformer

Xi

ΔHi

cysteine cysteine cysteine cysteine cysteine cysteine methionine methionine methionine methionine methionine

1 2 3 4 5 6 1 2 3 4 5

0.40 0.04 0.22 0.10 0.16 0.08 0.63 0.02 0.08 0.16 0.11

−96.1 −94.6 −95.0 −94.8 −94.8 −94.5 −104.3 −102.1 −103.3 −103.4 −103.4

the additivity approximation is a judicious approach to cut the computational expense without the loss of desired accuracy, even when dealing with multiple conformations as well. 3.5. Basis Set Dependence of Conformational Ordering. An interesting observation when multiple conformations are taken into account is that the size of the basis set used plays an important role in determining the conformational ordering. This is primarily due to the fact that energy differences between different conformations can be very small, i.e., many conformations can occur within an energy window of 1−2 kcal/mol. Thus, while the enthalpy of formation of a single conformer can be computed at the CBH-2 rung within an uncertainty of ±1−2 kcal/mol with modest basis sets, larger basis sets such as 6-311++G(3df,2p) may be needed to get the correct conformational ordering. For instance, even with the highly accurate CCSD(T) method, a small double-ζ basis set such as 6-31+G(d,p) reverses the conformational ordering sufficiently that the well-established lowest energy conformation of cysteine ceases to be the most stable conformer (Table 10). However, the overall enthalpies of formation are only affected slightly (by about 1 kcal/mol or less) by including multiple low energy conformations.

a

The conformers are arranged in the same order as their co-ordinates are given in the SI. The geometries were obtained at the B3LYP/631G(2df,p) level of theory. The energies were extrapolated to the 6311++G(3df,2p) basis-sets.

given in Table 8. For methionine, our weighted ΔH values (Table 8) with density functionals as well as post-HF methods agree with the recommended value of −102.8 ± 2.4 kcal/mol. For cysteine, the values computed here (barring the HF values, Table 8) are in the same range (±1−2 kcal/mol) as the G3 and G3MP2 results of Roux et al.7 but differ by 2−3 kcal/mol from their empirically corrected G3MP2 results. Finally, as in the case of dealing with only one conformer (vide supra), it is noticed with multiple conformers as well that similar results are obtained with B3LYP or MP2 geometries (Tables 8 and 9 also see SI for the individual weights and enthalpies of formation for all the conformers considered in this study at the MP2 geometry). For the CCSD(T) results in Table 8, the enthalpies of formation of cysteine and methionine taking multiple conformations were obtained by an additivity approximation. Without any approximation, the direct CCSD(T)/6-311+ +G(3df,2p) enthalpies of formations of cysteine and methionine taking multiple conformations into account are −95.3 and −103.7 kcal/mol respectively (compared with very similar numbers obtained via additivity, see Table 8). Therefore

Table 10. Illustrating How with a Small Basis Set Such as 631+G(d,p), Even the Highly Accurate CCSD(T) Method Gives a Wrong Conformational Orderinga molecule

conformer

Xi

ΔHi

cysteine cysteine cysteine cysteine cysteine cysteine

1 2 3 4 5 6

0.19 0.02 0.24 0.12 0.34 0.09

−96.1 −94.6 −95.5 −95.3 −95.0 −95.0

a

The table contains weights and individual enthalpies of formation (kcal/mol) for all the conformers of cysteine considered in this study with the CCSD(T) method at CBH-2. [The geometries were obtained at the B3LYP/6-31G(2df,p) level of theory. 6-31+G(d,p) basis set was used.] The conformers are arranged in the same order as their coordinates are given in the SI.

Table 8. Enthalpies of Formation (kcal/mol) of Cysteine and Methionine Taking into Account Multiple Conformations, at the CBH-2 Runga molecule

CCSD(T)

MP2

BPW91

BMK

B3LYP

M05-2X

M06-2X

TPSSh

B2PLYP

HF

cysteine methionine

−95.3 −104.0

−95.6 −103.8

−95.1 −101.2

−94.2 −101.8

−93.6 −100.9

−94.7 −102.30

−94.8 −102.4

−95.5 −101.6

−94.2 −102.0

−92.1 −100.5

a

The geometries were obtained at the B3LYP/6-31G(2df,p) level of theory. 6-311++G(3df,2p) basis set was used throughout. CCSD(T)/6-311+ +G(3df,2p) energies obtained using an additivity approximation. The individual weights and enthalpies of formations for individual conformers are provided in the Supporting Information in a detailed manner. F



Article

3.6. Isoatomic Scheme Is Better than Isodesmic Bond Separation Scheme. A recurring theme in all of the discussions so far has been that the isoatomic scheme (CBH2) provides sufficient error-cancellation with regularly used density functionals or the MP2 method not only for simple organic molecules as noted in our previous works12,13 but also for small to medium sized biomonomers such as cysteine and methionine containing as many as four different heavy atoms (C, O, N, and S) in close proximity. The isoatomic scheme can hence be used with common density functionals36 or the MP2 method in organic and biomolecular applications to yield better results than with the simple isodesmic bond separation scheme.37,38

(6) The isodesmic bond separation scheme is currently widely used in the thermochemistry community to benefit from error cancellations. In conjunction with the CBH work on simple organic molecules, we recommend the use of the isoatomic scheme (CBH-2) for the thermochemistry of small/ medium sized biomonomers as well. This offers much superior error cancellation over the isodesmic bond separation scheme so that less expensive density functionals or the MP2 method with modest sized basis sets can be used for accurate thermochemistry. Overall, this work establishes a proper computational protocol to apply CBH to interesting biomonomers and thus opens the door to advancing the use of error-canceling thermochemical hierarchies for applications to larger biomolecules.

■

4. CONCLUSION In conclusion, the most salient results from this study can be summarized as follows: (1) Our computed enthalpy of formation of methionine in the gas phase at room temperature using the isoatomic scheme, taking into account multiple conformations, is in agreement with the value of −102.8 ± 2.4 kcal/mol recommended by Roux et al. as well as the values provided by Dixon et al. The enthalpy of formation provided in the NIST WebBook, however, differs by 4−5 kcal/mol from each of our independent studies. We thus suggest the use of the value recommended by Roux et al. over the NIST value. (2) For cysteine, the weighted enthalpy of formation we obtain using the isoatomic scheme with commonly used density-functionals and post-HF methods is in the range of −93 to −96 kcal/mol, within 2 kcal/mol of Dixon’s values as well as the G3 and G3(MP2) values of Roux et al. However, the computed values deviate significantly from the available experimental value of −91.4 ± 0.4 kcal/mol.7 It will be interesting to consider this further in future studies. (3) The method employed in the geometry optimization of cysteine and methionine (MP2 vs B3LYP) does not change the computed enthalpies of formation. (4) A relatively small basis set such as 6-31+G(d,p) might be sufficient to accurately compute the enthalpy of formations of small/medium biomonomers such as the amino acids using the CBH-2 scheme. However, to properly account for the correct conformational ordering, relatively larger basis sets such as 6311++G(3df,2p) (which are still modest relative to the routinely used basis sets in theoretical thermochemistry) are necessary. (5) Putting together the enthalpies of formations obtained in this study, that by Dixon and co-workers using only the lowest energy conformation and Roux et al.’s G3 and G3(MP2) results, a striking observation is that, as long as the different conformations are closely spaced the enthalpy of formation of the lowest energy conformer is usually within 1−2 kcal/mol of the enthalpy of formation taking multiple conformers into account. The higher energy conformers can only make the calculated heat of formation more positive compared to the lowest energy conformation. It is likely that the multiple conformations may play a bigger role in the free energy than in the enthalpy due to the sampling contribution to TΔS. Therefore, it will be very interesting to perform a careful study in the future on a larger test set of molecules to see whether the lowest energy conformation is sufficient to achieve an accuracy of 1−2 kcal/mol in the computed enthalpies of formations.

ASSOCIATED CONTENT

* Supporting Information S

Several tables of additional data, the CBH-1 and CBH-2 reaction schemes for cysteine and methionine, the structures and the optimized geometries of all of the conformers used. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The authors acknowledge support from NSF Grant CHE0911454 at Indiana University. REFERENCES

(1) Nelson, D, L.; Cox, M. M. Lehninger Principles of Biochemistry; Palgrave-Macmillan: New York, 2008. (2) Pena, I.; Eugenia Sanz, M.; Lopez, J. C.; Alonso, J. L. Preferred Conformers of Protoionogenic Glutamic Acid. J. Am. Chem. Soc. 2012, 134, 2305−2312. (3) Cocinero, E. J.; Lasarri, A.; Grabow, J.-U.; Lopez, J. C.; Alonso, J. L. The Shape of Leucine in the Gas Phase. ChemPhysChem 2007, 8, 599−604. (4) Blanco, S.; Sanz, E.; Lopez, J. C.; Alonso, J. L. Revealing the Multiple Structures of Serine. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 20183−20188. (5) Sanz, M. E.; Blanco, S.; Lopez, J. C.; Alonso, J. L. Rotational Probes of Six Conformers of Neutral Cysteine. Angew. Chem., Int. Ed. 2008, 47, 6216−6220. (6) Perez, C.; Mata, S.; Blanco, S.; Lopez, J. C.; Alonso, J. L. JetCooled Rotational Spectrum of Laser-Ablated Phenylalanine. J. Phys. Chem. A 2011, 115, 9653−9657. (7) Roux, M. V.; Foces-Foces, C.; Notario, R.; da Silva, M. A. V. R.; da Silva, M. D. M. C.; Santos, A. F. L. O. M .; Juaristi, E. Experimental and Computational Thermochemical Study of Sulfur-Containing Amino Acids: L-cysteine, L-Cystine, and L-Cysteine-Derived Radicals. S−S, S−H, and C−S Bond Dissociation Enthalpies. J. Phys. Chem. B 2010, 114, 10530−10540. (8) Roux, M. V.; Notario, R.; Segura, M.; Chickos, J. S.; Liebman, J. F. The Enthalpy of Formation of Methionine Revisited. J. Phys. Org. Chem. 2012, 25, 916−924. (9) Dixon and coworkers have computed the heats of formations of all the naturally occurring amino acids: Stover, M. L.; Jackson, V. E.; Matus, M. H.; Adams, M. A; Cassady, C. J.; Dixon, D. A. Fundamental Thermochemical Properties of Amino Acids: Gas-Phase and Aqueous G



Article

Acidities and Gas-Phase Heats of Formation. J. Phys. Chem. B 2012, 116, 2905−2916. (10) For a recent review on the Gn methods, see: Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gn theory. WIREs Comput. Mol. Sci. 2011, 1, 810−825. (11) Radom, L.; Lathan, W, A.; Hehre, W, J.; Pople, J. A. Molecular orbital theory of the electronic structure of organic compounds. VIII. Geometries, energies, and polarities of C3 hydrocarbons. J. Am. Chem. Soc. 1971, 93, 5339−5342. (12) Ramabhadran, R. O.; Raghavachari, K. Theoretical Thermochemistry for Organic Molecules: Development of the Generalized Connectivity-Based Hierarchy. J. Chem. Theory. Comput. 2011, 7, 2094−2103. (13) Ramabhadran, R. O.; Raghavachari, K. Connectivity-Based Hierarchy for Theoretical Thermochemistry: Assessment Using Wave Function-Based Methods. J. Phys. Chem. A 2012, 116, 7531−7537. (14) For an isodesmic bond separation (IBS) based protocol similar in computational expense to the G3 and G4 methods which has been successfully applied to organic and small inorganic molecules, see Bakowies’s ATOMIC protocol (a) Bakowies, D. Ab initio Thermochemistry using Optimal-Balance Models with Isodesmic Corrections: The ATOMIC protocol. J. Chem. Phys. 2009, 130, 144113. (b) Bakowies, D. Ab initio Thermochemistry with high- level Isodesmic Corrections: Validation of the ATOMIC Protocol for a Large Set of Compounds with First Row Atoms (H, C, N, O, F). J. Phys. Chem. A 2009, 113, 11517−11534. (c) Bakowies, D. Assessment of Density Functional Theory for Approaches based on Bond Separation Reactions. J. Phys. Chem. A 2013, 117, 228−243. (15) For Fishtik’s group additivity methods for organic molecules, see: (a) Fishtik, I. Group Additivity Methods Without Group Values. J. Phys. Chem. A 2006, 110, 13264−13269. (b) Fishtik, I. Unique Stoichiometric Representation for Computational Thermochemistry. J. Phys. Chem. A 2012, 116, 1854−1863. (16) The thermochemical reaction scheme obtained at the CBH-2 rung is called the isoatomic scheme because it preserves the immediate chemical-environments of all the atoms in an organic molecule. We therefore use the terms “CBH-2 rung” and the “isoatomic scheme” interchangeably in this work. See refs 12 and 15 for more details. (17) Wheeler et al. first pointed out the inconsistencies in the homodesmotic scheme for hydrocarbons. They derived specific schemes applicable for hydrocarbons that corrected for the inconsistencies. The corrected homodesmotic schemes for hydrocarbons were called “hypohomodesmotic”, “hyperhomodesmotic”, etc. When generalizing Wheeler et al.’s work to all organic molecules (ref 12), we suggested that a more illustrative name for the hypohomodesmotic scheme applicable to hydrocarbons as well as non-hydrocarbons is “isoatomic scheme”. For more details on Wheeler et al.’s important contribution to hydrocarbons, see: Wheeler, S. E.; Houk, K. N.; Schleyer, P. v. R.; Allen, W. D. A hierarchy of Homodesmotic Reactions for Thermochemistry. J. Am. Chem. Soc. 2009, 131, 2547−2560. (18) Frisch, M. J.; Trucks, G. W.; Schlegel. H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalemani, G.; Barone, V.; Mennucci, B.; Pettersson, G. A.; et al. Gaussian 09, revision h08; Gaussian, Inc.: Wallington, CT, 2009. (19) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B 1992, 46, 6671−6687. (20) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Erratum: Atoms, Molecules, Solids and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B 1993, 48, 4978. (21) Perdew, J. P.; Wang, Y. Generalized Gradient Approximation for the Exchange-Correlation Hole of a Many-Electron System. Phys. Rev. B 1996, 54, 16533−16539.

(22) Boese, A. D.; Martin, J. M. L. Development of Density Functionals for Thermochemical Kinetics. J. Chem. Phys. 2004, 121, 3405−3416. (23) Becke, A. D. Density-functional thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (24) Lee, C. T.; Yang, W. T.; Parr, R. G. Development of the ColleSalvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785−789. (25) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Design of Density Functionals by Combining the Method of Constraint Satisifaction with Parametrization for Thermochemistry, Thermochemical kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2006, 2, 364−382. (26) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited states, and Transition elements: Two New Functionals and Systematic Testing of Four M06-Class Funtionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (27) Tao, J, M.; Perdew, J. P.; Starvoverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical MetaGeneralized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401−146404. (28) Grimme, S. Semiempirical Hybrid Density Functional with Perturbative Second-Order Correlation. J. Chem. Phys. 2006, 124, 034108. (29) (a) For the CCSD(T) method: Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A Fifth-Order Perturbation Comparison of Electron Correlation Theories. Chem. Phys. Lett. 1989, 157, 479−483. (b) For the MP2 method: Moller, C.; Plesset, M. S. Note on Approximate Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618−622. (c) For the HF method: Fock, V. Näherungsmethode zur Lösung des quantenmechanischen Mehrkörperproblems. Z. Phys. 1930, 61, 126−148. (30) Wilke, J, J.; Lind, M. C.; Schaefer, H. F.; Csaszar, A. G.; Allen, W. D. Conformers of Gaseous Cysteine. J. Chem. Theory. Comput. 2009, 5, 1511−1523. (31) Pedley, J. B.; Naylor, R. D.; Kirby, S. P. Thermochemical Data of Organic compounds, 2nd ed.; Chapman and Hall: London, 1986. (32) Anantharaman, B.; Melius, C. F. Bond Additivity corrections for G3B3 and G3MP2B3 Quantum Chemistry Methods. J. Phys. Chem. A 2005, 109, 1734−1747. (33) Jensen, F. Introduction to Computational Chemistry; John Wiley & Sons: Chichester, UK, 2007. (34) Hobza, P.; Muller-Dethlefs, K. Non-covalent Interactions: Theory and Experiment; RCS Publishing: London, 2009. (35) For Grimme’s D3 correction, see: Grimme, S.; Antony, J.; Ehlrich, S.; Krieg, H. A Consistent and Accurate Ab initio Parameterization for Density Functional Dispersion Correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (36) For a recent, interesting application of how the isoatomic scheme can provide insights on the DFT errors based on 1, 3 interactions in organic molecules (alkanes), see: Johnson, E. R.; Contreras-Garcia, J.; Yang, W. Density-Funtional Errors in Alkanes: A Real-Space Perspective. J. Chem. Theory Comput. 2012, 8, 2676−2681. (37) The isodesmic bond separation scheme as well as the isoatomic bond separation schemes can however be inadequate in the cases of molecules possessing certain special features such as aromaticity. For more details see refs 12 and 13. (38) It is worthwhile to point out here that there is unfortunately some confusion in the literature concerning the term “isodesmic”. As originally defined by Pople and co-workers in the 1970s, an isodesmic reaction is one in which the number of each type of bond is conserved. A specific type of isodesmic reaction is the bond separation reaction but while the bond separation reaction is always isodesmic, an arbitrarily constructed isodesmic reaction is not necessarily a bond separation reaction. A major objective of the CBH is to provide solutions to such confusions arising due to definition-based inconsistancies.

H


Application of the Generalized Connectivity-Based Hierarchy to

Recommend Documents