Introduction

Density-functional theory (DFT [1]) is currently the method of choice for the study of structural, electronic, and dynamical properties in materials science. The reason for this lies both in the ability of DFT schemes to treat models of relatively large sizes compared to other first-principles methods, and in the degree of accuracy that can be achieved for diverse systems in a systematic fashion. In particular, molecular dynamics (MD) based on DFT (Car–Parrinello MD (CPMD) [2]) has added flexibility to the standard implementations and has made it possible to study dynamical processes such as structural transformations of covalently bonded materials as well as chemical reactions in condensed phases.

There are several reasons for introducing DFT-based simulations in biology: The principal motivation is a need for an ab initio description of the interatomic interactions whenever standard force fields do not yield reasonable results. This is true for problems such as the binding of ligand–protein complexes, the structure of metalloproteins, or the mechanisms driving enzymatic reactions that involve bond redistribution. Thus, DFT-based MD seems in principle to offer an ideal method that combines the accuracy of DFT and the power of MD necessary for biological systems which are intrinsically dynamical. However, in practice, the applicability and usefulness of CPMD simulations in the biological field are still unclear.

First of all, it must be emphasized that biological materials systems are much more complex than traditional materials systems, in which the performance of current DFT implementations has been tested or at least has long been under control. In the description of biological materials systems, all of the problems encountered by current DFT schemes, such as accounting correctly for Van der Waals and hydrophobic interactions, are heightened because of the crucial role played by diverse interactions in determining both stability and functionality. Moreover, the relevant energy scales in biochemical materials are at most of the order of kT and thus require a higher accuracy than for traditional materials: This is indeed a difficulty that DFT shares with all other nonempirical methods. Compared to classical MD, DFT–MD also presents practical disadvantages, including limitations in the size range (up to 1000 atoms) and the shorter time scales for simulations (1–10 ps). The size limit has increased rapidly over the past decade (by two orders of magnitude), primarily because of the advent of powerful parallel computing [3]. The time-scale limitation is more difficult to overcome, even by using the algorithms common to classical MD that were devised specifically to increase either the time step or the configuration sampling [4]. Moreover, in principle, activated processes such as enzymatic reactions can be simulated by imposing specific constraints [5], but this possibility remains highly dependent on the specific system and is thus uneven. On the other hand, a “brute-force” DFT–MD simulation is currently impractical for most biological systems and, if confined to oversimplified models and/or irrelevant time scales, may lose its relevance to experiment. However, there are ways in which one can exploit all benefits of this method: either by selecting case studies for which intelligent (but size-limited) models can provide guidance to the physical understanding of more complex systems, or by applying it to the simpler biological models usually investigated in laboratory experiments (in the crystal phase). This track presents useful tests of the computational schemes by extracting important electronic parameters for use in force fields and that are transferable to larger systems, or by combining DFT and classical methods as a special case of coupled quantum-mechanical (QM) and molecular-mechanics (MM) methodologies.

This paper reports on recent attempts to use DFT and DFT–MD as indicated above, with special emphasis on our own projects, both past and current. In particular, details are given about recent calculations, the results of which we present here for the first time. Therefore, in some respects, this paper may appear a little unbalanced. Indeed, it is not our intention either to present an exhaustive review or to describe a single application in all possible detail. As a basis for our discussion, we have selected a number of papers that well exemplify the state of the art of the calculations. Also, we attempt to provide a critical survey of the field in order to put our own work into its proper context. Finally, future directions for this type of simulation are also outlined.

Structure and dynamics with DFT calculations

Static calculations (local geometry optimizations) on key elements of biomolecules such as DNA base pairs “in vacuo” have been performed in a routine manner for some time. Generally, the primary purpose of these calculations was to examine the performance of DFT methods with gradient-corrected exchange-correlation functionals (GGA) versus post-Hartree–Fock methods [6] in accounting for structure, dipole moments, and energetics. Regrettably, calculations of this kind have long suffered from the limited accuracy of the computational approach used (e.g., small basis sets for the electronic wave functions) and thus were unable to provide reliable reference data for useful comparisons. Only very recently has an accurate and systematic DFT study of Watson–Crick base pairs appeared [7] which has allowed one to assess the degree of convergence of different basis sets using localized wave functions and the performance of different (gradient-corrected) exchange-correlation functionals, as well as to make a critical comparison with X-ray crystallographic data. The latter established that the molecular environment present in the crystal, especially the counter ions and water molecules, has an important influence on its structural properties.

