Catalyzed Reaction of Isocyanates (RNCO) with Water

The reactions between substituted isocyanates (RNCO) and other small molecules (e.g. water, alcohols, and amines) are of signiﬁcant industrial importance, particularly for the development of novel polyurethanes and other useful polymers. We present very high-level ab initio computations on the HNCO + H 2 O reaction, with results targeting the CCSDT(Q)/CBS//CCSD(T)/cc-pVQZ level of theory. Our results aﬃrm that hydrolysis can occur across both the N −− C and C −− O bonds of HNCO via concerted mechanisms to form carbamate or imidic acid with ∆ H 0K barrier heights of 38.5 and 47.5 kcal mol − 1 . A total of 24 substituted RNCO + H 2 O reactions were studied. Geometries obtained with a composite method and reﬁned with CCSD(T)/CBS single point energies determine that substituted RNCO species have a signiﬁcant inﬂuence on these barrier heights, with an extreme case like ﬂuorine lowering both barriers by close to 20 kcal mol − 1 and most common alkyl substituents lowering both by approximately 4 kcal mol − 1 . Natural Bond Orbital (NBO) analysis provides evidence that the predicted barrier heights are strongly associated with the occupation of the in-plane C − O* orbital of the RNCO reactant. Key autocatalytic mechanisms are considered in the presence of excess water and RNCO species. Additional waters (one or two) are predicted to lower both barriers signiﬁcantly at the CCSD(T)/aug-cc-pV(T+d)Z level of theory with strongly electron withdrawing RNCO substituents also increasing these eﬀects, similar to the uncatalyzed case. The 298 K Gibbs energies are only marginally lowered by a second catalyst water molecule, indicating that the decreasing ∆ H 0K barriers are oﬀset by loss of translational entropy with more than one catalyst water. Two-step 2RNCO + H 2 O mechanisms are characterized for the formation of carbamate and imidic acid. The second step of these two pathways exhibits the largest barrier and presents no clear pattern with respect to substituent choice. Our results indicate that an additional RNCO molecule might catalyze imidic acid formation but have less inﬂuence on the eﬃciency of carbamate formation. We expect that these results lay a ﬁrm foundation for the experimental study of substituted isocyanates and their relationship to the energetic pathways of related systems

