Abstract

A procedure for the theoretical study of chemical reactions in solution by means of molecular dynamics simulations of aqueous solution at infinite dilution is described using ab initio solute-solvent potentials and TIP3P water model to describe the interactions. The procedure is applied to the study of neutral hydrolysis of various molecules (HCONH2, HNCO, HCNHNH2, and HCOOCH3) via concerted and water-assisted mechanisms. We used the solvent as a reaction coordinate and the free energy curves for the calculation of the properties related with the reaction mechanism, namely, reaction and activation energies.

1. Introduction

The study, at microscopic level, of chemical reactions in solution is nowadays one of the topics of greatest interest in quantum chemistry [18]. One of the major difficulties in this type of study is the great number of interactions that take place, so that a full quantum mechanical treatment of both solute and solvent (using, for example, the CPMD method [9]) is costly, and recourse is made to approximate methods. Although there exist other procedures with an acceptable response to study the solvation of chemical systems (such as those which incorporate the effect of the solvent as a dielectric [10] or those which employed a mixed quantum/mechanical computational model, the QM/MM methods [11]), we shall use the classical method of molecular mechanics [12] using ab initio potentials to describe the solute-solvent interaction and the TIPnP potentials to describe the explicit solvent, allowing water molecules to take part in the reaction mechanism as takes place in a large number of the reactions in solution. The description of the solute-solvent interaction, however, will be based on the Lennard-Jones 12-6-1 analytical function whose molecular parameters are taken from fitting solute-solvent interaction energies calculated at the ab initio level [1317], instead of using the geometric-mean combining rules or parameter tables for model molecules. The solute charges are derived from a fit of the electrostatic component (ES) of the bimolecular solute-water interaction energy, which will henceforth be denoted as ESIE. With respect to the molecular parameters of the van der Waals terms in LJ-(12-6) interaction potential, we employed a fitting procedure based on using different components of the interaction energies in particular, the repulsion-exchange (EX), polarization (PL), and dispersion (DIS) components to describe the repulsive and attractive contributions.

Knowledge of the energy curves that guide fluctuations of the solvent is particularly important in the study of chemical processes in solution because it allows one to calculate the activation energy of the process without the necessity of determining the reaction path. This is of special interest as reactions in solution do not follow simple reaction mechanisms as occur in gas phase, but more complex mechanisms where the molecules of the solvent are part of the intermediates and transition structures. Outstanding among the different theories proposed is that which describes the free-energy curves of the reactant () and product () states as quadratic functions of the solvent coordinate [1821]. The transition state is associated with the point of intersection of the two curves, and the reaction energy is obtained directly from relationships between points of these two curves (see Figure 1). One of the most commonly used reaction coordinates is the difference in the solute-solvent interaction energy of a given set of solvent molecules in the presence of the reactant and product structures [22], for which one only needs the potential function that suitably describes this interaction.

The present work takes part of a series of studies focused on proposing a solute-solvent interaction potential, in this case determined from the EX-PL-DIS-ES components, to describe chemical reactions in aqueous solution. The principal objective of the work is the application of the methodology based in the construction of free-energy curves to hydrolysis reactions of organic molecules in which the concerted mechanism may be assisted by the aqueous solvent. We shall also study how the activation barrier is lowered when a second water molecule, part of the solvent, is included in the transition state.

2. Formalism and Calculation Details

About a thousand values of the SCF and MP2 solute-solvent interaction energy were calculated with the 6-311++ basis [29, 30]. In order to try to appropriately describe the attractive, repulsive, and long-range interactions, the grid of points was generated by placing the water molecule at different positions relative to the solute. The 6-311++ basis set used contains polarization and diffuse functions for all of the atoms in order to improve the description of the outermost orbitals, and hence of the reactant and transition state energies and geometries.

These energies were used to obtain the , , and interaction parameters of the potential function:

The net charges on each solute atom were obtained with the ESIE procedure [13, 16], fitting the values of the Coulomb component of the interaction energy using the variational scheme of Morokuma and coworkers [31, 32]: where the charges of the solvent water are preassigned as the TIP3P charges [33].

The Lennard-Jones parameters and were obtained in a similar way to , but now the energies used in the fits were those that describe the exchange (EX) and polarization (PL) components of the interaction energy at the SCF level, and the dispersion (DIS) component related to the MP2 correlation energy [13, 16]:

To construct curves of the free energy G, we used as reaction coordinate the solvent fluctuation, that is, the difference in the interaction energy of a given set of solvent molecules in the presence of the reactant and transition state structures [22], for which one only needs the potential function that suitably describes this interaction.

Thus, to obtain the free-energy curve associated with the reactant simulation, one can use the differences in the solute-water interaction energies, (), between the diabatic states of the solute in its reactant () and transition state (TS) structures for a broad set of configurations of solvent molecules around the solute in a molecular dynamics simulation of the reactant solvation: where denotes the solute-solvent interaction , at the reactant structure and the solvent configuration W.

