Pair Interaction Energy Decomposition Analysis for Density Functional

Subscriber access provided by Gothenburg University Library

Article

Pair Interaction Energy Decomposition Analysis for Density Functional Theory and Density-Functional Tight-Binding with an Evaluation of Energy Fluctuations in Molecular Dynamics Dmitri G. Fedorov, and Kazuo Kitaura J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b12000 • Publication Date (Web): 16 Jan 2018 Downloaded from http://pubs.acs.org on January 16, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Pair Interaction Energy Decomposition Analysis for Density Functional Theory and Density-Functional Tight-Binding with an Evaluation of Energy Fluctuations in Molecular Dynamics Dmitri G. Fedorov1,*, Kazuo Kitaura2,3 1

Research Center for Computational Design of Advanced Functional Materials (CD-FMat),

National Institute of Advanced Industrial Science and Technology (AIST), Central 2, Umezono 1-1-1, Tsukuba, 305-8568, Japan. 2

Advanced Institute for Computational Science (AICS), RIKEN, 7-1-26 Minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan.

3

Fukui Institute for Fundamental Chemistry, Kyoto University, Takano-Nishihiraki-cho 34-4, Sakyou-ku, Kyoto 606-8103, Japan.

* AUTHOR EMAIL ADDRES S: [email protected]

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 61

ABSTRACT. Pair interaction energy decomposition analysis in the fragment molecular orbital (FMO) method is extended to treat density functional theory (DFT) and density-functional tightbinding (DFTB). Fluctuations of energy contributions are obtained from molecular dynamics simulations. Interactions at the DFT and DFTB levels are compared to the values obtained with Hartree-Fock, second order Møller-Plesset (MP2) and coupled cluster methods. Hydrogen bonding in water clusters is analyzed. 200 ps NVT molecular dynamics simulations are performed with FMO for two ligands bound to the Trp-cage miniprotein (PDB: 1L2Y); the fluctuations of fragment energies and interactions are analyzed.

1. Introduction The fragment-based1,2,3,4,5,6,7,8,9,10,11,12,13,14 and low scaling 15,16,17,18,19,20,21,22 methods form a vibrant field of modern research23,24 taking advantage of the locality of the electron density distribution in large molecular systems. These methods can provide useful insights into the properties of fragments and interactions between them. The fragment molecular orbital method (FMO) 25,26,27,28,29 is a fragmentation approach, in which the electronic state of fragments and their pairs is determined in the electrostatic embedding of the whole system, and the total properties are computed from those of fragments, their pairs, and, optionally trimers or even tetramers. 30 FMO has been successfully applied to a large number of various molecular systems, such as organic polymers, 31 sugars, inorganic materials.

36

32

DNA,

33

proteins,

34

molecular clusters and crystals 35 and

Recently, molecular dynamics simulations (MD) using FMO are

becoming accessible. 37,38,39,40 FMO is suitable to study properties of subsystems and interactions.


2

Page 3 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


41,42

Various interaction analyses have been designed to analyze pair interaction energies (PIE).

43,44,45,46

Fluctuations of PIEs in MD have been studied by analyzing snapshots. 47,48

PIEs from FMO have been used in many applications, in particular, in drug discovery49 and structure activity relationship (SAR) studies. 50 , 51 , 52 PIEs can be used as descriptors of protein-ligand binding, which correlate with experimentally measured binding.34 FMO is also used in fragment-based drug discovery.53,54

2. Methodology There are various energy decomposition analyses (EDA), 55 , 56 , 57 many of which are derived from the original Kitaura-Morokuma scheme.

58 , 59

EDA developed in the FMO

framework is called the pair interaction energy decomposition analysis (PIEDA). 60 , 61 Pair interaction energies can be decomposed into 5 components: electrostatics (ES), exchange repulsion (EX), charge transfer and higher order mix terms (CT+mix), dispersion (DI) and solvent screening (SOLV). ∆E IJ = ∆E IJES + ∆E IJEX + ∆E IJCT + mix + ∆E IJDI + ∆E IJSOLV

(1)

PIEDA has been developed for restricted closed60 and open 62 shell Hartree-Fock and second order Møller-Plesset (MP2) perturbation theory. Here, PIEDA is formulated for two other methods, density functional theory (DFT) 63 and density-functional tight-binding (DFTB), 64,65 interfaced with the polarized continuum model (PCM) to treat solvent effects.66,67 In this work, PIEDA is also applied for the first time to FMO interfaced68 with coupled cluster (CC) with singles, doubles and perturbative triples (CCSD(T)). FMO-DFTB has been used in some biochemical applications; 69 PIE computed with DFTB highly correlate with the values from DFT67 and MP2. 70 It has also been shown70 that


3


Page 4 of 61

atomic charges in DFTB and RHF correlate. Geometry optimizations with FMO-DFTB have delivered reasonable structures of cellulose 71 and proteins in solution,67 as well as radial distribution functions for several systems.38,39 Many biological and other materials are flexible, and the question of what structure to use for an interaction analysis is both difficult and important. One can use a selection of structures from a classical MD to obtain some averaged PIEs.72 Alternatively, a statistical correction model based on a parameter has been suggested. 73 By performing FMO-MD simulations, one can obtain the lowest energy structure (referred to as MDmin), which can be compared to an energy minimum obtained by a normal geometry optimization at 0K (0Kmin). The internal energy U is defined 74 as an ensemble average of the total energy, given by the sum of the potential (pot) Epot and kinetic (kin) Ekin energies:

U = E = E pot + E kin

(2)

QM energies E pot are huge, for example, about 50000 kcal/mol for 1 water molecule. It is useful to subtract the bulk of the energy from its fluctuations.

E = E 0 + ∆E = E kin,0 + E pot,0 + ∆E kin + ∆E pot

(3)

Here, E 0 is the energy of a reference structure, E 0 = E kin,0 + E pot,0 . It is convenient to define it as the structure with the lowest potential energy E pot in MD, but other choices are also possible. The trajectory average (MDaver) is computed over M points representing a time grid i, ∆E =

1 M

M

∑ ∆E

i

(4)

i =1

The potential energy in FMO is computed as,


4

Page 5 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


~ E pot = ∑ E I + ∑ ∆E IJ I

(5)

I >J

~ where E I is the energy of fragment I, divided into the solute internal energy E I′′ and the ~ solvation energy ∆E Isolv , E I = E I′′ + ∆E Isolv . For monomers (X=I) or dimers (X=IJ), V X is the

electrostatic potential due to other fragments and D X is the electron density. The internal solute energy E ′X′ is computed from the total solute energy E ′X of X excluding the embedding energy, E ′X′ = E ′X − E VX

(6)

(

)

where E VX is the embedding energy, E VX = Tr D X V X and for DFTB it is given elsewhere.38, 71 The total PIE is computed as ∆E IJ = (E IJ′′ − E I′′ − E ′J′ ) + ∆E VX + ∆E IJSOLV

(

where ∆D IJ = D IJ − D I ⊕ D J

)

(7)

is the density transfer matrix and ∆E IJSOLV is the solvent

(

)

contribution defined elsewhere.61 Except for DFTB, ∆E VX = Tr ∆D IJ V IJ and for DFTB it is given elsewhere.38, 71 In PIEDA, prior to this work developed for RHF and MP2, PIEs ∆E IJ are decomposed according to eq 1. In this work, PIEDA is extended to treat DFT and DFTB, which requires a reformulation as follows. The electrostatic component of the pair interaction energy between fragments I and J is computed in PIEDA following the far-separated dimer approximation.77,65 In DFT, it is

(

)

(

) ∑ ∑D

∆E IJES = Tr D I u I ( J ) + Tr D J u J ( I ) +

I

µν

J D ρσ (µν ρσ ) + ∆E IJNR

(8)

µν ∈I ρσ ∈J


5


Page 6 of 61

where u I ( J ) is the matrix of one-electron Coulomb integrals for the nuclei of fragment J and the basis functions (labeled by µ, ν, ρ and σ) of I. ∆E IJNR is the nuclear repulsion (NR) between nuclei in I and J. In DFTB3, a point charge model is used, 71

(

)

1   ∆E IJES = ∑ ∑ γ AB ∆q AI ∆q BJ + ΓAB ∆q AI + ΓBA ∆q BJ ∆q AI ∆q BJ  3  A∈I B∈J  where γ AB

(9)

and ΓAB are the damping functions, related to chemical hardness and inverse

distances. In DFTB2, the ΓAB terms are neglected ( ΓAB = ΓBA = 0 ). ∆q AI is the Mulliken charge of atom A in fragment I. In DFT, the exchange-repulsion between occupied orbitals of two fragments due to the Pauli exclusion principle, can be computed as ∆E IJEX = E IJ′′′ EX − E I′′′ − E J′′′ − ∆E IJES

(10)

where E ′X′′ excludes the DFT correlation energy E Xc (but includes the DFT exchange energy): E ′X′′ = E ′X′ − E Xc

(11)

and EX in E IJ′′′ EX indicates that the energy of dimer IJ is obtained by doing a calculation while allowing mixing only between occupied molecular orbitals of fragments I and J (in practice, it is obtained by orthogonalizing occupied orbitals of I and J and performing one SCF iteration60). For DFTB, ∆E IJEX = 0 . The charge transfer and mix term (CT+MIX) is obtained in DFT as ∆E IJCT + MIX = ∆E IJ − ∆E IJES − ∆E IJEX − ∆E IJc − ∆E IJSOLV

(12)

The correlation functional energy E Xc in DFT is used to evaluate the correlation functional contribution to the interactions,57


6

Page 7 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∆E IJc = E IJc − E Ic − E Jc

(13)

In RHF and DFT(B), the dispersion interaction is computed using the empirical dispersion energy E XDI of dimers and monomers, 75 ∆E IJDI = E IJDI − E IDI − E JDI

(14)

The contribution of the electron correlation (CORR) to the interaction between fragments has the dispersion (DI) and non-dispersion components. The latter term, labeled remainder correlation (RC), represents the coupling of the electron correlation to the charge transfer between fragments (the extent to which the charge transfer affects the electron correlation). Some discussion on the relation between the MP2 correlation energy and dispersion can be found elsewhere. 76 The total contribution of the electron correlation to the interaction between fragments is DI+RC. However, the DFT correlation functional (which explicitly defines RC) can in principle contain some dispersion interaction. It should be noted that for simplicity, double hybrid functionals are not considered, although the method can be generalized to include them. The DI part of DI+RC in DFT is defined in eq 14 and the total correlation energy in DFT is, ∆E IJCORR = ∆E IJDI + RC = ∆E IJDI + ∆E IJc

(15)

In RHF or DFTB, there is only the dispersion contribution, ∆E IJDI + RC = ∆E IJDI . In MP2 and CC, the correlation contribution to the interaction energy, computed from the correlation energies of monomers and dimers, defines this term, ∆E IJDI + RC = ∆E IJCORR = E IJCORR − E ICORR − E JCORR

(16)

In this work, eq 1 is generalized: instead of the dispersion DI, DI+RC is used. ∆E IJ = ∆E IJES + ∆E IJEX + ∆E IJCT + mix + ∆E IJDI+ RC + ∆E IJSOLV

(17)

Decomposing the internal energy38,71 difference in eq 7 for DFTB, one obtains


7


Page 8 of 61

