A modified QM/MM Hamiltonian with the Self-Consistent-Charge Density-Functional-Tight-Binding Theory for highly charged QM regions

Hybrid Quantum Mechanical/Molecular Mechanical Methods For Studying Energy Transduction in Biomolecular Machines

Hybrid quantum mechanical/molecular mechanical (QM/MM) methods have become indispensable tools for the study of biomolecules. In this article, we briefly review the basic methodological details of QM/MM approaches and discuss their applications to various energy transduction problems in biomolecular machines, such as long-range proton transports, fast electron transfers, and mechanochemical coupling. We highlight the particular importance for these applications of balancing computational efficiency and accuracy. Using several recent examples, we illustrate the value and limitations of QM/MM methodologies for both ground and excited states, as well as strategies for calibrating them in specific applications. We conclude with brief comments on several areas that can benefit from further efforts to make QM/MM analyses more quantitative and applicable to increasingly complex biological problems.

Accelerating Ab Initio Quantum Mechanical and Molecular Mechanical (QM/MM) Molecular Dynamics Simulations with Multiple Time Step Integration and a Recalibrated Semiempirical QM/MM Hamiltonian

Molecular dynamics (MD) simulations employing ab initio quantum mechanical and molecular mechanical (ai-QM/MM) potentials are considered to be the state of the art, but the high computational cost associated with the ai-QM calculations remains a theoretical challenge for their routine application. Here, we present a modified protocol of the multiple time step (MTS) method for accelerating ai-QM/MM MD simulations of condensed-phase reactions. Within a previous MTS protocol [Nam J. Chem. Theory Comput. 2014, 10, 4175], reference forces are evaluated using a low-level (semiempirical QM/MM) Hamiltonian and employed at inner time steps to propagate the nuclear motions. Correction forces, which arise from the force differences between high-level (ai-QM/MM) and low-level Hamiltonians, are applied at outer time steps, where the MTS algorithm allows the time-reversible integration of the correction forces. To increase the outer step size, which is bound by the highest-frequency component in the correction forces, the semiempirical QM Hamiltonian is recalibrated in this work to minimize the magnitude of the correction forces. The remaining high-frequency modes, which are mainly bond stretches involving hydrogen atoms, are then removed from the correction forces. When combined with a Langevin or SIN(R) thermostat, the modified MTS-QM/MM scheme remains robust with an up to 8 (with Langevin) or 10 fs (with SIN(R)) outer time step (with 1 fs inner time steps) for the chorismate mutase system. This leads to an over 5-fold speedup over standard ai-QM/MM simulations, without sacrificing the accuracy in the predicted free energy profile of the reaction.

A simplified charge projection scheme for long-range electrostatics in ab initio QM/MM calculations

In a previous work [Pan et al., Molecules 23, 2500 (2018)], a charge projection scheme was reported, where outer molecular mechanical (MM) charges [>10 Å from the quantum mechanical (QM) region] were projected onto the electrostatic potential (ESP) grid of the QM region to accurately and efficiently capture long-range electrostatics in ab initio QM/MM calculations. Here, a further simplification to the model is proposed, where the outer MM charges are projected onto inner MM atom positions (instead of ESP grid positions). This enables a representation of the long-range MM electrostatic potential via augmentary charges (AC) on inner MM atoms. Combined with the long-range electrostatic correction function from Cisneros et al. [J. Chem. Phys. 143, 044103 (2015)] to smoothly switch between inner and outer MM regions, this new QM/MM-AC electrostatic model yields accurate and continuous ab initio QM/MM electrostatic energies with a 10 Å cutoff between inner and outer MM regions. This model enables efficient QM/MM cluster calculations with a large number of MM atoms as well as QM/MM calculations with periodic boundary conditions.

Benchmarking quantum mechanical methods for calculating reaction energies of reactions catalyzed by enzymes