One operates in a similar way to simulate the solvation of the transition state and to obtain the values with respect to the other system, considering the same direction in all the cases, that is, .

The difference fluctuates during the S-simulation, and its values are collected as a histogram of the number of times that a particular value of the macroscopic variable appears in the simulation. The probability of finding the system in a given configuration can be expressed in terms of the delta function , where represent the fluctuations of the energy in the state and of the number of equally spaced steps in the simulation:

This allows us to compute the free energy :

The values of free energy obtained presents some dispersion in areas far from equilibrium, so it is advisable to make a search for the polynomial function that best fits these free energies , and the result is plotted. In all cases, the separation between the two minima representing equilibrium structures, when the curve intersects the curve at its lowest point , is calculated as with and being the most probable values of in the free-energy curves and , respectively, and a, b, and are the coefficients of the polynomial fit to the curve .

Finally, some details of the simulations merit mention. Molecular dynamics simulations of an NVT ensemble of a solute molecule in an aqueous environment represented by about 200 water molecules were carried out at 298 K using the AMBER program [34]. The time considered for the simulations was 2000 ps with time steps of 0.1 fs. The first 1000 ps were used to ensure that equilibrium was reached completely, and the last 1000 ps were to store the configurations of the water molecules required for the determination of the thermodynamic properties studied in this work. Those water molecules initially located at distances less than 1.6 Å from any solute atom were eliminated from the simulations. Long-range electrostatic interactions were treated by the Ewald method [35], and the solutes were kept rigid using the shake algorithm [36]. A cut-off of 7 Å was applied to the water-water interactions to simplify the calculations, and periodic boundary conditions were imposed to describe the liquid state. The grid of points used to fit the interaction potential to the Lennard-Jones 12-6-1 function was obtained with SCF and MP2 energies using the Gaussian/98 package [37], and the decomposition of the interaction energies was performed with the Gamess program [38]. The calculations were made on a QS16-2500C-X64Q QuantumCube multiprocessor computer at Extremadura University.

3. Reactions and Mechanism Studied

The selected reactions are typical examples of hydrolysis reactions. These reactions are of great interest because of their presence in numerous processes that occur in living organisms. Thus, the hydrolysis of formamide is a reaction of great importance because it can be considered as a model for the study of peptidic bonds reactivity. The hydrolysis of isocyanate is a reaction of interest for the polymer industry, the pharmaceutical industry, and agriculture. The hydrolysis of the amidine derivatives and carboxylic esters is interesting because of their biochemical and medical importance, being relevant to biology and pharmacology, as well as to study the hydrolysis of formamidine and methyl formate.

In general, these reactions are of the type A+ nH2O C + D + () H2O, that it to say, the addition of water to isocyanate, formamide, formamidine, and methyl formate molecules, taking into account the possibility of the participation of water molecules in the reaction mechanism (Figure 2). The water molecules can act in an intramolecular fashion, taking part in the transition state formation and helping to stabilize strained rings (hydrogen isocyanate TS1, formamide TS1, formamidine TS1), or they can act as a solvent, stabilizing the transition state through hydrogen bonds (methyl formate TS3 and formamidine TS2). In this case, the position of water molecules was optimized in the transition state and kept fixed during the molecular dynamics simulation.

In the studied mechanisms involving one or more water molecules in the reactant and transition system, the reaction takes place when the oxygen of a water molecule attacks nucleophilically the reactant, and there is a simultaneous proton transfer from the water hydrogen to the next atom of the reactant (i.e., TS1 in Figure 2 for the formamide).

In the decomposition reactions of formamidine and hydrogen isocyanate, the mechanism is somewhat different. Thus, the water molecule involved in the reaction assists the proton transfer process by transferring a hydrogen atom from the reactant to the water molecule and, simultaneously, a hydrogen from water is transferred to another zone of the reactant molecule, forming a six-membered ring, much more stable than the four-membered ring that would be formed without the presence of water (i.e. TS2 in Figure 2 for the formamidine).

Theoretical and experimental studies of these reactions have been carried out by several authors. Thus, in the study of the hydrolysis reaction of formamide, a number of theoretical studies have been devoted to the neutral hydrolysis, among which it is worth highlighting the work by Almerindo and Pliego [39], Antonczak et al. [40] and the studies by Estiu and Merz [41] and Gorb et al. [42]. As for the isocyanate hydrolysis, we can find numerous experimental and theoretical studies [28, 4347]. The reaction of the formamidine has also been studied both experimentally [48] and theoretically [24, 4952] by several authors. Both the gas-phase [49, 50, 52] and aqueous phase reactions, assisted by one or more molecules of water [24, 25], have been studied. The methyl formate molecule, due to its small size, has been employed as a test case in both experimental and theoretical investigations. Its hydrolysis has been extensively studied experimentally [26, 52, 53] and computationally [5456].