E IJ′′ − E I′′ − E ′J′ = ∆E IJ0 + ∆E IJREP + ∆E IJP + ∆E IJ′ ES + ∆E IJDI

(18)

where ∆E IJ0 = E IJ0 − E I0 − E J0 = Tr (∆D IJ H 0, IJ ) is the coupling of interfragment charge transfer ∆D IJ to the unperturbed electronic state (of atoms and diatomics). H 0, X is the zeroth order Hamiltonian

of

X

and

E X0 = Tr ( D X H 0, X )

.

The

repulsive

(REP)

energy

is

∆E IJREP = E IJREP − E IREP − E JREP . The difference in the projection operator38 P X energies is ∆E IJP = E IJP − E IP − E JP , where E XP = Tr ( D X P X ) . The difference in the internal electrostatic

energies is ∆E IJ′ ES = E IJ′ ES − E I′ ES − E J′ ES , where E ′XES =

1 E AES( X ) B ( X ) ∑ 2 A, B∈X

E AES( X ) B (Y ) = γ AB ∆q AX ∆q YB +

(19)

(

)

1 ΓAB ∆q AX + ΓBA ∆q BY ∆q AX ∆q BY 3

(20)

Using the symmetry of γ AB = γ BA and E AES( X ) B (Y ) = E BES(Y ) A( X ) , it can be shown that ∆E IJ′ ES =

1 ES, IJ ∆E AB + ∑ E AES( I ) B ( J ) = ∆E IJ′ CT⋅ES + ∆E IJES ∑ 2 A, B∈IJ A∈I , B∈J

(21)

∆E IJ′ CT⋅ES quantifies the effect of charge transfer on the electrostatics, computed as the sum of pairwise terms for atoms A, B ∈ IJ ES, IJ ∆E AB = E AES( IJ ) B ( IJ ) − E AES( I ′) B ( J ′)

(22)

where I’ and J’ are the fragments to which atoms A and B belong, respectively (when either A or B at the fragment boundary between I and J is shared by both fragments, eq 22 should be modified to have a subtraction of all relevant monomer values). In practice, ∆E IJ′CT⋅ES is calculated from eq 21 as the difference ∆E IJ′ ES − ∆E IJES . Using eq 7, one arrives at the final PIEDA expression for DFTB,


8

Page 9 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


~ ∆E IJ = ∆E IJES + ∆E IJ0 + ∆E IJCT⋅ES + ∆E IJDI + ∆E IJSOLV

(23)

~

where ∆E IJ0 = ∆E IJ0 + ∆E IJREP + ∆E IJP . The coupling of charge transfer to electrostatics is ∆E IJCT⋅ES = ∆E IJ′ CT⋅ES + ∆E IJV , where ∆E IJ′CT⋅ES and ∆E IJV is the energy of charge transfer in dimer IJ

coupled to the internal (charges in IJ) and external (charges outside IJ) embedding, respectively. The repulsive potential in DFTB is very short-ranged and, except for covalently connected fragments, ∆E IJREP is usually negligibly small. The projection energies E XP ≈ 0 and ∆E IJP ≈ 0 are also negligible60 because the electron density is such that Tr ( D X P X ) ≈ 0 (in test

calculations of polypeptides, the largest value of ∆E IJP was about 0.03 kcal/mol, for molecular

~ ~ clusters E XP = 0 ). Thus, ∆E IJ0 ≈ ∆E IJ0 . The sum ∆E IJ0 + ∆E IJCT⋅ES in DFTB formally corresponds ~

to ∆E IJEX + ∆E IJCT + MIX + ∆E IJc in DFT (compare eqs 12 and 23). ∆E IJ0 and ∆E IJCT⋅ES + ∆E IJES are non-

~ electrostatic (excluding dispersion) and electrostatic interactions in the solute, respectively. ∆E IJ0 is the interaction energy at the first-order DFTB (DFTB1) level (excluding the dispersion and solvent contributions). The solvent screening term61 in PCM is computed in all methods as ∆E IJSOLV = ∆E IJes + ∆E IJdisp + ∆E IJrep + ∆E IJCT ⋅es

(24)

where solvent-solute electrostatic (es), dispersion (disp) and repulsion (rep) interaction energies (Y=es, disp or rep) are obtained as ∆E IJY = ∆E IY( J ) + ∆E JY( I ) . ∆EIY( J ) is the Y-type interaction energy of solute fragment I with the solvent cavity section surrounding fragment J (see Figure 1). Solute-solute interactions are capitalized, for example, DI, whereas solute-solvent screening interactions are in small letters (e.g., disp). The last term in eq 24 describes the solute charge transfer coupling to the solvent electrostatic embedding,


9


(

1 ∆E IJCT⋅es = Tr ∆D IJ W IJ 2

Page 10 of 61

)

(25)

where W IJ is the matrix of the one-electron Coulomb integrals for basis functions of IJ due to the induced solvent charges.61 Only ∆E IJes and ∆E IJCT⋅es are method dependent (e.g., differ for DFT and DFTB), other components ∆E IJdisp and ∆E IJrep are parametrized and do not depend on the method. A summary of terms is given in Table 1. Bond detached atom (BDA) corrections60 can be used as in RHF; they are omitted from the above equations for simplicity. In this work, the fully polarized state (PL)

60

is used in PIEDA, and it is possible to reformulate equations for the

partially polarized state PL0. In the electrostatic dimer (ES-DIM) approximation77 in FMO, far separated dimers contribute only to the ES component. Introducing reference energies and fluctuations in eq 5, one obtains

~ ~ E pot = ∑ E I0 + ∑ ∆E IJ0 + ∑ ∆E I + ∑ ∆∆E IJ I

I >J

I

(26)

I >J

~ ~ ~ The values of ∆E I = E I − E I0 and ∆∆E IJ = ∆E IJ − ∆E IJ0

are the fluctuations in fragment

energies and PIEs, respectively. Another measure of fluctuations is given by the minimum and

{ }

{ }

~ ~ ~ ~ maximum values in MD, i.e., ∆E Imin = min i ∆E Ii , ∆E Imax = max i ∆E Ii , ∆E IJmin = min i {∆E IJi }

and ∆E IJmax = max i {∆E IJi }, where i numbers trajectory points in MD. In the comparative study of various methods for water clusters (H2O)n, averaged components are defined for Y=ES, EX, CT+mix, DI+RC, SOLV and REM values as ∆E Y =

1 HB ∑ ∆E IJY m I >J

(27)


10

Page 11 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where HB indicates that only dimers connected by a hydrogen bond (HB) are included in the sum (m is the number of hydrogen bonds). The physical meaning of fluctuations in monomer energies should be clarified. They encompass two kinds of effects: deformation and polarization; the two are coupled. For a fixed structure, monomer energies define the destabilization component of the polarization.60 In other words, depending on the environment surrounding a fragment, the fragment polarization changes. Secondly, fluctuations in the structure bring about a change in the deformation energies, that is, how the energy of the fragment changes with deformation of its own structure.

~ Thus, ∆E I

shows the degree of polarization and deformation of fragment I relative to

~ ~ its reference state (i.e., E I0 ). The values of ∆E I

for individual fragments may be both

repulsive and attractive. A separation of the deformation and polarization components is accomplished in the subsystem analysis;42 here, however, the total fluctuation in fragment energies describing the combined effect of the polarization and deformation is discussed. It should be noted that interaction energies in FMO are not the same as binding energies. The difference is the reference: for interaction energies the reference is the polarized state of fragments in the field of each other (in eq 7, the polarized fragment energies E I′′ and densities

D I (in ∆D IJ ) are subtracted); for binding energies the reference is the state of isolated noninteracting fragments (so that for binding energies one would subtract isolated fragment energies from the dimer energy). The two are related by the polarization energy and desolvation components and both can be evaluated in the subsystem analysis.42 Finally, the kinetic energy is discussed. The averaged kinetic energy is E kin = E kin,0 + ∆E kin = E kin,0 +

1 M

∑ (E M

kin, i

− E kin, 0

)

(28)

i =1


11


Page 12 of 61

where E kin,i and E kin,0 is the kinetic energy at time step i and for the reference structure, respectively. The total kinetic energy can be divided into a sum over fragments as

1 ∑ m A v A2 ,i = ∑I E Ikin,i 2 A

E kin,i =

(29)

where v A,i is the velocity of nucleus A with mass m A at the trajectory point i and

1 m A v A2 ,i ∑ 2 A∈I

E Ikin,i =

(30)

Averaging the fluctuation ∆E Ikin,i of the kinetic energy of fragment I, one obtains ∆E Ikin . ∆E Ikin,i = E Ikin,i − E Ikin, 0

(31)

Analogously to the definition of the local dielectric constant for a fragment in FMO,61 one can define an effective (“local”) temperature TI of fragment I, using the relation E Ikin = 3 N Iat kTI / 2 , where E Ikin = E Ikin,0 + ∆E Ikin , N Iat is the number of atoms in I and k is the Boltzmann constant: TI =

2 E Ikin 3 N Iat k

(32)

Here, TI is simply the normalized kinetic energy of fragment I per atom, having the dimension of temperature (K). The global temperature is given by

T=

2 E kin , at ∑ I 3N k I

(33)

where N at is the total number of atoms, so that the global temperature is the weighted average effective temperature of fragments.

T=

1 N at

∑N

at I

TI

(34)

I


12

Page 13 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The kinetic energy of fragments may be used to discuss dynamic aspects of proteinligand binding, showing how the kinetic energy of the ligand changes during MD simulations, which may be related to the dynamic behavior of drugs in vivo and drug retention and residence time, which are a hot topic in modern drug discovery.78 The concept of the effective temperature of a fragment can be used to identify “cold” and “hot” fragments, which may be related to their dynamic behavior relevant to biochemical activity. As a related quantity, B-factors (also called temperature factors) of atoms are used in X-ray crystallography. A higher rotational and vibrational freedom of a hotter fragment may be correlated to its larger entropy. In practical applications of FMO/MD to protein-ligand binding, it may be of interest to decompose the binding energy of a complex AB. In order to do this, one has to apply eq 2 to each term in the binding energy, ∆E AB = E AB − E A − E B

(35)

The internal binding energy is pot kin ∆U AB = E AB − E A − E B = ∆ E AB = ∆ E AB + ∆ E AB

(36)

where pot pot,0 pot ∆ E AB = ∆E AB + ∆ ∆E AB

(37)

kin kin,0 kin ∆ E AB = ∆E AB + ∆ ∆E AB

(38)

pot kin pot,0 kin,0 + ∆ ∆E AB ∆E AB + ∆E AB is the binding energy for the reference structures, and ∆ ∆E AB is

the difference in the fluctuations of the binding energy due to complexation. The contribution to the potential binding energy for the reference structures is (and, likewise, the kinetic contribution), pot,0 pot,0 ∆E AB = E AB − E Apot,0 − E Bpot,0

(39)


13


Page 14 of 61

The fluctuation in the potential energy contribution to the binding energy can be written using eq 26 (and, likewise, the kinetic fluctuation is defined), pot pot ∆ ∆E AB = ∆E AB − ∆E Apot − ∆E Bpot =

∑

I ∈AB

~ ∆E IAB +

∑

∆∆E IJAB