To assess the accuracy of different quantum mechanical methods for biochemical modeling, the reaction energies of 20 small model reactions (chosen to represent chemical steps catalyzed by commonly studied enzymes) were calculated. The methods tested included several popular Density Functional Theory (DFT) functionals, second-order Møller Plesset perturbation theory (MP2) and its spin-component scaled variant (SCS-MP2), and coupled cluster singles and doubles and perturbative triples (CCSD(T)). Different basis sets were tested. CCSD(T)/aug-cc-pVTZ results for all 20 reactions were used to benchmark the other methods. It was found that MP2 and SCS-MP2 reaction energy calculation results are similar in quality to CCSD(T) (mean absolute error (MAE) of 1.2 and 1.3 kcal mol−1, respectively). MP2 calculations gave a large error in one case, and are more subject to basis set effects, so in general SCS-MP2 calculations are a good choice when CCSD(T) calculations are not feasible. Results with different DFT functionals were of reasonably good quality (MAEs of 2.5–5.1 kcal mol−1), whereas popular semi-empirical methods (AM1, PM3, SCC-DFTB) gave much larger errors (MAEs of 11.6–14.6 kcal mol−1). These results should be useful in guiding methodological choices and assessing the accuracy of QM/MM calculations on enzyme-catalyzed reactions.

Quantum Mechanical/Molecular Mechanical Analysis of the Catalytic Mechanism of Phosphoserine Phosphatase

Phosphoserine phosphatase (PSP), a member of the haloacid dehalogenase (HAD) superfamily that comprises the vast majority of phosphotransferases, is likely a steady-state regulator of the level of d-serine in the brain. The proposed catalytic cycle of PSP consists of a two-step mechanism: formation of a phospho-enzyme intermediate by phosphate transfer to Asp11 and its subsequent hydrolysis. Our combined quantum mechanical/molecular mechanical (QM/MM) calculations of the reaction pathways favour a dissociative mechanism of nucleophilic substitution via a trigonal-planar metaphosphate-like configuration for both steps, associated with proton transfer to the leaving group or from the nucleophile. This proton transfer is facilitated by active site residue Asp13 that acts as both a general base and a general acid. Free energy calculation on the reaction pathways further support the structural role of the enzymatic environment and the active site architecture. The choice of a proper reaction coordinate along which to bias the free energy calculations can be guided by a projection of the canonical reaction coordinate obtained from a chain-of-state optimisation onto important internal coordinates.

Rationalising pKa shifts in Bacillus circulans xylanase with computational studies

Bacillus circulans xylanase (BcX), a family 11 glycoside hydrolase, catalyses the hydrolysis of xylose polymers with a net retention of stereochemistry. Glu172 in BcX is believed to act as a general acid by protonating the aglycone during glycosylation, and then as a general base to facilitate the deglycosylation step. The key to the dual role of this general acid/base lies in its protonation states, which depend on its intrinsic pK_a value and the specific environment which it resides within. To fully understand the detailed molecular features in BcX to establish the dual role of Glu172, we present a combined study based on both atomistic simulations and empirical models to calculate pK_a shifts for the general acid/base Glu172 in BcX at different functional states. Its pK_a values and those of nearby residues, obtained based on QM/MM free energy calculations, MCCE and PROPKA, show a good agreement with available experimental data. Additionally, our study provides additional insights into the effects of structural and electrostatic perturbations caused by mutations and chemical modifications, suggesting that the local solvation environment and mutagenesis of the residues adjacent to Glu172 establish its dual role during hydrolysis. The strengths and limitations of various methods for calculating pK_as and pK_a shifts have also been discussed.

Calculations of the absolute binding free energies for Ralstonia solanacearum lectins bound with methyl-α-l-fucoside at molecular mechanical and quantum mechanical/molecular mechanical levels

In this work, both a molecular mechanical (MM) method and a hybrid quantum mechanical/molecular mechanical (QM/MM) method have been applied in the study of the binding affinities of methyl-α-l-fucoside to Ralstonia solanacearum lectins.

Multiscale Quantum Mechanics/Molecular Mechanics Simulations with Neural Networks

Molecular dynamics simulation with multiscale quantum mechanics/molecular mechanics (QM/MM) methods is a very powerful tool for understanding the mechanism of chemical and biological processes in solution or enzymes. However, its computational cost can be too high for many biochemical systems because of the large number of ab initio QM calculations. Semiempirical QM/MM simulations have much higher efficiency. Its accuracy can be improved with a correction to reach the ab initio QM/MM level. The computational cost on the ab initio calculation for the correction determines the efficiency. In this paper we developed a neural network method for QM/MM calculation as an extension of the neural-network representation reported by Behler and Parrinello. With this approach, the potential energy of any configuration along the reaction path for a given QM/MM system can be predicted at the ab initio QM/MM level based on the semiempirical QM/MM simulations. We further applied this method to three reactions in water to calculate the free energy changes. The free-energy profile obtained from the semiempirical QM/MM simulation is corrected to the ab initio QM/MM level with the potential energies predicted with the constructed neural network. The results are in excellent accordance with the reference data that are obtained from the ab initio QM/MM molecular dynamics simulation or corrected with direct ab initio QM/MM potential energies. Compared with the correction using direct ab initio QM/MM potential energies, our method shows a speed-up of 1 or 2 orders of magnitude. It demonstrates that the neural network method combined with the semiempirical QM/MM calculation can be an efficient and reliable strategy for chemical reaction simulations.

