Discussion

The present results show that Car-Parrinello MD simulation is a powerful tool to assess solvent effects on a simple S$_\mathrm{N}$2 reaction. Semi-quantitative results for the energetics and a detailed picture of the structural aspects of the reaction are obtained. However, the study also shows that at present time the calculations are not without limitations. The most important ones, the relatively small box size and the short trajectories that can be calculated, cannot easily be removed because they are related to the high computational expense. For example, the time necessary to calculate one point from figure 3.5 (say, 6 ps total time for equilibration and sampling) was about 11 days on 6 IBM-SP nodes (power2sc processor) or 29 days on a cluster of 16 IBM RS6000 43P workstations.

The hysteresis found in the force of constraint is due to a memory effect in the solvent, which is directly related to the fact that the last configuration of every constrained AIMD run is used to construct the initial configuration for the simulation at the next point of the reaction coordinate, in combination with the limited time scale accessible to an AIMD simulation. Apparently, the rearrangements that have to be made by the solvent, to follow the change in the reaction coordinate from point to point, occur on a much longer time scale than the rotations and vibrations and diffusion of water molecules that are well described by our 10 ps molecular dynamics trajectories.

Simulating for longer times is the trivial solution to avoid memory effects. However, this could be as impractical as waiting for the reaction to occur spontaneously because there are indications that the adaptation of the aqueous solvent to the changing solute configurations is an activated process instead of a slowly diffusive process. In the 9 ps transition state simulation, for instance, a drift in the force of constraint would be expected, which should asymptotically decay towards zero, in the case of a slowly diffusively adapting solvation shell, but was not observed. Note also that, although sampling for extremely long times for each point $\xi^\prime$ will in principle remedy the hysteresis in the free energy profile, a too low barrier is found as these activated solvent adaptations are missing in the free energy profile. For comparison with the experimental rate constant, this is corrected by calculation of the pre-exponential factor, which includes the transmission coefficient $\kappa$, but the computational demand of this exercise increases dramatically with the amount of the underestimation of the free energy barrier (see further appendix 3.5). Secondly, we note that the force-field Monte Carlo simulations of the same S$_\mathrm{N}$2 reaction in 250 water molecules by Chandrasekhar et al.[65], using the umbrella sampling technique, also shows a significant hysteresis. The calculated probability distribution $P(\xi=\xi_{\mathrm{TS}})$ (figure 3.5 in ref ChSmJo) is not symmetric: the second peak at the reactant side next to the transition state is absent at the product side. This is also reflected in the radial distribution functions obtained from the MC simulation of the transition state.

It has become clear that certain necessary changes in the solvent as the reaction takes place are rare events on the typical AIMD time scales. One way to handle the problem of hysteresis is therefore to include the required solvent rearrangements into the reaction coordinate. Unfortunately, due to the large number of molecules involved and their complex rearrangements, finding a proper reaction coordinate is virtually impossible. The need to include solvent degrees of freedom into the reaction coordinate was also concluded recently by Geissler, Dellago and Chandler from their molecular simulations of the dissociation of NaCl in water[125]. In this work, they found additional free energy barriers in the solvent that have to be overcome as the reaction (i.e. the dissociation) occurs. For instance, the addition of a water molecule in the first solvation shell of Na$^+$ to bring the coordination from the initial five-fold to the final six-fold, required an amount of work equal to 1.7 $k_{\mathrm{B}}T$ ($\approx$1 kcal/mol)). During the S$_\mathrm{N}$2 reaction in the present work, a Cl$^-$ solvation shell has to lose two water molecules to reach the transition state, while at the same time the CH$_3$Cl chloride has to bond with three water molecules which have to be taken from the water network.