Combined Quantum Mechanical and Molecular

Combined Quantum Mechanical and Molecular

(Parte 1 de 4)

Combined Quantum Mechanical and Molecular

Mechanical Methods for Calculating Potential Energy Surfaces: Tuned and Balanced Redistributed-Charge Algorithm

Bo Wang and Donald G. Truhlar*

Department of Chemistry and Supercomputing Institute, UniVersity of Minnesota, 207 Pleasant Street South East, Minneapolis, Minnesota 55455-0431

Received July 15, 2009

Abstract: The combined quantum mechanical and molecular mechanical (QM/M) method is one of the most powerful approaches for including correlation and polarization effects in simulations of large and complex systems, and the present article is concerned with the systematics of treating a QM/M boundary that passes through a covalent bond, especially a polar covalent bond. In this study, we develop a new algorithm to treat such boundaries; the new method is called the balanced redistributed charge (balanced RC or BRC) scheme with a tuned fluorinelink atom. The M point charge on the M boundaryatom is modifiedto conserve the total charge of the entire system, and the modified charge is redistributed to the midpoints of the bonds between an M boundary atom and its neighboring M atoms. A pseudopotential is added to the fluorine link atom to reproduce the partial charge of the uncapped portion of the QM subsystem. We select proton affinities as the property used to validate the new method because the energy change associated with the addition of an entire charge (proton) to the QM system is very sensitive to the treatment of electrostatics at the boundary; we apply the new method to calculate proton affinities of 25 molecules with 13 different kinds of bonds being cut. The average proton affinity in the test set is 373 kcal/mol, and the test set provides a more challenging test than those usually used for testing QM/M methods. For this challenging test set, common unbalanced schemes give a mean unsigned error (MUE) of 15-21 kcal/mol for H link atoms or 16-24 kcal/mol for F link atoms, much larger than the 5 kcal/mol obtained by simply omitting the M region with either kind of link atom. Balancing the charges reduces the error to 5-7 kcal/mol for H link atoms and 4-6 kcal/mol for F link atoms. Balancing the charges and also tuning an F link atom lowers the MUE to 1.3-4 kcal/mol, with the best result for the balanced RC scheme. We conclude that properly tuning the link atom and correctly treating the point charges near the QM/M boundary significantly improves the accuracy of the calculated proton affinities.

1. Introduction

The application of quantum chemistry to large and complex systems is one of the most challenging areas of current computational chemistry and also one that is seeing the most progress.1 An important tool for such applications is the combined quantum mechanical/molecular mechanical (QM/ M) method for calculating potential energy surfaces and interatomic forces; the reader is directed to several reviews and overviews for background information.2-23

A stubborn issue in QM/M calculations is the treatment of the boundary between the QM and M regions when it

* Corresponding author phone: (612) 624-7555; fax: (612) 624- 9390; e-mail:

J. Chem. Theory Comput. X, x, 0 A

10.1021/ct900366m C: $40.75 X American Chemical Society passes through a bond, which is practically unavoidable in the treatment of many solids, polymers, and complex systems. In general the QM region is capped to saturate dangling valences caused by the cut. Three different kinds of methods have been proposed to deal with capping the QM boundary atom. The first one is the link atom approach (LA).24,25 The dangling bond of the QM region is capped with an additional atom (usually a hydrogen atom) and the QM calculations are performed on this capped system. The second method is localized orbitals.26-28 The dangling bond is saturated by orbitals rather than by an atom. Examples of this approach are the local self-consistent field (LSCF) method26 and the generalized hybrid orbitals (GHO) scheme.27,28 The third kind of methodinvolvesa pseudobond or an effective core potential (ECP). In this approach, a parametrized atom, modified to mimic the behavior of the original M boundary atoms or groups, is used to cap the QM system; examples of this approach are tuned capping atoms,29-31 adjustedconnectionatoms,32 a pseudobond,3-35 an effectivegrouppotential,36 a quantumcappedpotential,37-39 and a variationally optimized effective atom-centered potential.40 This third class of methods may be considered to be a second-generation link-atom method in which the link atom is optimized or tuned.

Though much progresshas been made, there are still many problems in the treatment of QM-M boundaries that pass through a bond. Most attention has been devoted to the cutting of C-C bonds, especially for modeling enzymatic bindingand reactions,but some proceduresare more general. The methods that have been developedexhibit a wide variety of differences in the precise way in which they have been implemented.

formulated. It should be generalParticular procedures
for particular moleculesshould be avoided.”41 If tests of