A few static DFT calculations [8] were also performed on crystal phases with the gradient-corrected exchange-correlation functionals (GGA) using plane-wave basis sets (e.g., [9–11]). Structural optimization generally begins with the X-ray refinement, as the results for the structure are primarily a test of the computational scheme, but they can also include some complementary features such as the hydrogen pattern. However, the new information provided by calculations is derived from the electronic structure. For instance, in the example reported in Reference [10], a platinum-modified nucleobase pair was studied in the crystal phase. Detailed comparison with calculations in vacuo in the absence of the platinum complex allowed simultaneous capture of the effects of metalation on the electronic structure, the chemical bonding, and the geometry of the nucleobase pair in the complex.

Applications of several DFT schemes (GGA and hybrid HF–DFT) to simple models of metal enzymatic centers have recently been reviewed in Reference [12]. Interestingly, comparing these model calculations to experimental data, one could in some cases show that treating the protein explicitly was not necessary. In fact, its effects could either be indirectly ignored, as in the case of a metal center in the blue-copper protein complex, or represented as a dielectric medium, as in the example of electron transfer in the photosynthetic reaction center of Rhodobacter sphaeroides. It would be interesting to understand why and to what extent this is true in other cases as well.

The important case of the Fe(II)- and Fe(III)-cytochrome c heme-protein complexes involving the nature of the bonding of iron to the protein prosthetic groups has recently been studied using model complexes, namely FeP-Im-DMS (P = porphyrin, Im = imidazole, DMS = dimethylsulfide) [13], in the plane-wave scheme of CPMD [8]. The effect of the environment on the iron–sulfur bond was also investigated by calculating only the electrostatic potential of the protein, using effective point charges obtained from the AMBER parameterization scheme [14]. Surprisingly, this effect was found to be very weak, a result contrary to chemical intuition. No explanation was provided, so that it is still unclear whether this was due to some subtle chemical reason or to an artifact of the modeling.

The possibility of explicitly treating proteins a few hundred atoms in size as well as studying their electronic properties and chemistry exhaustively has recently been realized. In particular, refinement of protein structures, which nowadays relies on ad hoc parameterized (empirical) force fields, can be complemented with DFT structural studies. We now describe a case we have investigated in detail in the plane-wave framework of CPMD: a relatively small protein, 1PNH, a scorpion toxin (PO₅–NH₂) analog with high affinity for apamin-sensitive potassium channels. As shown in Figure 1, it consists of 31 amino acid residues (approximately 500 atoms) and has a well-defined globular structure that has been resolved by nuclear magnetic resonance (NMR) [15]. By using a local optimization procedure, the structure of the protein was first refined computationally in the gas phase [16], starting from available NMR measurements [15], and derived with the help of simulations using a classical force field. In this first calculation, the protein was kept in its standard charge state (+6). The calculated root mean square displacement (RMSD) of the computationally refined structure from the experimental structure was only 0.181 Å/atom, demonstrating that the experimentally observed minimum is stable also in our gas-phase model. Moreover, partitioning the total RMSD by atom type, we obtain RMSD(C) = 0.126 Å/atom, RMSD(N) = 0.116 Å/atom, RMSD(O) = 0.164 Å/atom, RMSD(H) = 0.199 Å/atom, and RMSD(S) = 0.567 Å/atom. Clearly the two structures are very similar, with the exception of the positions of the sulfur atoms, as is evident from the comparison in Figure 2. The reason for the RMSD(S) discrepancy is not related to any inadequacy of the computational scheme, but lies in the incorrect values of the original structural refinement [15]. This can be ascribed to a lack of accuracy in the parameterization of the force field regarding the sulfur atoms engaged in –S–S– bridges, which are notoriously difficult to handle in a classical framework. From this result, the usefulness of our ab initio procedure for the structural refinement becomes readily apparent as a tool in correcting errors of this general type.

Figure 1 Figure 2

In a second stage, the charge state of the protein was modified; i.e., the global charge was neutralized and its structure reoptimized locally. The calculated RMSD with respect to the original refinement from experiment was 0.141 Å/atom, very close to the one found for the structure in vacuum. This indicates that this (indeed large) change of the charge state of the protein does not destabilize the local energy minimum. Extended to a general protein, this result can be quite important in view of the fact that charge states of proteins cannot easily be determined a priori. However, it must be stressed that, although it may play only a minor role in determining the details of a local energy minimum, the specific charge state could still be quite important in determining the global folding pathway.

