J. Phys. Chem. 1995,99, 12711-12717
12711
Calculation of the Proton Transfer Rate Using Density Matrix Evolution and Molecular Dynamics Simulations: Inclusion of the Proton Excited States Janez Mavrit and Herman J. C. Berendsen” BIOSON Research Institute, Department of Biophysical Chemistry, the University of Groningen, Nijenborgh 4, 9747 AG Groningen, The Netherlands Received: March 7 , 1995; In Final Form: June 20, 1995@
The methodology for treatment of proton transfer processes by density matrix evolution (DME) with inclusion of many excited states is presented. The DME method (Berendsen, H. J. C.; Mavri, J. J . Phys. Chem. 1993, 97, 13464) that simulates the dynamics of quantum systems embedded in a classical environment is extended to nonorthogonal basis sets. The method is applied to calculations of the rate of proton transfer in the doublewell intramolecular hydrogen bond in the hydrogen malonate ion in aqueous solution. Five optimized basis functions were applied, giving rise to five vibrational levels of the proton. A nearly 30-fold increase in the proton transfer rate relative to the earlier treatment using only two basis functions was observed, indicating the importance of flexible, optimized basis sets in the proton transfer calculations. Possible applications of the method to biologically interesting processes such as proton transfer in enzymatic reactions are discussed.
1. Introduction Proton transfer reactions are an important class of processes in chemistry and biochemistry. Since proton transfer usually takes place through a hydrogen bond, treatment of hydrogen bonding is of prime importance.’ Due to the relatively small mass of the proton, its thermal length at room temperature is 0.3 A, implying the breakdown of the Bom-Oppenheimer approximation. Therefore, one of the key assumptions for the validity of classical molecular dynamics (MD) simulations is violated.2 A complete quantum dynamical treatment of large systems like a protein in aqueous solution is not possible with present computer power, nor will it be possible in the near future. A possible and sufficiently accurate solution is the blending of classical and quantum degrees of freedom such that those degrees of freedom which have an essential quantum mechanical character are described by the laws of quantum mechanics, while the rest of the system is described classically. A hydrogenbonded system in solution would in this sense consist of a quantum subsystem comprising a proton, a proton donor, and a proton acceptor and a classical subsystem comprising the rest of the solute atoms and solvent molecules. In the last few years significant progress in the MD methodology for the quantum mechanical treatment of electronic degrees of f r e e d ~ m ~has -~ been recorded. Most methods are restricted to the electronic ground state, which is, in general, a very good approximation for electronic degrees of freedom when a large gap exists between the ground and first excited state. A key for understanding molecular activation-deactivation phenomena and therewith associated reaction dynamics in general is hidden in the dynamics of nuclear degrees of freedom.’O A proper quantum dynamical description of nuclear degrees of freedom, such as those of protons, requires the inclusion of excited states. The quantum dynamical treatment of the quantum dynamics of transitions between them is mandatory since the typical energy gaps between ground and excited states are comparable to ~ B T . Thus a nonadiabatic approach is required.
* To whom correspondence should be addressed.
’ Permanent
address: National Institute of Chemistry, P.O.B. 30, Hajdrihova 19, 61 115, Ljubljana, Slovenia. Abstract published in Advance ACS Abstracts, August 1, 1995. @
Several methods for the treatment of nonadiabatic systems have been developed. The wave packet method is able to describe quantum dynamics.”.’2 The path integral (PI) approach is an attractive formulation of quantum mechanic^;'^ a traditional view of the PI method is based on the isomorphism of the quantum particle with a necklace of classical particles connected by harmonic springs that provides correct ensemble averages but does not provide dynamical quantities. Recently, impressive progress was reported in the introduction of time dependence into the PI calculation^.'^-'^ Surface hopping is an efficient semiclassical method for treatment of the nonadiabatic transitions. For a recent formulation see ref 19. The literature on nonadiabatic transitions is too broad to be entirely included in this article; for a recent extensive review of the methods for treatment of nonadiabatic transitions, see Coker20and references therein. Recently, we developed the hybrid method of density matrix evolution (DME)*’ and applied it to the treatment of nonadiabatic transitions, ranging from a perturbed quantum harmonic oscillator to proton tunneling in a bistable hydrogen bond. The DME method consists of coupling the Liouville-von Neumann equation (evolution of the density matrix for the quantum subsystem) with the classical equations of motion (for classical degrees of freedom) in such a way that the basic conservation laws (e.g. of total energy) are obeyed. Together with the mathematical formulation of the DME simulation method a numerical test was performed.21 The system studied consisted of a two-level quantum harmonic oscillator (QHO) and a single, collinear, noble gas atom. The QHO was fixed at the coordinate origin and was initially entirely in its ground state. The classical atom was placed at some distance from the QHO with an initial velocity directed toward the QHO. An elastic collision was observed, and the total energy was conserved within the accuracy of the numerical integration of the equations of motion. We also simulated a three-level quantum harmonic oscillator immersed in a classical noble gas bath.22 The simulated populations of the fibrational levels of the QHO immersed in a noble gas bath appeared to be close to the Boltzmann distributions To our knowledge this is the first simulation of the equilibrium distribution over quantum level populations using
0022-365419512099-12711$09.0010 0 1995 American Chemical Society
Mavri and Berendsen
12712 J. Phys. Chem., Vol. 99, No. 34, 1995
this flexible basis set both nonadiabatic and adiabatic proton nonadiabatic quantum molecular dynamics. From the averaged transfer channels are taken into account simultaneously. We Boltzmann factor the free energy of transfer between ground restrict our treatment to initial rates with the proton still in the and excited vibrational levels was calculated and was found to reactant well; this allows us to neglect the feedback of the be larger than the corresponding in-vacuo value, as a result of quantum system to the classical system. There are two reasons the perturbation of the oscillator by the collisions. We applied that justify this restriction: first, the impossibility to simulate the DME method to study the inelastic collisions of a classical long enough to observe a statistically significant number of particle with a five-level quantum harmonic o ~ c i l l a t o r .The ~~ proton transitions, and second, we wish to make a valid results were compared with the exact quantum calculation^^^ comparison with the previous two-level calculation, based on and with the results of classical mechanics. The DME results the same MD run. are between the results of the full quantum treatment and classical mechanics, but closer to the latter. We demonstrated The organization of this article is as follows: Section 2 briefly the time reversibility of the DME-MD equations by reversing describes the DME method with an emphasis on introduction the velocity of the colliding particle and changing all complex of a nonorthogonal basis set. Requirements for the orthogonalization are due to the absence of dynamical coupling between quantities to their complex conjugates on the right side of the the quantum and classical subsystems. Section 3 describes the Liouville-von Neumann equation. It is essential for the time simulation of the hydrogen malonate ion in aqueous solution. reversibility of the process that the quantum phase information In section 4 results are given of the proton transfer rate stored in the off-diagonal elements of the density matrix be calculation that include five states produced by an optimized retained. Dephasing of the quantum subsystem that leads to basis set. The approach is discussed in section 5 , and possible irreversible collisions can be related to an increase in entropy. applications to biological systems are considered. Very recently,26we improved the numerical integration of the equations of motion for a collision between a multilevel quantum 2. Density Matrix Evolution harmonic oscillator and a classical particle. The Liouvillevon Neumann equation describing the quantum subsystem was For a detailed description of the DME method the reader is rewritten in the interaction repre~entation,~~ eliminating thus referred to ref 21, 22, and 23. Since the DME method is the frequencies of the unperturbed oscillator, allowing for a relatively new, we briefly present the formalism, but now larger time step in the numerical integration. including nonorthogonal basis functions since these occur in Calculations of the proton transfer rate constants in hydrogenthe description of the majority of realistic systems. bonded systems embedded in a fluctuating polar environment It is computationally convenient to calculate the matrix are applications of the methods of quantum molecular dynamics. elements over the original (nonorthogonal) basis set and use In the past separate calculations for adiabati~*~-~O and nonthe proper transformations to convert them to an orthogonal adiabatic proton transfer3' have been performed. One of the basis. Computationally this is more convenient than propagating first applications of the DME method was the calculation of he DME equations of motion in a nonorthogonal basis. the proton transfer rate in a double-well hydrogen bond in a ollowing the notation of ref 41, we denote the quantities over hydrogen malonate ion in aqueous solution.32 We performed the orthogonal basis set by a prime, while the quantities over the calculation with only two relatively narrow simple Gaussians the original, nonorthogonal basis set are denoted without the as basis functions: one centered in the reactant well and the prime. The DME procedure is represented here in an algorithother centered in the product well. This basis set was not mic way in order for it to be easily coded. optimized and was far from complete. Furthermore, coupling In the DME method the quantum subsystem is described in between classical and quantum subsystems was not realized a Hilbert space of basis functions. properly; that is, feedback from changes in the quantum system Step 1. The wave function of the quantum subsystem is to the classical subsystem was not incorporated. The role of expanded on a properly chosen orthonormal set of basis classical MD simulation was to provide time-dependent proton functions 4'. potentials for the proton restrained in the reactant well. The small basis set used to calculate the rate of proton transfer M was restricted to the nonadiabatic channel. Two basis functions giving rise to a two-level system were considered: one was localized in the reactant well and the other in the product well.32 For tunneling between states of almost equal energy, the transfer where 6 denotes the coordinate(s) of the quantum subsystem. is nonadiabatic. The poor basis set used does not allow a proper The time propagation of the quantum subsystem is described treatment of both adiabatic and nonadiabatic transfer. Therefore, by the time dependence of the coefficients in the linear extending the basis set in the DME method would present a combination. The M x M density matrix e' is defined as unifying approach (taking into account both adiabatic and = c:ch* where M is the number of basis functions applied in nonadiabatic channels including channels via excited states) for the linear combination; e' is complex but hermitian. Diagonal the calculation of proton transfer rates. A unifying approach elements represent the populations of the levels, while offfor the proton transfer rate calculations was offered by using diagonal elements contain phase information. For the time path integral techniques and transition state t h e ~ r y . ~The ~ - ~ ~ evolution we follow closely a formalism that is routinely used description of the quantum dynamics of the proton embedded in the dynamics of spin systems.27 This approach was recently in a classical environment using the method of wave packet followed for treatments of optical transitions rates49 and propagation represents a unifying a p p r o a ~ h . ~Very ~ . ~recently ~ vibrational relaxation50in the Redfield approximation:* which Hammes-Schiffer and T ~ l l y reported ~ ~ . ~ calculation ~ of the implies first-order rate constants for the time dependence of proton transfer rates using a surface-hopping method which is diagonal elements. In our approach such an approximation is also a unifying approach. not applied. In this article we recalculated the rate of intramolecular proton Step 2. Calculate Hamiltonian matrix elements over the transfer in aqueous solution of hydrogen malonate by using five original, nonorthogonal basis set in the absence of external optimized Gaussians, giving rise to five quantum levels. With perturbation:
L
Calculation of the Proton Transfer Rate
J. Phys. Chem., Vol. 99, No. 34, 1995 12713
where VO refers to potential energy of the unperturbed static system. The unperturbed system can be arbitrarily chosen, e.g. as the average potential experienced by the proton. Step 3. Calculate overlap matrix elements:
Step 4. Transform Ho to the orthogonal set of basis functions 4‘:
u Figure 1. Structure of the hydrogen malonate ion.
where X is obtained from the overlap matrix using Lowdin or canonical orthog~nalization.~’In our work the canonical procedure was used. The canonical orthogonalization is described in the Appendix. The matrix Ho’ is almost diagonal: the off-diagonal elements are nonzero because of the incompleteness of the basis set. Step 5. Diagonalize Ho’ to obtain the coefficient matrix C‘ over the orthonormal basis set. Step 6. Transform the coefficient matrix to the original, nonorthogonal basis set:4’
c = XC‘
(5)
A
Step 7. TOcalculate matrix elements AY of some operator (external perturbation, generalized Hellmann-Feynman force, population of the product state, ...) over the orthogonal basis, it is straightforward to show that M
M
I= I
k= I
I=lk=l
(7) and convert them to the orthogonal basis using the equation given in step 7. Matrix y, results. In our case V refers to the deviation from the averaged proton potential. Step 9. Add the matrices Ho’ and V’:
+ V’
(8)
The energy of the quantum subsystem can be calculated:
E, = tr(p’H’)
Note that in the treatment presented in this article quantum dynamics calculations were performed a posteriori, and the quantum forces calculation was not required. Step 11. Integrate simultaneously the equations of motion for classical and quantum subsystems. Note that in our case simultaneous integration was not performed, since we applied classical MD simulation only to provide a time course of proton potentials that were used in the DME equations. The density matrix evolves with time2’ according to the Liouville-von Neumann equation:
d
where the basis set is orthogonalized taking into account the shape of the unperturbed potential of the quantum subsystem in step 4. This procedure has the advantage that the ground state coincides with the reactant state so that simple initial conditions for the density matrix can be applied. Note that application of an orthogonalization procedure that takes into account knowledge about the average potential experienced by the proton is very convenient in the case where there is no feedback from dynamics of the quantum system to the classical system and the rate constant is calculated perturbatively (vide infra). Canonical orthogonalization,which does not take into account knowledge about the potential, would produce a symmetric ground state, and any pure state would predict an expectation value of the proton corresponding to the center of the hydrogen bond. Such an initial condition is unphysical and cannot correspond to the reactant state. Step 8. Calculate the external perturbation matrix elements,
H’ = Ho‘
quantum subsystem as well as the interaction energy between the quantum and classical subsystem. Step 10. If required, calculate the force acting on a classical particle from the quantum subsystem. The force on a kth classical particle reads
(9)
Note that this term includes kinetic and potential energy of the
i
$’ = @’H’ - H’p’) Note that the sign in this equation is a matter of convention and is not relevant for the observables since it only affects the quantum phase. If one applies proper DME coupling, the Newtonian equations of motion are integrated for classical particles, applying the force acting from the quantum subsystem on the classical subsystem in addition to classical particle to particle forces. Note that the steps 1-6 are performed only at the initialization, while the remaining steps are performed repeatedly when integrating the equations of motion. With e’ the expectation values can be calculated at every time instant also for any other operators (displacement, population of the product well, ...) using the transformation in step 7. Note that the transformation in step 7 is computationally more complex than canonical or Lowdin orthogonalization,but it takes into account the knowledge about the unperturbed potential. By the orthogonalization procedure defined in step 7, we defined a basis set on which the I$ is diagonal.
3. Computation of Proton Potentials The internal hydrogen bond in the hydrogen malonate ion is one of the strongest hydrogen bonds k n ~ w n . ~Its~structure .~~ is shown in Figure 1. Ab-initio calculations4 predict on all but the lowest levels of theory the proton potential to be a symmetric double well, but they are not conclusive about the exact height of the barrier. .We modeled the in-vacuo proton potential by 12-1 terms. A double-well potential with a classical barrier height of 2.52 kcal mol-’ and an equilibrium 0-0 distance of 2.54 8, resulted. All other nonbonding and bonding
Mavri and Berendsen
12714 J. Phys. Chem., Vol. 99, No. 34, 1995 parameters were taken from the GROMOS force field.45 The methylene group was treated as a united atom. Water was described in terms of the SPC (simple point charge) having charges of -0.82 and 0.41 on oxygen and hydrogen atoms, respectively, and a 12-6 term located on the oxygen atom. A truncated octahedron periodic boundary condition with a twincutoff method, with cutoff radii of 8 and 10 A, respectively, was used. Proton potentials were calculated by adding the electrostatic contribution of water to the in-vacuo value. For all details of the calculation of time-dependent proton potentials in ion hydrogen malonate the reader is referred to ref 32. Charges in hydrogen malonate were calculated with the procedure that preserves the dipole momene7 for hydrogen malonate in CI and CzVsymmetry on the HFDZP level of theory. Since the charges for the products and for the reactants are the same as a result of symmetry, the charges were available at three proton positions. A parabolic interpolation was used for other positions. In this work we address the calculation of the proton transfer rate in the very initial stage of the transfer where the population of the product state is still negligible. The proton was reduced to one dimension. Proton transfer via solvent molecules was not considered in this study. For the very initial stage of the proton transfer process one may assume that the quantum subsystem can develop autonomously without having a serious impact on the dynamics of the classical subsystem. By considering the proton fixed in the reactant well, we explored the solvent coordinate that is lssential in the proton transfer processes. Classical MD simulation was performed for a proton represented by a point charge that was held fixed in the reactant well. During the 100 ps of classical MD of hydrogen malonate ion in aqueous solution the quantum subsystem was therefore not evolving. Moreover, the proton, considered to be in its ground state, was represented by a point charge.32 The role of classical MD was therefore to provide a swarm of proton potentials. We stored the trajectory every 2 fs, and therefore 50 000 proton potentials resulted. The averaged proton potential was asymmetric due to solvent-induced asymmetry. The timeaveraged proton potential was chosen for the calculation of HO. This choice of the unperturbed potential is the most natural and is very convenient for the perturbational calculation of the proton transfer rate (vide infra).
1
9 (A) Figure 2. Time-averaged proton potential (thick solid line) of hydrogen
malonate (in kcaymol). Applied basis functions (thin lines) were located at 0.9, 1.0, 1.25, 1.5, and 1.6 A with the exponents of 45, 29, 14.837, 29, and 45 A-2. Accurate eigenvalues obtained with the Numerov method for this case are 3.894, 8.083, 11.102, and 15.2513 kcaymol, while with variational solving using the named basis set we obtained the eigenvalues of 3.892, 8.059, 11.81, and 19.40 kcal/mol.
reactant well is stabilized. Calculation of the rate constant would not be possible. Note that in the case of proper coupling between the classical and quantum subsystem the orthogonalization is arbitrary. We assumed that initially only the ground state corresponding to reactants was populated. We started each quantum dynamics run by a proton placed in its ground state of the averaged proton potential. Time-dependent deviations from the averaged proton potential were treated as extemal perturbations. For each proton potential matrix elements were calculated. Matrix elements were calculated numerically using an accurate 30-point GaussLegendre procedure. Values for the time between 2 fs were interpolated using cubic spline interpolation. DME equations of motions were integrated with a time step of 0.1 fs using a Runge-Kutta fourth-order integrator. The quantity
4. Quantum Dynamics by Density Matrix Evolution We performed quantum dynamical calculations a posteriori with five displaced simple Gaussian basis functions, the widths of which were optimized to match the eigenvalues of the averaged proton potential calculated by the Numerov method. Basis functions are shown in Figure 2. Five basis functions gave rise to five levels. Since the DME calculation was done a posteriori, calculation of the quantum force in step 10 was not necessary. On the other hand owing to the fact that in this work we address the problem of quantum subsystem without feedback to the classical system, the orthogonalizationprocedure as described in step 7 is very convenient, since it describes the reactant state as the pure ground state. We calculated the rate of proton transfer by performing quantum dynamics calculations for a proton in a fluctuating potential, averaged over many time origins, and comparing the increase of the product population by a phenomenological equation for the first-order reaction rate. The proton was initially confined to its ground state. If one would choose for the unperturbed proton potential a symmetric double-well potential, the population of the product well would be half and would even decrease with time. The explanation is simple, since owing to solvent-induced asymmetry, the
where ROis half of the proton donor-proton acceptor distance, corresponds to the population of the product state. Formally, 0 is the expectation value of the Heaviside operator H(( Ro). We took 2000 time origins, and the DME equations were integrated for 1.5 ps. The averaged time course of the product state populations is shown in Figure 3. From the coarse-grained slope a rate constant for the proton transfer of (4.0 f 0.9) x lo8 s-I is determined. The fact that there is a nonvanishing averaged initial population of the product state results from the quantum nature of the proton and the choice of the ground state as the initial reactant state. We used the expression for the phenomenological first-order, reversible rate process in the initial stage, in which the population of the product state is still negligible. The calculated rate constant is about 30 times larger than the value obtained using only two Gaussians with nonoptimized
5. Discussion We calculated the rate of proton transfer in a hydrogen malonate ion in aqueous solution using the method of density
Calculation of the Proton Transfer Rate
0013
00
J. Phys. Chem., Vol. 99, No. 34, I995 12715
I
d
io00 0 t (IS)
Figure 3. Time course of population of product state calculated by
the density matrix evolution method obtained by averaging solutions over 2000 time origins. matrix evolution that is able to combine the quantum and the classical degrees of freedom in molecular dynamics simulations. The quantum rate constant was determined perturbatively ; that is, the proton was fixed in the reactant well during the classical MD simulation. The dynamics of the quantum subsystem was calculated a posteriori, i.e. without feedback of the quantum subsystem to the classical subsystem. The rate constant was determined from the coarse-grained slope of the population of the product well in the very initial stage of the reaction, where the population of the product well is still negligible. The observed rate constant implies that the proton transfer occurs on a nanosecond time scale. To calculate the rate constant from the properly coupled quantum classical MD, a DME-MD simulation over several hundred nanoseconds would be required. Since this is at present not feasible, we approximate the rate constant by the initial population decay. The calculated rate constant corresponds to the quantum transition state (disregarding barrier recrossing events), rather than to the real quantum rate constant. Methodology for calculation of the latter has been developed5' and successfully applied in the calculation of a d i a b a t i ~ ~and ~ . ~adiabatic~ nonadiabatic proton transfer rate constants39in model systems. It is worthy to stress that the reactive flux correlation function approach requires MD simulations with several hundreds (or at least tens) of recrossings of the (quantum) barrier. The quantum barrier in the model system^*^-^^ was relatively low and thus allowed several recrossings in the nanosecond time scale. Inspection of the rate constant for the proton transfer in the ion hydrogen malonate in aqueous solution reveals that a single transition of the proton would occur in about 2 ns. Having in mind the necessity for the simulation of several recrossings, it becomes evident that it is not yet possible to calculate directly and rigorously the rate constant for proton transfer in the studied system. The transmission coefficient that is on the order of unity ensures that the estimated rate constant is useful for most of the molecular modeling purposes. Since our system is more similar to realistic biological systems (hydrogen-bonded systems in aqueous or polar macromolecular environment) than the model systems (A-H-A in dipolar ~ o l v e n t ) , * a~ similar - ~ ~ unfavorable time behavior is to be expected there. The presently developed methodology can be easily applied to the systems of (bio)chemical importance by providing the time course of the proton potentials for the proton fixed in the reactant well during the classical MD, exploring thus the solvent coordinate, rather than the protonproton donor distance. Ando and H y n e ~very ~ ~ recently
performed a related study of the ionization of hydrochloric acid in aqueous solution. They stressed the sensitivity of proton transfer calculations to the level of theory and the importance of inclusion of the first layer of water molecules in the abinitio quantum calculations. Our study showed the importance of a flexible basis set in the calculation of proton transfer rate constants. The previous study3*using only two basis functions underestimated the rate constant severly (i) by taking only one channel into account, and (ii) even that value was underestimated since the Gaussians were too narrow to produce correct coupling matrix elements. It is important to note that our method simultaneously includes adiabatic and nonadiabatic proton transfer because the flexible basis set allows a proper description of both the ground state and the relevant excited states. The recent study of HammesSchiffer and T ~ l l showed y ~ ~ that the assumption of adiabatic proton transfer overestimates the overall rate constant. The reason lies in the assumption of adiabatic s t ~ d i e s *that ~ . ~the ~ time-dependent proton wave function adapts instantaneously to the time-dependent proton potential, which is not physically reasonable. In the method used here for proton transfer rate calculation the classical system did not feel the changes in the quantum subsystem. Therefore, this version of the DME method falls into the category of perturbational treatments. The main deficiency of the DME method, i.e. that the classical subsystem feels the averaged quantum systems,24is of no consequence in the present treatment of proton transfer. We would like to emphasize that in the case of proper quantum classical coupling this might present a problem. A known example is the failure of the time-dependent self-consistent-field method53 when applied to a double-well quantum system strongly coupled to a harmonic bath.54 Despite the outstanding progress in the experimental proton potential d e t e r m i n a t i ~ n , ~ ~ab-initio -~* calculations remain a valuable source of information about the proton potentials in hydrogen-bonded systems. Since the ab-initio and the semiempirical calculations are limited by the size of systems, there is always a question of reliability of data calculated from truncated systems (in most cases in-vacuo) when these are used in liquid simulations. Combining a solvent reaction field method with ab-initio calculations is one solution. For a recent review see ref 59. A more elegant (and more reliable) method to obtain reliable proton potentials consists of combining an ab-initio molecular orbital method with the classical environmenLm Adopting such a method, it is possible to provide proton potentials for processes in solutions or in such complex environments as enzymes; following the procedure described above, it is possible to overcome the Bom-Oppenheimer approximation for the proton. Methods for hiearchical treatment of complex systems are developed and ready to be ~ s e d . ~ - ~ HodoSEek and Brooks7 very recently coupled an ab-initio MO and a molecular mechanics method in order to perform quantum classical MD that can be used in the calculation of timedependent proton potentials. If one opts for traditional parametrization of the molecular mechanics force field based on quantum methods, ab-initio methods are a valuable source of information for nonbonding parameters (for example, atomic charges) that may vary significantly with transfering the proton.32.61
.62
Since the proton transfer processes are electronically adiabatic, the electronic degrees of freedom follow the proton dynamics. Therefore, application of an electronically polarizable environment, rather than usage of averaged polarizable (effective) force fields, as is now quite often the case, would be desirable in the
12716 J. Phys. Chem., Vol. 99, No. 34, 1995 simulations of such processes. Development of a polarizable force field based on the polarizable shell model is in progress in our The vibration associated with the proton donor-proton acceptor distance was in our study treated classically. It was demonstrated elsewhere3I that quantum treatment of this mode, closely related to the coupling matrix element, results in a higher rate of proton tunneling. Treatment of this mode quantum dynamically using DME would be possible but would result in a many quantum body problem, giving rise to a somewhat more computationally demanding procedure for solving the DME equations. Future work will be directed mainly toward applications in more complex systems such as treatment of reaction dynamics for proton transfer in aqueous solutions and enzymatic environments using the approach applied in this study. An efficient approach to examine high quantum banier systems would include biased sampling, which was already used in our previous treatment.32 Inclusion of constraint or restraint for the given value of the reaction coordinate to the DME equations is possible and results in a potential of mean force for the proton. The procedure is somewhat related to the centroid dynamics in the path integral calculations.64 The transmission coefficient could be determined by starting activated trajectories at the top of the quantum barrier. We would like to stress here again the weak point of the DME method, i.e. that the classical subsystem feels the averaged quantum subsystem; this might be deficient for calculations of proton transfer based on the reactive flux correlation function. One of the future activities will be calculation of the rate constants for proton transfer processes using DME based on the reactive flux correlation functions, allowing thus several recrossings of the barrier. For such calculations it is essential to examine for which realistic systems the DME approximation breaks down and for which systems and what desired properties it becomes necessary to switch to advanced and computationally intensive methods for blending quantum and classical degrees of freedom, such as surface-hopping methods.I9 DME might be a promising method in the computational support of vibrational spectroscopy. Due to the fact that it is relatively easy to include an external optical field into the DME equations, simulations of pulse-echo experiments are quite possible.
6. Conclusions The density matrix evolution (DME) method is applied to calculate the rate of proton transfer in the double-well intramolecular hydrogen bond in a hydrogen malonate ion in aqueous solution. The proton was held fixed in the reactant well during the classical molecular dynamics (MD) simulation. Using proton potentials from the classical MD, the quantum dynamics of the proton was examined and the rate of proton transfer was calculated from the time-averaged coarse-grained increase of the product state population. With a flexible basis set the adiabatic and nonadiabatic channels for proton transfer were included. A nearly 30-fold increase in the proton transfer rate relative to the early treatment using only two basis functions was observed, indicating the importance of a flexible, optimized basis set in proton transfer calculations. We hope to report applications of the method on biologically interesting processes, such as proton transfer in enzymatic reactions, in the near future.
Acknowledgment. This work was supported by COST grant D-3. One of us (J.M.) is grateful for a long-term fellowship from the Human Frontier Science Program and to the University
Mavri and Berendsen of Groningen for the hospitality during his stay in The Netherlands. We are grateful Prof. Dugan Hadii (National Institute of Chemistry, Slovenia) and Marc F. Lensink (University of Groningen) for critical reading of the manuscript. We are grateful to Prof. James T. Hynes (University of Colorado at Boulder) for many stimulating discussions. We thank Dr. Daniel Borgis (Universite Pierre et Marie Curie) for many stimulating discussions and for providing us with the subroutine for solving the one-dimensional Schrodinger equation using the Numerov method.
Appendix Canonical orthogonalization of the basis set'" can be described as follows. Canonical orthogonalization uses the unitary matrix U,the columns of which are eigenvectors of HO. The overlap matrix S is diagonalized, and the matrix X is calculated from its diagonal form.
where indices i and j extend over M basis functions. The matrix X is used in step 4. An important advantage of the canonical orthogonalization over the Lowdin procedure is4' that in the case of linear dependence of the basis set, i.e. when eigenvalues of the overlap matrix approach 0, the matrix X can still be applied by eliminating the critical basis functions.
References and Notes Hadti, D. J. Mol. Struct. 1988, 177, 1. Chandler, D. J. Star Phys. 1986, 42, 49. Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471. Selloni, A.; Carnevali, P.; Car, R.; Parrinello, M. Phys. Rev. Lett. 1987, 59, 823. ( 5 ) Laasonen, K.; Sprik, M.; Paninello, M. J . Chem. Phys. 1993, 99, 9080. (6) Curioni, A.; Andreoni, W.; Hutter, J.; Schiffer, H.; Parrinello, M. J . A m . Chem. SOC. 1994, 116, 11251. (7) HodoSEek, M.; Brooks, B. R. To be published. (8) Field, M. J.; Bash, P. A.; Karplus, M. J . Comput. Chem. 1990, 11, 700. Thery, V.; Rinaldi, D.; Rivail, J. L.; Maigret, B.; Ferenczy, G. C. lbid. 1994, 15, 269. (9) Stanton, R. V.; Hartsough, D. S.; Merz, K. M., Jr. J . Comput. Chem. 1995, I , 113. (10) Hanggi, P.; Talkner, P.; Borkovec, M. Rev. Mod. Phys. 1990, 62, 251. (11) Heller, E. J. J . Chem. Phys. 1975, 62, 1544. (12) Truong, N. T.; Tanner, J. J.; Bala, P.; McCammon, J. A.; Kuori, D. J.; Lesyng, B.; Hoffman, D. K. J . Chem. Phys. 1992, 96, 2077. (13) Feynmann, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (14) Makri, N. Chem. Phys. Lett. 1992, 193, 435. (15) Makarov, D. E.; Makri, N. Phys. Rev. A 1993, 48, 3626. (16) Topaler, M.; Makri, N. Chem. Phys. Left. 1993, 210, 285. (17) Mak, C. H. Phys. Rev. Lett. 1992, 68, 899. (18) Mak, C. H.; Gehlen, J. N. Chem. Phys. Lett. 1993, 206, 130. (19) Tully, J. C. J . Chem. Phys. 1990, 93, 1061. (20) Coker, D. F. Computer Simulation Methods for Nonadiabatic Dynamics in Condensed Systems. In Computer Simulation in Chemical Physics; Allen, M. P.; Tildesley, D. S., Eds.; Kluwer Academic Publishers: Boston, 1993; p 315. (21) Berendsen, H. J. C.; Mavri, J. J . Phys. Chem. 1993, 97, 13464. (22) Mavri, J.; Berendsen, H. J. C. Phys. Rev. E 1994, 50, 198. (23) Mavri, J.; Berendsen, H. J. C. J . Mol. Struct. 1994, 322, 1. (24) Mavri, J.; Lensink, M.; Berendsen, H. J. C. Mol. Phys. 1994, 82, 1249. (25) Secrest, D.; Johnson, B. R. J . Chem. Phys. 1966, 45, 4556. (26) Lensink, M.; Mavri, J.; Berendsen, H. J. C. J . Comput. Chem., submitted. (27) Emst, R. R.; Bodenhausen, G.; Wokaun, A. Principles of Nuclear Magnetic Resonance in One and Two Dimensions; Clarendon Press: Oxford, 1987. (28) Borgis. D.; Tarjus, G.;Azzouz, H. J. Phys. Chem. 1992,96,3188. (29) Borgis, D.; Tarjus, G.; Azzouz, H. J . Chem. Phys. 1992.97, 1390. (30) Laria, D.; Ciccotti, G.; Ferrario, M.; Kapral, R. J . Chem. Phys. 1992, 97, 378. (31) Borgis, D.; Hynes, J. T. J . Chem. Phys. 1991, 94, 3619. (1) (2) (3) (4)
J. Phys. Chem., Vol. 99, No. 34, 1995 12717
Calculation of the Proton Transfer Rate (32) Mavri, J.; Berendsen, H. J. C.; van Gunsteren, W. F. J . Phys. Chem. 1993, 97, 13469.
(33) Lobaugh, J.; Voth, G. A. Chem. Phys. Lett. 1992, 198, 311. (34) Lobaugh, J.; Voth, G. A. J . Chem. Phys. 1994, 100, 3039. (35) Azzouz, H.; Borgis, D. J . Mol. Liq. 1995, 63, 89. (36) Pomes, R.; Roux, B. Chem. Phys. Lett. 1995, 234, 416. (37) Bala, P.; Lesyng, B.; McCammon, J. A. Chem. Phys. 1994, 180, 271. (38) Bala, P.; Lesyng, B.; Truong, T. N.; McCammon, J. A. Ab Initio Studies and Quantum-Classical Molecular Dynamics Simulations for Proton Transfer Processes in Model Systems and in Enzymes. In Molecular Aspects of Biotechnology: Computational Models and Theories; Bertran, J., Ed.; Kluwer Academic Publishers: Boston, 1992; p 299. (39) Hammes-Schiffer, S.; Tully, J. C. J . Chem. Phys. 1994, 101,4657. (40) Hammes-Schiffer, S.; Tully, J. C. J . Phys. Chem. 1995, 99, 5793. (41) Szabo, A,; Ostlund, N. S. Modern Quantum Chemistry; Macmillan: New York, 1982; p 144. (42) Djinovit, K.; Golit, L.; Hadii, D.; Orel, B. Croat. Chim. Acta 1988, 61, 405. (43) Kidrit, J.; Mavri, J.; Podobnik, M.; HadZi, D. J . Mol. Struct. 1990, 237, 265. (44) Mavri, J.; HodoEek, M.; Hadii, D. J . Mol. Struct. (THEOCHEM) 1990, 209, 421. (45) van Gunsteren, W. F.; Berendsen, H. J. C. GROMOS, Groningen Molecular Simulation Package; Biomos, b. v. Nijenborgh 4,9747 AG Groningen: The Netherlands, 1987. (46) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hemans, J. In lntermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; p 331. (47) Thole, B. T.; van Duijnen, P. Th. Theor. Chem. Acta 1980, 55, 307.
(48) Redfield, A. G. Adv. Magn. Reson. 1965, I, I. (49) Jean, J. M.; Friesner, R. A.; Fleming, G. R. J . Chem. Phys. 1992, 96, 5827. (50) Figueirido, F. E.; Levy, R. M. J . Chem. Phys. 1992, 97, 703. (51) Yamamoto, T. J . Chem. Phys. 1960, 33, 281. (52) Ando, K.; Hynes, J. T. In Structure and Reactivity in Aqueous Solution; Cramer, C. J., Truhlar, D. G., Eds.; American Chemical Society: Washington, DC, 1994; p 143. Ando, K.; Hynes, J. T. J . Am. Chem. SOC., accepted for publication. (53) Heller, E. J. J . Chem. Phys. 1976, 65, 4979. (54) Makri, N.; Miller, W. H. J . Chem. Phys. 1987, 86, 1451. ( 5 5 ) Stockli, A.; Meier, B. H.; Emst, R. R. J . Chem. Phys. 1990, 93, 1502. (56) Kearley, G. J.; Fillaux, F.; Baron, M. H.; Bennigton, S.; Tomkinson, J. Science 1994, 264, 1285. (57) Perrin, C. L. Science 1994, 266, 1665. (58) DolinSek, J.; Zalar, B.; Blinc, R. Phys. Rev. B 1994, 50, 805. (59) Tomasi, J.; Persico, M. Chem. Rev. 1994, 94, 2027. (60) De Vries, A. H.; Van Duijnen, P. Th.; Juffer, A. H.; Rullman, J. A. C., Dijkman, J. P.; Merenga, H.; Thole, B. T. J . Compur. Chem. 1995, 16, 37. (61) Florian, J.; Scheiner, S. J . Comput. Chem. 1994, 15, 553. (62) Florian, J.; Hrouda, V.; Hobza, P. J . Am. Chem. SOC.1994, 116, 1457. (63) Jordan, P. C.; van Maaren, P. J.; Mavri, J.; van der Spoel, D.; Berendsen, H. J. C. J . Chem. Phys., accepted for publication. (64) Voth, G. A.; Chandler, D.; Miller, W. H. J . Chem. Phys. 1989, 91, 7749.
Jp9507206