Additional theoretical research has helped further understanding of substituted isocyanate species. In 1994, McAllister and Tidwell performed comprehensive structural and isodesmic reaction analyses on a large selection of substituted isocyanates, with the goal of theoretically motivating the discovery of novel isocyanate derivatives. 59 Their work highlights many insightful features of a diverse range of RNCO species, but is not connected to any particular reaction mechanisms. In the past decade, a series of studies 45,[60][61][62][63][64][65][66][67] by Konovalov and coworkers predicted the quantum mechanical features of the R 1 NCO + R 2 OH (R 1 , R 2 =Ph, CH 3 , H) reactions in various conditions, including solvent effects and cooperative catalysis with the B3LYP/6-311++G(df,p) method. However, their small sample size precludes a deep understanding of how substituents generally influence the reaction energetics. Our preliminary computations also indicate water + HNCO reaction energies differs by more than 1 kcal mol −1 between the CCSD(T)/cc-pVQZ and B3LYP/cc-pVQZ levels of theory computed on the same geometries. This discrepancy between DFT and CCSD(T) in the simplest of cases exemplifies the need for rigorous ab initio methods to draw meaningful conclusions from theoretical study of more complex H 2 O + RNCO reactions.
In 2013 Wagner and coworkers predicted the barrier heights to the formation of carbamates given a small variety of substituted isocyanates (RNCO, R=CH 3 , CH 2 F, C 6 H 5 , SiH 3 , SiH 3 CH 2 ) with the B3LYP/6-31G* method. 44 They determined that greater electron withdrawing substituents on the RNCO molecule lowered the barrier to carbamate formation. However, Wagner and coworkers only consider the reaction across the N− −C double bond and exclude any mention of how substituents influence reaction across the C− −O double bond. They also neglect to explain why electron withdrawing effects influence the reaction barriers. The results of Wagner and coworkers also demonstrate significant discrepancies based on predictions made at the MP2/6-31G* and B3LYP/6-311++G(d,p) levels of theory employed in their study (e.g. a 5 kcal mol −1 disagreement for the CH 3 NCO + CH 3 OH transition state). More recently, Zhao and Suppes predicted the enthalpy of reaction for various isocyanates reacting with increasingly larger alcohols to form carbamates with the B3LYP/6-31G(d,p) method 68 . Their work elucidates trends concerning an array of aryl isocyanates and predicts that increasing the size of the reacting alcohol has little influence on the enthalpy of the urethane product relative to the reactants. In addition to these works, there are plenty of examples in the literature of specific R 1 NCO + R 2 OH reactions studied at various levels of theory. 43,50,51,69,70 All of this previous research emphasizes the importance of clearly understanding how substituted RNCO species influence the electronic structure and energetic landscape of the RNCO + H 2 O reactions. Despite its significant importance to industrial chemistry, the literature lacks a comprehensive and reliable theoretical benchmark for this system and a detailed analysis of the relationship between the isocyanate substituents and the electronic structure features of these reactions.
Our research builds on this body of research and improves previous characterization of the RNCO + H 2 O reactions in four important ways. First, we thoroughly consider the H 2 O reaction across the isocyanate C− −O bond which is often ignored and assumed to be high enough in energy to be unimportant. Even if that is generally the case, it would be helpful to test this assumption for many substituents and understand what conditions might result in exceptions. It is a possibility that imidic acid formation could be favored over carbamate depending on the energetic barriers of each process. Second, our work rigorously characterizes the energetic landscape of the parent HNCO + H 2 O reaction and uses these results to benchmark the levels of theory used to analyze the substituent trends. Previous studies used differing theoretical methods, making it difficult to compare results across the literature and determine consistent trends. Our reliable ab initio results will lay a firm and consistent foundation for any future work on isocyanate alcoholysis reactions. Third, we systematically study how substituents influence the electronic structure of important cooperative catalytic pathways (i.e. multiple water or RNCO molecules) that previous research [50][51][52] has suggested as very important for these systems. Finally, our sophisticated ab initio predictions are analyzed with Natural Bond Order(NBO) analyses to provide detailed understanding of the electronic structure manifest in each process in order to characterize relationships that could generalize to larger or more complex isocyanate reactions.
The RNCO + H 2 O reactions are characterized considering substituents grouped into two categories: the first group is a fundamental collection of carbon groups, pnictogens, chalcogens, and halogens that demonstrate periodic trends (R=CH 3 , SiH 3 , GeH 3 ; NH 2 , PH 2 , AsH 2 ; OH, SH, SeH; F, Cl, Br) while (R=CH 2 CH 3 , CH(CH 3 ) 2 , CH 2 CH 3 CH 3 , C(CH 3 ) 3 , C 6 H 5 , CHCH 2 , COH,CF 3 , COOH, SO 2 Cl, and CN) consists of larger substituents that are more relevant to industrial applications. 3,7,8 Three important autocatalyitc mechanisms (RNCO + 2 H 2 O, RNCO + 3 H 2 O, and H 2 O + 2 RNCO) are characterized considering a smaller subset of the aforementioned substituents. The influence of substituents on the efficiency of these catalyzed reactions is systematically studied for the first time. Our work should motivate and ground future study of these and related systems and, our findings can be extended in conjunction with other advanced analyses (e.g. kinetic models, solvent effects, etc.).