Some details of the computations are useful for obtaining a clearer idea of the computational effort they involve. On an IBM RS/6000* SP consisting of 128 200-MHz POWER3 processors, optimization of the electronic wave functions required around 300 seconds per step, while the geometry required 1.5 hours per step, and approximately 40 of these steps were necessary to obtain a converged result for the geometry. These values are reasonable enough to render this technique now feasible and routinely usable in the near future.

Static calculations are far from being exhaustive. Moreover, when testing a method such as DFT vs. experiment, local energy minimization starting from an educated guess is not sufficient and does not guarantee the absence of spurious minima and/or unrealistic conformations. More significantly, the study of the dynamics is vital for understanding the physical and chemical behavior of biomolecules [17] and is certainly essential for a correct description of the influence of water on their conformations and interactions. With respect to liquid water itself, CPMD simulations have shown that it can be treated at a reasonably accurate level within some of the DFT–GGA schemes [18–20]; structural and mobility properties as well as the solvation of diverse ions have been simulated successfully (see for example [21] and [22]). The ability to represent liquid water with a reasonable degree of accuracy is an important achievement, because if use is made of classical force fields, many difficulties are encountered in properly representing the intrinsic properties of water and specific interactions (see for example [23]); it opens up several new possibilities such as determining “structures in solution” by fully ab initio methods. This has recently been partially achieved for the monostrand major cisplatin– DNA¹ adduct [22], a system of intense investigation given the well-known antitumoral activity of cisplatin and the interest in the design of less toxic cisplatin derivatives. Here we report only few results from Reference [22], where the interested reader can find more details. As the initial model of the adduct, its experimentally known structure in the crystal phase was considered, and immersed in a water solution. First this system was equilibrated with classical MD simulations, keeping the cisplatin geometry fixed in order to avoid the expected dramatic effects of classical potentials. Subsequently, CPMD simulations [8] were run; these had previously been shown (also in Reference [22]) to reproduce well a number of molecular properties of cisplatin as well as the dynamics of its hydrolysis. Although very short (about 2 ps), the CPMD simulations of the cisplatin–DNA adduct were able to show the “immediate” effects of the aqueous solution, some of which can be argued from the results reported in Table 1, where the average sugar-backbone conformation (as depicted in Figure 3) found in this way (under the heading MD) is compared with the initial conformation (under the heading Start). Indeed, the changes are important and compare well with NMR data measured in aqueous solution for another, very similar, cisplatin adduct (see under the heading NMR in Table 1) that differs from the model investigated only in its lack of the 5' phosphate moiety. There is a noticeable discrepancy in the angles describing the nucleobase/sugar orientation between the experimental NMR-derived (GpG) results and the simulation in water that, apart from the time limitation of the CPMD runs, may also be ascribed to the known conformational flexibility of this molecule and to its sensitivity to the lack of the 5' phosphate moiety in the NMR structure. Also, an apparent effect of the interaction with the solvent is the reduction of the degree of destacking of the two guanines, which could reasonably be ascribed to the loss of hydrophobic interactions between the nucleobases of different complexes present in the crystal phase. In fact, the dihedral angle decreases from 81 to 61(9) degrees, a tendency also confirmed by the NMR value of 53 degrees measured for the analog complex in water.

Figure 3

Table 1	Sugar-backbone conformation of the cisplatin dpGpG adduct (Start = starting configuration; MD = result of CPMD simulation) compared to that of the [Pt(NH3)2]d(GpG) (from NMR) (see text). Notation of the torsion angles as in Figure 3. Reprinted with permission from [22], © 2000 American Chemical Society.

[Pt(NH3)2]–dpGpG Start   MD NMR*   5'-Guanidine –143   –147(15) ~–115 ß 165   132(12)   37   50(10) 52 100   90(10) 87(3) –133   –170(12) –162(4) 64   –70(11) –60   3'-Guanidine –127   –147(9) ~–115 –56   –68(9) –65 ß –159   –170(10) 180(3) 41   60(13) 50–60 143   129(31) 137