4. Results and Discussion

4.1. Energies

Ground state geometries of the reactant and transition state in solution that appear in Figure 2 were optimized at the MP2 level with the 6-311++G** basis starting from standard geometries and using the PCM model. All transition-state structures were confirmed by checking that there was one negative eigenvalue in the diagonalized hessian corresponding to the movement along the reaction path.

To calculate the reaction and activation free energies, we constructed the , TS, and diabatic free energy curves following a procedure similar to that employed by Ando and Hynes for acid ionization in water [57, 58]. The relative positions of the curves were corrected in accordance with the work of Tachiya [59], so that the curve is vertically shifted at the point of minimum energy of the transition state curve . The procedure is similar for the activation energy of the inverse reaction except that now the curve that is shifted is . The results from these curves are given in Table 1 together with the available data.

We can see that in the majority of the cases the mechanism assisted by solvent water molecules shows lower activation barriers. The reason is the formation of a ring in the transition state, which is more stable than when no water molecules participate in the reaction (TS1 versus TS2 case in formamide). In cases where the transition state with a water molecule forms a quite stable ring, the presence of another water molecule does not reduce noticeably the value of the activation energy (formamidine TS2 case), since this water molecule is not part of the transition state structure and it only assists the reaction, by forming hydrogen bonds. Comparing the TS1 structures of the reactions of formamide and formamidine assisted by a single water molecule, we find that the free energy of activation is significantly lower in the case of the formamidine (22.60 kcal/mol versus 51.28 kcal/mol), due to the lower tension in the ring formed in the transition state. In the case of the formamide reaction this tension is further reduced by introducing a second water molecule, which lowers the activation energy to a value closer to that of the formamidine

The participation of the second molecule of water in the mechanism of the formamide reaction decreases the energy barrier, as shown in Figure 3, leading to a result in good agreement with the experimental value of 31 kcal/mol.

As for hydrolysis of the formamidine (Figure 4), bibliographic values [24, 25] for the mechanisms assisted by a water molecule [24, 25] are available, showing a good agreement with those obtained using our methodology. Similar conclusions can be achieved when comparing the results for the mechanism assisted by two water molecules (TS2), although in this case the activation barrier remains practically constant when the free energy curves are used, since this second water molecules is not part of the structure and only assist the reaction by forming hydrogen bonds (Figure 4).

Similar results can be found for the reaction of isocyanate and methyl formate, finding values for the free energies of activation close to those from other authors when water molecules are considered to assist the reaction mechanism [2628]. In these cases, we considered 2 and 3 molecules of water, respectively. The presence of more water molecules has been observed that does not significantly changes the free energy of activation because they are not involved directly in the reaction mechanism.

To check the validity of our method MD/ESIE, we also calculated the activation energies for some of the previously discussed reactions taking the solvent as a continuum media, for which we used the PCM model implemented in the Gaussian98 package [37], considering the transition state with the explicit water molecules that take part in the reaction mechanism and taking account the nonelectrostatic components of the energy. The results obtained at level MP2 are shown in Table 2, together with the activation energies calculated using the proposed model (MD).

We find that the activation energy depends on the model used in the description of the solvent. For example, in the case of formamidine, assuming that the reaction is assisted by water molecules and considering the solvent as a continuum, the barrier to reaction is 38.70 kcal/mol. This barrier is lowered to 22.66 kcal/mol using ESIE charges and free energy curves, according to the proposed methodology. The later result of 22.6 kcal/mol is in agreement with that obtained by Nagaoka et al. [25] of 15.69 kcal/mol obtained from a chemical reaction dynamics in solution study.

Similar conclusions can be reached upon the study of the methyl formate reaction. The energy barrier of 37.73 kcal/mol from gas-phase MP2/6-311++ calculations can be improved appreciably by using our proposed MD procedure, that leads to free energies of activation of 28.17 kcal/mol in good agreement with the experimental values [26, 27] of 25.9 kcal/mol.

Finally, in the case of isocyanate, we calculated the activation energy of the direct and reverse reaction to check the ability of the method to calculate free energies of reaction. In this case, we obtained a value of −11.10 kcal/mol for the reaction energy, which can be compared to the value provided by Raspoet et al. of  kcal/mol.

The clear conclusion of this work is that the free energy curves of the species involved in the reaction provide a good description of the reaction thermodynamics in an aqueous medium. These curves, which are constructed based on the solute-solvent interaction energies with the solvent fluctuation being chosen as reaction coordinate, respond acceptably to the activation barrier of the process. Also, to obtain reasonable results for these reactions using the geometries in solution, one must take into account a mechanism assisted with water molecules, since one observes that the activation barrier for the process decreases appreciably when water molecules are involved in the reaction.

Acknowledgment

This research was sponsored by the Consejería de infraestructuras y Desarrollo Tecnolçogico de la Junta de Extremadura (Project GRU10036).