Methods
The geometries for the benchmark H 2 O + HNCO reaction are obtained using the CCSD(T) 71-74 method via the CFOUR 2.0 75 software package. The Dunning cc-pVQZ basis set is utilized for all stationary points except for the van der Waals complexes which are optimized with the aug-cc-pVTZ basis set to properly describe long range interactions. [76][77][78] Harmonic vibrational frequencies are obtained for each stationary point using the same level of theory as the geometry optimization and verify each geometry as a minimum or first-order saddle point (transition state). The connectivity of each stationary point is verified by performing qualitative intrinsic reaction coordinate (IRC) scans. Relative energies of each stationary point are further refined utilizing the focal point method of Allen and coworkers. 22,[79][80][81] The SCF/cc-pVXZ (X=Q,5,6) energy is extrapolated using the three point extrapolation formula of Feller 82 (E HF (X) = E ∞ HF + ae −bX ) and the correlation energy up to CCSD(T)/cc-pVXZ (X=Q,5) is extrapolated using the two point formula of Helgaker 83 (E corr (X) = E ∞ corr + aX −3 ) towards the complete basis set (CBS) limit .
A series of additive corrections are obtained to further refine the energy predictions and justify the approximations used in the geometry optimization. The zero-point vibrational energies (ZPVE) are obtained from the harmonic vibrational frequency computations producing 0 Kelvin enthalpies (H 0K ). Higher order CCSDT and CCSDT(Q) additive corrections [84][85][86][87] are computed with the cc-pVDZ basis set to capture as much electron correlation as possible and demonstrate convergence towards the FCI limit. A frozen core correction (∆ FC ) is obtained using the cc-pCVQZ basis set 88 to account for the difference in the CCSD(T) energy correlating all electrons and the CCSD(T) energy freezing the core. To account for scalar relativistic effects, the spin-free exact two-component one-electron (SF-X2c-1e) method along with the decontracted cc-pCVTZ basis set is used to compute a scalar relativistic correction (∆ rel ). [89][90][91][92][93][94] The correction accounts for the difference in energy with and without the (SF-X2c-1e) method turned on, while correlating all electrons. Finally, a diagonal Born-Oppenheimer correction (∆ DBOC ) is included as a diagnostic at the Hartree-Fock/cc-pVQZ level of theory to ensure that each stationary point is not influenced by any nearby conical intersection. 95,96 For the van der Waals complexes, corrections were computed with the appropriate augmented basis functions, which are detailed in Table S14. All corrections are obtained using the following software packages: CFOUR 2.0 75 , Molpro 2010 97 , and PSI4. 98 Geometries and harmonic frequencies for all species in the H 2 O + RNCO reactions are obtained using the MP2[TZ,QZ] + ∆CCSD(T)/DZ composite method recently highlighted 99 by Sherrill and coworkers and implemented in PSI4. The MP2[TZ,QZ] term refers to a density fitted second-order Møller-Plesset (MP2) gradient extrapolated to the CBS limit using the cc-pV(X+d)Z (X=T,Q) basis sets and the two-point extrapolation formula of Helgaker. The ∆CCSD(T)/DZ term corrects this extrapolated gradient with an additive density fitted CCSD(T)/cc-pV(D+d)Z gradient, which showed strong correlation with CCSD(T)/cc-pVQZ results, detailed in SI Section 1. The energy of each stationary point is further refined with CCSD(T)/cc-pV(X+d)Z (X=T,Q) single points extrapolated to the CBS limit. For the analysis considering the influence of multiple catalytic H 2 O or RNCO molecules, the CCSD(T)/augcc-pV(T+d)Z//MP2/jul-cc-pV[TZ,QZ]Z + ∆CCSD(T)/6-31+G** level of theory 100-106 is utilized to properly describe the noncovalent interactions involved in catalysis of the transition state with augmented basis functions while maintaining a feasible computational cost.
Natural Bond Orbital (NBO) analyses 107,108 are performed where appropriate to elucidate the electronic structure features associated with our predicted results. The NBO 6.0 program is interfaced with QCHEM 109 using a def2-QZVP basis set 110 and the B3LYP functional 111 to describe the exchange.