I > J ∈AB

  ~ ~ − ∑ ∆E IA + ∑ ∆∆E IJA + ∑ ∆E IB + ∑ ∆∆E IJB  I > J ∈A I ∈B I > J ∈B  I ∈A  ~ ~ = ∑ ∆ ∆E I + ∑ ∆ ∆∆E IJ + ∑ ∆ ∆E I + ∑ ∆ ∆∆E IJ + I ∈A

I > J ∈A

I ∈B

I > J ∈B

(40)

∑

∆∆E IJ

I ∈A, J ∈B

~ ~ ~ = ∆ ∆E A + ∆ ∆E B + ∆∆E AB where subscripts AB and A(B) indicates that the value is taken from the calculation of complex AB or isolated A(B), respectively. pot The fluctuation contribution to the potential binding energy, ∆ ∆E AB

in eq 40, is

~ ~ decomposed into the polarization of the fluctuation of the energy A(B), ∆ ∆E A ( ∆ ∆E B ) and ~ the fluctuation of the interaction between A and B, ∆∆E AB . ~ ∆ ∆E A

is called the polarization of the fluctuation of A, because it is the difference in

the fluctuation of A in the complex AB and in isolated A (i.e., the change in the fluctuation of A from its isolated state due to the complex formation, due to the polarization by B), computed as

~ ~ ∆ ∆E A = ∑ ∆ ∆E I + ∆ I∈A

∑

∆∆E IJ

(41)

I > J ∈A

where for I , J ∈ A , the polarization of the fluctuation of the energy of I and I,J interactions are, respectively,

~ ~ ~ ∆ ∆E I = ∆E IAB − ∆E IA

(42)

~ ~ ~ ∆ ∆∆E IJ = ∆∆E IJAB − ∆∆E IJA

(43)


14

Page 15 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


