J. Phys. Chem. B 2007, 111, 5243-5248
5243
Car-Parrinello Molecular Dynamics Study of Anharmonic Systems: A Mannich Base in Solution Aneta Jezierska,† Jarosław Panek,† Urban Borsˇtnik,‡ Janez Mavri,‡ and Dusˇanka Janezˇ icˇ *,‡ UniVersity of Wrocław, Faculty of Chemistry, 14 F. Joliot-Curie, 50-383 Wrocław, Poland, and National Institute of Chemistry, 19 HajdrihoVa, SI-1001 Ljubljana, SloVenia ReceiVed: December 18, 2006; In Final Form: February 23, 2007
A Car-Parrinello molecular dynamics study was performed for 4,5-dimethyl-2-(N,N-dimethylaminomethyl)phenol, a Mannich base, to investigate the vibrational properties in solution of its intramolecular hydrogen bond. The dynamic behavior of this hydrogen-bonded system was investigated using an explicit solvent model. Addition of a nonpolar solvent permitted inclusion of delicate environmental effects on the strongly anharmonic system which was studied from first principles. Molecular dynamics and a posteriori quantization of the O-H motion were applied to reproduce the vibrational features of the O-H stretching mode. Consistent application of Car-Parrinello dynamics based on the density functional theory with subsequent solution of the vibrational Schro¨dinger equation for the O-H stretching motion offers an effective method for strongly anharmonic systems, and this is supported by the comparison of the results with experimental spectra. As a further element of the intramolecular hydrogen bond study, the effects of deuteration were taken into account and a successful application of the O-H stretching mode quantization technique to the liquid phase is demonstrated. This provides a valuable computational methodology for investigations incorporating nuclear quantum effects in the liquid phase and enzyme active centers and can be used to investigate numerous systems that are not readily susceptible to experimental analysis.
Introduction Proton dynamics play a crucial role in many chemical, biochemical, and industrial processes1-4 because most of the properties of the water, which is the solvent in these systems, can be traced to intermolecular hydrogen bonds.5-10 The chemical and biochemical implications of intramolecular hydrogen bonds, however, are no less relevant to an understanding of atomic scale mechanisms11-14 such as proton transfer, an essential component of many enzymatic reactions.15-19 Hydrogen bonding, a ubiquitous chemical phenomenon, exhibits great variability of parameters with respect to the participating atoms, the environment, and the thermodynamic state of the system.20-22 One such example is the influence of the π-electron coupling between donor and acceptor centers. This can play a role in the intramolecular hydrogen bonds formed by moieties linked to an aromatic ring. Such coupling is seen in 3-hydroxyflavones, benzoxazoles, and Schiff bases.23-27 Attenuation of this π-electron coupling by inert spacers leads to Mannich bases such as o-hydroxy-N-alkylaminophenols.28,29 Such moderately strong and easily modified hydrogen bonds have potential industrial and material chemistry implications, and investigations of the properties of this class of compounds, seen from the perspective of basic physicochemical research, can provide data on the physical and chemical nature of the hydrogen bonds. The o-Mannich base, 4,5-dimethyl-2-(N,Ndimethylaminomethyl)phenol (DMM), was chosen for this study.28 This is a member of a class of compounds with intramolecular hydrogen bonds in which the hydrogen donor * Corresponding author. Phone: +386 1 476 0321. Fax: +386 1 476 0300. E-mail:
[email protected]. † University of Wrocław. ‡ National Institute of Chemistry.
and acceptor are separated by a methylene bridge and thus are only weakly coupled electronically.28-32 Computational investigations of the proton dynamics are complicated because of the anharmonic potentials involved and the low mass of the proton. This report describes calculations combining Car-Parrinello Molecular Dynamics (CPMD) with a posteriori quantum treatment of the proton motion in order to reproduce the vibrational features of an anharmonic hydrogenbonded system. In CPMD, electronic degrees of freedom are introduced as dynamical variables, leading to a system of coupled equations of motion for both ions and electrons. In this way, the two difficult obstacles mentioned above can be ameliorated. The method based on the density functional theory generates a time history of the system which incorporates environment fluctuations, and second, quantum corrections which are included in the computational model allow the possibility of generating vibrational spectra or studying kinetic isotope effects with relatively low computational cost. Accurate representation of the proton dynamics (e.g., vibrational features of the hydrogen bond) requires accurate potential energy functions, adequate sampling of the molecular potential energy surface (PES), proper inclusion of the environment and, in many cases, quantum treatment of the nuclear motion, at least for light nuclei such as protons. These requirements, respectively, can be met by expensive correlated ab initio approaches, systematic PES scanning or statistical methods such as molecular dynamics,33-35 Monte Carlo sampling, and nuclear quantum corrections using path integral techniques36,37 or a posteriori corrections.38 Several recent applications of molecular dynamics in hydrogen bridge investigations by MD sampling of proton motion39-41 promise to replace systematic, but tedious, sampling of the PES with the recording of time-dependent proton potentials deter-
10.1021/jp068676p CCC: $37.00 © 2007 American Chemical Society Published on Web 04/21/2007
5244 J. Phys. Chem. B, Vol. 111, No. 19, 2007 mined by molecular dynamics. This approach is a very efficient means by which both intramolecular fluctuations and environmental effects can be described. Proper treatment of proton motion requires its description in quantum mechanical terms. Methods of mixed quantum/classical dynamics have been developed (density matrix evolution (DME), surface-hopping, centroid molecular dynamics),42-47,36 but they are too computationally demanding to be used routinely in the context of CPMD. The computational methodology employed here can be briefly defined as the extraction of snapshots of proton potentials from the CPMD trajectory followed by solution of the onedimensional Schro¨dinger equation describing the O-H motion. Solving the vibrational Schro¨dinger equation for potential snapshots in the context of the Car-Parrinello scheme beyond the harmonic approximation represents our novel aspect. In this way the O-H stretching envelope is obtained as a bar spectrum corresponding to the 0f1 transition. Thus, by considering the quantum nature of the motion, this leads to the envelope of a selected mode of a strongly anharmonic system at a cost which is comparable to that of a plain CPMD simulation. In summary, this study is an application of a general theoretical model to the investigation of proton dynamics in a fluctuating environment. The system investigated in this study, a Mannich base in solution in CCl4, was selected for several reasons. The Mannich base under consideration has a moderately strong intramolecular hydrogen bond, but there is no experimental evidence for continuous proton shuttling between the N and the O atoms. Consequently, the anharmonicity will be substantial, and while quantum effects associated with the proton motion may be visible, they should not dominate the system. Moreover, the liquid environment, with its inherent flexibility and disorder is a good object for study by molecular dynamics techniques. The available experimental data for the Mannich base studied was obtained in solution in CCl4, permitting direct comparison with calculated results. Computational Methodology Car-Parrinello Molecular Dynamics48 based on density functional theory49 was used to reproduce the molecular properties of 4,5-dimethyl-2-(N,N-dimethylaminomethyl)phenol in a relatively nonpolar solvent, CCl4.28 Version 3.9.2 of the CPMD program50 was used, and calculations were performed on computer clusters and multiprocessors to complete the timeconsuming calculations in the liquid phase.51-53 The exchangecorrelation functional proposed by Perdew et al. coupled with a plane-wave basis set was applied during the study.54 The current formulation of DFT does not explicitly include dispersion forces, and the amount of dispersion recovered by implicit parametrization is disputed. However, we are interested in a spectroscopic description of the solute, in particular the O-H mode, which is dominated by a relatively strong intramolecular hydrogen bond. This turns dispersion interactions with the solvent into secondary interactions in our model. Therefore, we chose to retain the PBE functional model, which has been shown to perform reasonably closely to the reference MP2 and coupled cluster methods for hydrogen-bonded dimers.55 The ultrasoft pseudopotentials based on Vanderbilt’s scheme56 were applied to replace the core electrons in all simulations, the time step was consistently set to 3 au and initial structure optimization was carried out with the Schlegel Hessian.57 Subsequently, onthe-fly sampling of the PES and the small polarizing effect on the solute hydrogen bond of the continually moving, electron-
Jezierska et al. rich CCl4 molecules were incorporated by means of molecular dynamics. The MD calculations were performed at 300 K with allowed variations of (50 K. The liquid-phase MD simulation used a cubic box with a side of 16.35 Å and 24 CCl4 molecules. The fictitious electron mass (EMASS) parameter was set to 400 a.u., while the plane wave cutoff was 30 Ry. Up to eight neighboring cells (TESR ) 8) in each direction were included in the electrostatic interaction evaluation. In the initial step of the MD run, the box was equilibrated, and a further 5.2 ps of molecular dynamics was used for data collection. Two approaches were applied to the analysis of the proton position and the vibrational behavior of the hydrogen bridge. First, a standard autocorrelation of atomic velocity study provided direct information from the molecular dynamics data. The use of the ultrasoft Vanderbilt pseudopotentials to replace the core electrons prohibits subsequent dipole moment collection during the MD run. Berry phase algorithms for dipole moment calculations are not implemented in this case. Fortunately, although we lose information on the infrared intensities, the band positions are conserved in the atomic velocity analysis. However, there is also a fundamental reason for our choice of the atomic velocity power spectrum method; only in this way can we capture the dynamic nature of processes at the molecular level. Static calculation of the IR or Raman spectrum would lose, for example, different frequencies associated with the bridge proton motion due to the modulation of the O-H potential by varying the O‚‚‚N distance. Second, the nuclear effects of proton motion were included based on a consistent approach using Car- Parrinello dynamics for potential energy surface sampling and subsequent a posteriori nuclear motion corrections by solving the Schro¨dinger equation. In this approach, the trajectory obtained from the MD simulations was sampled at 250-fs intervals. The resulting snapshots were used as the starting points for the generation of one-dimensional proton potential functions for the motion of the hydrogen atom from the donor to the acceptor along a circular arc.58 This approach eliminates those arc motion end points with excessive energy. The potential function is then used for the solution of the nuclear Schro¨dinger equation to give the vibrational energy levels,59 and the set of generated vibrational frequencies is used to define the envelope of the O-H stretching mode. The deuterium isotope effect was also calculated.60,61 The procedure of one-dimensional nuclear motion quantization was repeated with the particle mass adjusted for deuterium, thus both O-H and O-D envelopes were generated. Molecular visualizations were performed with the VMD program.62 Results and Discussion The molecule under study has an intramolecular hydrogen bond and belongs to a class of compounds (Mannich and Schiff bases) well-known for their vibrational features,63 which are strongly influenced by their environment. Mannich bases seem to be ideal molecules with which to test a methodology proposed for accurate representation of the dynamics of strongly correlated systems, such as proton motion in the hydrogen bridge. From a classical dynamics trajectory analysis it is possible to proceed to a quantum description of the nuclear motion, accounting simultaneously for the influence of the solvent. Figure 1 shows the molecular structure of the Mannich base studied and the liquid-phase model used during the theoretical investigations. We found the run length of 5.2 ps for data collection to be sufficient for our purpose. It is known experimentally that the Mannich base under study does not exhibit proton-transfer phenomena. Additionally, the relatively inert solvent molecules
CPMD Study of Anharmonic Systems
J. Phys. Chem. B, Vol. 111, No. 19, 2007 5245
Figure 1. 4,5-Dimethyl-2-(N,N-dimethylaminomethyl)phenol (DMM) and the computational model used for the explicit solvent CPMD calculations.
Figure 2. Time evolution of the interatomic distances in the hydrogen bridge of 4,5-dimethyl-2-(N,N-dimethylaminomethyl)phenol (DMM) during the CPMD liquid-phase calculations.
will not form hydrogen bonds or any other directional interactions with the solute. On the other hand, the simulation time covers hundreds of O-H stretching cycles and several tens of slow O‚‚‚N bridge motions. In case of more polar or active solvents (e.g., water, methanol), longer simulation would indeed be desirable for ensuring adequate sampling of solute-solvent configurational space with possible intermolecular hydrogen bond breaking/forming events. As the first step in data analysis, the interatomic distances in the hydrogen bridge were analyzed as a function of time, and
the graphical representation of this time series is presented in Figure 2. No proton transfer was observed; the proton remained at the donor site during the entire computational study. During the simulation, the classical amplitude of the donor-proton (OH) motion was 0.20 Å, with strong variations on the short time scale. There is a period, from 3.6 to 4.1 ps, of larger O-H motions, with a subsequent phase of very small vibrations. These variations do not seem to be correlated with large-amplitude N‚‚‚H and N‚‚‚O motions although the latter appear to follow the same pattern. Interestingly, the N‚‚‚O motion has a smaller
5246 J. Phys. Chem. B, Vol. 111, No. 19, 2007
Jezierska et al.
Figure 3. The in-plane and out-of-plane lateral motions of the hydrogen bridge proton in the DMM compound studied. The plane is defined by atoms C10, O1, and N3.
Figure 4. The power spectra of atomic velocities for all atoms of the Mannich base and the hydrogen bridge proton. The intensities are in arbitrary units, while the wavenumbers correspond to the actual vibrational frequencies of the compound.
amplitude (0.728 Å, from 2.407 to 3.135 Å) than the N‚‚‚H motion (1.0034 Å, from 1.3531 to 2.3565 Å). This signifies that the proton can exhibit lateral motions (see Figure 3) in or out of the plane of the six-membered ring formed by the O1H2‚‚‚N3-C4-C9-C10 atoms. The in-plane O-H‚‚‚N angle deviates up to 20° from 180° as is often the case for
intramolecular, strained hydrogen bonds forming pseudo rings. The out-of-plane distortion is an additional form of relief of the strain induced by nonlinearity of the hydrogen bond. The power spectrum of the atomic velocities provides information on the frequencies in the system during the CPMD simulation. These frequencies, presented in Figure 4 and Table
CPMD Study of Anharmonic Systems
J. Phys. Chem. B, Vol. 111, No. 19, 2007 5247
TABLE 1: Comparison between Experimental28 and Calculated Vibrational Parameters for the O-H Stretching Mode in the 4,5-Dimethyl-2-(N,N-dimethylaminomethyl)phenol (DMM) Molecule in Solution band range [1/cm] mode
experimental
O-H stretch O-D stretch
2200-3400 1700-2600
atomic velocity power spectrum
quantized O-H stretching
2300-3400
1600-3400 1200-2400
1, correspond to the vibrational features of the system. A continuous slope, which vanishes at 500 cm-1, underlies the low-frequency part of Figure 4 and is due to the slow motion of the solvent molecules. The proton stretching region exhibits some sharp peaks associated with the C-H motions, but the O-H band is substantially broadened and is located between 2300 and 3400 cm-1 with the most intense O-H motions between 2600 and 2800 cm-1. Thus, even in the nonpolar solvent the proton potential function is strongly modulated, primarily by the N‚‚‚O distance and secondarily by the environment and the residual intramolecular oscillation. The quantum treatment of the O-H stretching mode, performed according to the snapshot extraction methodology, permits calculation of both the O-H and the O-D band envelopes using the same trajectory in order to estimate the isotope effect in the IR spectrum. The envelopes obtained are shown in Figure 5. From the complete comparison of the computed O-H stretching position with the available experimental data shown in Table 1, the spectral ranges obtained via the applied theoretical methodologies can be seen to be in good agreement with experimental liquid-phase O-H assignment.28
based on a plane-wave theoretical background, evidently provides an accurate potential energy surface of the system. It can therefore be used to generate instantaneous one-dimensional potential functions for selected atoms. These time-dependent functions may be used to solve the Schro¨dinger equation, and they provide an efficient method for incorporation of quantum corrections. This approach is not limited to hydrogen bond investigations of small molecules, and other possible applications include strongly anharmonic motions in the condensed phase and the gas phase (floppy molecules). Extension of the methodology to two or more dimensions is conceptually straightforward, although computationally demanding. The method applied here could be useful for nuclear tunneling investigations in enzymes either by CPMD simulation or by using one of the established QM/MM approaches. The approach used in our manuscript is intended to be a relatively inexpensive way of introducing a posteriori quantum corrections to the classical CPMD trajectory. It is thus in some respect akin to the quantum correction factor (QCF) technique. However, there are significant differences between the two methods. Generally, the modes that are being corrected by QCF are intramolecular or low-frequency rotational features in systems such as liquid neon or diatomic molecules.64-67 If quantum corrections are used for the hydrogen-bonded systems, either their use is restricted to the low-frequency non-stretching modes (for example, the far IR spectrum of water, up to 1200 cm-1)68 or a special theoretical treatment is required, thus limiting the generality of the model.69 The approach used in our manuscript restricts quantization to one mode of particular interest, but it is general and tuned for the needs of the vibrational spectroscopy. Therefore, we consider this approach to complement, rather than replace, the QCF technique.
Perspectives The method used in this study is very effective for the description of the dynamics of strongly correlated systems such as hydrogen bonds. The Car- Parrinello Molecular Dynamics,
Summary Hydrogen bridge proton dynamics in a Mannich base was studied by combined techniques of CPMD and nuclear motion
Figure 5. The simulated O-H and O-D stretching modes based on the snapshot method.
5248 J. Phys. Chem. B, Vol. 111, No. 19, 2007 quantization. The vibrational properties of the hydrogen bridge in this compound were accurately reproduced by this extension of the conventional CPMD scheme, which retains access to all the dynamical data. With an additional expense of less than 10% of the CPU time spent on the CPMD run, the method provides an alternative and attractive route to the estimation of the influence of nuclear quantum effects on the vibrational properties of a system. The successful application of the method to this solvated system implies that it can be used with gasphase, solid-state, and macrobiological systems, once an accurate molecular dynamics potential energy surface is generated using, for example, one of several standard QM/MM schemes. This will greatly facilitate studies of hydrogen bond dynamics in biological systems. Acknowledgment. A.J. and J.J.P. gratefully acknowledge the Wroclaw Centre for Networking and Supercomputing (WCSS) and Poznan´ Supercomputing and Networking Center for providing computer time and facilities. In addition, these calculations were partially performed using the CROW dual Opteron cluster at the National Institute of Chemistry, Ljubljana, Slovenia.51 The authors acknowledge the financial support from the state budget by the Slovenian Research Agency under Grant No. P1-0002. References and Notes (1) Desiraju, G. R.; Steiner, T. The Weak Hydrogen Bond. In Structural Chemistry and Biology; Oxford University Press: NewYork, 2001. (2) Computational Molecular Biology; Leszczynski, J., Ed.; Elsevier Science: Amsterdam, 1999. (3) Supramolecular Assembly Via Hydrogen Bonds I; Mingos, D. M. P., Ed.; Springer-Verlag: Berlin, 2004. (4) Supramolecular Assembly Via Hydrogen Bonds II; Mingos, D. M. P., Ed.; Springer-Verlag: Berlin, 2004. (5) Sauer, J.; Doebler, J. ChemPhysChem. 2005, 6, 1706-1710. (6) Kuo, I.-F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I.; VandeVondele, J.; Sprik, M.; Hutter, J.; Chen, B.; Klein, M. L.; Mohamed, F.; Krack, M.; Parrinello, M. J. Phys. Chem. B 2004, 108, 12990-12998. (7) Silvestrelli, P. L.; Bernasconi, M.; Parrinello, M. Chem. Phys. Lett. 1997, 277, 478-482. (8) Geissler, P. L.; Dellago, C.; Chandler, D.; Hutter, J.; Parrinello, M. Chem. Phys. Lett. 2000, 321, 225-230. (9) Vuilleumier, R.; Sprik, M.; Alavi, A. J. Mol. Struct. THEOCHEM 2000, 506, 343-353. (10) Hydrogen-Transfer Reactions; Hynes, J. T., Klinman, J. P., Limbach, H.-H., Schowen, R. L., Eds.; Wiley-VCH: New York, 2006. (11) Lee, Y. S.; Hodosˇcˇek, M.; Brooks, B. R.; Kador, P. F. Biophys. Chem. 1998, 70, 203-216. (12) Wong, S.; Bernacki, K.; Jacobson, M. P. J. Phys. Chem. B 2005, 109, 5249-5258. (13) Matanovic´, I.; Dosˇlic´, N. J. Phys. Chem. A 2005, 109, 4185-4194. (14) Bako´, I.; Hutter, J.; Pa´linka´s, G. J. Phys. Chem. A 2006, 110, 21882194. (15) Olsson, M. H. M.; Mavri, J.; Warshel, A. Philos. Trans. R. Soc., B 2006, 361, 1417-1432. (16) Pantano, S.; Alber, F.; Carloni, P. J. Mol. Struct. THEOCHEM 2000, 530, 177-181. (17) Frigyes, D.; Alber, F.; Pongor, S.; Carloni, P. J. Mol. Struct. THEOCHEM 2001, 574, 39-45. (18) Voth, G. A. Acc. Chem. Res. 2006, 39, 143-150. (19) Bała, P.; Grochowski, P.; Lesyng, B.; McCammon, J. A. J. Phys. Chem. 1996, 100, 2535-2545. (20) Rozas, I.; Alkorta, I.; Elguero, J. J. Phys. Chem. A 1997, 101, 42364244. (21) Shubina, E. S.; Belkova, N. V.; Epstein, L. M. J. Organomet. Chem. 1997, 536-537, 17-29. (22) Chamberlain, A. K.; Bowie, J. U. J. Mol. Biol. 2002, 322, 497503. (23) McMorrow, D.; Kasha, M. J. Phys. Chem. 1984, 88, 2235-2243.
Jezierska et al. (24) Mordzin´ski, A.; Grabowska, A.; Ku¨hnle, W.; Kro´wczyn´ski, A. Chem. Phys. Lett. 1983, 101, 291-296. (25) Mordzin´ski, A.; Grabowska, A. Chem. Phys. Lett. 1982, 90, 122127. (26) Grabowska, A.; Kownacki, K.; Karpiuk, J.; Dobrin, S.; Kaczmarek, L. Chem. Phys. Lett. 1997, 267, 132-140. (27) Sekikawa, T.; Kobayashi, T.; Inabe, T. J. Phys. Chem. A 1997, 101, 644-649. (28) Filarowski, A.; Szemik-Hojniak, A.; Głowiak, T.; Koll, A. J. Mol. Struct. 1997, 404, 67-74. (29) Koll, A.; Parasuk, V.; Parasuk, W.; Karpfen, A.; Wolschann, P. J. Mol. Struct. 2004, 690, 165-174. (30) Filarowski, A.; Koll, A. Vib. Spectrosc. 1996, 12, 15-24. (31) Koll, A.; Wolschann, P. Monatsh. Chem. 1999, 130, 9831001. (32) Koll, A.; Melikova, S. M.; Karpfen, A.; Wolschann, P. J. Mol. Struct. 2001, 559, 127-145. (33) Handgraaf, J.-W.; van Erp, T. S.; Meijer, E. J. Chem. Phys. Lett. 2003, 367, 617-624. (34) Morrone, J. A.; Haslinger, K. E.; Tuckerman, M. E. J. Phys. Chem. B 2006, 110, 3712-3720. (35) Sillanpa¨a¨, A. J.; Simon, C.; Klein, M. L.; Laasonen, K. J. Phys. Chem. B 2002, 106, 11315-11322. (36) Marx, D.; Tuckerman, M. E.; Martyna, G. J. Comput. Phys. Commun. 1999, 118, 166-184. (37) Benoit, M.; Marx, D. ChemPhysChem 2005, 6, 1738-1741. (38) Ramirez, R.; Lopez-Ciudad, T.; Kumar, P.; Marx, D. J. Chem. Phys. 2004, 121, 3973-3983. (39) Pejov, L.; Spångberg, D.; Hermansson, K. J. Phys. Chem. A 2005, 109, 5144-5152. (40) Bourˇ, P. Chem. Phys. Lett. 2002, 365, 82-88. (41) Wolf, K; Mikenda, W.; Nusterer, E.; Schwarz, K. J. Mol. Struct. 1998, 448, 201-207. (42) Berendsen, H. J. C.; Mavri, J. J. Phys. Chem. 1993, 97, 1346413468. (43) Tully, J. C.; Preston, R. K. J. Chem. Phys. 1971, 55, 562-572. (44) Cao, J.; Voth, G. A. J. Chem. Phys. 1993, 99, 10070-10073. (45) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 100, 5106-5117. (46) Martyna, G. J. J. Chem. Phys. 1996, 104, 2018-2027. (47) Cao, J.; Martyna, G. J. J. Chem. Phys. 1996, 104, 2028-2035. (48) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471-2474. (49) Hohenberg, P.; Kohn, W. Phys. ReV. 1964, 136, B864-B871. (50) CPMD Copyright IBM Corp. 1990-2004, Copyright MPI fuer Festkoerperforschung Stuttgart 1997-2001. (51) Borsˇtnik, U.; Hodosˇcˇek, M.; Janezˇicˇ, D. J. Chem. Inf. Comput. Sci. 2004, 44, 359-364. (52) Andreoni, W.; Curioni, A. Parallel Computing 2000, 26, 819842. (53) Hutter, J.; Curioni, A. ChemPhysChem 2005, 6, 1788-1793. (54) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865-3868. (55) Ireta, J.; Neugebauer, J.; Scheffler, M. J. Phys. Chem. A 2004, 108, 5692-5698. (56) Vanderbilt, D. Phys. ReV. B 1990, 41, 7892-7895. (57) Schlegel, H. B. Theor. Chim. Acta 1984, 66, 333-340. (58) Panek, J.; Stare, J.; Hadzˇi, D. J. Phys. Chem. A 2004, 108, 74177423. (59) Stare, J.; Mavri, J. Comput. Phys. Commun. 2002, 143, 222-240. (60) Isotope Effects in Chemistry and Biology; Kohen, A., Limbach, H.-H., Eds.; CRC Press: New York, 2005. (61) Pavon, J. A.; Fitzpatrick, P. F. J. Am. Chem. Soc. 2005, 127, 16414-16415. (62) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33-38. (63) Steiner, T.; Majerz, I.; Wilson, C. C. Angew. Chem., Int. Ed. 2001, 40, 2651-2654. (64) Egorov, S. A.; Skinner, J. L. J. Chem. Phys. 1996, 105, 70477058. (65) Egorov, S. A.; Rabani, E.; Berne, B. J. J. Chem. Phys. 1997, 108, 1407-1422. (66) Egorov, S. A.; Rabani, E.; Berne, B. J. J. Chem. Phys. 1998, 110, 5238-5248. (67) Lawrence, C. P.; Nakayama, A.; Makri, N.; Skinner, J. L. J. Chem. Phys. 2004, 120, 6621-6624. (68) Poulsen, J. A.; Nyman, G.; Rossky, P. J. PNAS 2005, 102, 67096714. (69) Rey, R.; Hynes, J. T. J. Chem. Phys. 1996, 104, 2356-2368.