High-Level HNCO + H 2 O Reaction
The energetic corrections included in the focal point analysis (Tables S12-S14) behave uniformly across all stationary points and indicate no anomalous features of this system. The ZPVE corrections are the largest and are no greater than 5.5 kcal mol −1 with transition state ZPVE corrections generally closer to 1.5 kcal mol −1 . The frozen core and scalar relativistic corrections are both small and never larger than 0.24 kcal mol −1 . The DBOC corrections are negligible and indicate that none of the stationary points is in the vicinity of a conical intersection or surface crossing. In all cases, excellent convergence is exhibited with respect to basis set and higher order coupled cluster terms. These results affirm that our predictions are well within the bounds of chemical accuracy (i.e. one kcal mol −1 ) with respect to our electronic energies. It should be noted that our enthalpy results correspond to gas phase results at 0 K. Any Gibbs Free Energy results correspond to water vapor at 298 K and 1 atm of pressure which overestimates the entropic effects manifest in the liquid phase.
Water and HNCO can react via two different concerted mechanisms, with the water O−H bond breaking across the N− −C double bond to form carbamate or across the C− −O double bond to form imidic acid ( Figure 1). Both mechanisms proceed through transition states (TS1 and TS2, respectively) that form fourmembered rings between the two reactants. The predicted TS1 barrier (38.5 kcal mol −1 ) is much lower than the TS2 barrier (47.5 kcal mol −1 ). This is qualitatively in agreement with previous research 46,48,52 but we predict TS1 to be somewhat higher than the recent work of Nicolle and coworkers who predicted a relative TS1 ∆H 0K of 35.3 kcal mol −1 at the CCSD(T)//M06-2X/6-311++G(d,p) level of theory. 48 IRC computations confirm the concerted nature of these pathways and indicate that the reactants begin separated and form the depicted products. The possible van der Waal complexes formed between HNCO and water (VDW1 and VDW2) are included in Figure 1, but are not definitively part of the reaction pathway and have small enough binding energies (1.7 and 4.3 kcal mol −1 , respectively) that they are mostly excluded from our discussion and would certainly be decreasingly relevant at high temperatures. Figure S7 depicts the 298 K Gibbs free energy surface with the loss of translational en- tropy leading to an average increase of relative energies by about 9 kcal mol −1 for all stationary points and confirms our assumptions concerning the van der Waal complexes.
The TS1 stationary point exhibits an imaginary mode of 1726i as one of the water hydrogens begins to form a bond with the nitrogen atom, and the water oxygen begins to bond with the carbon. Natural Resonance Theory (NRT) predicts the bond order between the nitrogen and the water hydrogen to be 0.49 and the bond order of the nitrogen and carbon atoms to be 1.48 at TS1. Likewise, the carbon and water oxygen interaction exhibits a bond order of 0.56. The NBO second-order perturbation (E (2) ) analysis indicates that the most significant delocalization corresponds to electron density donated from the water lone pair into one of the isocyanate C−O* orbitals (E (2) =144.2 kcal mol −1 ), which is depicted in Figure 2. This transition state leads to the highly exothermic formation (-15.6 kcal mol −1 ) of carbamate (M1). M1 is nearly planar, with the NH 2 moiety slightly displaced from the plane of the molecule. The NBO E (2) analysis additionally predicts that significant delocalization of the carbonyl lone pair electron density into the N−C* and O−C* orbitals might be one of the contributing factors to the favorable product energy, affirming the findings of Bharatam and coworkers concerning the importance of delocalization in ureas and related species. M1 can proceed through TS4 with a small barrier of 9.5 kcal mol −1 corresponding to the rotation of the hydroxyl group towards the amino group. The resulting product, M3, is moderately higher in energy than M1 but is a necessary intermediate for TS5 and the ultimate dissociation into NH 3 and CO 2 . The reaction barrier for this process is quite high at 35.0 kcal mol −1 and is characterized by an imaginary mode of 1799i. Dissociation to CO 2 and NH 3 is an exothermic process (9.4 kcal mol −1 lower than M3) and results in the most energetically favorable stationary point on the surface. One noteworthy feature of TS5 is that the IRC path in the direction of M3 terminates at a new first-order saddle point. This is because TS5 possesses C s symmetry and the IRC is unable to break symmetry. Slightly projecting in either direction along this mode and following the resulting IRC leads to M3 and confirms the connectivity of M4 and NH 3 + CO 2 .
Water and HNCO can alternatively react through TS2 which exhibits a 1718i imaginary mode as the water oxygen begins to form a bond with the isocyanate carbon and one of the water hydrogens begins to bond with the isocyanate oxygen. The new O HNCO -H water bond is predicted to have a NRT bond order of 0.32 while the carbon and water oxygen nearly exhibit a formal single bond order (0.90). The dominant E (2) term in the TS2 structure (90.2 kcal mol −1 ) corresponds to delocalization of the isocyanate oxygen lone pair electron density into the most proximal water O-H* orbital as represented in Figure 2. The magnitude of this interaction is much smaller than the delocalization predicted in the TS1 structure and the preclusion of delocalization could be a partial explanation for why TS2 is much higher in energy. The product M2 of reaction through TS2 is the planar hydroxylated imidic acid molecule which lies only 2.8 kcal mol −1 higher in energy than water + HNCO. NBO computations indicate that M2 also exhibits significant delocalization as electron density from both hydroxyl groups is donated into the out of plane C−N* orbital. M2 also has two other conformers which correspond to rotations of the hydroxyl groups and are predicted to have CCSDT(Q)/CBS energies of 2.0 and 5.1 kcal mol −1 relative to M2.
M2 is connected to M1 via TS3, which corresponds to a hydrogen transfer from one of the hydroxyl groups to the nitrogen with an associated imaginary mode of 2012i. The barrier for this process is 30.7 kcal mol −1 . Unsurprisingly, M2 lies 18.4 kcal mol −1 higher in energy than M1, confirming that the reaction through TS1 is both kinetically and thermodynamically favored. The formation of M1 via the step-wise procession through TS2 and TS3 has been debated in the literature 51 and our results seem to indicate that formation of M1 directly though TS1 is the most likely option absent of any other effects such as catalysts.