Charge-dependent many-body exchange and dispersion interactions in combined QM/MM simulations

Accurate modeling of the molecular environment is critical in condensed phase simulations of chemical reactions. Conventional quantum mechanical/molecular mechanical (QM/MM) simulations traditionally model non-electrostatic non-bonded interactions through an empirical Lennard-Jones (LJ) potential which, in violation of intuitive chemical principles, is bereft of any explicit coupling to an atom's local electronic structure. This oversight results in a model whereby short-ranged exchange-repulsion and long-ranged dispersion interactions are invariant to changes in the local atomic charge, leading to accuracy limitations for chemical reactions where significant atomic charge transfer can occur along the reaction coordinate. The present work presents a variational, charge-dependent exchange-repulsion and dispersion model, referred to as the charge-dependent exchange and dispersion (QXD) model, for hybrid QM/MM simulations. Analytic expressions for the energy and gradients are provided, as well as a description of the integration of the model into existing QM/MM frameworks, allowing QXD to replace traditional LJ interactions in simulations of reactive condensed phase systems. After initial validation against QM data, the method is demonstrated by capturing the solvation free energies of a series of small, chlorine-containing compounds that have varying charge on the chlorine atom. The model is further tested on the SN2 attack of a chloride anion on methylchloride. Results suggest that the QXD model, unlike the traditional LJ model, is able to simultaneously obtain accurate solvation free energies for a range of compounds while at the same time closely reproducing the experimental reaction free energy barrier. The QXD interaction model allows explicit coupling of atomic charge with many-body exchange and dispersion interactions that are related to atomic size and provides a more accurate and robust representation of non-electrostatic non-bonded QM/MM interactions.

QM/MM free energy simulations: recent progress and challenges

Due to the higher computational cost relative to pure molecular mechanical (MM) simulations, hybrid quantum mechanical/molecular mechanical (QM/MM) free energy simulations particularly require a careful consideration of balancing computational cost and accuracy. Here we review several recent developments in free energy methods most relevant to QM/MM simulations and discuss several topics motivated by these developments using simple but informative examples that involve processes in water. For chemical reactions, we highlight the value of invoking enhanced sampling technique (e.g., replica-exchange) in umbrella sampling calculations and the value of including collective environmental variables (e.g., hydration level) in metadynamics simulations; we also illustrate the sensitivity of string calculations, especially free energy along the path, to various parameters in the computation. Alchemical free energy simulations with a specific thermodynamic cycle are used to probe the effect of including the first solvation shell into the QM region when computing solvation free energies. For cases where high-level QM/MM potential functions are needed, we analyze two different approaches: the QM/MM-MFEP method of Yang and co-workers and perturbative correction to low-level QM/MM free energy results. For the examples analyzed here, both approaches seem productive although care needs to be exercised when analyzing the perturbative corrections.

QM/MM Analysis of Transition States and Transition State Analogues in Metalloenzymes

Enzymology is approaching an era where many problems can benefit from computational studies. While ample challenges remain in quantitatively predicting behavior for many enzyme systems, the insights that often come from computations are an important asset for the enzymology community. Here we provide a primer for enzymologists on the types of calculations that are most useful for mechanistic problems in enzymology. In particular, we emphasize the integration of models that range from small active-site motifs to fully solvated enzyme systems for cross-validation and dissection of specific contributions from the enzyme environment. We then use a case study of the enzyme alkaline phosphatase to illustrate specific application of the methods. The case study involves examination of the binding modes of putative transition state analogues (tungstate and vanadate) to the enzyme. The computations predict covalent binding of these ions to the enzymatic nucleophile and that they adopt the trigonal bipyramidal geometry of the expected transition state. By comparing these structures with transition states found through free energy simulations, we assess the degree to which the transition state analogues mimic the true transition states. Technical issues worth treating with care as well as several remaining challenges to quantitative analysis of metalloenzymes are also highlighted during the discussion.