*	J. H. J. den Hartog, C. Altona, G. C. Chottard, G. P. Girault, J. Y. Lallemand, F. A. A. M. de Leeuw, A. T. M. Marcelis, and J. Reedijk, Nucl. Acids. Res. 10, 4715 (1982).

In the study of the 1PNH protein illustrated above, the effect of the aqueous solution on the local structural minimum was also investigated. Simulating a sufficient number of water molecules to represent a realistic water solvation shell around the protein was computationally too demanding for full DFT calculations. To obtain qualitative information, a much simpler strategy was adopted. First, the protein was fixed at the structure resulting from a DFT refinement and solvated with water treated classically using the simple point charge (SPC) water model [24] and the GROMOS96 force field [25]. The water was then thermalized using the GROMOS code.² Finally, only the 85 water molecules contained in a spherical shell having a radius of 2.5 Å around the protein were retained in the model, and the structure of the protein was reoptimized ab initio, keeping the water molecules fixed. Preliminary results yielded an RMSD value of 0.120 Å/atom with respect to the gas-phase geometry, confirming again the relative stability of this structure. Note that the CPU time required for such a calculation is nearly double that for the gas-phase protein.

We note that this type of calculation is not meant to be restricted to proteins of affordable size, but can also be applied to protein fragments. An additional interesting result of these calculations was the determination of effective atomic charges from the best fitting of the electrostatic potential (electrostatic potential procedure (ESP) [26]) on the entire protein. Atomic charges are important because they represent one of the simplest parameters for describing the electrostatic potential generated by a protein—a quantity of fundamental importance in determining the functioning of the protein itself. Often these charges are extracted from calculations on small fragments (e.g., amino acids) or from the Coulomb terms of the force field (with or without charge-equilibration procedures [27]), and assumed to be transferable. In particular, when charge-equilibration procedures are not used, the charge associated with an atom of a certain type, belonging to a certain amino acid, is assumed to be the same irrespective of its position. Our ab initio calculations show that this transferability hypothesis is not justified. In fact, the values we find for the atomic charges are such that the charge associated with a carbon atom on the same type of residue—but in a different location in the protein—can vary by more than one charge unit! Therefore, one can reasonably expect relatively large errors in the electrostatic potential calculated using charges that are not derived from an electronic structure calculation on the entire protein (or a sufficiently large protein fragment). It would be very interesting to explore the impact of these findings on the accuracy of the docking models widely used in drug design.

DFT–MD/MM approach

The “modern” approach combining DFT and MM is a special case of QM/MM methods for which the concepts were introduced in the late 1970s in a well-known publication by Warshel and Levitt [28]. Aimed especially at the simulation of enzymatic reactions, the QM/MM approach consists in partitioning the enzyme into two subsystems in order to treat the portion directly involved in the catalytic reaction at the QM level (the QM subsystem); the remainder, where events such as bond breaking and bond formation do not take place, is represented with an MM force field (the MM subsystem). Although the basic philosophy is simple, and practical implementation requires no special computational skills, the practical use of QM/MM methods requires an ingenious treatment of the interface between the two subsystems [29]. Apart from the specific force field used for the MM subsystem and the specific method (theoretical as well as computational) used for the QM subsystem, the model chosen for the description of the interface and the definition of coupling as well as of the information to be exchanged comprise the distinguishing features of the various applications deployed so far [30]. A particularly comprehensive and critical review has recently been published by Monard and Merz [31].

Several attempts have been made to combine the DFT–MD and MM methods as a natural extension of the CP method, using the CPMD code for the QM part interfaced with classical force-field codes (see for example [32]). An interesting example of applications to biologically relevant systems is that of exploiting the possibility of comparing, within the same framework, the activity of a real enzymatic system treated with the hybrid approach to that of low-molecular-weight analogs (biomimetic compounds) treated fully at the DFT level. In the case of the copper enzyme galactose oxidase, this type of investigation [33] has led to the conclusion that small synthetic analogs can be devised that reproduce the activity of the real system. This was based on the identification of a structural difference between a real synthetic analog that renders it much less active than the real enzyme, and also on the fact that in these calculations the structural and electronic properties of the QM portion (~80 atoms) appear to be insensitive to the behavior of the MM region. Note that the interface was represented by the scaled position method for link atoms, as in [32].