What does this decomposition provide? The potential binding energy in eq 37 is pot pot,0 ) and dynamic ( ∆ ∆E AB ) contributions. The former is the value decomposed into static ( ∆E AB

for the minimum energy structures in MD (MDmin); the latter is the MD-averaged deviation pot (fluctuation) from it. ∆ ∆E AB has a clear physical picture – it is the polarizing effect on the

~ ~ fluctuation in the internal energy of A by B, ∆ ∆E A , and vice versa, ∆ ∆E B , and the ~ fluctuation in the interaction energy between A and B, ∆∆E AB . As shown in the static subsystem analysis,

42

the sum of pair interactions is an

overestimate of the binding energy, and it is necessary to add the contributions of polarizations of A and B, which also describes the desolvation and deformation effects. In the dynamic analysis in this work, applied to MD, one can evaluate how the presence of the ligand (protein) affects the fluctuations in the energies of fragment residues (ligand), and residue-residue interactions. 3. Computational details For all calculations, FMO as implemented 79 in GAMESS 80 was used in parallel, 81 modified to have the above analysis automatically performed in an MD simulation. Hybrid orbital projection operators82 were used to describe fragment boundaries. The solvent (water) was treated using C-PCM,83 FIXPVA tesselation84 and 60 tesserae per atom. Water clusters were fragmented as one molecule per fragment. The protein was divided as 1 residue per fragment at Cα atoms so that residue fragments differ from residues by a CO group; ligand was treated as a separate fragment. PIEDA is applied at the RHF, DFT, DFTB, MP2 and CC levels of theory to three water clusters (H2O)n, n=3 (triangular), 4 (square) and 8 (cubic), whose geometries, taken from the


15


Page 16 of 61

previous work,60 are shown in Figure 2. There are 3, 4 and 12 hydrogen bonds for n equal to 3, 4 and 8, respectively. The aug-cc-pVTZ basis set is used in FMO with the point charge embedding. 85 The following DFT functionals are used: PBE, B3LYP, CAM-B3LYP and M11, chosen to test the effect of long-range corrections (CAM-B3LYP vs B3LYP ) and to test dispersion in M11.86 PBE is selected because DFTB is typically parametrized with it. Throughout this work, the third order DFTB was used (DFTB3) with 3ob parameters.87 For empirical dispersion, the D3(BJ) method75 was employed. To study fluctuations, two complexes of the Trp-cage protein (PDB: 1L2Y) with ligands (neutral or deprotonated p-phenolic acid) were computed with DFTB3. For comparison, RHF, DFT and MP2 calculations were done using the 6-311G** basis set for the minimum energy structure (0Kmin from DFTB). The initial structure of the complexes was obtained with docking simulations88 refined by a consequent optimization.41 In this work, each complex was equilibrated for 10 ps, followed by a production run of 200 ps (NVT MD at 298 K) with RATTLE constraints89 and a step size of 1 fs. Thus, M=200000 energies were averaged in MD (eqs 26 and 28). On a dual 12-core 2.3 GHz Xeon node, 200 ps FMO-DFTB/MD for the anionic ligand took about 6 days (wall-clock time), i.e., one single point gradient took about 2.6 s on the average.

4. Results and Discussion 4.1. Comparison of RHF, DFTB and DFT The results for (H2O)8 are shown in Table 2 (in this Table, the empirical dispersion model is not used). It should be noted that not all 12 hydrogen bonds are similar; they are divided into


16

Page 17 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


three groups of 4 bonds each (Figure 2). For CCSD(T), the interaction per hydrogen bond is for the three groups -7.3 (strong), -6.3 (medium) and -5.5 (weak) kcal/mol. The reason why the trigonal cluster has hydrogen bonds weaker by about 1 kcal/mol than the tetragonal one is attributed to the strain of hydrogen bonding caused by the angles; triangular geometry does not allow for the optimum orientation to maximize the hydrogen bonding energy. The geometry of (H2O)8 resembles a slightly distorted cube. The distinction between bonds is related to their geometry; the strongest group has a shorter HB distance than the other two groups. There are two kinds of facets of the cube: made of 2 weak, 1 medium and 1 strong bonds or 2 strong and 2 medium bonds. In Table 2, the electrostatic component ∆E ES is quite similar for all methods (about -10 kcal/mol), with the exception of DFTB, whose ES component is much smaller. The ES value computed with PBE/STO-3G is -4.98 kcal/mol; the small ES in DFTB is therefore attributed to the use of a minimum basis set, which cannot describe polarization accurately. The exchangerepulsion ∆E EX energy of 8-10 kcal/mol (not computed in DFTB) shows more variation than ES. DFT has a larger EX than HF, by 1-2 kcal/mol. There is a negligible difference (0.1 kcal/mol) between ES components in B3LYP and CAM-B3LYP, perhaps, because this system is neutral. It should be noted that ∆E ES is the electrostatic interaction between polarized fragments and it can be divided into the electrostatic interaction between isolated fragments plus the stabilization part of the polarization energy.60 The charge transfer and mix component ∆E CT+ MIX is from -1.6 to -3.0 kcal/mol for RHF and DFT. The DFT correlation energies ∆E c are between -1.6 and -2.5 kcal/mol. The solvent


17


Page 18 of 61

screening ∆E SOLV is -0.3 for DFTB and otherwise between -0.44 and -0.49 kcal/mol. The averaged HB energies ∆ E are similar: -4.3 (RHF) and from -5.1 to -5.8 kcal/mol for DFT(B).

~ The ∆E IJ0 + ∆E IJCT⋅ES term in DFTB conceptually corresponds to ∆E IJEX + ∆E IJCT + MIX + ∆E IJc in DFT; however, the calculated values are not similar (neither with aug-cc-pVTZ as shown in Table 2, nor with STO-3G, not shown). It was observed that ∆E IJ′ CT⋅ES and ∆E IJV tend to have the opposite sign, and their sum ∆E IJCT⋅ES is typically small (for hydrogen bonds in (H2O)8, about ±0.2

~ kcal/mol). For hydrogen bonds in (H2O)8, ∆E IJ0 is about -0.6 kcal/mol and ∆E IJREP ≤0.0002 ~ kcal/mol. Both ∆E IJ0 and ∆E IJCT⋅ES are short-ranged and have substantial values only for very close fragments. When comparing DFT and DFTB, one should bear in mind that DFTB uses a minimal Slater basis set, compared to Gaussian-type basis sets used in DFT in this work, and DFTB treats only valence electrons as particles. Another point is that the electrostatics in DFTB is not a simple product of two charges divided by distance; instead, a more complex form (eq 9) is used with γ and Γ functions compensating for deficiency of the polarizable point charge model, which in principle, can bring in effects other than pure electrostatics, although in the case of interfragment interaction, excluding covalently connected pairs, the distance is sufficiently large so that the DFTB2 electrostatics may be essentially the Coulomb interaction of point charges. Summarizing, the total value of PIE in DFTB (-5.4 kcal/mol) for water is very similar to DFT; however, individual components deviate much. The ES component in DFTB may be improved with split valence basis sets,90 or, perhaps, with the use of multipoles instead of a charge model.91


18

Page 19 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Next, the DI+RC energy in DFT(B) and other methods is compared (Table 3). In this Table, the empirical dispersion model D3(BJ) is used in RHF and DFT(B). RHF, MP2 and M11 are close to CCSD(T); other DFT functionals show some overestimation. In DFTB, some part of the electron correlation is assigned to the REM component, which perhaps explains why the DI+RC values for DFTB are much underestimated compared to other methods. The total pair interaction energies are shown in Table 4. It is remarkable that they are rather similar to each other, for all methods, and the deviation is 0.7 kcal/mol at most. All methods underestimate the total interaction compared to CCSD(T). DFTB underestimates the total interaction more than other methods. The total interaction in DFTB and PBE results differ by 0.4 kcal/mol at most. It should be noted that pair interaction energies are not equal to hydrogen bond energies; the difference between them is the reference state: for the former (PIE), it is the polarized (interacting) state of monomers, for the latter (binding energies), isolated noninteracting monomers, so the difference is polarization and also deformation.42 Next, the correlation between the interaction energies of DFT(B) and CCSD(T) is studied, by taking all pair interactions in all 3 water cluster, as one set of interactions representing water (3+6+28=37 points). In applications such as drug discovery, it is important to have a good correlation, and the absolute value agreement is not always necessary. The results are shown in Figure 3. All interactions can be clustered into two groups, hydrogen bonding and the rest. The former has larger values for each component than the latter. ES values have a less local clustering pattern because of the long-range, slow decaying (1/r) nature of the Coulomb interaction. For DFTB/D, the total values show a high correlation to CCSD(T), R2 of 0.9945. The ES and DI+RC correlation (R2) is 0.9770 and 0.9523, respectively. The slopes are very different


19


Page 20 of 61

from unity for both ES (0.42) and DI+RC (0.19), but the slope for the total values is close to 1 (0.93). CAM-B3LYP/D results show a very high correlation to CCSD(T) with the R2 of 0.9992 (DI+RC) to 0.9997 (ES). However, the slope of DI+RC is 1.34 and 0.97 for ES, showing that absolute values systematically deviate. These results, extended to include non-HB interactions, provide further support of the above conclusions, that the total interaction in DFTB well agrees with other methods in the absolute sense, and for components a good correlation is observed, but the slope is rather far from unity. For DFT, all components are reasonably close to CCSD(T) in the absolute sense (the electron correlation is overestimated by about 30%), and a high correlation is also observed.

4.2. Protein-ligand complexes The RMSD between the optimized (0Kmin) and MD (MDmin) structures is 0.4492 and 0.4229 Å for the complexes of Trp-cage with the neutral and anionic ligands, respectively. The binding pose of the ligand is not much affected by temperature, and the main difference in the structures is for side chains and residues far from the ligand. The potential energy profile during MD is shown in Figure 4. The MD results for monomer properties in protein-ligand complexes are shown in Figure 5. By its nature (a square of velocity), the kinetic energy for the reference structure (MDmin) E Ikin,0 is positive. Fluctuations in the kinetic energy ∆E Ikin

can be of either sign, but the

averaged kinetic energy E Ikin, 0 + ∆E Ikin is always positive. The neutral and anionic ligands have the averaged kinetic energy of 15.3 and 14.2 kcal/mol, respectively, which is relatively large (the kinetic energy is additive for all atoms and the ligand has more atoms than most residues though).

~ Fluctuations in the monomer potential energy ∆E I

are between -6.4 and 10.5 kcal/mol,


20

Page 21 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


showing changes in the polarization and solvation energy of fragments. Some charged residues have large polarization fluctuations of 10 kcal/mol: e.g., Asn-1 (neutral ligand), Arg-16 (anionic ligand). The effective temperatures TI of fragments are shown in Figure 6. The global temperatures T computed in MD are 297.997, 298.005 and 298.019 K for the protein by itself, neutral and anionic ligand complexes, respectively, which corresponds to the temperature set to 298 K in NVT MD (deviations show the efficiency of the thermostat). One can see that there is a significant variation in effective temperatures in Figure 6 that vary in the range of 280-323 K. For many fragments, the effective temperatures for the complexes of the two ligands resemble each other, but some substantial deviations are found: for example, for Asn-1 (strongly bound to the anionic ligand, not bound to the neutral). Some differences cannot be explained by the binding: Asp-9 is the “hottest” residue fragment in both complexes, it is strongly bound to the neutral ligand, but has a strong repulsion to the anionic ligand. It is also interesting to see how complexation affects effective temperatures. For some residue fragments, there is no substantial effect (Trp-6, Leu-7 and Lys-8). The ligands are “hotter” in the complex than the global temperature of 298 K: 320 K (neutral) and 316 K (anionic); the much larger protein by the virtue of its size (note the size-related weight in eq 34) is “cooled” to 296.8 and 297.1 K in the complex with the neutral and anionic ligand, respectively. Some residues are “cooled” down (Pro-18, Ile-4), whereas some are “heated” up (Pro-12, Asn-1). Thus, it was observed that in the two protein-ligand complexes the protein is “cooled” down and the ligand is “heated” up by the complexation, in the sense of their effective temperature, determined by the normalized kinetic energy.


21


Page 22 of 61

Pair properties for protein-ligand complexes are shown in Figure 7. The neutral ligand is strongly bound to Asp-9; the anionic ligand has an interesting bonding (Asn-1, Lys-8, Arg-16) – antibonding (Asp-9, Ser-20) pattern, which correlates well with the charges (opposite for bonding and same for antibonding). The MD and 0K minimum interactions are essentially the same. The fluctuations in the pair interactions can have a weakening effect due to thermal motions (Asp-9 bound to the neutral ligand). This can be described as kinetic pressure leading to a weaker bonding, compared to the minimum energy structure at 0K; in other words, minimum energy structure at 0K has a tendency to overbind compared to room temperature. The neutral ligand is bound on the surface of the protein, whereas the anionic ligand is bound in the pocket deep inside the protein. This explains the difference in the binding pattern (only one important residue for the former and many for the latter), and the larger effect of thermal fluctuations for the former. The bounds for interactions (minimum and maximum values) show that occasionally, interactions can be very different, especially for the neutral ligand: Tyr-3, Arg-16 and Pro-18 have large minimum values of PIE, which means that occasionally they can strongly bind the ligand, but on the average their interaction is small. For the anionic ligand, buried inside the protein, only Trp-6 has a similar property of occasionally binding during MD, but no binding on the average (it can be noted that Trp-cage is the key residue binding Trp-cage in its native conformation,26 so its atoms are too busy interacting with other residues to also bind a ligand); also, occasionally, the strong repulsion to Asp-9 is reduced. The trends for the maximum values is similar – occasionally in MD, some residues (Gln-5, Arg-16) can be repulsive for the neutral ligand. Asp-9 binds the neutral ligand on the average but occasionally in MD it can be repulsive and thermal motions weaken its binding to the ligand.


22

Page 23 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The total binding energies and their decomposition are shown in Table 5. The neutral and anionic ligands have endothermic (2.01) and exothermic (-9.90) energies (kcal/mol), respectively (the binding pose of the neutral ligand found in MD differs from that found in earlier docking simulations and refined with an FMO optimization88). Interestingly, the kinetic energy (the sum of the reference value and fluctuation) has a negligible contribution. The fluctuations in the potential energy contribute relatively little about (1 kcal/mol) for both ligands. More details for the binding can be found in Table 6 and Table 7 for the neutral and anionic ligand, respectively. For the neutral ligand, the total potential binding energy in Table 5 is 3.49, which is further decomposed in Table 6 into -11.87 and 15.36 contributions (kcal/mol) of monomers and dimers, respectively, which describe polarization and desolvation of the protein and ligand, and the interaction between them. Similarly, the fluctuation value of -1.37 in Table 5 is decomposed into -14.90 and 13.53 kcal/mol in Table 6. In many cases, fluctuations in the values have the opposite sign, for example, Asn-1 has the reference and fluctuation values of -16.06 and 11.46 kcal/mol (thermal weakening). This is in agreement with the classical statistical theory applied to scaling pair interactions in FMO.73 However, as was also reported for solvent antiscreening,61 there are cases when fluctuations have the same sign increasing the interaction, for example, Arg-16 (neutral ligand) or Asn-1 (anionic ligand). The anionic ligand, due to its charge, has much larger contributions in Table 7, compared to the neutral ligand in Table 6. In particular, its own monomer contribution of 89.35 is mainly determined by the desolvation penalty, because the ligand is buried in the protein and due to complexation it loses much of its very exothermic solute-solvent interactions. The protein polarization is very large, as seen in the large value of residue-residue contributions (125.32 kcal/mol), which shows the weakening of residue-residue interactions caused by the presence of


23


Page 24 of 61

the ligand. The total strong destabilization of the protein and ligand, (polarization and desolvation) is compensated by the very strong protein-ligand interactions, so that the total binding (-8.98 kcal/mol for the reference structure) is exothermic. A comparison of protein-ligand interactions at different levels of theory is shown in Table 8 and Figure 8. The exact agreement between RHF and MP2 for most residues is because at large separation between fragments, the ES-DIM approximation is used in both, resulting in the same values. CAM-B3LYP also uses this approximation, but it is applied with DFT densities, so that DFT and RHF values slightly differ. DFTB uses atomic charges for this approximation. Only Lys-8, Asp-9 and Gly-10 are computed without this approximation for the neutral ligand. There is a good agreement between all theories for the neutral ligand. The largest deviation is for Asp-9: DFTB predicts about -42 kcal/mol and other theories about -57±1 kcal/mol. This fragment residue is very close to the ligand and strongly binds to it. The bulk of the difference is in the electrostatic component, as in water clusters; the CAM-B3LYP and DFTB

values

are

-56.10

and

-33.64

kcal/mol,

respectively.

Interestingly,

the

∆E IJEX + ∆E IJCT + MIX + ∆E IJCORR value for CAM-B3LYP is 8.58 kcal/mol (driven by a large repulsion

~ ∆E IJEX ) but the DFTB value of ∆E IJ0 + ∆E IJCT ⋅ES + ∆E IJDI is -2.56 kcal/mol.

The anionic ligand is buried deeply into the protein, only residue fragments 8,10-15 are treated with the ES-DIM approximation. Overall, DFTB values are close to the other theories. There are two exceptions: for Leu-7, DFTB predicts no binding (-0.5 kcal/mol) but other theories a substantial binding (-7.5±0.5 kcal/mol). The binding of Pro-17 in kcal/mol is: -2.4 (MP2), -0.6 (RHF), -7.4 (DFT) and 3.1 (DFTB), that is, DFTB predicts repulsion and other theories attraction. The reason why DFTB differs for Leu-7 is electrostatics, it is -0.9 (DFTB) and -7.5 (other theories) kcal/mol. For Pro-17, somewhat surprisingly, DFT predicts a large charge transfer, both


24

Page 25 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in terms of charge (0.06 a.u.) and energy (-8.1 kcal/mol), compared to 0.01 a.u. and -2.3 kcal/mol,

~ respectively, in RHF (and MP2). For DFTB, the ∆EIJ0 + ∆EIJCT⋅ES term is small (0.0 kcal/mol) and the electrostatics is repulsive, 1.2 kcal/mol (in RHF and DFT the electrostatics is -0.2 and -1.5 kcal/mol, respectively). All theories predict a somewhat similar solvent screening of 2.4...3.0 kcal/mol; whereas the DI+RC term in DFTB is -0.9 and in other theories -0.5...-3.0 kcal/mol. For Pro-18, DFTB predicts a binding of about -9 and other theories -19±1 kcal/mol. Ligand-residue interactions in DFTB show more deviations in the absolute values from other theories, compared to water clusters. It may in part be because the latter are important test systems in fitting DFTB parameters, and strong interactions including polarization and charge transfer in ligands, and especially charged residues, may be more difficult to describe with a general set of DFTB parameters. Nevertheless, the correlation coefficients between full sets of residue-ligand interactions for all pairs of methods including DFTB are 0.99 or larger, so that overall all methods agree well with each other (Figure 9). However, the slope is substantially different from 1: comparing DFTB/D and MP2, the slope is 0.73 and 0.91 for the neutral and anionic ligand, respectively. This shows the trend of DFTB to systematically underestimate interactions (mainly, electrostatic). Some part of the difference may be due to the basis set superposition error.

5. Conclusions In this work, pair interactions between fragments have been decomposed into 5 components: electrostatics, exchange-repulsion, charge transfer plus mix terms, dispersion plus remainder correlation and solvent screening (DFT) or electrostatics, coupling of charge transfer to electrostatics and model Hamiltonian, dispersion and solvent screening (DFTB). In addition,


25


Page 26 of 61

FMO analysis has been integrated into molecular dynamics simulations to study fluctuations of fragment energies and pair interactions. The developed analysis has been applied to water clusters. It is found that DFTB has a much smaller electrostatic component than other methods, attributed to the use of a minimum basis set. In water clusters, the absolute magnitude of charge-transfer and exchange-repulsion components in DFT is larger than in RHF by 1-2 kcal/mol. This is attributed to a smaller gap between occupied and virtual orbitals in DFT compared to RHF,67 which promotes charge (density) transfer. A density transfer due to the Pauli principle increases the exchange interaction. The DI+RC energy per hydrogen bond, underestimated in DFTB (-0.5 kcal/mol), varies in RHF/D, MP2, DFT/D and CCSD(T) between about -2 and -3 kcal/mol. Among other methods, RHF/D, M11/D and MP2 are the closest to CCSD(T), whereas PBE/D, B3LYP/D and CAMB3LYP/D deviate more. The solvent screening is very similar in all methods except for DFTB, where it is underestimated by about 40%). The total interactions, however, are rather similar in all methods, with a deviation not exceeding 0.7 kcal/mol (about 12%). Comparing DFTB/D or CAM-B3LYP/D with CCSD(T), the correlation between PIEDA components for a set of all pair interactions in water clusters is very high: the correlation coefficient R2 is between 0.95 and 1.00; the slope, however, can be quite different from 1, which shows that the absolute values are different. For instance, for the electrostatic interaction in DFTB vs CCSD(T), the slope is 0.42 (a systematic underestimatation problem of DFTB) but the correlation coefficient R2 is 0.98. It may be expected that because DFTB is a series expansion64 with respect to charge fluctuations, the accuracy of DFTB may be worse for short contacts and large charge transfer.