Substrate-Assisted Catalysis in the Reaction Catalyzed by Salicylic Acid Binding Protein 2 (SABP2), a Potential Mechanism of Substrate Discrimination for Some Promiscuous Enzymes

Although one of an enzyme's hallmarks is the high specificity for their natural substrates, substrate promiscuity has been reported more frequently. It is known that promiscuous enzymes generally show different catalytic efficiencies to different substrates, but our understanding of the origin of such differences is still lacking. Here we report the results of quantum mechanical/molecular mechanical simulations and an experimental study of salicylic acid binding protein 2 (SABP2). SABP2 has promiscuous esterase activity toward a series of substrates but shows a high activity toward its natural substrate, methyl salicylate (MeSA). Our results demonstrate that this enzyme may use substrate-assisted catalysis involving the hydroxyl group from MeSA to enhance the activity and achieve substrate discrimination.

Molecular simulation of water and hydration effects in different environments: challenges and developments for DFTB based models

We discuss the description of water and hydration effects that employs an approximate density functional theory, DFTB3, in either a full QM or QM/MM framework. The goal is to explore, with the current formulation of DFTB3, the performance of this method for treating water in different chemical environments, the magnitude and nature of changes required to improve its performance, and factors that dictate its applicability to reactions in the condensed phase in a QM/MM framework. A relatively minor change (on the scale of kBT) in the O-H repulsive potential is observed to substantially improve the structural properties of bulk water under ambient conditions; modest improvements are also seen in dynamic properties of bulk water. This simple change also improves the description of protonated water clusters, a solvated proton, and to a more limited degree, a solvated hydroxide. By comparing results from DFTB3 models that differ in the description of water, we confirm that proton transfer energetics are adequately described by the standard DFTB3/3OB model for meaningful mechanistic analyses. For QM/MM applications, a robust parametrization of QM-MM interactions requires an explicit consideration of condensed phase properties, for which an efficient sampling technique was developed recently and is reviewed here. The discussions help make clear the value and limitations of DFTB3 based simulations, as well as the developments needed to further improve the accuracy and transferability of the methodology.

Density functional tight binding: values of semi-empirical methods in an ab initio era

Semi-empirical (SE) methods are derived from Hartree-Fock (HF) or Density Functional Theory (DFT) by neglect and approximation of electronic integrals. Thereby, parameters are introduced which have to be determined from reference calculations and/or by fitting to available experimental data. This leads to computational methods that are about 2-3 orders of magnitude faster than the standard HF/DFT methods using medium sized basis sets while being about 3 orders of magnitude slower than empirical force field methods (Molecular Mechanics: MM). Therefore, SE methods are most appropriate for a specific range of applications. These include the study of systems that contain a large number of atoms and therefore being too large for ab initio or DFT methods and also problems where dynamic or entropic effects are particularly important. In the latter case, the errors made by considering a very limited number of molecular structures or neglecting entropic contributions can be much larger than the accuracy lost due to the use of SE methods. Another area where SE methods are attractive concerns the analysis of systems for which reliable MM models are not readily available. Therefore, even in an era when rapid progress is being made in ab initio methods, there is considerable interest in further developing SE methods. We illustrate this point by focusing on the discussion of recent development and application of the Density Functional Tight Binding method.

Comparison of ab Initio, DFT, and Semiempirical QM/MM Approaches for Description of Catalytic Mechanism of Hairpin Ribozyme

We have analyzed the capability of state-of-the-art multiscale computational approaches to provide atomic-resolution electronic structure insights into possible catalytic scenarios of the hairpin ribozyme by evaluating potential and free energy surfaces of the reactions by various hybrid QM/MM methods. The hairpin ribozyme is a unique catalytic RNA that achieves rate acceleration similar to other small self-cleaving ribozymes but without direct metal ion participation. Guanine 8 (G8) and adenine 38 (A38) have been identified as the catalytically essential nucleobases. However, their exact catalytic roles are still being investigated. In line with the available experimental data, we considered two reaction scenarios involving protonated A38H(+) as a general acid which is further assisted by either canonical G8 or deprotonated G8(-) forms. We used the spin-component scaled Møller-Plesset (SCS-MP2) method at the complete basis set limit as the reference method. The semiempirical AM1/d-PhoT and SCC-DFTBPR methods provided acceptable activation barriers with respect to the SCS-MP2 data but predicted significantly different reaction pathways. DFT functionals (BLYP and MPW1K) yielded the same reaction pathway as the SCS-MP2 method. The activation barriers were slightly underestimated by the GGA BLYP functional, although with accuracy comparable to the semiempirical methods. The SCS-MP2 method and hybrid MPW1K functional gave activation barriers that were closest to those derived from experimentally measured rate constants.