RNCO Substituent Analysis
In most industrial applications of isocyanates, the parent HNCO species is often substituted with some much larger RNCO species. The most common -R groups are generally large aromatic rings or long polymers which can both be modified to possess a diverse array of electronic structure features. Therefore, it is necessary to understand how different substituents influence the electronic structure and the resulting energetics of the RNCO + H 2 O reactions. The ∆H 0K barrier heights and relative product enthalpies are predicted for a diverse set of 24 different substituted RNCO molecules proceeding through TS1 and TS2. In all cases, the isocyanate substituents are far enough away from the active site of each reaction mechanism to ensure that the reaction pathway remains qualitatively invariant to choice of -R group. All possible conformers of each RNCO species were searched for and considered in this analysis. They are denoted alphabetically with a subscript A corresponding to the lowest energy conformer if applicable. Conformers other than the lowest are only considered in the analysis if either barrier height is lowered by more Physical Chemistry Chemical Physics than 1 kcal mol −1 relative to conformer A. We admit that many of the substituents chosen might be unrealistic for industrial applications but justify the selections for two reasons. 1) The selected subset of substituents (listed in Figure 3) exhibits a wide array of electronic effects that allows for greater confidence that our predicted trends are robust and sufficiently span the effects normally exhibited by more standard substituents. 2) The chosen substituents are small enough to be well described by high-level methods and eliminate as much uncertainty as possible due to compromises in level of theory. The goal of our results is that the electronic structure trends are informative and reliable enough to guide the intuition of more practical isocyanate applications.  The parent case, R=H, has a predicted barrier of 37.8 kcal mol −1 which is within the range of chemical accuracy of our benchmark value targeting the CCSDT(Q)/CBS limit and further reinforces the reliability of this analysis. The barriers for TS1 are quite sensitive to substituent selection, spanning a range of 19.8 kcal mol −1 . The lowest barrier is predicted for R=F at 23.5 kcal mol −1 and the highest barrier corresponds to R=SiH 3 at 43.3 kcal mol −1 , significantly higher than the R=H parent case. The small main group substituents unsurprisingly have the largest effect on the barrier height, due to their extreme electron donating and withdrawing abilities. Generally speaking, most of the more standard hydrocarbon substituents lower the barrier to carbamate formation by about 3-4 kcal mol −1 . This tight range seems to indicate that the size of the substituent has only minor influence on the reaction barrier, which is not surprising as Zhao and Suppes came to similar conclusions considering the size of the reacting alcohol. 68 In TS1 and TS2, the substituents are opposite from the active site of the reaction and have no proximity to any moieties that might result in steric hindrance. The one exception to this might be the R=COH B case, as the substituent hydroxyl group is oriented towards the isocyanate oxygen, resulting in a favorable intramolecular interaction.

TS2 TS1
The predicted barrier heights for TS2 manifest trends that are qualitatively similar to the TS1 results, with the TS2 process always significantly less favorable. The TS2 barrier predictions are slightly more sensitive than the TS1 predictions and span a range from 21.2 kcal mol −1 with the lowest barrier predicted for R=F (30.1 kcal mol −1 ) and the highest R=SiH 3 (51.2 kcal mol −1 ). The most competitive case, R=F, still predicts TS1 to be the more favorable pathway by 6.5 kcal mol −1 . Our results confirm that it is highly unlikely that the TS2 barrier could be made competitive with the TS1 barrier by substituent selection alone in the uncatalyzed RNCO + H 2 O case, also noting that this conclusion is agnostic to other factors such as catalysts, solvent effects, etc. Figure  S1 presents a scatterplot between the substituted TS1 and TS2 barrier heights, which correlate excellently, indicating that both barriers exhibit similar relationships to substituent identity. The supplementary information (Tables S1-S6) also contains the 298 K Gibbs free energy correction for each reaction which is generally insensitive to substituent choice and is close to 11 kcal mol −1 for all species.  One of the primary goals of this analysis is to elucidate any features of the RNCO molecule that could predict a priori the barrier heights and provide a deeper understanding of the key electronic structure features of each mechanism. NBO predictions were utilized to generate results that might be related to the transition state barriers such as atom charges, orbital occupations, and NRT bond orders. We determined that one of the best predictors of barrier heights for the RNCO molecules is the occupation of the in-plane C−O* orbital. Figure 4 depicts a scatter-plot between the TS1 barrier and the C−O* orbital occupation, as well as a line of best fit that has a high R 2 value of 0.88. (The R=GeH 3 case is excluded as an outlier because it exhibits a C−O triple bond at the transition state and is not comparable to the other species.) The C−O* orbital occupation is an excellent proxy for the electron withdrawing ability of the substituents. Our results are in qualitative accord with the work of Wagner and coworkers 44 who determine, with little explanation, that electron withdrawing groups lower the barrier to reaction of alcohols across the N− −C double bond. Utilizing the C−O* orbital occupation as a predictive quantity is also appropriate, due to its direct relationship with the previously discussed dominant E (2) interaction characteristic of the TS1 structure. The association is so strong because less electron density present in the C−O* orbital allows for more facile delocalization of the water oxygen lone pair. Thus, the extent to which the C−O* orbital is filled in the isolated RNCO molecule will directly relate to energetic barrier of carbamate formation.