26

Page 27 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The binding energy of ligands to the Trp-cage protein have been analyzed in detail. The magnitude of the protein and ligand polarization and desolvation has been evaluated; for the charged ligand these effects are expectedly large. A weakening of contributions due to thermal fluctuations has been described and attributed to kinetic pressure, and rare occasions of a thermal strengthening have been observed. The weakening corresponds to the statistical theory suggested before73 although no single parameter describing the weakening appears to exist. Here, kinetic pressure refers to the kinetic energy of nuclei, which tends to destabilize the system and populate conformational space with higher potential energy. It should not be confused with the concept of kinetic pressure92 of electrons that was introduced to describe the concept of chemical bonding and electron sharing. The total binding energy between two subsystems A and B in MD is decomposed into the binding energy for a reference structure (e.g., the minimum of potential energy) and fluctuations of three kinds: polarization of the fluctuation energy of A and B, and fluctuation of the interaction energy between A and B. The polarization of the fluctuation energy of subsystem A or B is the change in the fluctuation energy due to the complex formation (the difference in the fluctuation energy in the complex and isolated state). Thermal fluctuations due to kinetic pressure substantially offset the binding energy computed for the potential energy minimum. The averaging of properties over the whole trajectory should be more reliable than considering only a small subset of snapshots. The distribution of effective temperatures of residue fragments identified hot and cold residues, and the differences between the two ligand complexes of Trp-cage. As a result of complexation, the protein is cooled down and the ligand is heated up (a transfer of kinetic energy


27


Page 28 of 61

from the protein to a ligand). The relation between effective temperatures of residues are biological activity is an interesting topic that may be investigated with the analysis in this work. It has been shown that MP2, RHF/D, DFT/D and DFTB/D predict similar, highly correlated interactions between residue fragments and ligands, although for DFTB the slope is different from one and the method appears to systematically underestimates interactions (mainly, electrostatics). Because DFTB is a much faster approach (by 3-4 orders of magnitude67,70) than the other three, it can be used for a rapid evaluation of interaction suitable for high throughput screening is drug discovery. The results of this work may be used to improve parametrization in DFTB. PIEDA/DFT, developed in this work, may be useful to study interactions in various systems, containing metals, where RHF or MP2 may not work well. The comparative study of components may be useful in designing force fields based on fragment-based methods. 93,94 The dynamical aspects of proteinligand binding may be useful to analyze drug retention and residence time.

ACKNOWLEDGMENT This work was supported by JSPS KAKENHI, Grant Number 16K05677 (D.G.F.). We thank Prof. Andreas Savin for critical remarks on this work.


28

Page 29 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure Captions Figure 1. Solvent screening of the interaction between fragments I and J, ∆E IJSOLV , is computed SOLV es disp rep CT ⋅es as ∆E ISOLV , where ∆E ISOLV is the sum of the ( J ) + ∆E J ( I ) ( J ) = ∆E I ( J ) + ∆ E I ( J ) + ∆E I ( J ) + ∆E I ( J )

electrostatic (es), dispersion (disp), repulsion (rep) and charge-transfer-electrostatic coupling (CT·es) interactions between solute fragment I and solvent cavity around fragment J). Figure 2. Structures of water clusters. Pair interaction energies are shown in kcal/mol at the CCSD(T)/aug-cc-pVTZ level. Figure 3. Correlation of ∆E IJES , ∆E IJDI + RC , and ∆E IJ values (kcal/mol) in water clusters: (a) DFTB3(3ob)/D3(BJ)

vs

CCSD(T)/aug-cc-pVDZ

and

(b)

CAM-B3LYP/D3(BJ)

vs

CCSD(T)/aug-cc-pVDZ. Figure 4. Potential energy of the Trp-cage complexes with the neutral (blue) and anionic (red) ligands during 200 ps NVT FMO-DFTB/MD. Figure 5. Monomer properties from 200 ps NVT FMO-DFTB/MD: MDmin fragment kinetic

~ energies E Ikin,0 , fluctuations in kinetic ∆E Ikin and potential ∆E I

energies: a) neutral and b)

anionic ligands. Figure 6. Effective temperatures TI of fragments in Trp-cage (isolated) and its complexes with the neutral and anionic ligands (200 ps NVT FMO-DFTB/MD at 298 K). The global temperature T is also drawn.


29


Page 30 of 61

Figure 7. Pair properties from 200 ps NVT FMO-DFTB/MD: MDmin interaction energies ∆E IJ0 , fluctuations ∆∆E IJ , minimum ∆E IJmin and maximum ∆E IJmax values, and interaction energies at the energy minimum 0Kmin ∆E IJ0K for the a) neutral and b) anionic ligands. Figure 8. Interactions of residue fragments in Trp-cage with a) neutral and b) anionic ligands (0Kmin). Figure 9. Correlation between residue-ligand interaction energies in Trp-cage complexes, computed with DFTB/D and other theories (MP2, RHF/D and CAM-B3LYP/D) for the a) neutral and b) anionic ligands. A linear regression line and corresponding coefficients are shown for MP2 vs DFTB/D.


30

Page 31 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47


TABLES Table 1. Correspondence between pair interaction energy (PIE) contributions. term

RHF

MP2, CC

DFT

DFTB

electrostatic ∆E IJES

energy interactions between same electron density and nuclei RHF

as

in similar to RHF, but using DFT interactions between density atomic charges

solvation ∆E IJSOLV

energy interactions of nuclei and same electron density with RHF solvent surface charges

as

in similar to RHF, but using DFT interaction of solute charges with solvent density surface charges

correlation energy (dispersion + remainder correlation) ∆E IJDI+ RC exchange-repulsion energy ∆E IJEX

dispersiona

electron correlation

Pauli repulsion between same electron densities RHF

as

charge transfer energy remainder (total PIE minus same + high order mix terms the other four components) RHF ∆E IJCT+ MIX

as

a

dispersion and DFT correlation dispersiona functional energy a

in similar to RHF, but using DFT coupling of charge density, plus a contribution transfer to the model ~ from DFT exchange Hamiltonian ∆E IJ0 and electrostatics in similar to RHF, but using DFT to CT⋅ES b ∆E IJ . PIE

~ The values of empirical dispersion are calculated using parameters that are method dependent. b ∆E IJ0 in DFTB may also be related to

∆E IJc in DFT.


31

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

Page 32 of 61

Table 2. Averaged interaction energies ∆ E for hydrogen bonds (kcal/mol) in (H2O)8, and their components: electrostatic ∆E ES , exchange-repulsion ∆E EX , charge transfer + mix terms ∆E CT+ MIX , correlation ∆E c and solvent screening ∆E SOLV energies.a method RHF

a

∆E ES

-10.18

∆E EX

7.89

∆E CT+ MIX

∆E c

∆E SOLV

∆E

-1.58

-0.44

-4.30

-0.57b

-0.30

-5.35

DFTB

-4.47

PBE

-9.63

9.77

-3.04

-2.15

-0.49

-5.53

B3LYP

-9.82

10.36

-2.67

-2.46

-0.48

-5.08

CAM-B3LYP

-9.92

9.59

-2.53

-2.46

-0.48

-5.81

M11

-9.80

8.84

-2.38

-1.59

-0.48

-5.41

DFTB3 with 3ob parameters; other methods with aug-cc-pVTZ. The empty entrees correspond to quantities not computed in a

~ particular method. b For DFTB, the averaged contribution of ∆E IJ0 + ∆E IJCT⋅ES is shown, which corresponds to ∆E IJEX + ∆E IJCT + MIX + ∆E IJc in other methods.


32

Page 33 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3. Averaged correlation energy ∆E DI+ RC (kcal/mol) per hydrogen bond.a method

(H2O)3

(H2O)4

(H2O)8

RHF

-1.91

-1.97

-1.87

DFTB3

-0.50

-0.50

-0.50

PBE

-2.45

-2.82

-2.58

B3LYP

-3.05

-3.35

-3.13

CAM-B3LYP

-2.69

-2.99

-2.78

M11

-2.08

-2.30

-2.12

MP2

-1.86

-2.14

-1.97

CCSD(T)

-1.98

-2.24

-2.08

a

DFTB3 with 3ob parameters; other methods with aug-cc-pVTZ. RHF and DFT(B) methods use D3(BJ).

Table 4. Averaged total interaction ∆E (kcal/mol) per hydrogen bond.a method

(H2O)3

(H2O)4

(H2O)8

RHF

-6.01

-6.72

-6.17

DFTB3

-5.34

-6.36

-5.84

PBE

-5.77

-6.52

-5.97

B3LYP

-5.55

-6.30

-5.76

CAM-B3LYP

-5.92

-6.71

-6.13

M11

-5.87

-6.45

-5.93

MP2

-5.96

-6.89

-6.27

CCSD(T)

-6.08

-6.99

-6.38

a

DFTB3 with 3ob parameters; other methods with aug-cc-pVTZ. RHF and DFT(B) methods use D3(BJ).


33


Page 34 of 61

Table 5. Decomposition of the protein-ligand binding energies (BE, in kcal/mol) for Trpcage at 298 K (at the FMO-DFTB3/3ob level). contribution

meaning

pot,0 ∆EAB

reference potential BE

pot ∆ ∆E AB

neutral ligand anionic ligand 3.49

-8.98

fluctuation of potential BE

-1.37

-0.82

kin,0 ∆EAB

reference kinetic BE

-1.19

-30.51

kin ∆ ∆E AB

fluctuation of kinetic BE

1.09

30.41

pot pot,0 pot ∆ E AB = ∆E AB + ∆ ∆E AB

potential binding BE

2.11

-9.80

kin kin,0 kin ∆ E AB = ∆E AB + ∆ ∆E AB

kinetic binding BE

-0.10

-0.10

pot pot ∆U = ∆ EAB + ∆ EAB

total BE

2.01

-9.90


34

Page 35 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 6. Contributions (kcal/mol) to the potential binding energy of the neutral ligand to ~ Trp-cage (FMO-DFT3/3ob, at 298 K): fragment reference energies ∆E I0 and their ~ fluctuations ∆ ∆E I , reference interaction energies ∆EIL0 and their fluctuations ∆∆ E IL . ~0 ~ ∆EIL0 ∆ ∆E I fragment ∆E I -16.06 11.46 2.58 Asn-1

a

∆∆ E IL

-1.35

Leu-2

-6.43

3.35

-0.29

-0.71

Tyr-3

0.06

6.70

-0.59

-0.87

Ile-4

6.93

-7.78

-0.46

0.01

Gln-5

-10.88

7.51

-1.51

-1.45

Trp-6

4.31

-4.37

-1.64

-0.75

Leu-7

16.87

-18.24

-1.33

0.46

Lys-8

3.91

-9.81

12.19

-6.83

Asp-9

-0.33

-0.29

-41.32

15.63

Gly-10

10.45

-8.00

0.18

-0.41