The advantage of using DFT and in particular the plane-wave scheme (as in CPMD) for the QM subsystem, however, comes to light only when much larger systems, i.e., of the order of a few hundred to a thousand atoms, are being treated; as mentioned above, this has become possible only very recently. In fact, increasing the size of the system to be modeled does not imply simply that one can explicitly treat a larger number of atoms around the active region at a higher level of theory, but rather that the problem of “how to treat the interface” becomes much less critical because the MM subsystem is farther away from the active region. Indeed, it should be emphasized that most of the problems encountered in QM/MM applications that are related to detailed modeling of the interface [34], such as the definition of the constraints that must be imposed on the frontier bonds to avoid artifacts in the chemical behavior of the system, arise from the limited size of the QM component, or at least are greatly amplified. Moreover, compared to a QM approach using localized functions, the approach based on plane waves avoids the accuracy problem related to the basis superposition error and allows a better description of the polarization of the electronic wave functions.

Within the DFT–MD/MM scheme we have recently considered BamHI, a prototype of Type II restriction enzymes [35], which are responsible for a particular cleavage of short DNA sequences at both strands, producing 5'-phosphate and 3'-hydroxyl groups. A number of experimental structural studies have elucidated the structure of BamHI–DNA complexes [36]. However, it is unclear how the catalytic reaction evolves, although a few hypotheses about the mechanism have been made [37]. For example, it is known that the catalytic reactions require the presence of divalent metal ions such as Mg, Mn, Co, Zn, or Cd. However, whereas the majority of them work as catalysts and differ only in the efficiency of their catalytic action, Ca acts as an inhibitor of the enzymatic reaction, and the reason for this is not yet understood. Clearly, the metal ions cannot simply be assumed as point charges, and we risk biasing the result if we make any unsubstantiated hypothesis on the difference of the cations. To answer these questions, we must rely on an ab initio approach and thus have decided to carry out a parallel investigation of BamHI with Ca and with Mg using DFT–MD/MM simulations.³

An important feature has emerged from experimental evidence that aided our modeling and simulations, namely that the catalytic reaction drives negligible conformational changes in the BamHI–DNA complex. In fact, the pre-reactive structure measured with the inhibiting calcium ion as well as the post-reactive one (measured in the presence of Mn) are very similar (RMS deviation of only 0.33 Å for the C positions). Thus, we chose the pre-reactive state of the complex with Ca as the initial configuration, as indicated in Figure 4, and locally optimized the structure within our scheme.

Figure 4

Out of a total number of about 6000 atoms, 297 (protein and DNA backbone, plus crystallographic water molecules and metal atoms) were included in the QM subsystem. The complex, maintained at a fixed geometry, was first solvated in water and equilibrated with GROMOS simulations using the SPC model [24]. The effect of the MM subsystem, considered as a geometrically fixed environment and surrounded by an additional ~2300 water molecules, was taken into account via the inclusion of electrostatic and Van der Waals interactions. The electrostatic potential generated by the MM atoms and calculated using the standard GROMOS96 force-field charges [25] was mapped over the tridimensional grid containing the QM atoms. Because there are no conformational changes during the reaction, this calculation could be performed only once. Van der Waals interactions between QM and MM atoms, on the other hand, were calculated at each optimization or MD step. These were represented by the GROMOS96 6–12 Lennard–Jones potential with a cutoff of 10 Å. The open bonds of the QM atoms at the QM/MM interface were saturated with hydrogen atoms [29] and also kept fixed during the entire simulation. This simple modeling of the interface was made possible by two considerations: on the one hand, the large size of the QM system (the atoms of the active region are at least 5 Å away from the interface), and on the other the a priori knowledge (mentioned above) that the structure of the MM region is quite insensitive to the reaction. In conclusion, in the above scheme, the electronic structure of the QM portion is optimized under the influence of the MM system, represented by the electrostatic potential, and the updating of the atomic coordinates takes into account the Van der Waals contribution of the interaction with the MM atoms. The optimized structure of the pre-reactive complex with Ca does not deviate significantly from the experimental one. The value of the RMSD evaluated over all of the QM atoms (water excluded) is 0.2 Å/atom, while the Ca–Ca distance remains at 4.3 Å and the Ca–O distances vary by about 0.1 Å. When Ca is replaced with Mg, the Mg–O distances, as expected, shrink by 0.3 Å on average, but the overall RMSD value, as calculated above, turns out to be approximately 0.2 Å/atom.