Occupation of RNCO C-O* Orbital
Our computations predict that a similar relationship also holds for the TS2 barriers. The line of best fit for these barriers is red in Figure 4 and also possesses an excellent R 2 of 0.86. In the unsubstituted TS2 structure, the dominant E (2) interaction involves the delocalization of the isocyanate oxygen lone pair into the closest water O−H* orbital. As electron density is removed from the C−O* orbital via substituent effects, the C−O bond shortens and the isocyanate oxygen lone pair can delocalize more easily. To summarize, electron withdrawing substituents on the isocyanate species lower the barrier to reaction via the pathways through both TS1 and TS2. This principle can almost certainly be generalized to RNCO reactions with different R-OH species and even related reactants such as R 2 N−H, because the reaction motif will remain rather consistent. However, more complicated molecules reacting with RNCO might convolute these well-behaved trends, introducing new factors such as sterics, dispersion interactions, and exotic electronic structure features. Nevertheless, we expect our trends to hold given similar reaction motifs and the fact that the isocyanate substituents are not involved in the active site of the reaction. It is also key to note that the same correlations hold and maintain their fidelity when using the 298 K Gibbs energies, which usually shift all estimated values by approximately 11 kcal mol −1 regardless of the substituent. Thus, our substituent analysis seems robust to higher temperature conditions.
The relative enthalpies of the products of each pathway are presented in Figure 5. The carbamate relative energies (shown in blue) are always exothermic and somewhat sensitive to substituent substitution. They range from -28.9 (R=F) to -9.7 kcal mol −1 (R=SiH 3 ). Highly electronegative substituents such as R=F, OH, and NH 2 make the carbamate more favorable than the reactants relative to the R=H case. This is likely associated with the aforementioned proclivity of the carbonyl to donate electron density into the neighboring C−N* and C−O* orbitals. The COH B carbamate product also stands out as be- ing particularly favorable (see Figure S2 for scatterplot between substituted M1 and M2 relative enthalpies), but this is partially due to an intramolecular interaction originating from the substituent hydroxyl group. The imidic acid products are much more sensitive to substituent selection, which is unsurprising as NBO analysis of the parent imidic acid predicts significantly more delocalization into the C−N* orbital compared to the carbamate. The imidic acid enthalpies relative to the reactants range from -18.3 (R=F) to 9.2 (R=SiH 3 ) kcal mol −1 and many of the alkyl substituted imidic acid products are predicted to be slightly less thermodynamically favored than the reactants. There is no case where the imidic acid arrangement is lower in energy than the carbamate analogue. The closest the two isomers get in energy is 10.1 kcal mol −1 (R=F), indicating that they are unlikely to be thermodynamically competitive in any case. Assuming no catalyst, our data supports that the H 2 O + RNCO reaction would favor proceeding through the TS1 pathway independent of substituent.

Multimolecular Mechanisms
Recent research has highlighted the importance of cooperative mechanisms leading to the same carbamate and imidic acid products in the presence of excess water or isocyanate molecules. [50][51][52] Wei and coworkers predicted at the MP2/6-311++G** level of theory that additional water molecules lower the barrier to reaction across both the N=C or C=O bonds. 52 We present an analysis of how substituents influence the barrier to carbamate and imidic acid formation considering one catalyst water molecule, two catalyst water molecules, and one catalyst RNCO molecule. Our CCSD(T)/aug-cc-pV(T+d)Z//MP2/jul-cc-pV([T,Q]+d)Z + ∆CCSD(T)/6-31+G** results are a significant improvement over previous computations and allows for confident determination of the trends that manifest in the autocatalyzed energetics due to substituent effects. Physical Chemistry Chemical Physics