Gly-11

6.17

-7.09

0.76

-0.33

Pro-12

4.88

-3.67

-0.73

0.28

Ser-13

3.11

-3.93

0.15

-0.06

Ser-14

-6.29

4.49

0.50

-0.38

Gly-15

-3.52

2.22

-0.05

-0.04

Arg-16

-1.61

-4.73

7.22

-1.89

Pro-17

-7.31

6.29

-0.45

-0.45

Pro-18

1.15

-3.38

-0.20

-3.50

Pro-19

-7.01

7.40

0.03

-0.36

Ser-20

-27.65

12.88

-2.09

0.69

other

17.39a

-5.91b

42.42c

15.84d

total

-11.87

-14.90

15.36

13.53

~ ligand value ∆E L0 ,

∑ ∆∆E

I > J ∈A

0 d IJ

b

~ ligand fluctuation ∆ ∆EL ,

sum of residue-residue interaction fluctuations

c

sum of residue-residue values

∑ ∆ ∆E

IJ

I > J ∈A


35


Page 36 of 61

Table 7. Contributions (kcal/mol) to the potential binding energy of the anionic ligand (L) ~ to Trp-cage (FMO-DFT3/3ob, at 298 K): fragment reference energies ∆E I0 and their ~ fluctuations ∆ ∆E I , reference interaction energies ∆EIL0 and their fluctuations ∆∆ E IL .

~ ~ fragment ∆E I0 ∆ ∆E I ∆EIL0 11.01 0.97 -69.99 Asn-1 6.21 -3.31 -34.17 Leu-2

a

∆∆ E IL 5.88

4.41

Tyr-3

2.90

-4.08

12.21

0.41

Ile-4

2.83

-4.45

5.78

0.05

Gln-5

-1.75

0.25

-31.77

4.57

Trp-6

4.50

-3.66

-6.07

-3.16

Leu-7

7.10

-7.97

-0.21

0.09

Lys-8

3.19

-9.10

-42.41

1.67

Asp-9

-11.06

3.46

35.75

-1.06

Gly-10

-1.35

2.23

-0.77

-0.78

Gly-11

5.97

-5.85

-2.82

-0.32

Pro-12

-2.68

3.01

4.30

0.93

Ser-13

6.68

-5.64

-0.94

-0.31

Ser-14

1.56

-0.23

-2.92

-0.10

Gly-15

-4.36

3.49

-2.01

1.98

Arg-16

-4.21

2.86

-71.29

2.74

Pro-17

-2.43

2.19

2.55

-1.53

Pro-18

2.77

-3.06

-10.18

2.58

Pro-19

-2.71

3.80

-9.39

0.81

Ser-20

-34.61

14.42

11.12

-1.14

other

89.35a

-11.17b 125.32c

3.31d

total

78.93

~ ligand value ∆E L0 ,

∑ ∆∆E

I > J ∈A

0 d IJ

-21.86 b

-87.91

21.04

~ ligand fluctuation ∆ ∆EL ,

sum of residue-residue interaction fluctuations

c

sum of residue-residue values

∑ ∆ ∆E

IJ

I > J ∈A


36

Page 37 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 8. Ligand (L) – residue (I) fragment interaction energies ∆E IL (kcal/mol).a Neutral ligand CAMRHF/D B3LYP/D

residue I

MP2

Asn-1

2.95

2.95

Leu-2

-0.54

Tyr-3

Anionic ligand CAMRHF/D B3LYP/D -77.49 -76.05

DFTB/D

MP2

DFTB/D

2.74

2.59

-79.24

-0.54

-0.58

-0.36

-37.33

-35.05

-37.01

-33.17

-0.67

-0.67

-0.61

-0.63

7.74

7.86

4.80

12.41

Ile-4

-0.58

-0.58

-0.63

-0.48

4.85

4.84

4.63

6.30

Gln-5

-1.88

-1.88

-1.95

-1.61

-38.70

-37.20

-33.63

-33.01

Trp-6

-2.69

-2.69

-2.27

-2.12

-6.70

-6.29

-9.63

-8.39

Leu-7

-1.23

-1.23

-1.14

-1.25

-7.08

-6.92

-7.90

-0.46

Lys-8

11.82

11.72

9.03

12.54

-47.05

-47.05

-46.12

-45.42

Asp-9

-58.42

-56.94

-55.58

-41.88

45.95

45.98

44.23

42.38

Gly-10

1.87

1.97

1.57

0.05

-4.07

-4.07

-4.37

-1.50

Gly-11

1.11

1.11

1.00

0.72

-2.84

-2.84

-2.70

-3.11

Pro-12

-0.38

-0.38

-0.31

-0.65

4.67

4.67

4.44

5.14

Ser-13

0.18

0.18

0.14

0.12

-1.82

-1.82

-1.68

-0.99

Ser-14

0.24

0.24

0.23

0.59

-1.60

-1.60

-1.51

-3.39

Gly-15

0.35

0.35

0.31

0.01

-3.15

-3.15

-2.88

-1.97

Arg-16

8.05

8.05

7.44

7.20

-95.66

-94.37

-94.71

-83.59

Pro-17

-0.43

-0.43

-0.38

-0.51

-2.36

-0.55

-7.37

3.06

Pro-18

-0.36

-0.36

-0.26

-0.19

-20.13

-17.76

-19.21

-8.87

Pro-19

0.25

0.25

0.31

0.08

-11.59

-9.68

-11.74

-12.13

Ser-20

-2.19

-2.19

-2.05

-2.10

7.46

8.03

6.30

11.92

-70.17

a

Using D3(BJ) as the dispersion (D) model in RHF and DFT, 6-311G** basis set for MP2, RHF and CAM-B3LYP , and 3ob for DFTB3. The same 0Kmin (DFTB) structure is used for all methods. Bold numbers are for values computed as SCF; other values are obtained with the separated dimer approximation in eqs 8-9.


37


Page 38 of 61

REFERENCES ( 1 ) Otto, P.; Ladik, J. Investigation of the Interaction Between Molecules at Medium Distances: I. SCF LCAO MO Supermolecule, Perturbational and Mutually Consistent Calculations for two Interacting HF and CH2O Molecules. Chem. Phys. 1975, 8, 192–200. (2) Gao, J. Toward a Molecular Orbital Derived Empirical Potential for Liquid Simulations. J. Phys. Chem. B 1997, 101, 657–663. (3) Fang, T.; Li, Y.; Li, S. Generalized Energy-Based Fragmentation Approach for Modeling Condensed Phase Systems. WIREs: Comp. Mol. Sc. 2017, 7, e1297. ( 4 ) Söderhjelm, P.; Ryde, U. How Accurate Can a Force Field Become? A Polarizable Multipole Model Combined with Fragment-Wise Quantum-Mechanical Calculations. J. Phys. Chem. A 2009, 113, 617–627. (5) Liu, J.; Zhang, J. Z. H.; He, X. Fragment Quantum Chemical Approach to Geometry Optimization and Vibrational Spectrum Calculation of Proteins. Phys. Chem. Chem. Phys. 2016, 18, 1864–1875. (6) Kobayashi, R.; Amos, R.; Collins, M. A. Microsolvation within the Systematic Molecular Fragmentation by Annihilation Approach. J. Phys. Chem. A 2017, 121, 334–441. (7) Yu, H.; Leverentz, H. R.; Bai, P.; Siepmann, J. I.; Truhlar, D. G. Water 26-Mers Drawn From Bulk Simulations: Benchmark Binding Energies for Unprecedentedly Large Water Clusters and Assessment of the Electrostatically Embedded Three-Body and Pairwise Additive Approximations. J. Phys. Chem. Lett. 2014, 5, 660–670.


38

Page 39 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(8) He, X.; Merz, K. M. Jr. Divide and Conquer Hartree−Fock Calculations on Proteins. J. Chem. Theory Comput. 2010, 6, 405–411. ( 9 ) Steinmann, C.; Fedorov, D. G.; Jensen, J. H. Effective Fragment Molecular Orbital Method: a Merger of The Effective Fragment Potential and Fragment Molecular Orbital Methods. J. Phys. Chem. A 2010, 114, 8705–8712. ( 10 ) Sahu, N.; Gadre, S. R. Vibrational Infrared and Raman Spectra of Polypeptides: Fragments-in-Fragments within Molecular Tailoring Approach. J. Chem. Phys. 2016, 144, 114113. (11) Gao, J.; Truhlar, D. G.; Wang, Y.; Mazack, M. J. M.; Löffler, P.; Provorse, M. R.; Rehak, P. Explicit Polarization: A Quantum Mechanical Framework for Developing Next Generation Force Fields. Acc. Chem. Res. 2014, 47, 2837–2845. (12) Jose, K. V. J.; Raghavachari, K. Evaluation of Energy Gradients and Infrared Vibrational Spectra through Molecules-in-Molecules Fragment-Based Approach. J. Chem. Theory Comput. 2015, 11, 950–961. ( 13 ) Liu, J.; Herbert, J. M. Pair–Pair Approximation to the Generalized Many-Body Expansion: An Alternative to the Four-Body Expansion for ab Initio Prediction of Protein Energetics via Molecular Fragmentation. J. Chem. Theory Comput. 2016, 12, 572–584. (14) Gurunathan P. K.; Acharya, A.; Ghosh, D.; Kosenkov, D.; Kaliman, I.; Shao, Y.; Krylov, A. I.; Slipchenko, L. V. Extension of the Fragment Potential Method to Macromolecules. J. Phys. Chem. B 2016, 120, 6562–6574.


39


Page 40 of 61

( 15 ) Nishizawa, H.; Nishimura, Y.; Kobayashi, M.; Irle, S.; Nakai, H. Three Pillars for Achieving Quantum Mechanical Molecular Dynamics Simulations of Huge Systems: Divideand-Conquer, Density-Functional Tight-Binding, and Massively Parallel Computation. J. Comput. Chem. 2016, 37, 1983–1992. (16) Friedrich, J.; Hanrath, M.; Dolg, M. Energy Screening for the Incremental Scheme: Application to Intermolecular Interactions. J. Phys. Chem. A 2007, 111, 9830–9837. (17) Jacob, C. R.; Neugebauer, J. Subsystem Density-Functional Theory. WIREs Comput. Mol. Sci. 2014, 4, 325–362. (18) Aoki, Y.; Gu, F. L. An Elongation Method for Large Systems toward Bio-Systems. Phys. Chem. Chem. Phys. 2012, 14,7640–7668. (19) Shimazaki, T.; Kitaura, K.; Fedorov, D. G.; Nakajima, T. Group Molecular Orbital Approach to Solve the Huzinaga Subsystem Self-Consistent-Field Equations. J. Chem. Phys. 2017, 146, 084109. (20) Sun, Q.; Chan, G. K. L. Quantum Embedding Theories. Acc. Chem. Res. 2016, 49, 2705– 2712. (21) Elsohly, A. M.; Shaw, C. L.; Guice, M. E.; Smith, B. D.; Tschumper, G. S. Analytic Gradients for the Multicentred Integrated QM:QM Method for Weakly Bound Clusters: Efficient and Accurate 2-Body : Many-Body Geometry Optimizations. Mol. Phys. 2007, 105, 2777–2782.


40

