Conclusions

Using the Car-Parrinello method, we studied the solvation effects on the identity reaction between Cl$^-$ and CH$_3$Cl in a dilute (1:32 mole ratio) aqueous solution of HCl. We compared the free energy profile of the reaction in the gas phase and in solution and analyzed the structure of the solution at different stages of the reaction. The free energy profile of the reaction in aqueous solution was obtained using the method of constrained molecular dynamics simulations in combination with thermodynamic integration. The calculated barrier yields 27 kcal/mol after application of certain corrections (see below). This corrected value is in good agreement with the experimental value of 26.6 kcal/mol.

There are two important error sources in these calculations. The DFT (Becke-Perdew functional) electronic structure calculations describe the solvent and the solvent-solute interactions, in particular the solvation of the Cl$^-$ ion, the formation of the initial ion-dipole complex of Cl$^-$ with CH$_3$Cl, and the solvation of the reaction system in the course of the reaction, with reasonable accuracy. However, although DFT-BP does exhibit an energy rise when going from the ion-dipole complex to the transition state, thus reproducing the well-known double well energy profile for the gas phase reaction, it underestimates the barrier height of ca. 13 kcal/mol by as much as 8 kcal/mol. This error, being related to an imbalance in the description of exchange and correlation in the transition state by present day generalized gradient approximations like Becke-Perdew, cf. ref gritsenko, will persist in the transition state in solution. A correction of 8 kcal/mol has therefore been applied to the calculated barrier height. In the second place a "hysteresis", or inequivalence of forward and reverse reaction, has been found in the force of constraint profile. The hysteresis is due to a slow adaptation of the solvent to the reacting solutes and has the effect of pulling the reactants back toward the initial state when the constrained reaction coordinate is moved toward the product side. Apparently, the rearrangements that have to be made by the solvent occur on a much longer time scale than the rotations and vibrations and self-diffusion of water molecules that are well described by our 10 ps molecular dynamics trajectories. On the time scale of MD simulations, the rearrangements needed in the solvation shells of the reactants are rare events by themselves. Including solvent degrees of freedom in the reaction coordinate is in principle the way to handle this type of the rare events.[125]

The Car-Parrinello MD simulation of the reactants in a box with 32 water molecules is an important improvement on the numerous micro-solvation and dielectric continuum calculations. It is a great advantage that the varying solvation during the course of the reaction, due to the shift of negative charge from the attacking chloride initially to the leaving chloride finally, does not have to be described by parameterized model potentials that cannot be optimal for all occurring situations along the reaction coordinate. From the good agreement with experiment - after the corrections described above - we may infer that both the energetic and the entropic effects of the solvation on the free energy barrier are obtained to quite satisfactory accuracy with the AIMD simulation. Further improvements of this type of calculation can be envisaged, of both technical and more fundamental nature. At the technical side, we note that the limited size of the total system implies that approximations are made on the long range interactions in the solvent. Also the fixed volume has a small effect on the forces because the changing size of the reaction complex in solution affects the pressure on the system. The extension of the present type of calculation to larger systems will certainly be possible in the future due the rapid increase in available computing power. The most important deficiencies of the present work are of a more fundamental nature. The first is the underestimation of the transition state energy by the GGA (Becke-Perdew) functional. This will have to be remedied by the development of more accurate functionals, which, in view of computation time, should not rely on the incorporation of exact exchange. In the second place, it is desirable that the dynamics methods are improved in order to keep the necessary simulation times down to a manageable length. One possibility would be the development of methods to handle the slow adaptation of the solvent structure to the changing reaction coordinate by incorporating the necessary change in solvent structure in the reaction coordinate.