Catalyst Waters
The H 2 O + HNCO products, carbamate and imidic acid, can be formed via a concerted reaction mechanism in the presence of one or two catalyst waters. In these cases, the catalyst waters form six or eight atom ring transition states, respectively including either the N− −C or C− −O bonds. The IRCs from these processes connect to pre and post reactive van der Waals complexes which are reasonably strong relative to the separated molecules (between 5-10 kcal mol −1 for the single water catalyzed case at the ωB97x-D3/def2-QZVPPD level of theory). These complexes are excluded from our analysis as they are decreasingly relevant at higher temperatures, and preliminary computations indicated that substituents have far less influence on the energy of the van der Waals complexes than on the transition state energies. In certain cases, the determined enthalpies of the catalyzed transition states are negative relative to the separated products, but are still referred to as "barrier heights". This should be understood as a consequence of selecting the separated reactants as a reference, rather than selecting the water trimer and HNCO molecule as a reference, for example. Figure 6 presents the ∆H 0K barrier heights for each process relative to the separated reactants. TSN and TSO refer to reaction across the N− −C and C− −O bonds and are depicted with blue and red bars, respectively. The number following this designation indicates the number of water molecules involved in the reaction. In all cases, additional water molecules significantly lower the enthalpy of the transition state, with the first catalyst water hav-ing the largest effect and lowering the transition state barrier by at least 20 kcal mol −1 . This is unsurprising and NBO E (2) results indicate significant delocalization between the additional water molecules, lowering the energy of the transition state ring structures. Catalytic water molecules have a greater impact on the TSO barriers, which results in the TSO3 and TSN3 barriers being reasonably competitive in some cases with a difference of less than 4 kcal mol −1 across all substituents. Figure S4 presents 298 K Gibbs energies for each transition state. One catalyst water molecule still significantly lowers the ∆G 298K barriers for each mechanism (TSN1 by at least 11 kcal mol −1 and TSO1 by at least 15 kcal mol −1 ), regardless of substituent . However, the energetic benefit of a second catalyst water is negligible as it is overcome by the loss of translational entropy. This is in qualitative agreement with the work of Wei and coworkers 52 who demonstrate convergence in the Gibbs energy of the reaction barriers with more than two catalyst waters.

RNCO Substituent Water Catalyzed Transition
The data in Figure 6 clearly indicate that the substituent selection plays a significant role in the energetics of the water catalyzed transition states. The previous relationship between substituent and barrier height in the H 2 O + RNCO case no longer holds with these different mechanisms, at least not for the same reasons. In the case of TSN2 and TSN3, there is a weaker correlation between the RNCO C−O* occupation and the barrier height, as the dominant E (2) interactions have changed with the new mechanism. We observe that in TSN2 and TSN3, the NBO method predicts a transition state where the nitrogen exhibits two lone pairs. The first lone pair is of almost entirely p-orbital character, while the other, from the broken N− −C bond, is significantly hybridized. The p-orbital character of this second lone pair correlates quite strongly with barrier heights (R 2 =0.93 and 0.92 for TSN2 and TSN3, respectively. See Figure S9.) and can be explained in terms of Bent's rules. Strongly electronegative substituents direct more p-orbital character from the nitrogen in the R-N bond, which increases the p-orbital character of the second nitrogen transition state lone pair and increases the proclivity of a N− −C double bond breakage. In the TSO2 and TSO3 cases, we cannot make the same argument, but we do note that the RNCO C−O* occupation correlates extremely well with the barrier heights (R 2 = 0.84 and 0.85 for TSO2 and TSO3. See Figure  S10.). There is not a clear single explanation for this relationship, but it does confirm that the electron withdrawing ability of the RNCO substituent is a reliable predictor for how much one or two water molecules catalyze imidic acid formation.