Page 41 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(22) Parrish, R. M.; Sherrill, C. D. Spatial Assignment of Symmetry Adapted Perturbation Theory Interaction Energy Components: the Atomic SAPT Partition. J. Chem. Phys. 2014, 141, 044115. (23) Gordon, M. S.; Pruitt, S. R.; Fedorov, D. G.; Slipchenko, L. V. Fragmentation Methods: a Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632–672. (24) Akimov, A. V.; Prezhdo, O. V. Large-Scale Computations in Chemistry: A Bird’s Eye View of a Vibrant Field. Chem. Rev. 2015, 115, 5797−5890. (25) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment Molecular Orbital Method: an Approximate Computational Method for Large Molecules. Chem. Phys. Lett. 1999, 313, 701–706. (26) Fedorov, D.G.; Kitaura, K. Extending the Power of Quantum Chemistry to Large Systems with the Fragment Molecular Orbital Method. J. Phys. Chem. A 2007, 111, 6904–6914. (27) Fedorov, D.G.; Nagata, T.; Kitaura, K. Exploring Chemistry with the Fragment Molecular Orbital Method. Phys. Chem. Chem. Phys. 2012, 14, 7562–7577. (28) Tanaka, S.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Fukuzawa, K. Electron-Correlated Fragment-Molecular-Orbital Calculations for Biomolecular and Nano Systems. Phys. Chem. Chem. Phys. 2014, 16, 10310–10344. (29) Fedorov, D. G. The Fragment Molecular Orbital Method: Theoretical Development, Implementation in GAMESS and Applications. WIREs: Comp. Mol. Sc., WIREs, 2017, 7, e1322.


41


Page 42 of 61

(30) Nakano, T.; Mochizuki, Y.; Yamashita, K.; Watanabe, C.; Fukuzawa, K.; Segawa, K.; Okiyama, Y.; Tsukamoto, T.; Tanaka, S. Development of the Four-Body Corrected Fragment Molecular Orbital (FMO4) Method. Chem. Phys. Lett. 2012, 523, 128–133. (31) Nakata, H.; Fedorov, D. G.; Yokojima, S.; Kitaura, K.; Nakamura, S. Derivatives of the Approximated Electrostatic Potentials in Unrestricted Hartree-Fock Based on the Fragment Molecular Orbital Method and an Application to Polymer Radicals. Theor. Chem. Acc. 2014, 133, 1477. (32) Sawada, T.; Fedorov, D. G.; Kitaura, K. Structural and Interaction Analysis of Helical Heparin Oligosaccharides with the Fragment Molecular Orbital Method. Int. J. Quant. Chem. 2009, 109, 2033–2045. ( 33 ) Watanabe, T.; Inadomi, Y.; Fukuzawa, K.; Nakano, T.; Tanaka, S.; Nilsson, L.; Nagashima, U. DNA and Estrogen Receptor Interaction Revealed by Fragment Molecular Orbital Calculations, J. Phys. Chem. B. 2007, 111, 9621–9627. (34) Heifetz, A.; Chudyk, E. I.; Gleave, L.; Aldeghi, M.; Cherezov, V.; Fedorov, D. G.; Biggin, P. C.; Bodkin, M. J. The Fragment Molecular Orbital Method Reveals New Insight into the Chemical Nature of GPCR-Ligand Interactions. J. Chem. Inf. Model. 2016, 56, 159–172. (35) Fukunaga, H.; Fedorov, D. G.; Chiba, M.; Nii, K.; Kitaura, K. Theoretical Analysis of the Intermolecular Interaction Effects on the Excitation Energy of Organic Pigments: Solid State Quinacridone. J. Phys. Chem. A 2008, 112, 10887–10894.


42

Page 43 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(36) Avramov, P. V.; Fedorov, D. G.; Sorokin, P. B.; Sakai, S.; Entani, S.; Ohtomo, M.; Matsumoto, Y.; Naramoto, H. Intrinsic Edge Asymmetry in Narrow Zigzag Hexagonal Heteroatomic Nanoribbons Causes their Subtle Uniform Curvature. J. Phys. Chem. Lett. 2012, 3, 2003–2008. (37) Komeiji, Y.; Nakano, T.; Fukuzawa, K.; Ueno, Y.; Inadomi, Y.; Nemoto, T.; Uebayasi, M.; Fedorov, D. G.; Kitaura, K. Fragment Molecular Orbital Method: Application to Molecular Dynamics Simulation, 'ab initio FMO-MD'. Chem. Phys. Lett. 2003, 372, 342–347. (38) Nishimoto, Y.; Nakata, H.; Fedorov, D.G.; Irle, S. Large-Scale Quantum-Mechanical Molecular Dynamics Simulations Using Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method. J. Phys. Chem. Lett. 2015, 6, 5034–5039. (39) Nishimoto, Y.; Fedorov, D.G. Three-Body Expansion of the Fragment Molecular Orbital Method Combined with Density-Functional Tight-Binding. J. Comput. Chem. 2017, 38, 406– 418. ( 40 ) Nebgen, B.; Prezhdo, O. V. Fragment Molecular Orbital Nonadiabatic Molecular Dynamics for Condensed Phase Systems. J. Phys. Chem. A 2016, 120, 7205–7212. (41) Fedorov, D. G.; Asada, N.; Nakanishi, I.; Kitaura, K. The Use of Many-Body Expansions and Geometry Optimizations in Fragment-Based Methods. Acc. Chem. Res. 2014, 47, 2846– 2856.


43


Page 44 of 61

(42) Fedorov, D. G.; Kitaura, K. Subsystem Analysis for the Fragment Molecular Orbital Method and its Application to Protein-Ligand Binding in Solution. J. Phys. Chem. A 2016, 120, 2218–2231. ( 43 ) Mochizuki, Y.; Fukuzawa, K.; Kato, A.; Tanaka, S.; Kitaura, K.; Nakano, T. A Configuration Analysis for Fragment Interaction. Chem. Phys. Lett. 2005, 410, 247–253. (44) Ishikawa, T.; Mochizuki, Y.; Amari, S.; Nakano, T.; Tokiwa, H.; Tanaka, S.; Tanaka, K. Fragment Interaction Analysis Based on Local MP2. Theor. Chem. Acc. 2007, 118, 937–945. (45) Asada, N.; Fedorov, D. G.; Kitaura, K.; Nakanishi, I.; Merz, K. M. An Efficient Method to Evaluate Intermolecular Interaction Energies in Large Systems Using Overlapping Multicenter ONIOM and the Fragment Molecular Orbital Method. J. Phys. Chem. Lett. 2012, 3, 2604–2610. ( 46 ) Watanabe, C.; Fukuzawa, K.; Okiyama, Y.; Tsukamoto, T.; Kato, A.; Tanaka, S.; Mochizuki, Y.; Nakano, T. Three- And Four-Body Corrected Fragment Molecular Orbital Calculations with a Novel Subdividing Fragmentation Method Applicable to Structure-Based Drug Design. J. Mol. Graphics Modell. 2013, 41, 31−42. (47) Ishikawa, T.; Ishikura, T.; Kuwata, K. Theoretical Study of the Prion Protein Based on the Fragment Molecular Orbital Method. J. Comput. Chem. 2009, 30, 2594–2601. (48) Okamoto, T.; Ishikawa, T.; Koyano, Y.; Yamamoto, N.; Kuwata, K.; Nagaoka, M. A Minimal Implementation of the AMBER-PAICS Interface for ab initio FMO-QM/MM-MD Simulation. Bull. Chem. Soc. Jap. 2013, 86, 210–222.


44

Page 45 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


( 49 ) Alexeev, Y.; Mazanetz, M. P.;

Ichihara, O.; Fedorov, D. G. GAMESS as a Free

Quantum-Mechanical Platform for Drug Research. Curr. Top. Med. Chem. 2012, 12, 2013–2033. (50) Mazanetz, M. P.; Ichihara, O.; Law, R. J.; Whittaker, M. Prediction of Cyclin-Dependent Kinase 2 Inhibitor Potency Using the Fragment Molecular Orbital Method. J. Cheminf. 2011, 3, 2. (51) Zhang, Q.; Yu, C.; Min, J.; Wang, Y.; He, J.; Yu, Z. Rational Questing For Potential Novel Inhibitors of FabK From Streptococcus Pneumoniae by Combining FMO Calculation, CoMFA 3D-QSAR Modeling and Virtual Screening. J. Mol. Model. 2011, 17, 1483–1492. (52) Hitaoka S.; Chuman H.; Yoshizawa K. A QSAR Study on the Inhibition Mechanism of Matrix Metalloproteinase-12 by Arylsulfone Analogs Based on Molecular Orbital Calculations. Org. Biomol. Chem. 2015, 13, 793–806. (53) Baker, M. Fragment-Based Lead Discovery Grows up. Nature Reviews Drug Discovery 2013, 12, 5-7 . (54) Ozawa, M.; Ozawa, T.; Ueda, K. Application of the Fragment Molecular Orbital Method Analysis to Fragment-Based Drug Discovery of BET (Bromodomain and Extra-Terminal Proteins) Inhibitors. J. Mol. Graph. Model. 2017, 74, 73-82. (55) von Hopffgarten, M.; Frenking, G. Energy Decomposition Analysis. WIREs: Comput. Mol. Sc. 2012, 2, 43–62.


45


Page 46 of 61

(56) Parrish, R. M.; Parker, T. M.; Sherrill, C. D. Chemical Assignment of Symmetry-Adapted Perturbation Theory Interaction Energy Components: The Functional-Group SAPT Partition. J. Chem. Theory Comput. 2014, 10, 4417–4431. (57) Su, P.; Li, H. Energy Decomposition Analysis of Covalent Bonds and Intermolecular Interactions. J. Chem. Phys. 2009, 131, 014102. ( 58 ) Kitaura, K.; Morokuma, K. A New Energy Decomposition Scheme for Molecular Interactions within the Hartree-Fock Approximation. Int. J. Quant. Chem. 1976, 10, 325–340. (59) Chen, W.; Gordon, M. S. Energy Decomposition Analyses for Many-Body Interaction and Applications to Water Complexes. J. Phys. Chem. 1996, 100, 14316–14328. (60) Fedorov, D. G.; Kitaura, K. Pair Interaction Energy Decomposition Analysis. J. Comput. Chem. 2007, 28, 222–237. (61) Fedorov, D. G.; Kitaura, K. Energy Decomposition Analysis in Solution Based on the Fragment Molecular Orbital Method. J. Phys. Chem. A 2012, 116, 704–719. (62) Green, M. C.; Fedorov, D. G.; Kitaura, K.; Francisco, J. S.; Slipchenko, L. V. Open-Shell Pair Interaction Energy Decomposition Analysis (PIEDA): Formulation and Application to the Hydrogen Abstraction in Tripeptides. J. Chem. Phys. 2013, 138, 074111. (63) Fedorov, D. G.; Kitaura, K. On the Accuracy of the 3-Body Fragment Molecular Orbital Method (FMO) Applied to Density Functional Theory. Chem. Phys. Lett. 2004, 389, 129–134. (64) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge DensityFunctional Tight-Binding Method (SCCDFTB). J. Chem. Theory Comput. 2011, 7, 931−948.


46