Pople has emphasized the importance of theoretical models, where a theoretical model is “an approximate but well-defined mathematical procedure for simulation...T he approximate mathematical treatment must be precisely the model against a broad data set are successful, the model is said to be validated. The goal of this article is to develop and validatea new method,in the spiritof a theoreticalmodel chemistry, for the treatment of a boundary between bonded atoms in QM/M simulations. It is precisely defined in a general way applicable to all systems and all kinds of single bonds, and it is tested against a data set of 25 systems in which13 differentkindsof bondsare cut, in particular(where the atom listed first is in the QM subsystem, and the one listed second is in the M subsystem): C-C, N-C, O-C, S-C, C-N, O-N, C-O, Al-O, Si-O, C-Si, O-Si, C-S, and S-S.

2. Methods

Our group has developed redistributed charge (RC) and redistributed charge and dipole (RCD) methods to treat the charges near a QM/M boundary that passes through a bond.42 These methods give good results even when large charges are present near the boundary. In the current work, we improve the RC and RCD methods by adding two new elements, a charge balancing step and a tuned link atom. In particular, the redistributed charges are used to conserve the charge of the entire system, and a tuned fluorine atom is used to saturate the free valence of the QM region and to reproduce the partial charge of the uncapped portion of the cappedQM subsystem.The improvedmethodis used to treat polar bonds betweenthe QM and M subsystemswith large partial atomic charges near the boundary.In order to describe the algorithm, we label the atoms according to “tiers”. The definition is the same as what is used in previous work;4,42 in particular, the M boundary atoms are denoted as M1 atoms, and the M atoms directly bonded to M1 atoms are denotedas M2 atoms.M3 atomsare the third-tierMM atoms. The QM boundary atoms are denoted as Q1 atoms. The QM atoms directly bonded to Q1 atoms are labeled Q2 atoms. Q3 atoms are those bonded to Q2 atoms and so forth for Q4, Q5, etc. The QM region is also called the primary subsystem (PS) in this study. The sum of all QM atoms and M atoms before the cutting and capping is called the originalentire system.The sum of the capped QM subsystem and the whole M subsystem after the charge redistribution is called the QM/M entire system.

In the QM/M calculations, we use an additive QM/M scheme to define the total energy of the system:23 where EQM is the quantum mechanical energy of the QM region, EM is the molecular mechanical energy of the M region, and EQM/M accounts for the interaction energy betweenthe QM and the M regions.EQM/Mis decomposed into three terms; EelQM/M represents electrostatic interactions,

EvdWQM/M represents van der Waals interactions, and Eval QM/M represents valence interactions. In this study, we will con- centrate on the electrostaticcoupling term EelQM/M, which is the most technicallyinvolved term. The EvdWQM/M and Eval QM/M terms will cancel out in the present work because we study fixed-geometryprotonaffinitiesto isolatetheelectrostaticterms, but these other QM/M terms will be studiedlater when we considerQM/M geometryoptimization. 2.A. Treatment of Boundary Charge. It has been found that it is important to conserve the total charge of the QM/ M entire system in QM/M calculations,43 that is, the sum of the M partial atomic charges of the M region and the QM charge of the capped QM region should equal the total charge of the original entire system, as shown in eq 3:

However, when the original entire system is divided into QM and M regions, the sum of M charges of the M region does not necessarily equal zero or an integer. If M charges are not modified, the total charge of the QM/M entire system is not conserved. Several workers have recognized that this causes inaccuracies and have suggested various methods to remedy this.3,43-45 Sherwood et al.4 adjusted the charge on an M1 atom to conserve the total charge of the QM/M entire system,

B J. Chem. Theory Comput., Vol. x, No. x, X Wang and Truhlar

and they redistributed the adjusted charge on the M1 atom evenly to M2 atoms; point dipoles were added at the M2 atoms to compensate the changes in the M1-M2 bond dipoles due to the movement of the charges. Zhang et al.3 zeroed the charges on all M atoms that are in the same group as the M1 atom. Das et al.45 used a double link atom approach combined with delocalized Gaussian M charges. Walker et al.43 added the charge difference to the nearest M2 atom or evenly to all the M atoms except the M1 atom.

In the previousRC scheme,42 the charge on each M1 atom is redistributedto the midpointsof M1-M2 bonds.However, the total charge of the QM/M entire system is not conserved when the sum of M charges of the M region is not zero or an integer. In the balanced RC scheme, introduced here, we first adjust the charge on the M1 atom so that

where q0 is the modified M1 charge, {qi} are the M point charges of other M atoms (except M1), qQM is the charge of the QM region (that is, of the capped QM subsystem), and qtotal is the charge of the original entire system. This charge balancing step conserves the charge of the QM/M entire system.

Then the balanced RC scheme redistributes the charge q0 evenly to the midpoints of all M1-M2 bonds, with each bond midpoint obtaining a charge qRC ) (q0/n), where n is the number of M1-M2 bonds. For the balancedredistributed charge and dipole (balanced RCD) method, we double the redistributed charges and adjust the charges qM2RCD on M2 atoms to conserve the total charge of the QM/M entire system, as shown in eqs 5 and 6:

These two schemes are illustrated in Figure 1.