Benchmark Study of the SCC-DFTB Approach for a Biomolecular Proton Channel

The self-consistent charge density functional tight binding (SCC-DFTB) method has been increasingly applied to study proton transport (PT) in biological environments. However, recent studies revealing some significant limitations of SCC-DFTB for proton and hydroxide solvation and transport in bulk aqueous systems call into question its accuracy for simulating PT in biological systems. The current work benchmarks the SCC-DFTB/MM method against more accurate DFT/MM by simulating PT in a synthetic leucine-serine channel (LS2), which emulates the structure and function of biomolecular proton channels. It is observed that SCC-DFTB/MM produces over-coordinated and less structured pore water, an over-coordinated excess proton, weak hydrogen bonds around the excess proton charge defect and qualitatively different PT dynamics. Similar issues are demonstrated for PT in a carbon nanotube, indicating that the inaccuracies found for SCC-DFTB are not due to the point charge based QM/MM electrostatic coupling scheme, but rather to the approximations of the semiempirical method itself. The results presented in this work highlight the limitations of the present form of the SCC-DFTB/MM approach for simulating PT processes in biological protein or channel-like environments, while providing benchmark results that may lead to an improvement of the underlying method.

Multiscale Multiphysics and Multidomain Models I: Basic Theory

This work extends our earlier two-domain formulation of a differential geometry based multiscale paradigm into a multidomain theory, which endows us the ability to simultaneously accommodate multiphysical descriptions of aqueous chemical, physical and biological systems, such as fuel cells, solar cells, nanofluidics, ion channels, viruses, RNA polymerases, molecular motors and large macromolecular complexes. The essential idea is to make use of the differential geometry theory of surfaces as a natural means to geometrically separate the macroscopic domain of solvent from the microscopic domain of solute, and dynamically couple continuum and discrete descriptions. Our main strategy is to construct energy functionals to put on an equal footing of multiphysics, including polar (i.e., electrostatic) solvation, nonpolar solvation, chemical potential, quantum mechanics, fluid mechanics, molecular mechanics, coarse grained dynamics and elastic dynamics. The variational principle is applied to the energy functionals to derive desirable governing equations, such as multidomain Laplace-Beltrami (LB) equations for macromolecular morphologies, multidomain Poisson-Boltzmann (PB) equation or Poisson equation for electrostatic potential, generalized Nernst-Planck (NP) equations for the dynamics of charged solvent species, generalized Navier-Stokes (NS) equation for fluid dynamics, generalized Newton's equations for molecular dynamics (MD) or coarse-grained dynamics and equation of motion for elastic dynamics. Unlike the classical PB equation, our PB equation is an integral-differential equation due to solvent-solute interactions. To illustrate the proposed formalism, we have explicitly constructed three models, a multidomain solvation model, a multidomain charge transport model and a multidomain chemo-electro-fluid-MD-elastic model. Each solute domain is equipped with distinct surface tension, pressure, dielectric function, and charge density distribution. In addition to long-range Coulombic interactions, various non-electrostatic solvent-solute interactions are considered in the present modeling. We demonstrate the consistency between the non-equilibrium charge transport model and the equilibrium solvation model by showing the systematical reduction of the former to the latter at equilibrium. This paper also offers a brief review of the field.

A hybrid approach to simulation of electron transfer in complex molecular systems

Electron transfer (ET) reactions in biomolecular systems represent an important class of processes at the interface of physics, chemistry and biology. The theoretical description of these reactions constitutes a huge challenge because extensive systems require a quantum-mechanical treatment and a broad range of time scales are involved. Thus, only small model systems may be investigated with the modern density functional theory techniques combined with non-adiabatic dynamics algorithms. On the other hand, model calculations based on Marcus's seminal theory describe the ET involving several assumptions that may not always be met. We review a multi-scale method that combines a non-adiabatic propagation scheme and a linear scaling quantum-chemical method with a molecular mechanics force field in such a way that an unbiased description of the dynamics of excess electron is achieved and the number of degrees of freedom is reduced effectively at the same time. ET reactions taking nanoseconds in systems with hundreds of quantum atoms can be simulated, bridging the gap between non-adiabatic ab initio simulations and model approaches such as the Marcus theory. A major recent application is hole transfer in DNA, which represents an archetypal ET reaction in a polarizable medium. Ongoing work focuses on hole transfer in proteins, peptides and organic semi-conductors.