Page 47 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(65) Nishimoto, Y.; Fedorov, D. G.; Irle, S. Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method. J. Chem. Theory Comput. 2014, 10, 4801–4812. ( 66 ) Fedorov, D. G.; Kitaura, K.; Li, H.; Jensen, J. H.; Gordon, M. S. The Polarizable Continuum Model (PCM) Interfaced with the Fragment Molecular Orbital Method (FMO). J. Comput. Chem. 2006, 27, 976–985. (67) Nishimoto, Y.; Fedorov, D. G. The Fragment Molecular Orbital Method Combined with Density-Functional Tight-Binding and the Polarizable Continuum Model. Phys. Chem. Chem. Phys. 2016, 18, 22047–22061. ( 68 ) Fedorov, D. G.; Kitaura, K. Coupled-Cluster Theory Based Upon the Fragment Molecular-Orbital Method. J. Chem. Phys. 2005, 123, 134103. (69) González, R.; Suárez, F.; Bohórquez, J.; Patarroyo, A.; Patarroyo, M. E. Semi-empirical quantum evaluation of peptide - MHC class II binding. Chem. Phys. Lett. 2017, 668, 29–34. (70) Morao, I.; Fedorov, D. G.; Robinson, R.; Southey, M.; Townsend-Nicholson, A.; Bodkin, M. J.; Heifetz, A. Rapid and Accurate Assessment of GPCR-Ligand Interactions Using the Fragment Molecular Orbital Based Density-Functional Tight-Binding (FMO-DFTB) Method. J. Comput. Chem. 2017, 38, 1987–1990. (71) Nishimoto, Y.; Fedorov, D. G.; Irle, S. Third-Order Density-Functional Tight-Binding Combined with the Fragment Molecular Orbital Method. Chem. Phys. Lett. 2015, 636, 90–96. (72) Ishikawa, T.; Kuwata, K. Interaction Analysis of the Native Structure of Prion Protein with Quantum Chemical Calculations. J. Chem. Theory Comput. 2010, 6, 538–547.


47


Page 48 of 61

(73) Tanaka, S.; Watanabe, C.; Okiyama, Y. Statistical Correction to Effective Interactions in the Fragment Molecular Orbital Method. Chem. Phys. Lett. 2013, 556, 272–277. (74) Atkins, P. W. Physical Chemistry, Oxford University Press, Oxford, 1994. (75) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456–1465. (76) He, X.; Fusti-Molnar, L.; Cui, G.; Merz, K. M. Jr. Importance of dispersion and electron correlation in ab initio protein folding. J. Phys. Chem. B 2009, 113, 5290–5300. (77) Nakano, T.; Kaminuma, T.; Sato, T.; Fukuzawa, K.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment Molecular Orbital Method: Use of Approximate Electrostatic Potential. Chem. Phys. Lett. 2002, 351, 475. (78) Bernetti, M.; Cavalli, A.; Mollic, L. Protein-Ligand (Un)Binding Kinetics as a New Paradigm for Drug Discovery at the Crossroad Between Experiments and Modelling. Med. Chem. Commun. 2017, 8, 534–550. (79) Fedorov, D. G.; Kitaura, K. The Importance of Three-Body Terms in the Fragment Molecular Orbital Method. J. Chem. Phys. 2004, 120, 6832–6840. (80) 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., et al. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347–1363.


48

Page 49 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(81) Fedorov, D. G.; Olson, R. M.; Kitaura, K.; Gordon, M. S.; Koseki, S. A new Hierarchical Parallelization Scheme: Generalized Distributed Data Interface (GDDI), and an Application to the Fragment Molecular Orbital Method (FMO). J. Comput. Chem. 2004, 25, 872–880. (82) Nakano, T.; Kaminuma, T.; Sato, T.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment Molecular Orbital Method: Application to Polypeptides. Chem. Phys. Lett. 2000, 318, 614–618. (83) Nagata, T.; Fedorov, D. G.; Li, H.; Kitaura, K. Analytic Gradient for Second Order Møller-Plesset Perturbation Theory with the Polarizable Continuum Model Based on the Fragment Molecular Orbital Method. J. Chem. Phys. 2012, 136, 204112. (84) Su, P. F.; Li, H. Continuous and Smooth Potential Energy Surface for Conductorlike Screening Solvation Model Using Fixed Points with Variable Areas. J. Chem. Phys. 2009, 130, 074109. (85) Fedorov, D. G.; Slipchenko, L. V.; Kitaura, K. Systematic Study of the Embedding Potential Description in the Fragment Molecular Orbital Method. J. Phys. Chem. A 2010, 114, 8742–8753. (86) Peverati, R.; Truhlar, D. G. Improving the Accuracy of Hybrid Meta-GGA Density Functionals by Range Separation. J. Phys. Chem. Lett. 2011, 2, 2810–2817. (87) Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. J. Chem. Theory Comput. 2013, 9, 338–354.


49


Page 50 of 61

(88) Fedorov, D. G.; Alexeev, Y.; Kitaura, K. Geometry Optimization of the Active Site of a Large System with the Fragment Molecular Orbital Method. J. Phys. Chem. Lett. 2011, 2, 282– 288. (89) Pruitt, S. R.; Brorsen, K. R.; Gordon, M. S. Ab initio Investigation of the Aqueous Solvation of the Nitrate Ion. Phys. Chem. Chem. Phys. 2015, 17, 27027–27034. (90) Goldman, N.; Goverapet Srinivasan, S.; Hamel, S.; Fried, L. E.; Gaus, M.; Elstner, M. Determination of a Density Functional Tight Binding Model with an Extended Basis Set and Three-Body Repulsion for Carbon under Extreme Pressures and Temperatures. J. Phys. Chem. C 2013, 117, 7885–7894. (91) Bodrog, Z.; Aradi, B. Possible Improvements to the Self-Consistent-Charges DensityFunctional Tight-Binding Method within the Second Order. Phys. Status Solidi B 2012, 249, 259–269. (92) Ruedenberg, K.; Schmidt, M. W. Why Does Electron Sharing Lead to Covalent Bonding? A Variational Analysis. J. Comput. Chem. 2007, 28, 391–410. (93) Ji, C. G.; Mei, Y.; Zhang, J. Z. H. Developing Polarized Protein-Specific Charges for Protein Dynamics: MD Free Energy Calculation of pKa shifts for Asp26/Asp27 in Thioredoxin. Biophys. J. 2008, 95, 1080–1088. (94) Okiyama, Y.; Watanabe, H.; Fukuzawa, K.; Nakano, T.; Mochizuki, Y.; Ishikawa, T.; Ebina, K.; Tanaka, S. Application of the Fragment Molecular Orbital Method for Determination


50

Page 51 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of Atomic Charges on Polypeptides. II. Towards an Improvement of Force Fields Used for Classical Molecular Dynamics Simulations. Chem. Phys. Lett. 2009, 467, 417–423.


51


Figure 1.


Page 52 of 61

Page 53 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2.



a) DI+RC, DFTB/D

-1 ES, DFTB/D

0.1

y = 0.4243 x - 0.2090 R² = 0.9770

-2 -3 -4 -5

0

y = 0.1873x - 0.1035 R² = 0.9523

0.0 -0.1

Total、DFTB/D

0

-0.2 -0.3 -0.4 -0.5

-15

-10 -5 ES, CCSD(T)

0

y = 0.9346x + 0.1935 R² = 0.9945

-2 -4 -6 -8

-0.6

-6

-3

-2 -1 DI+RC, CCSD(T)

-8

0

-6

-4

-2

0

Total, CCSD(T）

b) DI+RC, CAM-B3LYP/D

y = 0.9730 x + 0.0081 R² = 0.9997

-5

-10

-15 -15

-10 -5 ES, CCSD(T)

0

0

0

y = 1.3363x - 0.0027 R² = 0.9992

-1

Total, CAM-B3LYP/D

0 ES, CAM-B3LYP/D

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

Page 54 of 61

-2 -3 -4

y = 0.9673x + 0.0409 R² = 0.9996

-2 -4 -6 -8

-3

-2 -1 DI+RC, CCSD(T)

0

Figure 3.


-8

-6 -4 -2 Total, CCSD(T)

0

Page 55 of 61

-256,700 potential energy, kcal/mol

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46


-256,750 -256,800 -256,850 -256,900 -256,950 neutral

anionic

-257,000 0

50

100 time, ps

150

200

Figure 4.



b) 25

20

20 energy, kcal/mol

25

15 10 5 0

15 10 5 0 -5

-10

-10 Asn-1 Leu-2 Tyr-3 Ile-4 Gln-5 Trp-6 Leu-7 Lys-8 Asp-9 Gly-10 Gly-11 Pro-12 Ser-13 Ser-14 Gly-15 Arg-16 Pro-17 Pro-18 Pro-19 Ser-20 ligand

-5

~

E I

Ekin0 E I

kin,0

E I kin

Asn-1 Leu-2 Tyr-3 Ile-4 Gln-5 Trp-6 Leu-7 Lys-8 Asp-9 Gly-10 Gly-11 Pro-12 Ser-13 Ser-14 Gly-15 Arg-16 Pro-17 Pro-18 Pro-19 Ser-20 ligand

a)

energy, kcal/mol

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

Page 56 of 61

~

E I

Figure 5.


Ekin0 E I

kin,0

E I kin

Page 57 of 61

330 temperature TI, K

320 310 300 290 280

isolated

neutral complex

anionic complex

ligand

Ser-20

Pro-19

Pro-18

Pro-17

Arg-16

Gly-15

Ser-14

Ser-13

Pro-12

Gly-11

Gly-10

Asp-9

Lys-8

Leu-7

Trp-6

Gln-5

Ile-4

Tyr-3

Leu-2

270 Asn-1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46


global T

Figure 6.



b)

30

50

15

25

0

0

-15

-25

-30

-50

-45

-75

-60

-100

0K E IJ0K

E0 E IJ0

E IJ

Emin E IJmin

Emax E IJmax

Asn-1 Leu-2 Tyr-3 Ile-4 Gln-5 Trp-6 Leu-7 Lys-8 Asp-9 Gly-10 Gly-11 Pro-12 Ser-13 Ser-14 Gly-15 Arg-16 Pro-17 Pro-18 Pro-19 Ser-20

a)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

Page 58 of 61

0K  E IJ0K

Figure 7.


E0  E IJ0

E IJ

Emin E IJmin

Emax E IJmax

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

20

0

-20

-40

-60

MP2 RHF/D CAM-B3LYP/D DFTB/D


a)

Interaction energy with ligand, kcal/mol


Interaction energy with ligand, kcal/mol

Page 59 of 61 The Journal of Physical Chemistry

b) 50

25 0

-25

-50

-75

-100

MP2

Figure 8.


RHF/D CAM-B3LYP/D DFTB/D


a)

b) 50

20 y = 0.7358x + 0.1709 R² = 0.9898

0

other theory. kcal/mol

other theory, kcal/mol

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

Page 60 of 61

-20 -40 -60

y = 0.9083x + 1.8608 R² = 0.9884 0

-50

-100 -60

MP2

-40 -20 0 DFTB residue-ligand interaction, kcal/mol RHF/D

CAM-B3LYP/D

20

線形 (MP2) MP2

-100

-50 0 DFTB residue-ligand interaction, kcal/mol

MP2

Figure 9.


RHF/D

CAM-B3LYP/D

50

線形 (MP2) MP2

Page 61 of 61 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46


Table of Contents Graphic (4.75 x 8.5 cm)


Pair Interaction Energy Decomposition Analysis for Density Functional

Recommend Documents