Energy Decomposition Analysis in Solution Based on the Fragment

ARTICLE pubs.acs.org/JPCA

Energy Decomposition Analysis in Solution Based on the Fragment Molecular Orbital Method Dmitri G. Fedorov*,† and Kazuo Kitaura†,‡ † ‡

NRI, National Institute of Advanced Industrial Science and Technology (AIST), Central 2, Umezono 1-1-1, Tsukuba, 305-8568, Japan Graduate School of Pharmaceutical Sciences, Kyoto University, Sakyo-ku, Kyoto 606-8501, Japan ABSTRACT: We develop the pair interaction energy decomposition analysis (PIEDA) in solution by combining the fragment molecular orbital (FMO) method with the polarizable continuum model (PCM). The solvent screening of the electrostatic interaction and the desolvation penalty in complex formation are described by this approach from ab initio calculations of fragments and their pairs. The applications to the complex of solvated sodium and chlorine ions, as well as to lysine and aspartic acid, show how the analysis helps reveal the physical picture. The PIEDA/PCM method is also applied to a small protein chignolin (PDB: 1UAO), and the solvent screening of the pair interactions is discussed.

numerous fragment-based approaches1122 recently reviewed,23 as well as other schemes.24,25 In this work we employ the fragment molecular orbital (FMO) method proposed by Kitaura et al.2628 combined with the polarizable continuum model (PCM)2 in the FMO/PCM approach.29,30 Geometry optimizations3134 and molecular dynamics35,36 can be performed with FMO, which has also been applied to excited states,37,38 radicals,39 proteinligand binding,4042 protein folding,43 and drug design.44,45 The distinctive feature of FMO is the inclusion of the electrostatic (ES) potential (ESP) from the whole system into each individual fragment calculation and in using the many-body expansion46 to account for the interfragment interactions. In addition to PCM, the PoissonBoltzmann47 and effective fragment potential48 solvation models are also developed for FMO. Among various continuum solvation models, PCM is one of the most widely used, and its recent development includes the contributions by Li and co-workers4951 as well as other groups.5254 The energy decomposition analyses (EDA) are an active area of research,55 and by using some models, the solutesolvent charge transfer can be discussed.56 FMO and other fragment methods operate with fragments, facilitating their use in fragment-based drug design57 by providing pair interaction energies (PIE) between fragments (sometimes called interfragment interaction energies, IFIE). These energies deliver the means to elucidate the mechanism of binding and ligand recognition in large systems and are used extensively in many applications.40,42,44,58,59 PIEs can be decomposed into the

1. INTRODUCTION Many chemical processes occur in solution. The description of solvent in quantum mechanics (QM) has presented very considerable difficulties, driven by its temperature-dependent macroscopic nature. There are two groups of solvation methods, explicit solvation,1 where the solvent is treated as discrete molecules, and various continuum2 models; these groups are bridged by the statistical3 methods that retain some atomistic details in the averaged way, hybrid solvation4,5 (i.e., adding a few explicit solvent molecules to continuum models), and approaches relying on molecular dynamics (MD) to generate an averaged solvent electrostatic field.6 The explicit solvation is ultimately the best approach, which, however, requires a long time MD, in order to describe the realistic solvent rather than a frozen molecular cluster, finding whose global minimum is a major problem itself. On the other hand, continuum models from the beginning deal with macroscopic solvent and offer a much faster shortcut to describing the solvent, at the cost of relying on parametrization and losing atomistic details, which in some cases can be important. Typically, the solventsolute charge transfer is not explicitly built in continuum solvation models, and in particular for properties such as the electronic excitations that can be a problem. One faces the dilemma in practical applications, whether the additional cost of explicit methods justifies the improvement of the results, which by no means is guaranteed if the configurational space is not properly sampled. An immensely important group of processes in solution comes from biology. There, in addition to the problem of describing the solvent,7 the size of biochemical systems is typically very large, and the QM models are difficult to apply because of their large cost. There are QM methods applicable to the whole system,810 r 2011 American Chemical Society

Received: October 5, 2011 Revised: November 18, 2011 Published: November 18, 2011 704

dx.doi.org/10.1021/jp209579w | J. Phys. Chem. A 2012, 116, 704–719

The Journal of Physical Chemistry A

ARTICLE

electrostatic, exchange-repulsion, dispersion, and charge transfer, plus mixed terms in the PIE decomposition analysis (PIEDA),60 or complemented with basis set superposition error (BSSE) corrections.61,62 However, these pair interaction energies have so far been computed using QM energies and, thus, do not reflect the solvent screening.63 Perhaps, surprisingly, the same has been true for PIEs coming from FMO/PCM as well, despite the presence of solvent. Such bare QM-based PIEs have very large values for the electrostatic interaction between two charged fragments, often exceeding 100 kcal/mol, and this makes the results appear unbiophysical. It is shown in this work that in FMO/PCM the solvent-induced screening is accounted for, but it is hidden in the bulk total values of the solutesolvent interactions, and the screening contributions to individual pair interactions are not extracted. On the other hand, the total FMO/PCM binding energies have provided reasonable results.41,42 The purpose of this work is to define the explicit solvent corrections to individual pair interactions so that FMO can become much more useful in biochemical applications.

where 1 X X X ΔEes X ¼ ½TrðD W e Þ þ W N 3 q 2 ðW Xe Þμν ¼ ¼

ðW XN Þi ¼

2.1. Basic FMO/PCM Equations. The total energy in the twobody FMO2/PCM method is expressed in terms of the internal energies E0 of fragments (also called monomers) I and their pairs (dimers) IJ. Earlier publications describe the general formulation of FMO64 and give the details of FMO/PCM.29 EFMO2=PCM ¼

∑I

þ

∑

I>J

½ðE0IJ

E0I

E0J Þ

Cq ¼

þ TrðΔD V Þ IJ

IJ

þ ΔEcav þ ΔEdisp þ ΔErep

ð1Þ

~ ¼ W

∑I

∑

I>J

ð6Þ

ε1 ~ W ε

ð7Þ

N

∑I W~ I

þ

N

~ IJ W ~I W ~ JÞ ðW ∑ I>J α ∑ α ∈ X jR α R i j

Z

ð8Þ ð9Þ

ð10Þ

One can also define the internal PIE (denoted by “in”) 00 00 00 IJ IJ ΔEin IJ ¼ ðEIJ EI EJ Þ þ TrðΔD V Þ

ð11Þ

Still, the goal is not achieved, because the values of ΔEint IJ are typically within 1 kcal/mol equal to ΔEin IJ , and neither is properly screened. What is the reason for the failure of definitions in ΔEint IJ and ΔEIJ ? The answer is that they fail to include the

es es ðΔEes IJ ΔEI ΔEJ Þ

þ ΔEcav þ ΔEdisp þ ΔErep

Z

0 0 0 IJ IJ ΔEint IJ ¼ ðEIJ EI EJ Þ þ TrðΔD V Þ

ð2Þ

ΔEes I þ

α ∑ jR Rαj i α∈X

ð5Þ

and C is the NT NT square matrix describing the geometric ~ and W ~ X in eqs 7 and 8 are details of the tessera interaction, q, W, vectors of size NT and the definition of wi can be inferred from ~ X is the soluteeq 5. ε is the dielectric constant of the solvent. W induced electrostatic potential due to fragment I (X = I) or dimer X = IJ upon the tesserae i representing the solvent. There are several models of PCM in FMO:29 PCM[2] uses eq 8 directly and requires iterating q and FMO2 calculations until mutual consistency, PCM[1] neglects the second sum in the ~ and PCM[1(2)] uses eq 8, but performs only definition of W, one iteration to obtain q and one more FMO2 calculation with this q to obtain the energy. PCM[1(2)] is an economic and accurate method, often used in applications. 2.2. Definition of the Pair Interaction Energy in Solution. The straightforward definition of the pair interaction energy in solution following the gas phase is

∑I E00I þ I∑> J½ðE00IJ E00I E00J Þ þ TrðΔDIJ VIJ Þ þ

q

~ Xi ¼ TrðDX w i Þ þ W

so that the FMO/PCM energy becomes EFMO2=PCM ¼

NT

i jνæ, μν ∈ X Æμj ∑ jr Rij i¼1

where

ΔDIJ = DIJ DI x DJ is the difference between the dimer and the direct (block-wise) sum of the monomer densities; VIJ is the electrostatic potential acting on dimer IJ exerted by all fragments K 6¼ I,J. ΔEcav is the cavitation energy, ΔEdisp and ΔErep are the solutesolvent dispersion and repulsion energies, respectively. The solutesolvent electrostatic energy is included in E0 . All monomer sums are from 1 to N, where N is the number of fragments. The pair interaction energies (PIEs) in vacuum physically should have very different ES components compared to solution; while other components can also be affected but probably not as dramatically. However, PIEs in solution (individual terms in the second sum in eq 1) when computed with FMO/PCM still feature very large total values for the case when fragments have a charge, even though one expects that solvent should screen the interaction. But the screening of pair interactions does not come out of the calculations. Let us examine why and where the screening hides. One can define the internal solute energies E00 X by separating the solutesolvent electrostatic interaction ΔEes X (X = I or IJ) E00X ¼ E0X ΔEes X

NT

∑ qi wiμν i¼1

E0 X is the QM energy of X in solution, including the solute solvent electrostatic interaction and excluding the ESP from other fragments. E00 X is the internal solute energy of X. The solvent in PCM is represented by the cavity surface, divided into NT pieces called tesserae; the charges of tesserae form the vector q (of size NT); α runs over atoms, i, over tesserae. Zα are atomic nuclear charges. μ and ν represent atomic orbitals. Ri and Rα are the coordinates of tesserae and atoms, respectively. WXe and WXN are the solvent-induced electrostatic potentials due to the interaction of the solvent with the electrons and nuclei of the solute, respectively. Now we introduce the PCM equation to obtain the apparent surface charges (ASCs) on the tesserae

2. THEORY

E0I

ð4Þ

ð3Þ 705

dx.doi.org/10.1021/jp209579w |J. Phys. Chem. A 2012, 116, 704–719


ARTICLE

