J. Phys. Chem. B 2004, 108, 7963-7968
7963
A Variational Definition of Electrostatic Potential Derived Charges Alessandro Laio* and Francesco Luigi Gervasio Computational Science, Department of Chemistry and Applied Biosciences-ETH Zurich, USI Campus, Via Giuseppe Buffi 13, 6900 Lugano, Switzerland
Joost VandeVondele Department of Chemistry, Cambridge UniVersity, Lensfield Road, Cambridge, UK
Marialore Sulpizi and Ursula Rothlisberger LCBC, EPFL Lausanne, Switzerland ReceiVed: January 26, 2004
In a recent work [Laio, A., VandeVondele, J., Rothlisberger, U. J. Phys. Chem. B 2002 106, 7300] a novel method has been proposed to define dynamical electrostatic potential derived (D-RESP) charges for systems described within a quantum mechanics/molecular mechanics (QM/MM) scheme. Here, we derive the analytic dependence of these charges on the quantum charge density and on the atomic positions. This Variational property can be exploited for defining interaction potentials between the quantum and the classical subsystems that depend explicitly on the value of the D-RESP charges. Such potentials can be used for a multitude of different purposes, such as improving the computational efficiency of the electrostatic coupling between the QM and the MM subsystems and for defining a QM/MM analogue of the exclusion schemes commonly used in classical biomolecular force fields.
Atomic point charges provide a simple mean to describe the charge distribution of a system at an atomic level. Due to their intuitive nature, it is common to write Hamiltonians as a function of these charges. Important examples are the force fields for the simulation of large molecules by molecular dynamics,1-3 in which the electrostatic potential is represented as a superposition of Coulombic fields centered on the atoms (or on additional auxiliary sites), and the assigned charges are considered as fixed parameters during the dynamics. Obviously, such a description is not sufficient for describing events in which significant rearrangements of the electronic density occur, such as, e.g., chemical reactions or charge-transfer processes. Furthermore, the response of an atom to changes in the environment cannot be accounted for, and polarizable models have been developed in order to circumvent this limitation.4 One possible manner to define atomic point charges is to fit their value to the electrostatic potential. The charges obtained in this manner are called ESP (respectively RESP, if a restraining penalty function is used during the fitting procedure)5,6 and provide an accurate description of the electrostatic properties of the system. RESP-type charges are nowadays largely employed in force fields for classical molecular dynamics simulations.1 In their standard definition, the RESP charges are fitted in such a way as to reproduce the electrostatic field due to the charge density of a system determined by a quantum chemical electronic structure method. For this purpose, the field is evaluated on a grid or on a surface around the system. Recently, in the context of mixed quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations,8-10,17,18 a novel procedure to compute RESP charges has been introduced.7 In this approach, the charges are fitted to the electrostatic field due to the quantum charge density evaluated on the
classical (MM) atoms close to the quantum (QM) system. The RESP charges can be computed on the fly during the QM/MM dynamics and can be used as an indicator of the chemical state of the system. Due to the fact that they evolve during the simulation time, the charges computed with this procedure have been named dynamical RESP charges (D-RESP). In this work, we show that it is possible to derive the analytic dependence of the D-RESP charges on the coordinates of the system, i.e., on the position of the atoms and on the electronic density. To the best of our knowledge, this dependence has never been explicitly derived for the standard RESP charges, whose values depend on the electronic density in a highly nontrivial and nonlinear manner. A Variational definition for atomic point charges can be used for several purposes in which it is necessary to use interaction potentials that depend explicitly on the value of the charges. Examples are implicit solvation models, in which the solvent modifies the electronic structure of a quantum solute by a term that depends on charges of the quantum atoms.20 In this work, we demonstrate that the variational definition of the charges can also be exploited to improve the accuracy and the efficiency of QM/MM calculations. In particular, in section 2, we will show that the D-RESP charges can be used to couple the QM system to the MM atoms in a QM/MM simulation within a fully Hamiltonian framework that includes polarization of the electronic density at a very low computational price. Finally, in section 3 we will demonstrate that the D-RESP coupling scheme allows the selective exclusion of interactions between pairs of QM and MM atoms from the QM/MM interaction potential. This permits to define the QM/MM equivalent of the exclusion scheme already used in the classical force field approach and it allows the accurate description of the electronic structure of systems in which the QM atoms are
10.1021/jp0496405 CCC: $27.50 © 2004 American Chemical Society Published on Web 05/15/2004
7964 J. Phys. Chem. B, Vol. 108, No. 23, 2004
Laio et al.
at covalent bonding distance from charged MM atoms, greatly extending the possibility to perform reliable and robust QM/ MM calculations. Method The D-RESP charges are fitted in order to reproduce the electrostatic field due to the QM charge density during the dynamics minimizing, with respect to the set qDi , i∈QM, the norm
E)
(
qDi
∑ ∑ j∈NN i∈QM r
)
2
- Vj + w q
ij
(qDi - qHi )2 ∑ i∈QM
∫ dr F(r)Vj(|r - rj|)
∫ dr Fel(r)πi(r) - Zi
πi(r) )
qDi )
ti ) Hik )
VRESP(r) )
∑j (∑i Aji qDi - Tj)2
{
{
δVRESP δF
VRESP(r) )
τj(r) ) (5)
(6)
The j summation runs over the set TS (target set), with TS ) NN∪QM. For the sake of clarity, summations on lower and
(8)
)
∑i
∂VRESP δqDi ∂qDi
δF
j τj(r) ∑j qEFF
(9)
where j ) qEFF
and the response matrix A is given by
1/r j ∈ NN Aji ) w δij q ij j ∈ QM
∑j Aji Ajk
The functional derivatives of the D-RESP charges with respect to the density follows directly from the explicit expression of the D-RESP charges. From eqs 5-7, we have
where the target vector T is given by
Vj j ∈ NN Tj ) H wqqj j ∈ QM
∑j AjiTj
(4)
where Fati is the atomic (pseudo) valence charge density of the atom i and Zi ) ∫ dr Fati (r) is its valence.19 To derive the analytic dependence of the D-RESP charges on the atomic coordinates and on the density, we rewrite eq 1 in the compact form
(7)
The charge set qDi , i∈QM can be used to construct a potential VRESP explicitly depending on the D-RESP charges, that can be added to the QM/MM Hamiltonian. Two examples of possible choices for this potential will be provided in the following sections. The goal of this section is to give an explicit derivation of the potential on the electrons and the forces on the atoms due to VRESP. The potential on the electrons VRESP(r) is obtained by taking the functional derivative of VRESP with respect to the charge density F:
Fati (|r - ri|)
∑k
∑k H-1 ik tk
with
(3)
Fatk (|r - rk|)
)
The solution of this system of linear equations (one for every QM atom) is given by
(2)
is the field due to the QM charge density F(r) on the MM atom j, Vj(|r - rj|) is a Coulomb potential suitably modified at short range in order to avoid spurious overpolarization effects.18 The summation in the first term in eq 1 runs over the subset NN of MM atoms that are within a shell of thickness rMP around any atom of the quantum system (see also section 2 for a more rigorous definition). The second term in eq 1 is introduced in order to restrain the value of the charges qDi close to the corresponding Hirshfeld value qHi . The term wq is an adjustable weight parameter introduced in order to reduce the unphysical component of the charge fluctuations observed in an unrestrained ESP fit,6,18 and the Hirshfeld charge qHi of the QM atom i is given by
qHi )
∑j (∑k Ajk qDk - Tj)Aji ) 0
(1)
where rij is the distance between the QM atom i and the MM atom j,
Vj )
higher indexes are intended to run on the QM set and the TS set, respectively. The least-squares problem in the charge set qDi , i ∈QM reads (∂/∂qDi ) ) 0, i.e.,
∑ ik
{
∂VRESP j H-1 ik Ak D ∂qi
V (|r - rj|) j ∈ NN δTj ) j wqτj(r) δF j ∈ QM
(10)
(11)
Hence, the external field on the electrons due to VRESP can be expressed as a linear superposition of the (smoothed) Coulombic fields Vj(|r - rj|) centered on the NN atoms and the Hirshfeld probability densities πj(r) for each QM atom. The coefficients of the linear superposition have the dimension of a charge and j . are denoted as qEFF The forces on the atoms are obtained by taking the derivatives of VRESP with respect to the atomic positions. Hence, we have,
Electrostatic Potential Derived Charges
J. Phys. Chem. B, Vol. 108, No. 23, 2004 7965
using definition (7) of qDj and eq 10,
dVRESP
Fi ) -
)-
dri -
∂ri +
F(1) i
-
∂ri
∂VRESP
))
∂VRESP
F(2) i
∑ jkm
+
∑j
D-RESP Derived Electrostatic Coupling
∂VRESP ∂qDj ∂qDj
∂ri
m ∂VRESP d(H-1 jk Ak ) Tm ∂qj dri
dTj
j ∑j qEFF dr
i
F(3) i
(12)
(2) (3) where eq 12 is a definition for F(1) i , Fi , and Fi . (1) The term Fi ) - ∂VRESP/∂ri comes from the explicit dependence of VRESP on the atomic positions, and is the classical force due to a set of D-RESP charges located on the QM atoms. (3) The terms F(2) i and Fi would be zero if the D-RESP charges would not depend on the atomic positions and on the electronic density. They include, respectively, the derivatives with respect to the atomic position of the response matrix A and of the target vector T. The explicit expression for F(2) i is
{
fik] ∑k∈QM [(Υi -k Ti)φkkfik +k qkqiEFF D k - ∑k∈NN [(Υ - T )φkfi + qi qEFF fki ]
F(2) i )
i ∈ NN i ∈ QM
with
φk )
∑l
Υk )
∑l
fji ) and
{
Vcoulomb ) EC
dri
VEC )
) (rj - ri)/r3ij, j ∈ NN
∫ dr F(r) ∫
∑
∑
j j∈NNqEFF
i ∈ QM
where the derivatives with respect to ri are given by
dri
)
Fati ′(|r - ri|)
r - ri
∑k Fatk (|r - rk|)
|r - ri|
∑
qMM qDj i
i∈EC, j∈QM
i ∈ NN
∫
dπj(r)
(13)
with a potential explicitly depending on the D-RESP charges of the form
Akl ql
dAji
∫ dr F(r)Vi(|r - ri|) ∑ qMM i
i∈EC
dVi(|r - ri|) dri dπj(r) j dr Fel(r) + j∈QMqEFF dri dgi(|r - ri|) dr Vj(|r - rj|) dri
qiEFF
F(3) i )
∂VRESP H-1 lk ∂ql
The computation of the external field due to the MM charges on the real space grid on which the electronic density is represented is, within a plane wave based approach, rather expensive, since it requires a number of operations that scales linearly with the number of grid points and the number of MM atoms that are explicitly included in the interaction Hamiltonian. In our QM/MM approach, the MM atoms that are further away than a fixed distance rEC from any atom of the quantum system are coupled to the QM charge density by a suitable multipolar expansion of the full electrostatic interaction,18 and only the MM atoms within a distance rEC from any QM atom (NN atoms) are explicitly coupled to the quantum system by terms of the form of eq 2. This allows, in most cases, a very significant reduction of the computational burden without a loss of accuracy.18 Still, the multipolar expansion is only valid to represent the interaction with MM atoms that are at a distance from the geometric center of the QM system that is at least double the maximum distance between any QM atom and the geometric center. Hence, for a large QM system, the number of NN atoms can be very large, in particular if the QM system is significantly elongated in one direction. This implies, in some cases, a significant increase of computational cost with respect to a pure QM calculation. To reduce the computational cost, we replace, for a subset EC (esp coupled) of the NN atoms, the expensive Coulombic interaction term
(δij - πj(r))
dgi(|r - ri|) r - ri ) g′i(|r - ri|) dri |r - ri| r - ri dVi(|r - ri|) ) V′i(|r - ri|) dri |r - ri| where the prime indicates derivatives with respect to the argument.
rij
(14)
It should be noticed that, although this potential has the functional form of a point charge Coulombic interaction, it depends nontrivially on the electronic density and on the atomic positions, since the D-RESP charges qDj depend on the positions of all the QM atoms and all the NN atoms, and moreover, depend functionally on the density. To construct the set of atoms EC, a charge group based pair list is first computed. We denote by r(i) the minimum distance between the MM atoms belonging to the charge group CG(i) and all the QM atoms:
r(i) )
min
k∈QM,j∈CG(i)
rjk
Then, a new cutoff is specified, denoted by rNN, with rNN e rEC. If for the charge group i r(i) < rNN, all the atoms with nonzero charge belonging to CG(i) are assigned to the NN set. If rNN e r(i) < rEC, the atoms belonging to CG(i) are assigned to the EC set. The atoms belonging to charge groups further away than rEC from the QM system are coupled by the multipolar interaction Hamiltonian described in ref 18. To provide a converged value for the multipolar expansion, rEC must be at least twice the maximum distance between any QM atom and the geometric center of the QM system. This can easily be afforded since we introduced the intermediate layer of atoms that are coupled using the D-RESP charges, and the computational cost of this intermediate layer is totally negligible.
7966 J. Phys. Chem. B, Vol. 108, No. 23, 2004
Laio et al.
Figure 1. Schematic view of the structure of the guanine/cytosine QM/MM system. The QM subsystem is represented in red, the MM subsystem in black, while the capping hydrogens (Hc) are in green.
As we have shown in section 1, the forces and the potential deriving from an Hamiltonian that is a function of the D-RESP charges are fully specified by the partial derivatives with respect j and of φk) to the charges (appearing in the definition of qEFF and with respect to the atomic positions (appearing in F(1) i ). For a potential of the form of eq 14, we have (∂VEC/∂qDi ) ) ∑k∈EC(qMM k /rik), and hence j qEFF
qMM k j ) H-1 mi Ai m,i∈QM,k∈EC rmk
∑
(15)
The external field VEC(r) deriving from VEC can be expressed, by eq 9, as a superposition of Coulombic fields centered on the NN atoms and Hirshfeld probability distributions, with efficacious charges given by eq 15. The field VEC(r) has the important property to be almost equal, on the position of the QM atoms, to the electrostatic field due to the MM atoms belonging to the set EC. In fact, if ri is the position of the QM atom i, we have
VEC(ri) )
qMM k j j H-1 ml Alτ (ri) mlj k∈EC rmk
∑∑
i ∈ QM
Since τj(ri) = Aji (see eqs 11 and 6), we have
qMM qMM k k -1 j j Hml Al Ai ) VEC(ri) = mlj k∈EC rmk k∈EC rik
∑∑
∑
(16)
The computational cost for evaluating the electrostatic field due to atoms belonging to the intermediate layer EC is totally negligible, since it can be computed together with the Coulombic interaction with the NN atoms redefining the charges on the NN atoms as qMM + qEFF i i . Also the efficacious field proportional to the Hirshfeld probability densities πi(r) can be computed with no significant computational overload with respect to the evaluation of the Hirshfeld charges. It should be anyway noticed that the quality of the D-RESP set depends on the number of NN atoms, since the D-RESP charges are fitted exactly in order to reproduce the field (eq 2) evaluated on the
NN atoms. Hence the two cutoffs have to be chosen carefully in order to achieve a good compromise between accuracy and computational efficiency. We test this electrostatic coupling scheme on a QM/MM simulation of a fully hydrated double strand DNA decamer, namely d(GpCpGpCpGpCpGpCpGpCp) whose crystal structure is experimentally known11 (see Figure 1). Our choice was motivated by the fact that the atomic and electronic structure of this laboratory-realizable molecule have already been thoroughly characterized via full ab initio optimizations,12 giving us the opportunity to compare the QM/MM results to the full ab initio ones on a large system of biophysical relevance. The elementary cell contains overall 654 heavy atoms and 540 hydrogen atoms. All the simulations used periodic boundary conditions. The calculations were performed in the framework of density functional theory using the BLYP functional14,15 and Martins-Troullier pseudopotentials for the core electrons.13 The Kohn-Sham orbitals were expanded up to an energy cutoff of 70 Ry. All the calculations were made with the CPMD code.16 For further technical details, see ref 12. Within the QM/MM framework, we describe a single guanine/cytosine base pair at the QM level, with a cut passing between the QM and MM system between the nitrogen of the base N and the C1′ of the sugar ring. This cut has been chosen because it preserves the aromatic nature of the base and furthermore the quantum moiety corresponds to a molecule stable in the gas phase whose reference properties are experimentally known. The valence of the nitrogen atom at the QM/ MM border is saturated with an extra (dummy) hydrogen atom. The classical region includes all of the remaining DNA atoms and the solvent and is treated with the Amber force field.1 To compare the performance of the D-RESP-based coupling scheme with the one presented in ref 18, we compute, for a snapshot configuration of this system, the total energy (Figure 2a) and the dipole moment (Figure 2b) of the QM subsystem as a function of different cutoffs. If the intermediate layer of D-RESP coupled atoms is not included, the energy is converged within ≈ 1 kcal/mol for rNN ≈ 8 Å, whereas the dipole continues to vary significantly up to rNN ≈ 10 Å. If the D-RESP coupling scheme is applied, a similar accuracy is obtained for rEC ) 15 Å and rNN ) 5 Å. This implies a very significant reduction of
Electrostatic Potential Derived Charges
J. Phys. Chem. B, Vol. 108, No. 23, 2004 7967
Vexcl ) -
qiMM qDje e
∑
(ie,je)
rieje
(17)
to the total Hamiltonian of the system, where the summation runs over the pair interactions that one wants to exclude from the Hamiltonian (ie and je are assumed to run over a subset of MM and QM atoms, respectively). Once again, although functionally identical to a Coulomb interaction, the potential in eq 17 is instead a complicated many-body expression, since the D-RESP charges on the QM atoms depend on the position of the QM and the MM atoms and on the electronic density. The potential Vexcl(r) due to this extra term will have the form of eq 9, with
j )qEFF
qiMM e Hj-1k Ajk rieje e
∑∑
(ie,je) k
(18)
The term Vexcl(r) has the property to be very small on any QM atom not involved in any exclusion, and to be approximately equal to the field due to the interactions that have to be excluded on the other quantum atoms. This property is a consequence of eq 1, i.e., of the choice to fit the D-RESP charges on the field evaluated on the atomic positions, and holds for every reasonable choice of the NN set and of the restraining constant wi. For wi ) 0 it holds exactly. In fact, using eqs 9 and 18, we have, at the position ri of the QM atom i,
Vexcl(ri) ) -
∑∑
(ie, je) jk
)-
the computational burden, since the number of atoms that are explicitly coupled to the QM system scales cubically with rNN. Exclusion of Selected Coulombic Pair Interactions from the QM/MM Hamiltonian A rather severe problem in QM/MM simulations concerns the proper treatment of electrostatic interactions between QM and MM atoms separated by only one or two covalent bonds (the so-called “frontier atoms”). All the classical interactions involving the dummy atoms are excluded,17 while the only nontrivial spurious term is the electrostatic one. An exclusion scheme, analogous to that used in the classical force field, can be introduced for these spurious electrostatic terms. The D-RESP coupling scheme provides, in fact, an elegant manner to avoid these artifacts. The idea is to exclude the unwanted pair interactions from the QM/MM electrostatic potential by adding a term of the form
qiMM e Hj-1k Hki rieje e
∑∑
(ie, je) k
Figure 2. Total energy and dipole as a function of the cutoff radii for two different coupling schemes for the QM/MM system depicted in Figure 1. The zero of the energy is defined by the value obtained for rEC ) rNN ) 15 Å. Red line with squares: the coupling scheme of ref 18. All the MM charge groups at a distance up to the cutoff radius are included in the NN set. Blue line with diamonds: D-RESP-based coupling scheme. The charge groups at a distance between the cutoff radius and 15 Å are coupled to the QM system by a potential of the form of eq 14. The MM charge groups further than 15 Å are coupled with the QM subsystem using the scheme of ref 18.
qiMM e Hj-1k AjkVj(|rij|) rieje e
qiMM e )δjei (ie, je) rieje
∑
where we used Vj(|rij|) = (1/rij) ) Aij and the definition (eq 8) of Hki. Such an exclusion scheme is particularly important if the cut between the QM and the MM subsystems runs in a region in which the atomic charges are large. Consider, e.g., the system depicted in Figure 1. The capping hydrogen atoms (Hc in Figure 1) do not exist in the real system, and are introduced only in order to saturate the valence of the nitrogen atom N of the bases. These hydrogen atoms are extremely close to the MM atoms of the sugar and experiences their unscreened electrostatic field. Their equilibrium position, e.g., is only 0.4 Å away from the C1′ carbon atom of the sugar ring, whose effective atomic charge is -0.012 in the case of cytosine and 0.036 in the case of guanine, and is 1.4 Å from the oxygen O4′ of the sugar, whose charge is -0.37. The electronic structure of this QM/MM system is seriously affected by the spurious interactions between the capping hydrogen and the close-by classical atoms. The qualitative description of the LUMO is wrong and large fluctuations are observed: the LUMO should be localized on the cytosine ring and should have a π* symmetry,12 while it is localized on the cytosine oxygen or alternatively on the phosphate backbone (see Figure 3a,c). Moreover, the system is forming an unphysical
7968 J. Phys. Chem. B, Vol. 108, No. 23, 2004
Figure 3. Isosurfaces of the HOMO (in blue) and LUMO (in red) of the cytosine/guanine base pair with full electrostatic coupling (A,C) and “exclusion” (B,D). For the system in panels C and D, the positions of the dummy hydrogens are relaxed. At the isosurfaces represented, the charge density has a value of 10-2 electrons Å-3. The MM atoms surrounding the system are not represented except in panel C, in which all the MM atoms closer than 10 Å from the LUMO isosurface are shown (thinner lines). The LUMO in this case is localized close to a MM phosphate oxygen.
hydrogen bond between the capping hydrogen and the O4′ of the sugar ring. The HOMO-LUMO gap, that for the same system described at the full QM level is approximately 2.6 eV, is lowered to less than 1 eV, and this remarkable difference is expected to significantly affect the chemical properties of the system. A mere inspection of the structural properties upon geometry optimization of the QM/MM system did not indicate this problem, since the structure of the system across the cut is held together by classical bond, angular, and dihedral terms that are not affected by problems in the electronic structure. A few ps of molecular dynamics, however, lead to the breaking of the N-Hc bond, causing an instability of the system. The most simple practical remedy to this problem would be to set the values of the charges of the MM atoms close to the capping QM atoms to zero. Unfortunately, this simple solution would significantly change the electrostatic properties of the MM system in a region critically close to the quantum system. In the example of Figure 1, the two sugar rings connected to the QM system would have to be described with ad hoc charges with undesirable consequences on the properties of the system. Setting, e.g., the charge on the sugar oxygen O4′ to zero would change the hydrogen bonding pattern to nearby QM water molecules. In our DNA model, the introduction of the “exclusion” scheme resolves the problems due to the capping hydrogen previously described. In Figure 3b and d, the isosurfaces representing the HOMO and the LUMO of the QM subsystem are depicted. The geometry is exactly the same as in Figure 3
Laio et al. a and c, but for the position of the capping hydrogen that was allowed to relax to its equilibrium position. With the new electrostatic treatment the qualitative description of the LUMO is now correct. Moreover, the HOMO-LUMO gap is now = 2.2/2.3 eV, much closer to the full quantum gap in the same geometry, and there are no more instabilities due to breaking of the N-Hc bond. In conclusion, the exclusion scheme we presented here allows for an accurate QM/MM treatment of a complex system like DNA. Although in the example we presented the QM subsystem is really minimal (only a single base pair), all the main features of the electronic structure of the full quantum system (e.g., HOMO-LUMO gap and nature of the frontier orbitals) are correctly reproduced, and the electrostatic coupling among the QM and the MM subsystems is correctly kept into account. This represents a major improvement with respect to other available schemes. References and Notes (1) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (2) CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations; Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. The Energy Function and Its Parameterization with an Overview of the Program; MacKerell, A. D., Jr.; Brooks, B.; Brooks, C. L.; Nilsson, L., III; Roux, B.; Won, Y.; Karplus, M. In The Encyclopedia of Computational Chemistry; Schleyer, P. v. R. et al., Eds.; John Wiley and Sons: Chichester; 1998; Vol. 1, pp 271-277. (3) Scott, W. R. P.; Hu¨nenberger, P. H.; Tironi, I. G.; Mark, A. E.; Billeter, S. R.; Fennen, J.; Torda, A. E.; Huber, T.; Kru¨ger, P.; van Gunsteren, W. F. J. Phys. Chem. A 1999, 103, 3596. (4) Palel, S.; Brooks, C. L. J. Comput. Chem. 2004, 25, 1. Kaminski, G. A.; Stern, H. A.; Berne, B. J.; Friesner, R. A.; Cao, Y. X.; Murphy, R. B.; Zhou, R.; Halgren, T. A. J. Comput. Chem. 2002, 23, 1515. Liu, Y.-P.; Kim, K.; Berne, B. J.; Friesner, R. A.; Rick, S. W. J. Chem. Phys. 1998, 108, 4739. (5) Momany, F. J. J. Phys. Chem. 1978, 82, 592. (6) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (7) Laio, A.; VandeVondele, J.; Rothlisberger, U. J. Phys. Chem. B 2002, 106, 7300. (8) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (9) Singh, U. C.; Kollman, P. A. J. Comput.Chem. 1986, 7, 718. (10) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (11) Ban, C.; Ramakrishnan, B.; Sundaralingam, M. Biophys. J. 1996, 71, 1215-1221. (12) Gervasio, F. L.; Carloni, P.; Parrinello, M. Phys. ReV. Lett. 2002, 89, 108102. (13) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993. (14) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (15) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (16) CPMD V3.5 (Copyright IBM Corp 1990-2001, Copyright MPI fuer Festkoerperforschung Stuttgart, 1997-2001). (17) Sherwood, P. In Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; NIC Series, Ju¨lich, Germany, 2000; Vol. 1, p 257. (18) Laio, A.; VandeVondele, J.; Rothlisberger, U. J. Chem. Phys. 2002, 116, 6941. (19) Hirshfeld, F. L. Theor. Chim. Acta 1977, 44, 129. (20) Tianhai, Z.; Jiabo, L.; Hawkins, G., D.; Cramer C., J.; Truhlar, D., G. J. Chem. Phys. 1998, 109, 9117.