Catalyst RNCO
An additional RNCO molecule can also catalyze the reaction between water and RNCO to form carbamate or imidic acid through two-step mechanisms, as previously described by Cheikh and coworkers 50 (see Figure S8). The mechanisms towards severing the N− −C and C− −O bonds both begin with six-membered ring transition states, with the catalyst RNCO nitrogen abstracting a water hydrogen and its carbon attacking the other RNCO nitrogen (TSN A ) or oxygen (TSN A ), respectively. The IRCs for these processes again trace to initial van der Waals complexes which need not be explicitly considered. Proceeding through these tran- The second barrier for the carbamate formation process is particularly high (54.4 kcal mol −1 ), due to the energetic favorability of the INT-N species which seems to indicate that this multimolecular mechanism is likely not a catalyst for carbamate formation, as this barrier is higher than the uncatalyzed barrier through TS1 in the H 2 O + HNCO case. However, the barrier between TSO B and INT-O is only 39.5 kcal mol −1 , which might indicate that the imidic acid route is slightly favored with an additional isocyanate in the unsubstituted case. Additionally, the fact that the RNCO catalyzed reactions must proceed through two significant barriers would likely inhibit the catalytic efficiency of these pathways.
Substituents significantly influence the energetic landscape of the 2 RNCO + H 2 O reaction. The barrier heights (0 K enthalpies and 298 K Gibbs energies) are presented in Figures S5 and S6 for each substituent case. Generally speaking, the barrier to form the TSN A and TSO A six-membered ring transition states are decreased by electron withdrawing groups and increased by electron donating groups, similar to the water catalyzed cases. In the carbamate formation mechanism, the barrier through TSN A is significantly less than the barrier through TSN B for all substituents. The lowest barrier through TSN B is the case of R=F The RNCO catalyzed imidic acid formation pathway exhibits many similarities to the carbamate formation pathway. The first barrier (TSO A ) follows a very similar pattern with respect to substituent choice, except that the barriers are generally higher than TSN A , in each case by about 5 kcal mol −1 . The key difference in the imidic acid formation pathway is that the second barrier through TSO B is much smaller and is predicted to be between 34 and 42 kcal mol −1 across all the substituents considered. In every case, the energetic barriers from the respective intermediate structures through TSO B are much lower than the barri-ers through TSN B . This suggests that catalysis via an additional RNCO molecule could favor the imidic acid product over the carbamate product. This is certainly a consequence of the INT-O structure lying so much higher in energy than INT-N. Again, there are too many factors to provide a definite single explanation for how the substituents influence the TSO B barriers, but many of the same reasons mentioned in the previous paragraph are factors here as well.
Even though a mechanism exists for the RNCO + H 2 O reaction with an additional catalyst RNCO species, the energetic barriers likely remain too high to be significant without the involvement of other factors. The necessity of overcoming two substantial barriers, partially due to the favorable intermediate structures, implies that these pathways are not nearly as competitive as the water catalyzed pathways. Considering the 298 K Gibbs energy makes matters worse as the TSN A and TSO A barriers increase by about 20 kcal mol −1 due to the loss of translational entropy. With the Gibbs energy considered, the RNCO catalyzed pathway would hardly compete with the uncatalyzed RNCO + H 2 O reaction and would likely be far less efficient than catalysis with excess water based on our results. One should note however that Cheikh and coworkers studied the PhNCO + isopropylalcohol reaction and determined that both excess alcohol and PhNCO catalyzed the formation of the carbamate product 50 . They note that solvent effects can significantly influence the energy of the hydrogen transfer process through TSN B and TSO B which confirms that other factors beyond gas phase predictions are necessary for a full description of these systems. Nevertheless, our substituent analysis of the 2 RNCO + H 2 O reactions should be insightful for further research on these systems in more complicated environments.

Conclusions
Our work presents the highest level ab inito study to date of the important RNCO + H 2 O reactions, with the goal of guiding experimental progress of novel isocyanate containing reactions. We characterize the fundamental stationary points of the reactions of water across the HNCO N− −C bond to form carbamate and across the C− −O bond to form imidic acid, with energies targeting the CCSDT(Q)/CBS level of theory. Composite method geometry optimizations consisting of MP2 and CCSD(T) refined with large basis CCSD(T) single points, describe the influence of 24 RNCO substituents on the barriers to carbamate and imidic acid formation. NBO analysis reveals that the occupation of the in-plane C−O* orbital of RNCO is strongly associated with higher TS1 and TS2 barriers, due to the the electron delocalization motifs present in the transition states. The most extreme electron withdrawing substituent (R=F) lowered both barriers (∆H 0K ) by at least 13 kcal mol −1 and most alkyl substituents lowered the barriers by around 4 kcal mol −1 .
The catalytic influence of one extra water, two extra waters, and one extra RNCO species are predicted and discussed with respect to different substituted RNCO reactants. Electron withdrawing substituents significantly lower the barrier to carbamate and imidic acid formation in the water catalyzed cases. For the TSN2 and TSN3 barriers, the lower barriers are likely related to Bent's rule influencing the proclivity of the nitrogen to break the N− −C bond. We also find that the TSO2 and TSO3 barriers are highly associated with the C−O* occupation in the RNCO reactant. The 298 K Gibbs energies predict very similar barriers in the one and two water catalyst cases. Thus, the ability of the second catalyst water molecule to lower the enthalpy of reaction is offset by the loss of translation entropy manifest at higher temperatures and confirms the findings of Wei and coworkers 52 that many additional waters produce only marginal catalytic efficiency. A catalytic RNCO molecule results in two-step mechanisms towards carbamate and imidic acid formation. Our results indicate that the highly favorable intermediate in the carbamate formation pathway (INT-N) makes the second barrier large and an unlikely path for carbamate formation regardless of substituent. However, the second barrier in the RNCO catalyzed imidic acid formation pathway (through TSO B ) is significantly lower and might indicate more efficient imidic acid production in the presence of excess RNCO species. The complex factors involved in these RNCO catalyzed pathways did not reveal any clear trends with respect to substituent implying needs for additional work on a case by case basis for more complicated systems in order to more fully understand these RNCO catalyzed pathways. Our work lays a firm theoretical foundation for the RNCO + H 2 O class of reactions and provides insights that might guide future experimental work on these industrially relevant species.

Conflicts of interest
There are no conflicts to declare.