In an attempt to understand why Ca behaves differently from Mg, one might suggest that perhaps the propensity to donate two electrons causes the difference in catalytic activity. From the calculated electronic structure of the pre-reactive geometries, however, it can be clearly demonstrated that this is not the case. For both, two valence electrons are removed, and the localization of the transferred charge on the protein is essentially the same. This striking similarity can be seen by inspecting Figures 5(a) and 5(b), which illustrate the electrostatic potential mapped on an isodensity hypersurface for the two cases. From a methodological point of view, it is interesting to note the extent and location of the changes of the charge distribution induced in the QM subsystem by the interaction with the MM subsystem. The answer can be found in Figure 6, which depicts the difference of the electron density (also mapped on an isodensity surface) of the QM region when it is isolated and when it interacts with the MM region. The difference is indeed sizable, especially in the close environment of DNA, where the reaction will eventually take place. Therefore, a purely quantum-mechanical treatment of the ~300-atom subsystem would not be adequate.

Figure 5 Figure 6

We are currently performing simulations of the cleavage reaction that should yield useful information on the mechanism of the catalytic reaction and especially identify the reasons why Ca cannot be used as a cofactor.

Conclusions and perspective

There is no doubt that the value and importance of computer simulations in biology are becoming more and more recognized. It is also generally believed, however, that in order to make a real impact, progress will be necessary both in hardware and software. The latter will require not only coding but also (and especially) methodological advances. From the scientific point of view, it is clear that computational chemistry methods, and in particular the combination of quantum-mechanical and molecular-dynamics methods, could make a difference by providing the degree of accuracy, reliability, and transferability that force fields generally do not offer. We have briefly surveyed the current state of DFT and DFT–MD calculations, and have reported our own new developments and results, which are a faithful example of what can currently be achieved with CPMD simulations when combined with molecular mechanics. It is probably worth emphasizing that a DFT-based approach allows an accurate calculation of a number of physical and chemical properties, some of which can either be compared directly to experiment or are complementary to it. Beyond the structural characterization of a protein and its electron density distribution reflected in the electrostatic potential (of which we have provided examples), these calculations can in fact provide details on properties such as atomic chemical shifts, measured by nuclear magnetic resonance (see for example [38]), electron spin density distributions (see for example [39]), probed by electron spin resonance, and infrared spectra (see for example [40]). Also, energy barriers for chemical reactions in solution can be inferred from CPMD simulations (see [3] and references therein), which opens the road to a quantitative assessment of the mechanism of enzymatic reactions.

However, DFT–MD procedures are complex and relatively inefficient when applied to biomolecular systems. The route one has to pursue to render them as efficient as the traditional ones is full of difficulties, but certainly worth the challenge. Regarding an ab initio approach to biomolecular systems, there are still a number of fundamental issues to be addressed, in addition to the known problems encountered at the level of chemical accuracy. For example, regarding the DFT approach, the description of dispersion forces still is an unsolved and critical problem of current implementations. This is an obvious consequence of the nature of the approximations (local or semilocal) made thus far for the exchange-correlation density functionals [41], thus requiring novel and clever implementations [42]. Another problem pertains to the correct treatment of hydrogen-bonded systems and of proton-transfer events that may well require that the quantum nature of the nuclei be taken explicitly into account. An extension of DFT–MD that addresses these problems combines it with the path-integral method that has recently been introduced [43, 44]. Coupling classical and quantum approaches is mandatory, given the typical size of biologically relevant systems; however, performing such a coupling in a “correct” and general way is still a challenge in itself. Here, a systematic effort aimed at elucidating the failures and successes of a given scheme and/or at providing a thorough comparison of various schemes would be very useful indeed.

Acknowledgment

We wish to thank an anonymous reviewer for comments and suggestions that we feel have helped to improve our paper.

References and notes

Footnotes

¹ The term cisplatin denotes cis-diamminedichloroplatinum(II).
² GROMOS Software (C) Biomos bv, University of Groningen, The Netherlands, and ETH Zurich, Switzerland.
³ For the QM subsystem, a box of (22.95 × 25.69 × 20.75) ų was used, and a plane-wave cutoff of 70 Ry. Other details as in [16]. In particular, note that the QM subsystem is treated as an “embedded cluster”; i.e., it is treated as isolated (Hockney method), but the effect of the interaction with the MM system is taken into account.

Received September 20, 2000; accepted for publication January 24, 2001