In this study, we compare balanced RC and balanced

RCD to other methods that differ in how the redistributed charges are handled, e.g., to what location are they redistributed. These methods include: balanced straight electrostatic embedding (SEE), balanced RC2, Amber-1,43 balanced RC3, Amber-2,43 and balanced Shift.4 Amber-1 and Amber-2 are the options called adjust_q ) 1 and adjust_q ) 2i n AMBER 10. The distinction between these methods is in the position of the redistributed charges and whether the dipoles of the M1-M2 bonds are corrected.

In balanced SEE, the charge on the M1 atom is set to q0, and it is not moved. In balanced RC2, we distribute q0 evenly to all M2 atoms. In balanced RC3, we distribute q0 evenly to all M2 and M3 atoms. In Amber-1, we move q0 to the nearest M2 atom. In Amber-2, we distribute q0 evenly to all M atoms, except the M1 atom. (Note that

Amber-2 is the default option in revision 10 of AMBER, 46 whereas Amber-1 can be selected in AMBER 10 by specifying adjust_q ) 1.) In balanced Shift, the redistributed charges are placed at M2 atoms, and dipoles are added around M2 atoms to compensate the movement of the charges. A summary of these charge schemes is shown in Table 1.

We call the methods in Table 1 balanced methods because they all conserve the total charge of the QM/ M entire system. Five unbalanced methods, in which the total charge of the QM/M entire system is not necessarily conserved, are also tested, including SEE, Z1, Z2, Z3, and RC.42 SEE is straight electronic embedding that makes no change of the charges of M boundary atoms, Z1 denotes that the charge of the M1 atom is zeroed (this can be chosen by specifying adjust_q ) 0i n AMBER 10, and it is the default method in CHARMM47), Z2 denotes that the charges of M1 and M2 atoms are zeroed (Z2 is the default scheme in both Gaussion 0348 and Gaussian 0949), and Z3 denotes that all the charges of all M1, M2, and M3 atoms are zeroed. RC denotes that the charge on the M1 atom is redistributed to the midpoints of M1-M2 bonds without the balancing step. Balanced methods and unbalanced methods are compared to test the importance of conserving the charge of the QM/M entire system. To make a comparison, we also carry out calculations on the capped primary system (CPS), in which the whole M region is substituted by the link atom. 2.B. Link Atom. Another issue in the boundary treatment is the choice of the link atom. A hydrogen atom can be used as the link atom when a C-C bond is cut. However, a Q1-H bond may be a poor model for the cut bond when the M1 atom is electronegative, such as in a Si-Oo rC -O bond. Therefore, we use a tuned capping atom as the link atom to mimic a cut polar bond and to reproduce the electronic structure of the QM subsystem. Redondo et al.29 used a tuned hydrogen atom to replace a silicon atom. Koga et al.30 added a shift operator on the hydrogen atom to reproduce the effect of the substitution. Zhang et al.3 and Nasluzov et al.31 used tuned fluorine atoms and derived pseudopotentials for carbon and oxygen boundary atoms. Here, we provide a more general rule to tune a capping atom for boundary atoms. The capping atom is always a tuned F atom. We first replace the 1s2 core by a conventional pseudopotential U, and then a

Figure 1. QM/M boundary treatments in (a) the balanced RC scheme and (b) the balanced RCD scheme.

qM2RCD ) qM2 - qRC (6)

Combined QM/M Methods J. Chem. Theory Comput., Vol. x, No. x, X C

ventional pseudopotential used here is the CRENBL effective core potential (ECP) for a fluorine atom developed by Pacios and Christiansen. 50 The form of this potential is


r is the distance of an electron from the capping nucleus, and |lm〉 is a spherical harmonic. The parameters for this pseudopotential are listed in Table 2. The form of U0′(r)i s where C and r0 are parameters. The basis set used for the tuned F atom is the same as for a conventional F atom. For example, if the QM subsystem is treated by the 6-31G* basis set, then the tuned F atom has the 6-31G* basis set of a conventionalF atom. To find an appropriatepseudopotential, we set r0 equal to 1 bohr and tune the parameter C of the pseudopotential.

The next key decision is how to choose C. In order to reproduce the electronic structure of the QM subsystem, we require that the total charge of the uncapped portion of the QM subsystem in a QM/M calculation is equal to the total charge of the same subsystem in a QM calculation of the original entire system or, in practice, of a system that mimics the original entire system better than the capped QM subsystem does (see below for more details of this large system). Mulliken charges were used as the indicator. Because Mulliken charges become unphysical when large basis sets are used, we used small basis sets without diffuse functions for this tuning step, in particular, 6-31G* when M1 is from the second period (Li through F) and STO-3G otherwise. Since the STO-3G basis set is defined for the entire periodic table, the tuning step is well-defined for the entire periodic table.

(Parte 1 de 4)