Modeling catalytic promiscuity in the alkaline phosphatase superfamily

In recent years, it has become increasingly clear that promiscuity plays a key role in the evolution of new enzyme function. This finding has helped to elucidate fundamental aspects of molecular evolution. While there has been extensive experimental work on enzyme promiscuity, computational modeling of the chemical details of such promiscuity has traditionally fallen behind the advances in experimental studies, not least due to the nearly prohibitive computational cost involved in examining multiple substrates with multiple potential mechanisms and binding modes in atomic detail with a reasonable degree of accuracy. However, recent advances in both computational methodologies and power have allowed us to reach a stage in the field where we can start to overcome this problem, and molecular simulations can now provide accurate and efficient descriptions of complex biological systems with substantially less computational cost. This has led to significant advances in our understanding of enzyme function and evolution in a broader sense. Here, we will discuss currently available computational approaches that can allow us to probe the underlying molecular basis for enzyme specificity and selectivity, discussing the inherent strengths and weaknesses of each approach. As a case study, we will discuss recent computational work on different members of the alkaline phosphatase superfamily (AP) using a range of different approaches, showing the complementary insights they have provided. We have selected this particular superfamily, as it poses a number of significant challenges for theory, ranging from the complexity of the actual reaction mechanisms involved to the reliable modeling of the catalytic metal centers, as well as the very large system sizes. We will demonstrate that, through current advances in methodologies, computational tools can provide significant insight into the molecular basis for catalytic promiscuity, and, therefore, in turn, the mechanisms of protein functional evolution.

Stabilization of different types of transition states in a single enzyme active site: QM/MM analysis of enzymes in the alkaline phosphatase superfamily

The first step for the hydrolysis of a phosphate monoester (pNPP(2-)) in enzymes of the alkaline phosphatase (AP) superfamily, R166S AP and wild-type NPP, is studied using QM/MM simulations based on an approximate density functional theory (SCC-DFTBPR) and a recently introduced QM/MM interaction Hamiltonian. The calculations suggest that similar loose transition states are involved in both enzymes, despite the fact that phosphate monoesters are the cognate substrates for AP but promiscuous substrates for NPP. The computed loose transition states are clearly different from the more synchronous ones previously calculated for diester reactions in the same AP enzymes. Therefore, our results explicitly support the proposal that AP enzymes are able to recognize and stabilize different types of transition states in a single active site. Analysis of the structural features of computed transition states indicates that the plastic nature of the bimetallic site plays a minor role in accommodating multiple types of transition states and that the high degree of solvent accessibility of the AP active site also contributes to its ability to stabilize diverse transition-state structures without the need of causing large structural distortions of the bimetallic motif. The binding mode of the leaving group in the transition state highlights that vanadate may not always be an ideal transition state analog for loose phosphoryl transfer transition states.

Efficient algorithms for the simulation of non-adiabatic electron transfer in complex molecular systems: application to DNA

In this work, a fragment-orbital density functional theory-based method is combined with two different non-adiabatic schemes for the propagation of the electronic degrees of freedom. This allows us to perform unbiased simulations of electron transfer processes in complex media, and the computational scheme is applied to the transfer of a hole in solvated DNA. It turns out that the mean-field approach, where the wave function of the hole is driven into a superposition of adiabatic states, leads to over-delocalization of the hole charge. This problem is avoided using a surface hopping scheme, resulting in a smaller rate of hole transfer. The method is highly efficient due to the on-the-fly computation of the coarse-grained DFT Hamiltonian for the nucleobases, which is coupled to the environment using a QM/MM approach. The computational efficiency and partial parallel character of the methodology make it possible to simulate electron transfer in systems of relevant biochemical size on a nanosecond time scale. Since standard non-polarizable force fields are applied in the molecular-mechanics part of the calculation, a simple scaling scheme was introduced into the electrostatic potential in order to simulate the effect of electronic polarization. It is shown that electronic polarization has an important effect on the features of charge transfer. The methodology is applied to two kinds of DNA sequences, illustrating the features of transfer along a flat energy landscape as well as over an energy barrier. The performance and relative merit of the mean-field scheme and the surface hopping for this application are discussed.