by its own density (eq 13); in the dimer case, they are induced by the density of dimer IJ (eq 14); and in the full-solute case, they are induced by the whole system (eq 12) at the level of FMO/PCM (e.g., PCM[1]). In this work, we choose to use the full-solute definition in eq 12, which reflects the induced charges by the whole system, and is consistent with the electronic density determination for the electronic state of the solute (in FMO, the fragment densities are determined self-consistently polarized by the electrostatic potential of the whole system). The full-solute screening of the IJ pair interaction includes the effect of the charge induction by all fragments, not just I and J, which should be remembered when analyzing the results. In particular, some unexpected values of the screening can be obtained if the fragmentation is performed in such a way that many-body effects are large in the determination of ASCs. In addition, having to solve eq 13 or 14 for each monomer or dimer is very time-consuming and is not practical. ΔEes I is the ES solvation energy of monomer I due to the whole cavity, which can be decomposed into the contributions from individual monomer tesserae making the cavity.

Figure 1. Induced surface charges δ around fragments I and J, separated by RIJ.

contributions from the terms that, at first sight, do not arouse interest, ΔEes I , and we next show how and why they contain the screening. Let us consider the pair interactions of fragments I and J in solution. In the presence of a solute, charges are induced in the solvent, represented by ASCs in PCM. These charges surround each fragment, as shown in Figure 1. The induced charges are the microscopic representation of the screening. In PCM, the interaction of these charges with the solute is contained in the ΔEes X terms. The most direct effect is es es es in ΔEes I while ΔEIJ ΔEI ΔEJ are typically small (see eq 3). To uncover the screening, one should consider the interaction of fragment I with the induced charges around fragment J and, similarly, of fragment J with the induced charges around fragment I. The sum of these interactions gives the direct screening in FMO/PCM. Define the ES interaction between fragment I and the solvent tesserae i of fragment J. es ΔEes IðJÞ ΔEI jq ¼ qJ , ( 0, i ˇ J ε1 J ~ W , Cq ¼ qi ¼ qi , i ∈ J ε

es ΔEes I ¼ ΔEIðIÞ þ

¼ ΔEes ΔEes1 I IðIÞ

ð12Þ

es es ΔEes2 IJ ¼ ΔEIðJÞ þ ΔEJðIÞ

ð17Þ

In this work, we do not focus on the individual components es es and ΔEJ(I) , which describe the screening contributions ΔEI(J) made by each solute fragment in the pair, interacting with its partner’s portion of the solvent cavity, but it can be done. This may be of interest if only one fragment in the pair has a nonzero charge. Now we clearly separate the total sum as follows !

∑I ΔEesI ¼ ∑I ¼

ΔEes IðIÞ þ

ΔEes ∑ IðJÞ J 6¼ I

∑I ΔEesIðIÞ þ I∑> J ðΔEesIðJÞ þ ΔEesJðIÞÞ

ð18Þ

Let us now consider the meaning of the electrostatic dimer term in eq 3 and use the fact that the nuclear contribution is additive.

ΔEesI jq ¼ qJ ,

ΔEesIðJÞ ΔEesI jq ¼ qJ , ( 0, i ˇ J ε 1 IJ J ~ W , CqIJ ¼ qi ¼ qIJi , i ∈ J ε

ð16Þ

The pairwise screening is (the es2 notation denotes the two-body electrostatic solutesolvent interaction)

computed for the whole cavity with the full consideration of the many-body effects (i.e., the coupling between all monomers). A different picture to represent the screening would be obtained if one considered the charges induced only by each fragment independently; such a description is a simpler representation of the screening of pair interactions; however, it is not very practical as the coupling is strong and it seems better to base the screening analysis on the standard polarized state in PCM, that is, when ASCs are determined from the potential of the whole system and, afterward, are divided into fragment contributions a posteriori. In the case of a large coupling, one can expect that the pair IJ screening can at times be strongly affected by other fragments K ¼ 6 I,J implicitly via the ASC determination based on the expansion in eq 8. Alternatively to the above full-solute definition of ΔEes I(J), one can also suggest the monomer and dimer definitions given, respectively, by ( 0, i ˇ J ε1 J J ~ qi ¼ , C~ qJ ¼ W ~qJi , i ∈ J ε

ð15Þ

es The meaning of ΔEI(I) is the ES interaction between the solute fragment I and its own part of the total solute cavity. These terms correspond to “es1”, the ES solutesolvent interaction for monomers (the label 1 in es1 refers to one-body).

It is important to understand that the charge assignment in qJi and ΔEes I(J) is a posteriori. The apparent surface charges in PCM are

ΔEesIðJÞ

ΔEes ∑ IðJÞ J6¼ I

ð13Þ

1 es es es IJ IJ ð19Þ ΔEes3 IJ ¼ ΔEIJ ΔEI ΔEJ ¼ TrðΔD W e Þ 2 The term contains the electrostatic interaction of the solvent (its electronic part) coupled with the change in the solute density due to the charge transfer in dimer IJ. This is a higher order term than the previously considered es1 and es2, and it can be decomposed (through WIJe , similar to eq 15) into two- and three-body terms, just like ΔEes I has both one- and two-body contributions. We denote the sum of these contributions as “es3” (the label three here points out to three-body terms included in this quantity, which also includes some two-body terms), and

ð14Þ

In other words, the difference in the three definitions of ΔEes I(J) is the extent of including many-body effects in the PCM equation used to generate the surface charges determining the screening. In the monomer case, the charges on J are only induced 706



ARTICLE

choose not to decompose it into individual components at this time (the absolute values of “es3” are typically small, less than 1 kcal/mol). Now we can define the PIE made of the internal and screening components,

Compared to gas phase FMO, these energies include the destabilization polarization due to the solvent. We note that the polarization of fragments includes two contributions: (a) destabilizing, which shows the change in the internal energy of a fragment from its free state to a bound state in a larger system, and (b) stabilizing, which is the gain in the electrostatic interaction energy between the fragments polarized in a system; more details = E00 I + ΔEsolv is the energy of can be found elsewhere.60 EPCM I I fragment I in solution, including the solutesolvent interactions (electrostatics, dispersion and repulsion) between fragment I and the part of the cavity surrounding it (made of atoms in I). It also includes the cavitation energy of I, that is, mainly the entropy loss of between the solvent to create the cavity for I. The change in ΔEsolv I bound and isolated states during complex formation is the desolvation energy, as described in detail below. Next, consider

þ scr es2 00 00 00 ¼ ΔEint ΔEin IJ IJ þ ΔEIJ ¼ ðEIJ EI EJ Þ es3 þ TrðΔDIJ V IJ Þ þ ΔEes2 IJ þ ΔEIJ

ð20Þ

This definition contains the electrostatic screening es2 es3 ΔEscr IJ ¼ ΔEIJ þ ΔEIJ

ð21Þ

Let us examine the remaining PCM specific terms, Ecav + Edisp + Erep, and decompose them as ΔEcav ¼

∑I ΔEcav I

ΔEdisp ¼

∑I

ΔErep ¼

∑I

disp ΔEI

rep ΔEI

ð22Þ

∑I

¼

¼

disp ΔEIðIÞ

∑I

rep ΔEIðIÞ

∑

þ

J6¼ I

þ

∑

J6¼ I

¼ ðE00IJ E00I E00J Þ þ TrðΔDIJ V IJ Þ ΔEPCM IJ

! disp ΔEIðJÞ

disp

es3 þ ΔEes2 IJ þ ΔEIJ þ ΔEIJ

ð23Þ

ð24Þ

disp where ΔEI(J) describes the dispersion interaction between atoms in solute fragment I and solvent tesserae of fragment J, and analogously for the corresponding repulsion. These terms, including two fragment indices, correspond to the solutesolvent contributions and can also be included into PIEs. The cavitation energy70,71 is computed as a sum of the contributions derived from the atomic surfaces; it can be divided into fragment contributions of formally one-body nature and makes no contribution to PIEs in our formulation. Because the surface of a fragment is affected by its environment (i.e., other fragments), it is in fact a complex quantity to be divided. It is also possible to leave the cavitation energy as is without its decomposition into fragment contributions to avoid the ambiguities in its division. The dispersion and repulsion components can be summed as

disp

ΔEIJ

disp

disp

¼ ΔEIðJÞ þ ΔEJðIÞ

rep

rep

ð25Þ

rep

ΔEIJ ¼ ΔEIðJÞ þ ΔEJðIÞ

ð26Þ

Now we are ready to define PIE in PCM, ΔEPCM IJ

¼

þ scr ΔEin IJ

þ

disp ΔEIJ

þ

rep ΔEIJ

EX CT þ mix solv ¼ ΔEES þ ΔEDI ΔEPCM IJ IJ þ ΔEIJ þ ΔEIJ IJ þ ΔEIJ

ð27Þ

ð32Þ

Finally, we can rewrite the total PCM energy equation as EFMO2=PCM ¼

∑I

E00I þ

∑I

ΔEsolv þ I

∑ I>J

ΔEPCM IJ

disp

es2 es3 ΔEsolv IJ ¼ ΔEIJ þ ΔEIJ þ ΔEIJ

ð28Þ

disp

rep

ð29Þ

The PCM energy of monomer I is ¼ E00I þ ΔEsolv EPCM I I

rep

þ ΔEIJ

ð33Þ

ΔEES IJ is computed as the ES interaction between two fragments in FMO/PCM and it is not screened, while the screening es3 solv ΔEes2 IJ + ΔEIJ is included in ΔEIJ . We compute the exchangeEX repulsion ΔEIJ using the same algorithm as for the gas phase, by estimating the energy of the orthogonalized wave function for the dimer (for the PCM Hamiltonian, so there is a small explicit PCM contribution due to the electronic density coupled with the + mix is obtained from the internal PCM potential), and ΔECT IJ in ES DI DI solute PIE as ΔEIJ ΔEIJ ΔEEX IJ ΔEIJ . ΔEIJ is calculated as the difference in the correlation energy of dimer IJ and two

where ΔEsolv ¼ ΔEcav þ ΔEes I I IðIÞ þ ΔEIðIÞ þ ΔEIðIÞ

ð31Þ

In the above expression, there are six terms. The first term E00 IJ E00 I E00 J is the internal solute fragment interaction energy, in which the influence of the polarizing ESP exerted by other fragments is separated into the second term Tr(ΔDIJVIJ), describing the explicit charge transfer effect between the two monomers I and J coupled with ESP. Implicitly, the charge transfer energy is also present as a part of the first term, because E00 IJ is the dimer energy already reflecting the charge transfer. Tr(ΔDIJVIJ) is equal to zero if there are only two fragments, because VIJ iz zero. The third term ΔEes2 IJ is the explicit solventinduced electrostatic screening of PIE. The fourth term ΔEes3 IJ is the coupling between the solvent-induced electrostatic field (its electronic contribution) and the explicit charge transfer between is the solutesolvent dispersion I and J. The fifth term ΔEdisp IJ correction to PIE. The sixth term, ΔErep IJ , is the solutesolvent repulsion correction to PIE, that is, the part of the solutesolvent interaction energy excluding both the electrostatics and dispersion; it corresponds to the exchange-repulsion in HartreeFock (HF). In the following discussion of PIEDA, we assume that it is performed for the PL state,60 that is, the polarized state of fragments, obtained by their self-consistent mutual polarization, which is a part of the FMO calculation. At this point we do not extend our formalism to the PL0 state (which refers to the interaction analysis of free monomers, involving their single, i.e., not selfconsistent, mutual polarization from the free state). In PIEDA/ PCM, PIEs are divided into the following five contributions,

! rep ΔEIðJÞ

rep

þ ΔEIJ

ð30Þ

Let us give a physical description of all quantities in eq 28. E00 I is the internal solute energy of fragment I in FMO/PCM. 707



ARTICLE

monomers I and J in solution. It should be noted that the solventsolute electrostatic interaction implicitly affects all EX CT+mix , and ΔEDI terms ΔEES IJ , ΔEIJ , ΔEIJ IJ . Let us consider what happens to the three-body formulation46,65 of FMO with explicit trimer corrections, FMO3/PCM.

surface of each atom taken without contacts to other atoms. It is a somewhat cumbersome quantity to obtain for atoms in molecules, also complicated by the tessera generation algorithms in PCM, which reflect the intersections to other spheres, as well as the surface area scaling (AS) adopted in some cavitation schemes in PCM (such as the generating polyhedra (GEPOL)-AS), dependent upon the surrounding atoms. Therefore, we introduce the following atomic coverage based on the theoretical value of the sphere surface as a function of its radius rα, σα sα ¼ 100% ð40Þ 4πrα 2

EFMO3=PCM ¼ EFMO2=PCM

∑ ðEIJK EI EJ EK EIJ EIK EJK Þ I>J>K þ ∑ ½TrðΔDIJK V IJK Þ TrðΔDIJ V IJ Þ I>J>K 0

þ

0

0

0

0

0

0

TrðΔDIK V IK Þ TrðΔDJK V JK Þ

ð34Þ

00

We rewrite the first sum in terms of E X, 0

0

0

0

0

0

0

00

00

00

Such a definition of the coverage also reflects to a minor extent the deviation of the tessellated surface from the perfect spherical surface. Some tests we performed for sodium and chlorine indicate that this definition results in a normalization very close to 100% for atomic PCM calculations, so it should suffice for the purpose of the analysis. Next, we introduce the fragment coverage. Here, we took a somewhat different definition not derived from the atomic coverages, because the atomic spheres overlap, and this reduces their surface irrespective of how close they are to the solvent. σI ð41Þ sI ¼ 100% σ

EIJK EI EJ EK EIJ EIK EJK 00

00

00

00

¼ EIJK EI EJ EK EIJ EIK EJK þ ΔEes4 IJK 1 IJK IJK ΔEes4 IJK ¼ ½TrðΔΔD W e Þ 2

ð35Þ

where ΔΔDIJK DIJK (DIJ x DK) (DIK x DJ) (DJK x DI) + 2DI x DJ x DK and ΔDIJK DIJK DIx DJ x DK. For an explicit definition of the three-body effects, one can divide ΔEes3 IJ as well as ΔEes4 IJK into components, which we do not do, as they are very small and it is more convenient to operate with their total values. The physical meaning of ΔEes4 IJK is the explicit three-body charge transfer (i.e., the additional charge transfer inside trimer IJK, complementing the two-body charge transfers in the three dimers IJ, IK, and JK) coupled with the solvent-induced electrostatic potential (its electronic contribution). The label 4 in es4 refers to the four-body effects included in this term, in addition to some three-body effects. 2.3. Fragment Properties in Solution. In the further analysis, it is convenient to define the fragment surfaces areas σI, σI ¼

∑

α∈I

σα

where σ is the total cavity surface, σ = ∑IσI. The dimer coverage can be defined as a sum of the monomer values, sIJ ¼ sI þ sJ

As indicated above, there are several cavities and their corresponding surfaces, of which we focus on the electrostatic terms es σes α and σI . In addition to the surface, it is useful to define atomic and fragment charges for the cavity by summing the corresponding ASCs,

ð36Þ

qI ¼

in terms of the atomic surface areas σα, σα ¼

∑ σi

ð37Þ

qα ¼

i∈α

which are computed by summing the surface areas for tesserae i located on atom α. σI gives a convenient measure of the surface , one of fragment I exposed to solvent. In the discussion of ΔEPCM IJ can also take advantage of the dimer surface area, σIJ ¼ σI þ σJ

disp, H

¼

σI

qα

ð43Þ

∑ qi

ð44Þ

i∈α

∑i qi ¼ ∑α qα ¼ ∑I qI ¼

ð38Þ

ε1 Q ε

ð45Þ

where Q is the total charge of the system, obtained as the sum of the monomer charges QI, Q ¼

∑I QI

ð46Þ

In the case of ICOMP=0, there is a small deviation of ∑iqi from the indicated theoretical value, usually explained with the charge penetration effects (i.e., the solute charge distribution penetrating into the continuum beyond the cavity). The extent in which the many-body effects are important in qI (and thus the appropriateness of the formal charge division, leading to the screening in eq 45) can be analyzed by comparing their values to the theoretical estimate of ((ε 1)/ε)QI; a deviation from this value also includes the charge penetration effect.

disp, O

þ σI ð39Þ 2 It is of interest to obtain a measure of how much an atom (or a fragment) is “buried”. Such a measure is related to the full disp

∑

α∈I

The cavity fragment charges q are important for the discussion of the screening; they are either unnormalized (ICOMP=0) or normalized (ICOMP=2)66 to their proper value, that is,

In PCM, the cavitation, dispersion, and electrostatic interactions are computed using three different cavities. The repulsion component uses the same cavities as the dispersion. Thus it is necessary to introduce a superscript indicating what cavity is meant, e.g., the electrostatic one, σes I . Moreover, there is a separate dispersion cavity ; for example, for water, there are for each solvent atom type A, σdisp,A I two cavities representing the oxygen (A = O) and hydrogen atoms (A = H). Therefore, one has to specify the atom type in that case, and one can also define the average in order to simplify the discussion. σI

ð42Þ

708



ARTICLE

One can also define the effective dielectric constant ε for the pair interactions based on the extent in which the interactions are screened by the solvent, either with respect to vacuum or the internal solute PIE, for example, for the latter case in þ scr εIJ ¼ ΔEin IJ =ΔEIJ

gets all of the electrons of BDA except the one electron assigned to the right fragment, so that the latter fragment contains the covalent bond. The atomic spheres representing the two atoms forming a covalent bond overlap, and the corresponding tesserae in the intersection between the two atoms are removed during the cavity construction. This, in particular, means that the parts of the two spheres corresponding to the bond are largely removed, thus, all retained tesserae on BDA are properly assigned to the left fragment and none to the right. This is how we treat the issue of redundant BDAs, resulting in the condition that each tessera contributes to exactly one fragment. Typically, BDAs are backbone atoms not exposed to the surface and have a small charge on them (in biological systems, BDAs are Cα); all of this makes the contributions from them to PCM-related properties small and thus the effect in which their contributions are divided also insignificant.

ð47Þ ΔEPCM IJ

are not where the nonelectrostatic contributions in included. An alternative definition can be done in PIEDA explicitly for the electrostatic interaction, which can be computed in the gas phase or solution, and the former choice is the most strictly proper: ES, gas

εIJ ¼ ΔEIJ

PCM es3 =ðΔEES, þ ΔEes2 IJ IJ þ ΔEIJ Þ

ð48Þ

In biochemistry, the idea of a local ε is widely used,67 for example, in the solvation models50,68 or to discuss the interactions.69 However, the values of ε are often taken at will, based on the experience of how well they represent the physics in a given problem. Here, we suggest an ab initio derived mechanism for the evaluation of ε. For the perfect screening, which does not take place for partially solvated fragments, εIJ is the same as the universal dielectric constant of the solvent. Now we can show that the screening in the limit indeed gives the expected relation. Consider the case when the distance RIJ between I and J is very large. Also, assume that we have a single point charge QI and QJ representing the two fragments I and J, respectively, whose induced charges are qI and qJ. The distance between the charges QI and QJ and qI and qJ can be taken to be RIJ (in the following discussion we assume that QI and QJ are not zero). Then in gas phase the interaction energy is gas

ΔEIJ ¼

QI QJ RIJ

3. COMPUTATIONAL DETAILS There are many ways to perform PCM calculations.2 First of all, there are several polarizable continuum models, conductorbased C-PCM, dielectric-based D-PCM and integral equation formalism IEF-PCM; second, there are different tessellation schemes (GEPOL-GB, GEPOL-AS, FIXPVA, etc.); third, there are charge compensation schemes (ICOMP=0, 2, etc.); fourth, there are atomic radii; and fifth, the cavitation discretization, that is, the density of tesserae per atomic sphere. In FMO there are several levels to evaluate the charges q, such as PCM[1], PCM[1(2)], and PCM[2]. All of these factors play either a large or a small role in the exact values of the components in the analysis. In the tests below, we do not undertake a comprehensive study of the effect all of these parameters have on the results, but we do compare several choices that appear to us to be particularly interesting. Throughout this work, we use atomic units of charges. We specified 60 tesserae per cavity and employed the Claverie-Pierotti cavitation energies.70,71 Namely, we compare GEPOL-AS72 to FIXPVA,50 PCM[1]29 to PCM[1(2)],29 no charge compensation (ICOMP=0) versus doing it (ICOMP=2),66 and C-PCM (closely related to COSMO73)74 to IEF-PCM75 (related to surface and simulation of volume polarization for electrostatics, SS(V)PE76). The analysis proposed in this work was implemented into GAMESS package77 and parallelized with the generalized distributed data interface (GDDI).78 To eliminate possible numerical artifacts when studying the distance dependence of properties, we switched off all FMO approximations in the two-fragment systems, while in the small protein the default values were used. All calculations were performed on a single core Pentium4 3.2 GHz with 2 GB RAM node of the Soroban cluster, which survived the great earthquake of 2011.

ð49Þ

whereas in solution (sol) it is ΔEsol IJ ¼

QI Q J 1 QI qJ 1 q I QJ þ þ 2 RIJ 2 RIJ RIJ

ð50Þ

We note that 1/2 in eq 50 is directly related to 1/2 in eq 4. At a large distance, the charges qI on I are induced only by QI qI ¼

ε1 QI ε

ð51Þ

and the same for J, ε1 QJ ε By substituting these charges, we get qJ ¼

ΔEsol IJ ¼

QI Q J ε 1 QI QJ QI QJ ¼ ε RIJ RIJ εRIJ

ð52Þ

ð53Þ

which is the expected relation, so that gas

ΔEIJ =ΔEsol IJ ¼ ε

4. RESULTS AND DISCUSSION

ð54Þ

4.1. Prototypical Ionic Screening: Sodium Chloride. As a perhaps the most typical and simple (prototypical) model of the solvent induced screening, we study the solvation of sodium chloride. There is a wide variety of radii used in solvation models.49,79 One choice is given by the experimental radial distribution functions (RDF) of solvated Na+ and Cl, which have the first NaO peak of about 2.40 Å,80,81 in agreement with 2.40 Å in CPMD;82 and the ClH peak of 2.28 Å,83 in reasonable agreement with 2.22 Å in classical MD.84 Therefore, we set the

This relation is satisfied only at large distances as otherwise various assumptions do not hold. Also, in practical calculations, the atomic radii as well as other details of PCM strongly affect the deviation of the screening from the limiting value in eq 54. Finally, it is necessary to deal with the bond detachment issue. In FMO,28 for each detached bond in the fragment definition, one atom called the bond detached atom (BDA) has to be put into two fragments, denoted as left and right. The left fragment 709



ARTICLE

Table 1. Dependence of the Cavity Surfaces σes (Å2), Induced Charges q, and Effective Dielectric Constants εIJ on the RHF/PCM Modela basis set

a

model

ICOMP

radiib

tesselation

σes Na

qNa

σes Cl

qCl

εIJ

6-31G*

C-PCM

0

HLi

FIXPVA

40.72

0.9871

95.73

0.9863

75.20

6-31G* 6-31G*

C-PCM IEF-PCM

2 2

HLi HLi

FIXPVA FIXPVA

40.72 40.72

0.9870 0.9870

95.73 95.73

0.9870 0.9870

76.88 76.88

6-31G*

C-PCM

2

HLi

GEPOL-AS

40.72

0.9870

95.73

0.9870

76.88

6-31G*

C-PCM

2

UAHF

FIXPVA

40.28

0.9854

70.94

0.9854

68.36

6-31G*

C-PCM

2

RDF

FIXPVA

72.38

0.9843

65.33

0.9843

63.56

6-31G*

C-PCM

2

vdW

FIXPVA

96.56

0.9825

59.28

0.9825

56.99

6-31+G*

C-PCM

0

HLi

FIXPVA

40.72

0.9870

95.73

0.9799

60.42

6-31+G*

C-PCM

2

HLi

FIXPVA

40.72

0.9846

95.73

0.9846

65.10

6-311G* 6-311G*

C-PCM C-PCM

0 2

HLi HLi

FIXPVA FIXPVA

40.72 40.72

0.9872 0.9868

95.73 95.73

0.9856 0.9868

73.48 75.72

6-311+G*

C-PCM

0

HLi

FIXPVA

40.72

0.9871

95.73

0.9803

61.39

6-311+G*

C-PCM

2

HLi

FIXPVA

40.72

0.9849

95.73

0.9849

66.02

Calculated for Na+ and Cl separated by 100 Å. b See main text for the values and references defining the radii.

RDF radii of Na and Cl to be 2.0 and 1.9 Å, respectively, so that the cavities, which are defined with the radii multiplied by α = 1.2, resemble the experimental distribution of solvent. These radii can be compared with 2.31 (Na) and 1.81 Å (Cl), the default in GAMESS based on the van der Waals (vdW) radii. Li et al.49,50 used the values of 1.5 and 2.3 Å for Na+ and Cl, respectively (the quoted values of 1.8 and 2.76 Å include the multiplication factor α = 1.2), and we denote them as the HLi radii (not “Li” radii to avoid possible confusion with lithium). We also tested the popular UAHF85 set of radii, 1.492 Å85 for Na+ and 1.980 Å79 for Cl, which we use as all other radii with the scaling factor of α = 1.2, as recommended.85 For the tests, we used RHF and each ion was made a separate fragment. FMO2-RHF/PCM[2] for two fragments is identical to ab initio RHF/PCM. Only for the PL0 state and two fragments PIEDA is nearly same as the original EDA,86 thus the analysis in this work for two fragments should not be called EDA/PCM, because the former would need a definition of the isolated states of monomers, which is not done here. At large distances there are two separate cavities and the solventsolvent dispersion and repulsion models cannot be applied; therefore, we did not use them in this system (and, correspondingly, instead of MP2 for the solute, we used RHF). First, we investigate the limiting case, which we represent by the 100 Å separation between the two ions. This is long enough and yet the interaction energy is substantial, so that the inaccuracies in the calculations such as the SCF convergence do not affect the results. At this separation, PCM[2] is numerically the same as other models, so we do not vary this degree of freedom. The results are summarized in Table 1. We observe a fairly broad range of εIJ values from 56.99 to 76.88. The correct result is 78.390, which is the value built into all PCM models for water. For all results, the electrostatic component in PIEDA, ΔEES IJ is the same, 3.32064 kcal/mol. Note that the simple model (QI = +1, QJ = 1) in eq 49 gives 3.32064 kcal/mol, in perfect agreement. The model and basis set dependence thus comes in not as the QM interaction between two fragments, but rather as the screening only, which depends upon the locally induced charges. For comparison, we also performed a gas phase calculation, and obtained 3.32064 kcal/mol, thus, the unscreened ES interaction in gas phase and PCM

in this case is practically the same. Other components in PIEDA are zero (EX and CT+mix). Several trends for εIJ can be discerned from the results. First, increasing the basis set from double-ζ to triple-ζ gives a shift of about 1.5. ICOMP=2 relative to ICOMP=0 for the basis sets with (without) diffuse functions gives a shift of about +4.7 (+2.0). Thus, we observe that the former basis sets are more strongly affected, which is quite reasonable as the charge penetration is larger. Changing the model from C-PCM to IEF-PCM (D-PCM is not implemented for FMO) as well as the tessellation from FIXPVA to GEPOL-AS makes almost no effect upon εIJ. On the other hand, the radii have the strongest effect among all factors considered. The HLi, UAHF, RDF, and vdW radii give the εIJ values of 76.88, 68.36, 63.56 and 56.99, respectively, at the same level. The origin of the radii dependence is clarified below. For ICOMP=2 the charges are renormalized, so that for the case when there are only two fragments in the system and the total charge Q = 0 qI þ qJ ¼ 0

ð55Þ

or qNa = qCl in the case when each fragment is made of one 6 0, atom. Then assuming the relation in eqs 51 and 52, and QI = QJ ¼ one obtains that εI ¼

QJ QI ¼ ¼ εJ q I þ QI q J þ QJ

ð56Þ

which can also be related to the more general one in eq 48 (we note that ΔEes3 IJ = 0 at large distances). Also, without the assumption in eq 55, one can define the internal ε for each fragment I with a nonzero charge QI, εI ¼

QI q I þ QI

ð57Þ

The internal dielectric constant εI here is system-specific (because the charges q are determined from eq 7) but pair independent. On the other hand, the effective dielectric constant εIJ is both system and pair dependent. 710



ARTICLE

Then for both ICOMP=0 and ICOMP=2, one gets from eq 50, " #! QI QJ 1 εI 1 εJ 1 Q I Q J ¼ ð58Þ 1 þ 2 εI εJ RIJ εRIJ hence ε¼2

εI εJ εI þ ε J

ð59Þ

which for ICOMP=2 with 2 fragments reduces to a single ε = εI = εJ. The internal εI in eq 57 is important conceptually. It clearly shows that the effective dielectric constant is essentially determined by the properties of each fragment, provided that the separation is large, and the polarization and quantum effects are negligible. In other words, the screening itself always appears for pair interactions, but at large distances it is determined by the locally induced solvent charges. However, qI in the atomic (Na+) and diatomic (NaCl) calculations are not the same for ICOMP=2 even at large separation, because the charge renormalization in diatomics is performed for the total charge, not individual atomic charges. So for both ICOMP=0 (see eq 59) and ICOMP=2, the effective dielectric constant in dimer IJ at large separation differs from that determined in atomic calculations. If one fragment in a pair has a zero charge, then at large distances qIQJ = qJQI = QIQJ = 0 and any effective dielectric constant can be assumed. Usually, the induced fragment charges qI in molecules are not exactly 0 so that eq 57 gives a zero εI for neutral fragments (QI = 0). We note that the radii fitting scheme of Dupuis and coworkers,87,88 is based on a simple relation with few parameters, determining the radii on the fly as function of the atomic charges of the solute. Because these solute charges also determine in some way the induced solvent charges, the screening argument can be used to justify that approach. Although we note that the HLi radii happened to give the best εIJ, the problem of finding a good set of PCM radii is a hard one, and perhaps looking at εIJ can be one factor in determining them. One other suggestion arising from this work is the need to define independent radii for the cavitation, dispersion(repulsion) and electrostatic interactions. At present, different but linearly dependent radii are used, calculated by using various scaling factors α. One could introduce a new charge renormalization scheme, so that the individual induced charges qI obey the relation in eq 51, while in the present ICOMP=2 approach the sum of ∑IqI is renormalized (separately for the electronic and nuclear components). However, as shown below for chignolin, manybody effects make this idea dubious. It is clear that the charge renormalization has a very considerable effect upon the induced charges at large distances and thus the screening (at shorter distances the effect of the charge renormalization becomes less significant, because only at large distances a small variation in the charges introduced by the renormalization has a considerable effect). The individual charges qI reported in Table 1 for ICOMP=2 depend strongly upon the radii of both atoms (not just those in I), because of the nature of the charge renormalization applied to scale the charges with respect to the total charge (see eq 45).

Figure 2. Internal dielectric constant as a function of the atomic radius rα for α = Na and Cl, computed at the RHF/6-31G*/C-PCM level and for atoms Na+ and Cl separately with ICOMP=0 (ICOMP=2 automatically gives εI = 78.39 for atoms).

One should not be confused into thinking that ICOMP=2 completely removes the charge penetration problem from PCM. All it does is to rescale all charges with a single factor. Thus in the NaCl system, sodium is made to pay for the escaped charge on chlorine. There are other schemes such as ICOMP=4,89 which are not implemented in our program, and some improvement may be possibly obtained with them. It is instructive to perform atomic PCM calculations with ICOMP=0 and calculate the internal dielectric constant (eq 57). The results are plotted in Figure 2. Sodium has reasonable effective dielectric constants in the range of 77.6477.31 for the atomic radii of 1.4922.31 Å. Chlorine, on the other hand, has ε varying between 38.95 and 73.03 for the radii of 1.812.3 Å. The dielectric constant decreases with the atomic radius for sodium but increases for chlorine. Thus, one can see that in the diatomic system the results are determined almost entirely by the change in the radius of chlorine, which affected the observed dielectric constant. Of course, the ICOMP=2 charge renormalization automatically generates the right εI in atoms for any atomic radius, but it fails to do so in the diatomic system. The trends are somewhat obscured by the combined effect of the two radii; and it is possible to obtain more details upon the radii effect by fixing one radii and varying the other. However, in this system the radius of sodium makes a little difference. Especially in anions, small changes in the radii evoke a large effect in the induced charges and the dielectric constants. Next, we varied the NaCl distance in the range of 210 Å, and investigated what happens to the PIEDA components as well as the effective dielectric constant. In these calculations, we used C-PCM, ICOMP=2, and the HLi radii. The results are plotted in Figures 3 and 4. In Figure 3, we varied the tesselation scheme to see if the observed peak is due to some peculiarity in any of them. However, GEPOL-GB, GEPOL-AS, GEPOL-RT, and FIXPVA all exhibited a similar peak shape, differing in the peak position and height (all of them have the same limit of εIJ = 76.02 and εI = 75.91 at R = 10 Å). For FIXPVA (Figure 3d), one can see that the effective dielectric constant determined by the ratio in eq 48 has a maximum near RNaCl = 5.75 Å, with the effective εIJ = 115.28. The atomic charges (Figure 4a) show a fairly smooth dependence on the distance. Interestingly, the atomic coverages 711



ARTICLE

Figure 4. RHF/6-31G*/C-PCM/HLi/ICOMP=2 properties of solvated Na+ 3 3 3 Cl as a function of the interatomic distance R: (a) charges (CT is the charge transfer in fractions of electron between Cl and Na+), qNa = qCl is the induced solvent charge, QNa = QCl is the solute Mulliken atomic charge, (b) atomic coverages, sNa and sCl for the electrostatic cavity, and (c) PIEDA/PCM energy components. Figure 3. Effective (eff, εIJ) and internal (int, εI = εJ) dielectric constants of solvated Na+ 3 3 3 Cl for RHF/6-31G*/C-PCM/HLi/ ICOMP=2 as a function of the interatomic distance R: (a) GEPOLGB, (b) GEPOL-AS, (c) GEPOL-RT, and (d) FIXPVA.

reject the supposition that the peak is caused by FIXPVA. On the other hand, the sum of the two RDF peaks for Na+ and Cl is 4.68 Å, which when combined with 0.95 Å of the OH distance, gives the approximate maximum radius of the first solvation shell of 5.63 Å. Both experimentally and theoretically, solventseparated ion pairs have attracted considerable attention,90 corresponding to our region of the ionic separation. Thus, it is reasonable to expect that something may be happening with the effective dielectric constant at that separation but we do not have the confidence to expect that even a qualitatively correct description of this is delivered by the continuum model where all atomistic details are but shadows on the cavity. At 10 Å, εIJ = 76.02, and the dielectric constant slightly increases at larger distances, reaching at 100 Å the value of εIJ = 76.88, as discussed above (i.e., there is a very shallow minimum in this range of distances). Among other interesting results (Figure 4c), we note that the exchange interaction ΔEEX IJ is very short-ranged, as expected, and is quite the energy of charge transfer plus mixed terms ΔECT+mix IJ long-ranged in this system, probably related to the diffuse nature of the excessive charge on chlorine. Also, at the distances of 4 Å or

(Figure 4b) in the vicinity of 5.5 Å show the last signs of the sphere intersection (sNa = 95.0%) before they go to 100% at larger distances (>5.75 Å). The sum of the radii of the actual spheres defining the cavity is 1.8 + 2.76 = 4.56 Å. On the contrary, the internal εI (eq 57) hides this maximum and the dielectric constant has a much simpler behavior, because its definition is less sensitive to the actual dimer structure (in fact, this definition is derived in the assumption that the other monomer in the pair is at a large distance). It is somewhat surprising that both definitions, εI and εIJ, show results which agree well with each other except for the maximum peak. Now we focus on the FIXPVA results. This maximum of εIJ is probably related to the jump in the total energy, observed by Li and co-workers49 at about 5.8 Å, who suggested two possibilities for its cause: (a) an artifact of the FIXPVA tessellation (related to the area scaling) and (b) the physical nature of the restructuring of the first solvation shell. Our results in Figure 3 at least for εIJ 712



ARTICLE

less, a considerable portion of the Na cavity is lost because it is hidden by the larger sphere of Cl. To provide more physical insight into the details of the polarization of the ionic cavities, we define the induced charge angular distribution showing the competition between the polarizing electrostatic potentials of the opposite sign coming from the two ions, which are added at each tessera.

∑

1 qi q̅ α ðθa Þ ¼ Nθ i ∈ α, i ∈ fθa g

ð60Þ

where we average all induced charges on atom α for those tesserae i on α, whose spherical angle θi = arccos ((zi Rα,z)/rα) lies within the interval {θa}, θa Δ/2 e θi < θa + Δ/2 (in NaCl, atoms are put on the z-axis, zi, and Rα,z are the coordinates of tessera i and atomic center α, respectively; rα is the sphere radius). The whole range of θ, [0;π], is divided into equally spaced 60 intervals, so that Δ = π/60, and θa gives the medium point of the interval. The normalization constant Nθ is equal to the number of tesserae in the interval (eq 60 is applied only to those {θa}, for which Nθ is nonzero). Due to the nature of the discrete tesselation algorithm, the points are not distributed equally with respect to θ, some ranges of θ happen to have no tesserae in them; it is necessary to normalize for a proper comparison (out of 60 intervals only at most 8 have some tesserae in them). θ is defined in such a way that for each atom in NaCl, θ = 0 corresponds to the closest contact between atoms. The range of θ reflects that some part of the sphere is lost in close contact of the two atoms, so that there are no small θ values for Na+ at a close separation R. The results are plotted in Figure 5. We see that in both cases of Na+ and Cl, at large R (i.e., R = 5.7 Å) the angular distribution is flat, meaning that the induced charges on the sphere are almost uniformly distributed (little affected by the other ion). As the distance becomes shorter, the charges near the contact of the two spheres (small θ region) are strongly polarized, whereas those at the far points (θ approaching π) are almost unaffected in Na+ and shifted by less than 30% in Cl. We observe two other interesting points. The charge densities for both atoms at small θ for R = 5.0 Å (near the point where the two spheres touch each other, 4.56 Å) differ very considerably from the other values for R around 5.0 Å. Another is that some induced charges on Cl are negative (at small R), and some are nearly zero for both atoms. The former ~ in eq 7 is nearly zero (cancels out) is because the potential W or is dominated by that of Na+, resulting in the negative charges. This is clearly depicted in Figure 4a (decrease in q at small R) and Figure 5 (decrease in q at small θ); such phenomenon can be called the induced charge quenching effect caused in this system by the partial loss of the solvation shell and the potential pressure of the nearby ion. Finally, we remark here that, when small effects are sought in such calculations, it is probably advisible to use a higher density of tesserae. The deviations of the tessera distribution from the perfect sphere should be expected to have some effect on the properties like the screening. However, we verified by using 960 tesserae per sphere (vs 60 used elsewhere in this work) that the εIJ peak is little affected (its height decreased by about 8), thus, it is not an artifact of the tesselation, neither due to its discretization degree nor a particular scheme of defining the tessera; the conclusion is that the peak in εIJ is apparently brought about by the particular ratio of the atomic radii.

Figure 5. Induced charge angular distribution q(θ) for Na+ 3 3 3 Cl at the RHF/6-31G*/C-PCM/HLi/FIXPVA/ICOMP=2 level: (a) Na+, (b) Cl. Curves are shown for various separations R (in Å).

Figure 6. Components of the pair interaction energies of solvated Na+ 3 3 3 Cl at the RHF/6-31G*C-PCM/FIXPVA level as a function of the interatomic distance R. The HLi or UAHF radii are used, and ICOMP is shown after the slash.

It is important practically to probe how much PIEs depend upon the solvation model. The results are presented in Figure 6. The internal PIE, ΔEin IJ depends very weakly upon the PCM model through the polarization by the solute solvent electrostatic interaction, and the values for the three 713



ARTICLE

which is one reason for the accuracy loss of FMO, and also for the purpose of the analysis, as the molecular mechanics study91 with explicit solvation showed, whose conclusion was that the two residues forming a salt bridge individually appear to make a large contribution to the ligand binding, which however, is offset by the desolvation penalty, with the net result that the salt bridges as a whole contribute relatively little and are better analyzed as one unit. differs very largely from ΔEPCM, Finally, we note that ΔEPCM IJ even though only two fragments are present. The former is the pair interaction energy, which is a part of the latter. In addition to that, ΔEPCM also includes the monomer energies, which reflect the fragment polarization (by other fragments and solvent) and other solvation effects (electrostatic and nonelectrostatic solutesolvent interaction energy). To illustrate the last point, which is important, we study the following process: þ Naþ aq þ Claq f ðNa 3 3 3 Cl Þaq

ð61Þ

The energy for it is given by (I = 1, Na+; J = 2, Cl) ΔE ¼ EPCM ðNaClaq Þ EPCM ðNaþ Þ EPCM ðCl Þ ¼ ΔEPCM þ ΔEPLd þ ΔEPLd þ ΔEdesolv þ ΔEdesolv IJ I J I J

ð62Þ E00 NaCl I +

ΔE00 Na+ I

= is the destabilization polarizawhere tion energy of I (Na ) induced by fragment J (Cl), as well as by = ΔEsolv,NaCl ΔEsolv,Na+ is the the solvent, and ΔEdesolv I I I + desolvation penalty of I (Na ). The same type of these two terms for J (Cl) is defined. The monomer terms (the last four terms in eq 62) give the difference in the PCM energies (see eq 30) of Na+ and Cl in the complex and as isolated ions. The sum of the monomer terms constitutes the difference between the PIE and the actual energy of formation ΔE. As noted above, sodium is partially buried at close distances in the larger sphere of chlorine, so its desolvation penalty should be expected to be larger. In other words, PCM for such radii predicts that during the complex NaCl formation sodium looses a larger part of its solvation shell than chlorine, leading to a larger solvation energy loss (i.e., larger desolvation penalty). The stabilization polarization energy is not explicitly separated in PIEDA/PL and it forms a portion of ΔEES IJ . We plot the energies in eq 62 in Figure 7 for a single set of HLi radii and ICOMP=2 (the same level as used for Figure 4). We note that the purpose of these calculations is to show the conceptual meaning of the terms rather than to evaluate them accurately at a high level of theory. Moreover, we do not explicitly discuss the components of the desolvation penalty (see eq 29), although it can be easily done. The cavitation term is contributes only about 12 kcal/mol, so the bulk of ΔEdesolv I given by the change in the electrostatic component, ΔEes I(I) (and the dispersion and repulsion terms are not computed, because of the two cavities at longer distances). Figure 7 reveals that the total formation energy is quite small, on the order of 17 kcal/mol, which is obtained from several large contributions of the opposite sign. The PIE in this system alone has the largest value of about 89 kcal/mol, which is very different from ΔE, pointing to the failure of the idea of looking only at PIEs as a measure of the binding energy. The monomer contribution which is strongly positive (repulsive) cancels most of the PIE resulting in the small ΔE. ΔEPLd I

Figure 7. Energetics of solvated Na+ 3 3 3 Cl at the RHF/6-31G*/ C-PCM/HLi/ICOMP=2/FIXPVA level as a function of the interatomic distance R. I = 1 refers to Na+ and J = 2 to Cl. The insert shows the close-up at 4.56.5 Å.

models agree within 2 kcal/mol. For the interatomic distances of 26 Å, the solvation component of PIE, ΔEsolv IJ and solv = ΔEin as the result, the total PIE/PCM, ΔEPCM IJ IJ + ΔE IJ show quite considerable dependence (up to 16.2 kcal/mol) upon the radii, but are very weakly affected by the charge renormalization (less than 0.3 kcal/mol). For R > 5.5 Å, for all models differ by less than 1 kcal/mol. The large ΔEPCM IJ electrostatic interaction ΔEin IJ at 10 Å (33.207 kcal/mol at all three levels in Figure 6, in accord with the simple predic= QIQJ/RIJ = 33.206 kcal/mol) is effectition ΔEgas IJ values between 0.4 vely screened to yield the ΔEPCM IJ and 0.5 kcal/mol. Thus, we see that the charge renormalization is important only at large distances to get the small total PIE accurately, which determines the local dielectric constant, but as the value of PIE itself is small, for practical applications of PIEs, for example, to drug design, the charge compensation appears to be relatively unimportant. On the other hand, the radii can be very important at short distances; however, typically, charged groups in biomolecules do not come to each other as close as the region where the radii have a large effect, with the possible exception of salt bridges. The salt bridges, however, may be better to put into the same fragment, both for the purpose of reducing the charge transfer, 714



ARTICLE

Figure 8. Structures of the complex of aspartic acid and lysine (at 3 Å separation).

Interestingly, the sodium destabilization is negligible in this system, 1.0 kcal/mol or less, while that of chlorine is as large as 16.7 kcal/mol (both values are at 2 Å). This means that sodium (hard ion) is mainly polarized by the solvent, and the polarization change along the reaction coordinate is small. On the opposite, chlorine (soft ion) perhaps because of its diffusive anionic electron distribution is strongly destabilized by the presence of sodium at shorter distances of 3 Å or less; we also note that the destabilization polarization energy also reflects the change in the solvent-induced polarization. The desolvation penalties of both sodium and chlorine are large, with the former being the larger, determined by the larger extent in which sodium looses its solvent surface. We observe the same small energy barrier49 in the vicinity of 5.7 Å (clearly seen as a wiggle in ΔE in the insert of Figure 7), which may be related to the first solvation shell restructurization, as described above. 4.2. Prototypical Biochemical Screening: Aspartic Acid and Lysine. In the previous subsection, we investigated an inorganic system consisting of two ions. Now we consider the system of two interacting charged molecules, which have a complex nature of the charge delocalization and renormalization (by ICOMP=2); also, it is of interest to see the effect of hydrogen atoms. We apply the RHF/6-31G* level with several PCM types. The structures were generated as follows. First, we optimized the two NH2-CHR-COOH molecules separately at the MP2/631+G* level in gas phase: deprotonated aspartic acid, R = -CH2COO and protonated lysine, R= -(CH2)4-NH3+. Then we reoriented them so that the CC bond in CCOO and the CN bond in CNH3+ are aligned along the z-axis, and the two molecules face each other with the -COO and -NH3+ groups (Figure 8). The reaction coordinate R was defined as the distance between the carbon atom of -COO and the nitrogen atom of NH3+. R was varied between 3 and 10 Å. We compare the results for the van der Waals (vdW) and SUAHF radii using PCM[1], for ICOMP=0 and ICOMP=2. Also, we compare PCM[1] to PCM[1(2)] for one combination (vdW, ICOMP=2). The vdW atomic radii are 1.2, 1.7, 1.6, and 1.5 Å for H, C, N, and O atoms, respectively. For SUAHF, the respective radii are 0.01, 1.77, 1.68, and 1.59 Å (to be multiplied by α = 1.2 in both radii sets). C-PCM/FIXPVA is used with 60 tesserae per sphere. First, we investigate the properties at 100 Å. The effective dielectric constants εIJ according to eq 48 are 94.69, 36.93, 150.27, and 78.37 for vdW/ICOMP=2, vdW/ICOMP=0, SUAHF/ICOMP=2, and SUAHF/ICOMP=0, respectively. The values for PCM[1] and PCM[1(2)] calculated for vdW/ ICOMP=2 differed by 0.004. We see that ICOMP at large distances has a strong effect on εIJ, and so do the atomic radii,

Figure 9. Dielectric constants in the complex of aspartic acid and lysine, for van der Waals (vdW) and SUAHF atomic radii, with ICOMP=0 and ICOMP=2 (as indicated after slash on the legend). PCM[1] is used, except vdW/2a, where PCM[1(2)] is employed: (a) effective εIJ, (b) internal εI. For ICOMP=2, εI = εJ, otherwise, both values for lysine and aspartic acid are shown (Lys and Asp on the legend). The insets show the close-ups in the range of 35 Å.

while the level of PCM[n] does not. Next, we look at the monomer properties at these levels, shown in Table 2. The PCM[1] and PCM[1(2)] results were the same, so the latter are not shown. The internal εI agree quite well with the values of the effective dielectric constant, except for SUAHF/ICOMP=0. In comparing vdW and SUAHF results for ICOMP=2, we see that the induced charges and εI, as well as the surface area σes, are larger for the latter because the radii of heavy atoms (C, N, and O) are larger. For SUAHF/ICOMP=0, the internal εI for Lys and Asp differ tremendously. It is interesting to observe that, according to eq 59, the effective pair ε from εI = 415.93 and εJ = 36.25 is predicted to be 79.41, which agrees well with the value of 78.37 obtained from eq 48. This peculiar behavior must be related to the fact that the -NH3+ group does not have spheres around hydrogens in SUAHF and it is solvated by a tight sphere determined by the nitrogen radius. Thus, an overestimated screening appears to take place in charged functional groups with hydrogens. The lack of the charge renormalization appears to be a problem in this case with a negative dielectric constant, but ICOMP=2 results in some averaged ε closer to the expected value. On the other hand, the vdW radii show very similar charges qI and qJ, as well as εI and εJ for ICOMP=0 versus ICOMP=2, which indicates less severe charge renormalization problems. The SUAHF monomer solute solvent ES energies ΔEes I(I) for lysine are much larger than those for vdW, also apparently caused by the overestimated electrostatic 715



ARTICLE

2 Table 2. Cavity Surfaces σes I (Å ), Induced Charges qI, Internal Dielectric Constants εI, and Monomer Solvation Energy a Components ΔE (kcal/mol): The RHF/6-31G*/C-PCM/FIXPVA Model Applied to Lysine (LYS) and Aspartic Acid (ASP), Separated by 100 Å

radii

ICOMP

vdW

2

vdW

a

0

SUAHF

2

SUAHF

0

σes I

qI

LYS

135.1

0.9900

99.99

72.3

26.0

46.3

ASP

118.0

0.9900

99.99

74.3

18.4

55.8

fragment I

εI

ΔEes I(I)

ΔEcav I

ΔEsolv I

LYS

135.1

0.9733

37.39

71.0

26.0

45.0

ASP

118.0

0.9739

38.28

72.9

18.4

54.5

LYS ASP

170.1 138.0

0.9933 0.9933

149.28 149.28

87.2 75.1

20.1 16.4

67.1 58.7

LYS

170.1

1.0024

415.93

87.2

20.1

67.1

ASP

138.0

0.9724

36.25

73.9

16.4

57.5

Defined in the main text.

Table 3. Dependence of the RHF/6-31G* Pair Interaction Energy (kcal/mol) between Lysine and Aspartic Acid (ΔEPCM for PCM IJ a or ΔEint IJ for Gas Phase) upon the Separation R (Å) and Calculation Details

a

PCM[1] vdW/2b

PCM[1(2)] vdW/2b

PCM[1] vdW/0b

PCM[1] SUAHF/2b

PCM[1] SUAHF/0b

R

gas phase

3.0

120.0

64.1

66.0

64.8

56.2

56.6

4.0

89.3

31.4

31.9

32.2

22.4

22.8

5.0

67.9

9.1

9.1

9.9

6.2

6.6

6.0 7.0

55.6 47.2

2.2 1.5

2.2 1.5

2.9 2.2

0.6 0.5

0.9 0.7

8.0

41.1

1.2

1.2

1.7

0.4

0.6

9.0

36.4

1.0

1.0

1.5

0.3

0.5

10.0

32.7

0.8

0.8

1.3

0.3

0.5

C-PCM/FIXPVA is used. b Atomic radii/ICOMP.

constant differs by not more than 0.05. On the other hand, the εI and εIJ curves separate into two distinctive vdW and SUAHF groups for R < 5 Å (see the insets in Figure 9). At longer distances, the induced-charge derived εI for SUAHF/ICOMP=0 goes to positive infinity (see Figure 9b; for example, at R = 7 and 10 Å, the εI values are 764.55 and 3475.14, respectively, and the latter is not shown in Figure 9b), followed by negative values (see Table 2 for R = 100 Å, the εI value is 415.93) caused by the proximity of the induced charges to 1, so that a small variation in the induced charges makes a large effect upon the dielectric constant because the denominator in eq 57 approaches zero. A resonance of εI takes place if qI becomes equal to QI. However, the averaged value of ε according to eq 59 is reasonable for the separations with very large εI, as was the case at 100 Å (e.g., at 8 Å, εI = 1156.07, εJ = 39.73, and ε is 76.81). Probably, SUAHF with ICOMP=0 gave very good dielectric constants in this system by coincidence. Practically, the long-range behavior is not very important in terms of PIEs because their values are small. As the data in Table 3 show, PIEs depend on the atomic radii; with the difference of about 10 kcal/mol or less (about 15% of the PIE value), the dependence on the ICOMP method or the PCM[n] level is much smaller, 1.2 kcal/mol or less. 4.3. Biochemical Application: Chignolin. As a demonstrative biochemical application, we chose chignolin (PDB: 1UAO), divided into 10 fragments (Figure 10). The structure was obtained from the previous study,31 where it was optimized at the FMO/6-31G(+)*/PCM/TIP3P level (TIP3P was used for explicit solvent, water). We apply the analysis in this work at the

Figure 10. Chignolin (PDB: 1UAO) divided into 10 fragments (fragments differ from residues by the carboxyl group).

interaction because of the absence of hydrogen spheres. In any case, one can conclude that for both this system and NaCl, the dielectric constant in PCM is very radius- and ICOMP-dependent. Next, we varied the distance between lysine and aspartic acid and studied the dependence of the pair interaction energy (Table 3) and the effective dielectric constants (Figure 9). Comparing PCM[1] and PCM[1(2)], the effective dielectric 716



ARTICLE

Figure 11. (a) Components of the pair interaction energies in chignolin at the MP2/6-31G* with C-PCM/ICOMP=2/FIXPVA level of theory. On the x-axis, the I,J pairs are shown: gas phase, solvent, and PCM components are shown as diamonds, squares and triangles, respectively. (b) The difference (shown with circles) in the PCM pair interaction energies ΔEPCM for ICOMP=2 minus those with ICOMP=0. IJ

internal εI, as was the case in two-fragment systems with zero total charge, see eqs 55 and 56. The pair interaction energies in the gas phase and PCM are is not exactly equal to ΔEgas plotted in Figure 11 (note ΔEPCM IJ IJ + solv ΔEIJ ). The nearest neighbor interactions for the fragment pairs between which a bond is detached are not shown because they include the large fragmentation correction.60 The “loose ends flat middle” pattern of the PIEs seems to be reasonable. One can see that the charged residue fragments in this system, the N-terminus (I = 1), Asp-3, Glu-5, and the C-terminus (I = 10), have the large interaction energies in the gas phase, which are made small in solution, by the solvation component ΔEsolv IJ of the opposite sign. The solutesolvent dispersion and repulsion components of PIE in this system did not exceed 2.0 and 0.1 kcal/mol, respectively. An interesting PIE is found for the pair 6,3 (Thr-6 and Asp-3). The gas phase and PCM values are 15.2 and 22.5 kcal/mol, respectively. Thus, the solvent effects here lead to an increase in PIE, which is almost entirely due to the induced charges determining ΔEes2 IJ . The neutral Thr-6 has the induced charge of q6 = 0.3223, while the negatively charged Asp-3 has the charge q3 = 0.6431. The fact that a neutral fragment has a considerable induced charge is attributed to the many body effects in inducing charges (see eqs 7 and 8). Thus, the screening increases the PIE because of these two induced charges. The largest electrostatic interaction between charged fragments in solution is for the pair 1,3 (the N-terminus fragment, + NH3-CαH2 and Asp-3). Fragment 1 is rather small and some screening of it is counted in the interaction for its neighbor, fragment 2. That is, the charge of fragment 1 induces ASCs on the cavities of its own fragment and also on some atoms in fragment 2, for example, qα for atom 3 (the closest atom in fragment 2)

MP2/6-31G* level with C-PCM[1]/FIXPVA/vdW for the fragmentation of one residue per fragment This calculation of PIEDA/PCM, which is the longest in this work, took 16.8 h on a single core Pentium4 3.2 GHz node. Before investigating the results, it is necessary to introduce a mental picture. In all other tests there were only two fragments, so the surface charges on each are induced by its own, as well as its neighbor’s potential. Therefore, if the two fragments are well separated, then only each fragment itself determines its induced charges, and therefore, the internal dielectric constant model in eq 57 works well. However, in polyfragment systems, the same reasoning applies only if a fragment is well separated from the rest. But, typically, one has some fragments close and some far from a given fragment. Probably, in general, neutral fragments have a fairly weak effect upon the induced charges, so a more complicated behavior should be typified by the case of three charged fragments. The effective dielectric constant in this case is largely determined by the strong interactions at a short distance. If we consider two close fragments I and J, and another K which is far from both, then the interactions IK and JK are largely affected by the presence of the third fragment J and I, respectively, and the effective dielectric constant is not determined by the simple relation in eqs 57 and 59 for the pairs IK and JK. Even a neutral fragment J if close can exert a considerable effect on the effective dielectric constant for IK, because it also contributes to the total potential, determining the induced charges. Thus, we see that many-body effects in the induced charge determination complicate and obscure the relatively simple picture in a two-fragment system. In addition, the charge renormalization does not lead to the equality of all 717


The Journal of Physical Chemistry A is 0.028. This is, however, a small fraction of the induced charge on the N-terminus, q1 = 0.8399. Another reason is the manybody effects leading to the induced cavity charges qI. As mentioned above, the charge on Asp-3 (0.6431) is much smaller than the charge, which would be if Asp-3 was the only fragment in the system (i.e., according to eq 51, such charge is 0.9872). Therefore, the screening is incomplete, mainly because of Asp-3 and its environment. can be The effects of the charge renormalization upon ΔEPCM IJ seen in the lower part of Figure 11. There is a strong correlation (Figure 11a) and the difference between the values of ΔEsolv IJ for the ICOMP=2 and ICOMP=0 between the values of ΔEPCM IJ (Figure 11b). Thus, the charge compensation scheme in general made a small scaling difference to the pair interactions, with the absolute magnitude of less than 0.5 kcal/mol. The charge renormalization brings back the escaped charge, that is, electronic and nuclear charges are scaled separately with a constant slightly more than one (e.g., in this system they were 1.00893 and 1.00892, respectively). Therefore, a small proportional increase in the absolute magnitude of ΔEsolv IJ is typically observed for both signs of it. Summarizing, the escaped charge takes along a small part of the screening, and using ICOMP=2 increases the screening by a small amount.

ARTICLE

upon the screening of the long-range interactions can add up to considerable values because of their large number. The complex formation of NaCl in solution has been analyzed, and the observed induced charge distribution and the dielectric constants have been explained with the induced charge quenching effect caused by the partial loss of the solvation shell and the potential pressure of the nearby ion. Our findings, including the effect of using zero as the hydrogen radius, may be useful in a future refinement of atomic radii in PCM or other continuum solvation methods. It is our hope that with the screening FMO will be put into a new level of usefulness in biochemical applications.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT We thank the Next Generation SuperComputing Project, Nanoscience Program (MEXT, Japan) for the financial support. ’ REFERENCES

5. CONCLUSIONS We have developed an ab initio based method for the determination of the solvent screening of the interactions, which does not require any parameters (apart from the usual PCM parameters). In particular, it does not involve a dielectric constant for the solute, and only the universal dielectric constant for the solvent is used in the determination of the solute-induced solvent charges. It is clear that these induced charges on the solvent govern the screening and, thus, the interactions in the solvated systems. In this work we have introduced the mathematical foundations for the screening in PCM and developed the PIEDA/PCM approach. We have demonstrated that it does describe the effective reduction of the electrostatic interaction in solution, as determined in the model test calculations of the ionic and molecular systems. Also, the relevance of the monomer terms related to their polarization and solvation in the binding process has been stressed. Three types of the effective screening can be named: regular screening ε > 1, negative screening ε < 0, and antiscreening 0 < ε < 1. The regular screening decreases the electrostatic interaction, the negative screening changes its sign, while the antiscreening increases the interaction. ε = 1 corresponds to no screening, while ε = 0 is not well-defined. The origin of the antiscreening is the many-body nature of the determination of the induced charges in PCM. These charges as well as the charge renormalization problem can be thought to be the origins of the negative screening. The charge renormalization appears to have a shifting effect, that is, it shifts the total value of the screened interaction by a small amount on the order of 1 kcal/mol or less; this is not important at short distances, but has a large relative effect otherwise, leading to the effective dielectric constant very different from its proper value. On the other hand, the choice of the atomic radii is important for the whole range of separations. However, the effect of the particular model (such as ICOMP)

(1) Gordon, M. S.; Mullin, J. M.; Pruitt, S. R.; Roskop, L. B.; Slipchenko, L. V.; Boatz, J. A. J. Phys. Chem. B 2009, 113, 9646–9663. (2) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999–3093. (3) Casanova, D.; Gusarov, S.; Kovalenko, A.; Ziegler, T. J. Chem. Theory Comput. 2007, 3, 458–476. (4) Li, H.; Pomelli, C. S.; Jensen, J. H. Theor. Chim. Acta 2003, 109, 71–84. (5) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A 2006, 110, 2493–2499. (6) Sanchez, M. L.; Martín, M. E.; Aguilar, M.; Olivares del Valle, F. J. Comput. Chem. 2000, 21, 705–715. (7) Wong, S. E.; Lightstone, F. C. Exp. Opin. Drug Discovery 2011, 6, 65–74. (8) Scuseria, G. E. J. Phys. Chem. A 1999, 103, 4782–4790. (9) Nikitina, E.; Sulimov, V.; Zayets, V.; Zaitseva, N. Int. J. Quantum Chem. 2004, 97, 747–763. (10) Stewart, J. J. P. J. Mol. Model. 2009, 15, 765–805. (11) Mitani, M.; Aoki, Y.; Imamura, A. Int. J. Quantum Chem. 1997, 64, 301–323. (12) Gao, J. J. Phys. Chem. B 1997, 101, 657–663. (13) Mei, Y.; Ji, C. G.; Zhang, J. Z. H. J. Chem. Phys. 2006, 125, 094906. (14) Mahadevi, A. S.; Rahalkar, A. P.; Gadre, S. R.; Sastry., G. N. J. Chem. Phys. 2010, 133, 164308. (15) Sorkin, A.; Dahlke, E. E.; Truhlar, D. G. J. Chem. Theory Comput. 2008, 4, 683–688. (16) Hua, W. J.; Fang, T.; Li, W.; Yu, J. G.; Li, S. H. J. Phys. Chem. A 2008, 112, 10864–10872. (17) S€oderhjelm, P.; Ryde, U. J. Phys. Chem. A 2009, 113, 617–627. (18) Mata, R. A.; Stoll, H.; Cabral, B. J. C. J. Chem. Theory Comput. 2009, 5, 1829–1837. (19) Huang, L. L.; Massa, L. J. Mol. Struct.: THEOCHEM 2010, 962, 72–79. (20) Kobayashi, M.; Kunisada, T.; Akama, T.; Sakura, D.; Nakai, H. J. Chem. Phys. 2011, 134, 034105. (21) He, X.; Merz, K. M., Jr. J. Chem. Theory Comput. 2010, 6, 405–411. (22) Steinmann, C.; Fedorov, D. G.; Jensen, J. H. J. Phys. Chem. A 2010, 114, 8705–8712. 718


The Journal of Physical Chemistry A (23) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Chem. Rev. 2011, DOI: 10.1021/cr200093j. (24) Vreven, T.; Morokuma, K.; Farkas, O.; Schlegel, H. B.; Frisch, M. J. J. Comput. Chem. 2003, 24, 760–769. (25) Combined Quantum Mechanical and Molecular Mechanical Methods; Gao, J., Thompson, M. A., Eds.; ACS Symposium Series 712, Oxford University Press: U.K., 1998. (26) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Chem. Phys. Lett. 1999, 313, 701–706. (27) Fedorov, D. G.; Kitaura, K. J. Phys. Chem. A 2007, 111, 6904–6914. (28) The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems; Fedorov, D. G., Kitaura, K., Eds.; CRC Press: Boca Raton, FL, 2009. (29) Fedorov, D. G.; Kitaura, K.; Li, H.; Jensen, J. H.; Gordon, M. S. J. Comput. Chem. 2006, 27, 976–985. (30) Li, H.; Fedorov, D. G.; Nagata, T.; Kitaura, K.; Jensen, J. H.; Gordon, M. S. J. Comput. Chem. 2010, 31, 778–790. (31) Fedorov, D. G.; Ishida, T.; Uebayasi, M.; Kitaura, K. J. Phys. Chem. A 2007, 111, 2722–2732. (32) Ishikawa, T.; Yamamoto, N.; Kuwata, K. Chem. Phys. Lett. 2010, 500, 149–154. (33) Fedorov, D. G.; Alexeev, Y.; Kitaura, K. J. Phys. Chem. Lett. 2011, 2, 282–288. (34) Nagata, T.; Brorsen, K.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. J. Chem. Phys. 2011, 134, 124115. (35) Komeiji, Y.; Mochizuki, Y.; Nakano, T.; Fedorov, D. G. J. Mol. Struct.: THEOCHEM 2009, 898, 2–7. (36) Fujita, T.; Nakano, T.; Tanaka, S. Chem. Phys. Lett. 2011, 506, 112–116. (37) Milne, B. F.; Marques, M. A. L.; Nogueira, F. Phys. Chem. Chem. Phys. 2010, 12, 14285–14293. (38) Taguchi, N.; Mochizuki, Y.; Nakano, T. Chem. Phys. Lett. 2011, 504, 76–82. (39) Pruitt, S. R.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. J. Chem. Theory Comput. 2010, 6, 1–5. (40) Nakanishi, I.; Fedorov, D. G.; Kitaura, K. Proteins: Struct., Funct., Bioinf. 2007, 68, 145–158. (41) Sawada, T.; Fedorov, D. G.; Kitaura, K. J. Phys. Chem. B 2010, 114, 15700–15705. (42) Sawada, T.; Fedorov, D. G.; Kitaura, K. J. Am. Chem. Soc. 2010, 132, 1686–16872. (43) He, X.; Fusti-Molnar, L.; Cui, G.; Merz, K. M., Jr. J. Phys. Chem. B 2009, 113, 5290–5300. (44) Fujimura, K.; Sasabuchi, Y. Chem. Med. Chem. 2010, 5, 1254–1257. (45) Mazanetz, M. P.; Ichihara, O.; Law, R. J.; Whittaker, M. J. Cheminf. 2011, 3, 2. (46) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2004, 120, 6832–6840. (47) Watanabe, H.; Okiyama, Y.; Nakano, T.; Tanaka, S. Chem. Phys. Lett. 2010, 500, 116–119. (48) Nagata, T.; Fedorov, D. G.; Sawada, T.; Kitaura, K.; Gordon, M. S. J. Chem. Phys. 2011, 134, 034110. (49) Su, P.; Li, H. J. Chem. Phys. 2009, 130, 074109. (50) Si, D.; Li, H. J. Chem. Phys. 2009, 131, 044123. (51) Wang, Y.; Li, H. J. Chem. Phys. 2009, 131, 206101. (52) Weijo, V.; Randrianarivony, M.; Harbrecht, H.; Frediani, L. J. Comput. Chem. 2010, 31, 1469–1477. (53) Lipparini, F.; Scalmani, G.; Mennucci, B.; Cances, E.; Caricato, M.; Frisch, M. J. J. Chem. Phys. 2010, 133, 014106. (54) Lipparini, F.; Scalmani, G.; Mennucci, B.; Frisch, M. J. J. Chem. Theory Comput. 2011, 7, 610–617. (55) Su, P. F.; Li, H. J. Chem. Phys. 2009, 131, 014102, and the references therein. (56) Vaart, A.; van der; Merz, K. M., Jr. J. Am. Chem. Soc. 1999, 121, 9182–9190. (57) Ichihara, O.; Barker, J.; Law, R. J.; Whittaker, M. Mol. Inf. 2011, 30, 298–306. (58) Munei, Y.; Shimamoto, K.; Harada, M.; Yoshida, T.; Chuman, H. Bioorg. Med. Chem. Lett. 2011, 21, 141–144.

ARTICLE

(59) Mori, H.; Ueno-Noto, K. J. Phys. Chem. B 2011, 115, 4774–4780. (60) Fedorov, D. G.; Kitaura, K. J. Comput. Chem. 2007, 28, 222–237. (61) Ishikawa, T.; Ishikura, T.; Kuwata, K. J. Comput. Chem. 2009, 30, 2594–2601. (62) Okiyama, Y.; Fukuzawa, K.; Yamada, H.; Mochizuki, Y.; Nakano, T.; Tanaka, S. Chem. Phys. Lett. 2011, 509, 67–71. (63) Fedorov, M. V.; Kornyshev, A. A Mol. Phys. 2007, 105, 1–16. (64) Nagata, T.; Fedorov, D. G.; Kitaura., K. Mathematical Formulation of the fragment molecular orbital method. In Linear-Scaling Techniques in Computational Chemistry and Physics; Zalesny, R., Papadopoulos, M. G., Mezey, P. G., Leszczynski, J., Eds.; Springer: New York, 2011; pp 1764. (65) Fedorov, D. G.; Kitaura, K. Chem. Phys. Lett. 2006, 433, 182–187. (66) Cammi, R.; Tomasi, J. J. Comput. Chem. 2004, 16, 1449–1458. (67) Finkelstein, A. V.; Ptitsyn, O. B. Protein Physics: A Course of Lectures; Academic Press: Amsterdam, 2002. (68) Altman, M. D.; Bardhan, J. P.; White, J. K.; Tidor, B. J. Comput. Chem. 2009, 30, 132–153. (69) Teixeira, V. H.; Cunha, C. A.; Machuqueiro, M.; Oliveira, A. S. F.; Victor, B. L.; Soares, C. M.; Baptista, A. M. J. Phys. Chem. B 2005, 109, 14691–14706. (70) Pierotti, R. A. Chem. Rev. 1976, 76, 717–726. (71) Langlet, J.; Claverie, P.; Caillet, J.; Pullman, A. J. Phys. Chem. 1988, 92, 1617–1631. (72) Li, H.; Jensen, J. H. J. Comput. Chem. 2004, 25, 1449–1462. (73) Klamt, A. J. Phys. Chem. 1995, 99, 2224–2235. (74) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. J. Comput. Chem. 2003, 24, 669–681. (75) Cances, E.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 107, 3032–3041. (76) Chipman, D. M. J. Chem. Phys. 1997, 106, 10194–10206. (77) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., Jr. J. Comput. Chem. 1993, 14, 1347–1363. (78) Fedorov, D. G.; Olson, R. M.; Kitaura, K.; Gordon, M. S.; Koseki, S. J. Comput. Chem. 2004, 25, 872–880. (79) Mu, W.; Chasse, G. A.; Fang, D. Int. J. Quantum Chem. 2008, 108, 1422–1434. (80) Skipper, N. T.; Neilson, G. W. J. Phys.: Condens. Matter 1989, 1, 4141–4154. (81) Kameda, Y.; Sugawara, K.; Usuki, T.; Uemura, O. Bull. Chem. Soc. Jpn. 1998, 71, 2769–2776. (82) Ikeda, T.; Boero, M.; Terakura, K. J. Chem. Phys. 2007, 126, 034501. (83) Powell, D. H.; Neilson, G. W.; Enderby, J. E. J. Phys.: Condens. Matter 1993, 5, 5723–5770. (84) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757–3766. (85) Barone, V.; Cossi, M.; Tomasi, J. J. Chem. Phys. 1997, 107, 3210–3221. (86) Kitaura, K.; Morokuma, K. Int. J. Quantum Chem. 1976, 10, 325–340. (87) Camaioni, D. M.; Dupuis, M.; Bentley, J. J. Phys. Chem. A 2003, 107, 5778–5788. (88) Ginovska, B.; Camaioni, D. M.; Dupuis, M.; Schwerdtfeger, C. A.; Gil, Q. J. Phys. Chem. A 2008, 112, 10604–10613. (89) Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 106, 5151–5158. (90) Marcus, Y.; Hefter, G. Chem. Rev. 2006, 106, 4585–4621. (91) Murata, K.; Fedorov, D. G.; Nakanishi, I.; Kitaura, K. J. Phys. Chem. B 2009, 113, 809–817.

’ NOTE ADDED AFTER ASAP PUBLICATION This article was published ASAP on December 8, 2011, with an incomplete version of reference 23. The correct version was published ASAP on December 14, 2011. 719


Energy Decomposition Analysis in Solution Based on the Fragment

Recommend Documents