Efficient Geometry Optimization of Large Molecular Systems in

Subscriber access provided by The University of Melbourne Libraries

Article

Efficient Geometry Optimization of Large Molecular Systems in Solution Using the Fragment Molecular Orbital Method Hiroya Nakata, and Dmitri G. Fedorov J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.6b09743 • Publication Date (Web): 17 Nov 2016 Downloaded from http://pubs.acs.org on November 26, 2016

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 51

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

Efficient Geometry Optimization of Large Molecular Systems in Solution Using the Fragment Molecular Orbital Method Hiroya Nakata† and Dmitri G. Fedorov∗,‡ Department of Fundamental Technology Research, R and D Center Kagoshima, Kyocera, 1-4 Kokubu Yamashita-cho, Kirishima-shi, Kagoshima, 899-4312, Japan, and Research Center for Computational Design of Advanced Functional Materials (CD-FMat), National Institute of Advanced Industrial Science and Technology (AIST), 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan E-mail: [email protected]

Abstract The analytic gradient is derived for the frozen domain formulation of the fragment molecular orbital (FMO) method combined with the polarizable continuum model. The accuracy is tested in comparison to full FMO calculations for a representative set of systems in terms of the gradient accuracy, protein-ligand binding energies and optimized structures. The frozen domain method reproduced geometries optimized with full FMO within 0.03-0.09 Å in terms of reduced mean square deviations; whereas a single point gradient calculation is accelerated by the factor of 38 (Trp-cage protein in explicit solvent, PDB: 1L2Y) and 12 (crambin, PDB: ∗ To

whom correspondence should be addressed

† Department of Fundamental Technology Research,

R and D Center Kagoshima, Kyocera, 1-4 Kokubu Yamashitacho, Kirishima-shi, Kagoshima, 899-4312, Japan ‡ Research Center for Computational Design of Advanced Functional Materials (CD-FMat), National Institute of Advanced Industrial Science and Technology (AIST), 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan

1

ACS Paragon Plus Environment


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 51

1CRN). The method is applied to a geometry optimization of the K-Ras protein-ligand complex (4Q03) using two domain definitions, and the optimized structures are consistent with experiment. Pair interaction analysis is used to identify residues important in binding the ligand.

Introduction A number of efficient quantum-mechanical (QM) approaches have been developed, using linearscaling methods, 1–4 hybrid approaches 5,6 and fragmentation. 7–23 Although fully QM methods can be used in practice to estimate energies, for instance, protein-ligand binding energies, 24,25 the utility of fully QM methods for geometry optimizations and transition state search is much more limited, especially for systems as large as proteins, and often due to their cost one is forced to retort to refining geometries using other approaches, such as force fields. Biological systems are typically found in water, and one has to further increase the complexity of theory by adding a solvent model. Explicit models are in principle the most accurate but their added cost is substantial, which can be reduced with the use of QM-parametrized potentials used for solvent. 26 Alternatively, there are continuum solvent approaches such as the polarizable continuum model (PCM) 27 and the conductor-like screening model (COSMO), 28 as well as theories relying on a distribution function for solvent, such as the reference interaction site model (RISM). 29 Hybrid approaches 30–32 as well as fragmentation methods 33,34 can be combined with continuum models, such as PCM. The fragment molecular orbital (FMO) method 35–39 is a fragmentation approach, in which each fragment is calculated in the presence of the electrostatic potential (ESP) exerted by the rest of fragments, followed by calculations of fragment pairs. FMO has been applied to many large systems. 40–44 FMO has been combined with PCM for energy and gradient 45–48 and other solvation models. 49–51 Geometry optimizations 52 have been done with FMO using hybrid 53,54 approaches and full QM for small proteins 55 and inorganic nano rings. 56 Alternatively, one can use various partial geometry optimization schemes, 57,58 which benefit from fragmentation making calculations 2


Page 3 of 51

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


effective and accurate. The frozen domain (FD) approach to geometry optimizations 59–61 and spectra simulations 62 in FMO have been developed using the multilayer scheme 63 in gas phase. Although one can use explicit solvent, it is expensive to add several layers of water molecules around large molecular systems. In this study, energy and its analytic gradient are developed for the frozen domain (FD) approach and the dimer (D) approximation (FDD) and combined with PCM for an efficient treatment of solvent. The accuracy is evaluated on a set of test systems, comparing analytic and numerical gradients, optimized geometries and protein-ligand binding energies. Finally, the method is applied to a geometry optimization of K-Ras in the complex with GDP, for which an X-ray structure has been recently resolved with a high resolution of 1.2 Å(PDB: 4Q03). 64

Theory Introduction As shown in Figure 1, in FMO/FD the solute system S is divided fragmentwise into the frozen (F) and polarizable buffer (B) domains. B in turn is divided into the active domain A and the rest, b. Domain B is treated as the higher layer in the multilayer FMO. The solvent surface in PCM is also divided into frozen and polarizable buffer domains. Only some atoms in the active domain A are allowed to move during geometry optimizations. When they move, due to polarization effects the electronic state of other fragments also changes, and the approximation in FD is that the polarization effects are limited to B. FMO/FDD differs from FMO/FD in the neglect of interactions between fragments I and J in B, if neither I nor J are in A. As an example, one may want to optimize only ligand coordinates in a protein-ligand complex S. Then A will include only the ligand; b will be the binding pocket in the protein, F is the rest of the protein, and B will be the binding pocket and the ligand. A geometry optimization in FMO/FDD proceeds as follows (Figure 2). For the initial geometry, the whole system is calculated once at the FMO1 level, that is, fragments are calculated in the 3



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 4 of 51

embedding potential of the whole system, but no pair calculations are performed. At the end of this step, one obtains the polarized density of all fragments, and the induced apparent solvent charges (ASCs) on the whole solvent surface. These calculations are performed at the lower layer level (L1) for the whole system S. In the next second step, fragments in the polarizable buffer B are recalculated at the level of layer 2 (L2), in the embedding potential of both layer 2 (B) and layer 1 (F). ASCs for the solute surface surrounding B are recomputed in this step, and ASCs for the solute surface F are frozen. ASCs for B are determined using FMO/PCM equations based on the solute potential, which includes contributions from all fragments (both in B and F). In FDD only fragment pairs are calculated where at least one fragment is in A; that is, in terms of domain pairs, AF and AB. In FD, AF and BB pairs are computed. AF pairs are calculated using the electrostatic separated dimer approximation(ES-DIM), 65 whereas the pairs AB and BB are computed with self-consistent field (SCF) for close pairs and ES-DIM for the rest. At the end of this step, one obtains the energy and its gradient for fragments and pairs, and computes the total energy and gradient, as shown below. This gradient is then used to obtain the next geometry, and this is repeated until convergence.

Apparent solvent charges To obtain solvent charge for the whole system at the initial geometry, the following equation is used, which is the same equation as used in single layer FMO/PCM. Throughout this paper, FMO/PCM level is used. 46 In this method, ASCs are determined from fragment monomer densities neglecting pair corrections (1 in stands for one-body).

qtF = −

∑

∑

Cts−1VsK,L1 ,

(1)

s∈F+B K∈F+B

where C is the matrix related to distances between tesserae t and s; the scaling factor including the dielectric constant ε is included in C. In PCM, solvent surface is constructed as a union of atomic

4


Page 5 of 51

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


spheres, and each sphere is divided into pieces called tesserae, and a point charge is put into the center of each of them. The potential exerted by fragment K for layer 1 (L1) upon tessera t is VtK,L1 =

ZA t − ∑ DK,L1 µν wµν |R − R | t A µ ,ν A∈K

∑

(2)

where wtµν

1 ν . = µ |r − Rt |

(3)

A denotes atoms with the nuclear charges ZA . RA , Rt , and r are the coordinates of atoms, tessera centers, and electron, respectively. DK,L1 µν is the density matrix for fragment K at the level of L1. Solvent charges qtF are used in the Fock matrix as described elsewhere, 46 polarizing the solute; the polarization is mutual and the charge determination is done self-consistently (Figure 2). After ASCs and solute electron densities are obtained for the whole system at the level of layer L1, ASCs qB are recalculated in layer L2 for tessera t in B as follows,

∑ ∑ Cts−1VsK,L2,

(4)

ZA t − ∑ DK,L2 µν wµν µ ,ν A∈K |RA − Rt |

(5)

qtB = −

K∈B s∈B

where VtK,L2 =

∑

DK,L2 / B, qtB = 0. µν is the density matrix of K in L2. For t ∈ The final ASC q used in layer L2 is obtained as follows

q =qF + qB − qB(F)

5


(6)


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

B(F)

where qt

Page 6 of 51

= 0 if t ∈ / B, and for t ∈ B, B(F)

qt

=−

∑ ∑ Cts−1VsK,L1.

(7)

K∈B s∈B

Equation 6 for t ∈ S can be written as

qt = −

"

∑ Cts−1 ∑ VsK,L1 + ∑ Cts−1 ∑ s∈B

K∈F

s∈S

VsK,L2 −Vs

K∈B

K,L1

δt∈B

#

(8)

where δt∈B is 1 if t ∈ B and 0 otherwise. It is clear from eq 8 that the charge is determined by the total potential, with the explicit correction due to the second layer acting only on tesserae in B.

Energy definition The energy in FMO/FDD/PCM is E FDD = ∑ EI′L2 + I∈B

+

∑

∑

I>J I∈A,J∈B

′L2 ∆EIJ + Tr(∆DIJ,L2 VIJ )

′L2,L1 ∆EIJ + GL2 + Gcdr .

(9)

I∈A,J∈F ′L2 = E ′L2 − E ′L2 − E ′L2 . E ′L2 is the internal energy of monomers X = I or dimers X = IJ where ∆EIJ IJ I J X ′L2,L1 for layer L2. ∆EIJ is the interaction energy for pairs of fragments IJ, I in layer L2 and J in L1. IJ is the ESP acting on dimer ∆DIJ,L2 = DIJ,L2 − DI,L2 ⊕ DJ,L2 is the density transfer matrix and Vµν

IJ. GL2 is the electrostatic solute-solvent interaction in PCM. Gcdr is the cavitation (c), dispersion

(d) and repulsion (r) solute-solvent interaction terms, which are parametrized and independent on the electronic state.

6


Page 7 of 51

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


For FMO/FD/PCM, the expression is E FD = ∑ EI′L2 + I∈B

∑

+

∑

I>J I,J∈B

′L2 ∆EIJ + Tr(∆DIJ,L2 VIJ )

′L2,L1 ∆EIJ + GL2 + Gcdr .

(10)

I∈A,J∈F

Compared to eq 9, the difference is in the selection of dimers in the second sum. Below, FDD is assumed and FD equations can be constructed by replacing "I ∈ A, J ∈ B" with "I, J ∈ B". The internal fragment energy is EX′L2

=∑

µ ,ν

X DX,L2 µν hµν

1 1 X,L2 X,L2 + Dµν Dλ σ (µν |λ σ ) − (µλ |νσ ) , 2 µ ,ν∑ 2 ,λ ,σ

(11)

where hXµν and (µν |λ σ ) are one and two-electron integral terms, respectively. Throughout the paper, Greek (µ , ν etc) and Roman (i, j) letters denote atomic and molecular orbitals, respectively IJ is (however, t, s number tesserae). The ESP Vµν

IJ = Vµν

"

∑ ∑

K6=I,J A∈K

−ZA µ |r − R

A

ν + |

∑

λ ,σ ∈K

#

DλKσ (µν |λ σ ) .

(12)

The interaction energy between fragments in L2 and L1 is ′L2,L1 ∆EIJ

=

∑ ∑

µ ,ν ∈I A∈J

+

DI,L2 µν

∑ ∑

µ ,ν ∈I λ ,σ ∈J

−ZA −ZA J,L1 ν | ν + ∑ ∑ Dµν µ µ | r − RA |r − R | A µ ,ν ∈J A∈I

J,L1 DI,L2 µν Dλ σ ( µν |λ σ )

(13)

The electrostatic term GL2 in eq 9 is

GL2 =

1 T VI,L2 q + ∑ 2 I∈B

∑

T

∆VIJ,L2 q,

I>J I∈A,J∈B

7


(14)


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 8 of 51

where ∆VIJ,L2 = VIJ,L2 (DIJ,L2 ) − VI,L2 (DI,L2 ) − VJ,L2 (DJ,L2 )

(15)

V and q are vectors whose size is the number of tesserae. The electronic part of the solute potential acting on tessera t (the nuclear contribution cancels in pair corrections, see eq 5) is t VtX (DX,L2 ) = − ∑ DX,L2 µν wµν

(16)

µ ,ν

Analytic energy gradient The derivative of the total energy with respect to a nuclear coordinate a is

∂E ∂ EI′L2 =∑ + ∂ a I∈B ∂ a +

∑

∑

I>J I∈A,J∈B ′L2,L1 ∂ ∆EIJ

I∈A,J∈F

∂a

′L2 ∂ ∆EIJ + ∂a

+

∑

I>J I∈A,J∈B

IJ ∂ ∆DIJ,L2 µν Vµν ∑ ∂a µ ,ν

∂ GL2 ∂ Gcdr + . ∂a ∂a

(17)

The derivative of all terms except GL2 can be found in an earlier study. 66 One has to solve coupled-perturbed Hartree-Fock (CPHF) equations as in other FMO gradients, which are modified as given below. Using eq 6, GL2 can be divided into contributions from the domains B and F. i.e.,

GL2 =GB,B + GB,F GB,B = GB,F =

1 T VI,L2 qB + ∑ 2 I∈B

(18)

∑

T

∆VIJ,L2 qB ,

(19)

I>J I∈A,J∈B

1 T ∑ VI,L2 (qF − qB(F)) + 2 I∈B

∑

T

∆VIJ,L2 (qF − qB(F) ).

I>J I∈A,J∈B

8


(20)

Page 9 of 51

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 derivative of the first term in eq 18 is T

∂ GB,B ∂ VI,L2 B =∑ q + ∂a ∂a I∈B

T

∑

I>J I∈A,J∈B

∂ ∆VIJ,L2 B q ∂a

1 ∂ Cts B + ∑ ∑ qtB q + ∑ ∑ ∆qtB 2 I∈B t,s∈B ∂ a s K∈B t∈B

"

# ∂ Cts B ∂ VtK,L2 qs + , ∑ ∂a s∈B ∂ a

(21)

where ∆qtB = − ∑ Cts−1 s∈B

∑

∆VIJ,L2 . s

(22)

I>J I∈A,J∈B

Note that the PCM derivative for domain B contains only contributions that come from the ASCs of domain B. The derivative of the second term in eq 18 is

∂ GB,F ∂a

T

=VB (qF − qB(F) ) + VB

VB =

1 VI,L2 + ∑ 2 I∈B

∑

F − qB(F) ∂ q T

∂a

∆VIJ,L2 .

(23) (24)

I>J I∈A,J∈B

where VB is the potential used in the evaluation of G. In contrast to the B-B interaction derivatives

∂ GB,B /∂ a, the derivatives of B-F terms, ∂ GB,F /∂ a, run over all ASC in the system. However, the contribution to the gradient for atoms in the active domain is usually small, because the density matrix and apparent surface charges are frozen during geometry optimization. Thus, the last term T VB ∂ qF − qB(F) /∂ a in eq 23 is neglected in the implementation. The derivative of the solute potential in eq 21 is

all occ ∂ VtK,L2 ∂ ZA t,a a,K t =∑ − ∑ DK,L2 µν wµν − 4 ∑ ∑ Umi wmi ∂a m i µ ,ν A∈K ∂ a |RA − Rt |

9


(25)


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 10 of 51

t where wtmi is the MO-based wtµν , and wt,a mi is the derivative of wmi . The derivatives of ASCs are

# " L1,K ∂ qtF ∂ Cts−1 L1,K ∂ V s = ∑ ∑ − ∂ a Vs −Cts−1 ∂ a ∂ a s∈F+B K∈F+B # " B(F) L1,K ∂ qt ∂ Cts−1 L1,K ∂ V s . =∑ ∑ − Vs −Cts−1 ∂a ∂ a ∂ a s∈B K∈B

(26) (27)

Thus, the last term in eq 23 is neglected in the implementation, but its derivation is provided: BT

V

∂ qF ∂ qB(F) − ∂a ∂a

!

=

∑

B(B) ∂ Cts F qs

qt

∂a

t,s∈F+B

−

∑

t,s∈B

+

B(B) ∂ Cts B(F) qt qs −

∂a

∑

∑

L1,K B(B) ∂ Vt

qt

∂a

t∈F+B K∈F+B

∑∑

L1,K B(B) ∂ Vt qt ,

t∈B K∈B

∂a

(28)

where qB(B) = − C−1 VB .

(29)

CPHF equations The derivatives of total energy (eq 17) and the potential (eq 25) contain the unknown response a,I terms Umi , which can be obtained by solving the CPHF equations with the solvent potential.

The Fock matrix for molecular orbitals i and j ∈ I is occ

F˜iIj =hIi j + ∑{2(i j|kk) − (ik| jk)} +ViIj +WiIj k −qt I j , Wi j = ∑ i |r − Rt | t

10


(30) (31)

Page 11 of 51

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


Therefore, the derivative of the Fock matrix 67 is all occ all ∂ F˜iIj all a,I I I + Fia,I + = ∑ Umi Fjm + ∑ Uma,Ij Fim ∑ ∑ Umla,I A′i j,ml j ∂a m l m m all occ

−

∑ ∑∑

a,K I,K Uml Ai j,ml +

K6=I m∈K l∈K

∂ WiXj , ∂a

(32)

where A′i j,ml and AX,K i j,ml are the internal (same fragment) and external orbital Hessian terms for fragment ESP: A′i j,kl =4(i j|kl) − (ik| jl) − (il| jk)

(33)

AI,K i j,kl = − 4(i j|kl).

(34)

The last term in eq 32 is PCM specific (see supporting information for details), all all ∂ WiXj a,X X t,a X = ∑ Umi W jm + ∑ Uma,X j Wim − ∑ wi j qt ∂a m m t

+

∑

∂ Cst X,i j qt + qBs ∂a

t,s∈B X,i j

where qt

∂ VtI X,i j ∑ ∑ qt , t∈B I∈B ∂ a N

(35)

corresponds to the solvent charge induced by an interaction of a pair of orbitals i and j X,i j

qt

= ∑ Cts−1 wsij

(36)

s∈B

and the derivative of the solute potential is

∂ VtI ∂ =∑ ∂a A ∂a

all occ ZA − ∑ DIµν wt,a − µν ∑ ∑ 4Ukla,I wtkl |RA − Rt | µ ,ν k l

all occ

=Vta,I − ∑ ∑ 4Ukla,I wtkl k

(37)

l

Using the diagonality of the Fock matrix ∂ F˜iIj /∂ a = 0 for i 6= j, one obtains CPHF equations

11



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 12 of 51

for FMO/FDD/PCM vir occ

a,I ∑ ∑ AI,I i j,kl Ukl + k

l

∑

vir occ

vir occ

I,K,PCM a,K a,K Ukl =BI,a ∑ ∑ AI,K ij , i j,kl Ukl + ∑ ∑ ∑ Ai j,kl

K6=I k

K k

l

(38)

l

where a,I I ˜ a,I BI,a i j = − Si j ε j + Fi j −

occ 1 occ a,I ′ 1 occ a,I I,K,PCM 1 a,I I,K S A − S A − ∑ kl i j,kl 2 ∑ ∑ Skl Ai j,kl , kl i j,kl 2∑ 2 K6∑ K k,l k,l =I k,l

(39)

and F˜ia,I j is the integral derivative of the Fock matrix in eq 32. The orbital Hessian related terms are I I ′ AI,I i j,kl =δik δ jl ε j − εi − Ai j,kl ,

AiI,K,PCM =4 j,kl

∑

wti jCts−1 wskl

(40) (41)

t,s∈B

The PCM orbital Hessian terms in eq 41 are needed only for fragments in the B domain, thus the additional computational cost for solving eq 38 is comparable to gas phase. In practice, the Z-vector method is used 67 for adding response terms to the gradient.

Computational details FMO/FDD/PCM was implemented into a development version of GAMESS 68,69 parallelized with the generalized distributed data interface (GDDI). 70,71 The lower and higher layers were described with STO-3G and 6-31G(d), respectively. Hartree-Fock (HF) was used for the lower level and most of higher layer calculations, except that density functional theory (DFT) was used (LC-BLYP and LC-BOP functionals) for the higher layer where indicated. It has been shown that RHF/D works well for describing dispersion in proteins and other systems. 72,73 For accuracy tests, the default Lebedev grid was used, and in the K-Ras protein-ligand complex optimization, SG-1 grid was employed. In geometry optimizations all atoms in the active domain A were optimized. To evaluate the gradient accuracy, three systems were used (Figure 3), (a) an organic molecule 12


Page 13 of 51

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 a water cluster (divided as one molecule per fragment), 2,2,6,6-tetramethyl-cyclohexanone (TEMHEX), (b) an oligomer representing the cationic polymer 74 hexyl acrylate:N-(3-acryloyloxy)propyl)N,N,N-triprorpylammonium) mixed in the ratio of (6:1) and divided in FMO as 2 units per fragment, (c) a complex of the oligomer with an anionic ligand B(Phe)4. The domain definitions are shown in Figure 4. The TEMHEX molecule was assigned to the domain A, and the domain B also included water molecules within 5.2 Å from TEMHEX. In the cationic polymer, one fragment including a quaternary ammonium cation was selected as the domain A and fragments within 5.2 Å were assigned to B. For the oligomer complex, the ligand was selected as the domain A, and the domain B was assigned fragments within 5.2 Å from the ligand. The size of B (RB ) mentioned below refers to including all fragments within this distance from A. For testing the accuracy of optimized structures and measuring timings, the Trp-cage miniprotein construct (PDB: 1L2Y) immersed into 161 water molecules (total 787 atoms) and the crambin protein (PDB: 1CRN), 642 atoms, were used. The domain definitions are shown in Figure 5. For Trp-cage, Trp-6 was selected as the active domain A and all residues and water molecules within a set of distances (3.9, 5.2, and 6.5 Å ) from Trp-6 were chosen as domains B. For the crambin protein, Tyr-44 was selected as the active domain, and all residues within a set of distances (3.9, 5.8, and 7.2 Å ) from Tyr-44 were assigned to domain B. For these two systems, the effect of the size of the domain B on the optimized structure and timings was investigated. The accuracy in reproducing protein-ligand binding energies was evaluated for Trp-cage and two ligands, neutral 4-hydroxybenzoic acid and anionic 4-hydroxybenzoate. The active domain included the ligand and one residue: Pro-18 and Gln-5 for the neutral and anionic ligands, respectively. All residues within a set of distances (3.9, 5.2 and 6.5 Å ) from the ligand were included into domain B. The complex of K-Ras with its ligand 4-Bromophenyl (Figure 6) containing 2753 atoms was optimized as a demonstrative application. The Grimme dispersion model was used. 75 Two domain setups were used: (a) the ligand was assigned to A and all residues within 5.2 Å from it were

13



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 14 of 51

assigned to B and (b) the ligand and a part of the binding pocket of the protein (Cys-39, Ile-55, Val-7, and Thr-74) were assigned to A, and all residues within 5.2 Å from it were assigned to B. Proteins were fragmented as 1 residue per fragment, with the exceptions as follows. Sulphurbridged residues were merged into 1 fragment. In Trp-cage, Pro-18 and Gln-5 were merged with the ligand into one fragment for the neutral and anionic ligand, respectively. In K-Ras, the ligand is covalently bound to Cys39 by disulfide bond, and ligand is merged into 1 fragment with Cys39. The initial structures for proteins were taken from PDB and optimized with the AMBER force field 76 using AmberTools15. 77 The initial structures of other systems were constructed with JOCTA program 78 and optimized using COGNAC solver 79 with the general AMBER force field (GAFF). 80 To make FMO input files, Fragit 81 program was used for proteins using automatic fragmentation, and Facio 82 program was used for the oligomer using manual fragmentation. The effect of changing the tessera density was tested and found to be small due to the accuracy of the FIXPVA tesselation; 83 but to eliminate the tesselation as the source of error, all gradient accuracy tests for TEMHEX were performed with 240 tesserae per sphere; in other calculations 60 tesserae were used. C-PCM was used with the default (van-der-Waals) atomic radii. The cavitation (ICAV=1) and dispersion+repulsion (IDISP=1) Gcdr terms were added. The electrostatic point charge (ES-PTC) approximation 65 was used in the ESP with the unitless threshold of 3.5 except in the K-Ras optimization, where 3.0 was used to improve efficiency. The ESDIM approximation was used in all calculations with the unitless threshold of 2.0. The thresholds for geometry optimizations were 10−4 and 1/3 × 10−4 hartree/bohr for the maximum and RMS gradient values, respectively (for the geometry optimization of the ligand and binding pocket of K-Ras these thresholds were multiplied by 3). Numeric gradients were evaluated by double differencing with the step size of 0.0005 Å.

14


Page 15 of 51

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


Results and discussion Accuracy of the analytic energy gradient The accuracy of analytic energy gradient is shown in Table 1, computed as the difference between analytic and numerical gradients. For comparison, accuracy results for B in FDD taken to be the whole system S (i.e., no F), and FMO without FDD are also shown. It should be noted that FDD with B=S is not the same as full FMO. This is because in FDD, dimer contributions are ignored unless one of fragments is in A. The accuracy of FD and FDD gradients is quite similar. For TEMHEX a comparison of HF and LC-BLYP results is given, and the accuracy of DFT is very similar to HF. Also, either increasing the size of the domain B or switching to full FMO improves the gradient accuracy. This is because in the implementation of the analytic gradient of FMO/FDD/PCM, some small terms involving solvent were neglected, see eq 28 . The accuracy of the oligomer is similar to that of TEMHEX, and the complex has a higher accuracy (possibly, because the system adding an anionic ligand to the cationic polymer neutralizes it locally, although it is actually a zwitterion as the charges of the oligomer and ligand are 4 and -1, respectively).

Geometry optimizations Energy profiles for geometry optimizations of Trp-cage are shown in Figure 7. It can be seen that all curves are nearly parallel to each other and no oscillations are observed near convergence. The accuracy of optimized geometries evaluated for Trp-cage and crambin proteins, in comparison with full FMO2 is summarized in Table 2. Here, the effect of the size of the domain B is investigated. As expected, the accuracy increases with the size of B. In all cases RMSD values do not exceed 0.03 Å. Both Trp-cage and crambin show a residual error for B=S. This error comes from the neglect of dimers not including a fragment in A. The dependence of the results on the size of B is due to the polarization of the protein, affected by a partial geometry optimization of a subsystem. The protein polarization includes the change in 15



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 16 of 51

the interactions between residues and effects such as charge transfer within the protein. 84 The computational time for a single step increases with the size of B; taking B to be the whole system S makes the calculations about 4 times slower compared to B of 3.9 Å. Comparing FMO/FDD with the size of B equal to 3.9 Å to full FMO, the acceleration factor is 38 and 12 for Trp-cage and crambin, respectively. It can be seen that FMO/FDD does indeed make the calculations much faster. Comparing the timings of FD and FDD, the latter is faster about 2-6 times; the difference between Trp-cage (factor of 6) and crambin (factor of 2) lies in the fact that the former has many water molecules (small fragments). For the size of B equal to 3.9 in Trp-cage, FD was slightly less accurate than FDD. Summarizing, the size of B in FD and FDD should be at least 5-6 Å.

Protein-ligand binding energy The protein-ligand binding energy is computed using FMO/PCM subtracting the energies of the protein and ligand from the energy of the protein-ligand complex, E binding = Ecomplex − Eprotein − Eligand

(42)

The accuracy of the binding energies for Trp-cage is evaluated as a function of the size of B. The results are summarized in Table 3. Most presented results are for FDD so this method is discussed in detail. The protein polarization effect is substantial, and a difference of 2.02 kcal/mol is observed for the neutral ligand, comparing 3.9 Å and B=S. For the anionic ligand, the use of a small polarizable buffer B resulted in a large error of 21.15 kcal/mol in the binding energy. The cause of this was investigated. It was found that using a larger basis set for the second layer (i.e., for B), 6-31G*, does not substantially improve the results. The main reason for the large error is the fact that several residues assigned to the frozen domain F, most importantly, negatively charged Asp-9, are very strongly polarized by the ligand, which is properly accounted for when Asp-9 is placed in

16


Page 17 of 51

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


B (e.g., for the size of B of 5.2 Å). FMO/FDD was thus found to be prone to large errors in the binding energies for charged ligands when B is small. To separate the effect of geometry optimization from the binding energy at a given geometry, single point energies were evaluated with full FMO/PCM using geometries optimized with FD and FDD. The results, shown in parentheses of Table 3, are very accurate, the errors do not exceed 0.4 kcal/mol. FMO/FDD can be used for a reliable geometry optimizations even with small domains for charged ligands and produce good structures; however, the energies reflecting protein polarization may have to be refined with a single point FMO/PCM calculation. Using FD instead of FDD shows that although the error in the binding energy is improved a little for all energies evaluated with FD, the values refined with full FMO at the FD structures are somewhat worse than the corresponding FDD results, which indicates that 3.9 Å, especially for charged ligands, is not an appropriate size.

Application to a protein-ligand complex Two domain definitions were used, when only ligand was optimized and when ligand and a part of the binding pocket were optimized. The results are shown in Figure 8. For the smaller active domain A, the RMSD between the optimized and X-ray structures is 0.199 Å (calculated for optimized heavy atoms). The experimental structure was resolved with high resolution, and FMO/FDD/PCM geometry optimization reproduced it well. A small deviation is found for the benzene ring. For the larger active domain A, the RMSD between the optimized and X-ray structures is 0.241 Å (also calculated for optimized heavy atoms). It can be seen in Figure 8 that the structure optimized with FMO well reproduces the experimental one. Finally, details protein-ligand binding were clarified using pair interaction energies in FMO/FDD/PCM, defined as ′L2 DISP ∆EIJ = ∆EIJ + Tr(∆DIJ,L2 VIJ ) + ∆EIJ + ∆GSOLV IJ

17


(43)


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 18 of 51

DISP , evaluated with the Grimme’s method, 75 is explicitly separated. The where the dispersion ∆EIJ

solvent screening term ∆GSOLV is obtained with a decomposition of GL2 + Gcdr in eq 9 into fragIJ ment contributions. 85 Main interactions are plotted in Fig. 9 at the structure optimized with the larger A, including a part of the binding pocket. It can be seen that the largest attractive interaction is with residue fragments Leu-56, Ile-55, Met-72, Val-7, Lys-5 and Asp-57. Glu-76 has a large repulsive interac′L2 is about 11.1 kcal/mol, and tion with the ligand; looking at components, one can see that ∆EIJ

the explicit charge transfer Tr(∆DIJ,L2 VIJ ) is -4.5 kcal/mol. The ligand, which has a substantial dipole moment of 7.3 D, apparently perturbs the charge distribution in the charged residue fragment Glu-76, and the interaction is not electrostatic in nature, as verified by additional calculations (the electrostatic interaction between Glu-76 and ligand is about 0.4 kcal/mol). There was no substantial charge transfer between the ligand and Glu-76 (0.0015 a.u.). Apparently, due to a steric repulsion with the rest of the binding pocket, the ligand and Glu-76 cannot adjust their geometry to reduce the repulsion. In the attractive interactions, dispersion plays an important role; it is the main contributor for Leu-56. As can be seen in Figure 8, Leu-56 is arranged in a very special way, by having several hydrogen atoms face the benzene ring of the ligand; this results in a large gain of the attractive dispersion interaction between Leu-56 and ligand. Other interactions are mainly polar but some (e.g., Ile-55) have a substantial dispersion.

Conclusions The analytic gradient for FMO/PCM combined with the FD and FDD models has been derived and implemented into GAMESS (in the implementation some small terms are neglected). The accuracy of the method has been tested in terms of the gradient accuracy, optimized geometries and binding energies. Neutral ligands can be accurately calculated with a relatively small polarizable buffer, and

18


Page 19 of 51

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


charged ligands need a larger buffer, consistent with earlier observations. 59 The protein polarization effects in B, such as the interactions between residues are found to make noticeable contributions to the energies. In particular, charged ligands are prone to strongly polarize the protein and affect the interactions between residues in the protein. FMO/FD takes these protein polarization effects into account within B, whereas in FMO/FDD some effects in B are ignored; neither method accounts for polarization effects in F, which may be important if the size if B is small. However, interactions in F affect gradients for atoms in F, which are not relevant for geometry optimizations of atoms in A. It has been shown that FMO/FDD produces accurate geometries; but in general it is better to recalculate final energies with full FMO. The method has been applied to study a realistic protein-ligand complex K-Ras and the optimized structure is found to well reproduce the high resolution experimental structure. The interaction energy analysis identified the main residues binding the ligand, and one residue, which has a repulsive interaction. The developed FMO/FDD/PCM can be used to refine structures of large molecular systems, and make binding energy analysis 84–86 based on FMO more reliable and useful in drug discovery. 87,88 The method can be also used for optimizing proton positions, for instance, in pKa determination.

Acknowledgments A part of calculations were performed with supercomputer system for the Foundation for Computational Science (FOCUS), and computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID:hp160276)

Supporting Information Available Derivations of CPHF equations for FMO/FDD/PCM are provided. This material is available free of charge via the Internet at http://pubs.acs.org. 19



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 20 of 51

Table 1: Root mean square deviation (RMSD) and maximum error (Max error) in hartree/bohr of the analytic gradient vs numerical for systems shown in Figure 3. The size of the domain B, RB , is in Å (all means B=S). system method higher layera RB RMSD Max error TEMHEX FMO/FD HF 5.2 0.000203 0.000567 TEMHEX FMO/FDD HF 5.2 0.000203 0.000710 TEMHEX FMO/FDD LC-BLYP 5.2 0.000252 0.000752 TEMHEX FMO/FDD HF all 0.000044 0.000133 TEMHEX FMO (HF) - 0.000034 0.000084 oligomer FMO/FDD HF 5.2 0.000322 0.000971 oligomer complex FMO/FDD HF 5.2 0.000019 0.000059 a The lower layer is always HF; for FMO without frozen domain there is only one layer, HF. Table 2: Maximum and root mean square deviations (RMSD) between FMO/FDD and FMO optimized geometries in Å, as a function of size of the domain B, RB (Å); the timings T are given in minutes for a single point gradient on 16 CPU cores of Xeon E5-2630 v4 (2.2 GHz) used as one group in GDDI. RB method Max RMSD T Trp-cage+(H2 O)161 3.9 FMO/FD 0.08792 0.0172 155.0 3.9 FMO/FDD 0.06321 0.0210 26.5 5.2 FMO/FDD 0.05014 0.0158 59.2 6.5 FMO/FDD 0.05431 0.0173 60.7 all FMO/FDD 0.04278 0.0154 106.7 FMO 1018.8 Crambin 3.9 FMO/FD 0.0191 0.0037 107.8 3.9 FMO/FDD 0.0704 0.0299 63.5 5.8 FMO/FDD 0.0669 0.0277 75.7 7.2 FMO/FDD 0.0540 0.0259 88.2 all FMO/FDD 0.0354 0.0130 233.4 FMO 769.9

20


Page 21 of 51

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: Protein-ligand binding energies (kcal/mol) for Trp-cage and neutral and anionic ligands at the geometries optimized with the method shown, for several sizes of B, RB (in Å). a RB method neutral anionic 3.9 FMO/FD -11.07 (-11.44) 5.41 (-13.49) 3.9 FMO/FDD -11.20 (-11.96) 9.40 (-7.47) 5.2 FMO/FDD -10.26 (-12.25) -11.63 (-7.67) 6.5 FMO/FDD -10.49 (-12.06) -10.86 (-7.75) all FMO/FDD -9.18 (-12.06) -10.75 (-7.70) FMO (-12.37) (-7.77) a The energies are computed with the method shown; the values in parentheses are evaluated (refined) with FMO/PCM.

21



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 22 of 51

Figure captions.

Figure 1: Frozen F, polarizable buffer B and active A domains in FMO/FD are shown as purple, orange and red areas, respectively. Green ellipses indicate solvent.

22


Page 23 of 51

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: Computational procedure in FMO/PCM combined with FD or FDD.

23



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 24 of 51

(a) O

(b)

(c)

B

Figure 3: Systems used in accuracy tests: (a) 2,2,6,6-tetramethyl-cyclohexanone, (b) oligomer, (c) ligand forming complex with oligomer.

24


Page 25 of 51

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


(a)

(b)

(c)

Figure 4: Domain definitions for gradient accuracy tests: (a) 2,2,6,6-tetramethyl-cyclohexanone surrounded with water molecules, (b) oligomer, (c) oligomer complex. Frozen F, polarizable buffer B and active A domains in FMO/FD are shown in green, blue and red, respectively.

25



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 26 of 51

(a)

(b)

Figure 5: Domain definitions for geometry optimization tests: (shown for B with the radius of 5.2 Å): (a) Trp-cage protein in a cluster of explicit water, (b) Crambin protein. Frozen F, polarizable B and active A domains in FMO/FD are shown in green, blue and red, respectively.

26


Page 27 of 51

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 6: Frozen F, polarizable B and active A domains in FMO/FD are shown in green, blue and red, respectively, for the K-ras protein-ligand complex.

27



0 -1 Energy change

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 28 of 51

-2 -3 -4 -5 -6 -7 0

5

10 15 Optimization step

20

25

Figure 7: Energy change during geometry optimizations of the protein-ligand complex of Trp-cage and its neutral ligand. FMO/FDD/PCM results with the 5.2 Å size of B and full B=S are shown in red and blue, respectively. Full FMO/PCM results are in green.

28


Page 29 of 51

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


(a)

(b)

VAL-7 LEU-56 Ligand LYS-5

CYS-39 THR-74 Figure 8: Protein-ligand complex of K-Ras: (a) only ligand was optimized; experimental and theoretical structures are shown in red and blue, respectively (atoms not included in the optimization are shown in cyan) and (b) ligand and some residues in the binding pocket were optimized; experimental structure is shown in blue and theoretical structure is colored by atom (red: O, grey: C, dark red: Br, white: H and yellow: S); atoms included in the calculations, but only optimized atoms are shown.

29



∆Ε’IJ

∆ΕdispIJ

Tr(∆DIJVIJ)

∆GsolvIJ

total ∆ΕIJ 15 Pair interaction energy (kcal/mol)

10 5 0 -5 -10 -15 GLU-76

ASP-57

LYS-5

VAL-7

MET-72

ILE-55

LEU-56

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 30 of 51

Figure 9: Important interactions of the ligand with the residue fragments in the protein-ligand complex of K-Ras.

30


Page 31 of 51

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


References (1) Mezey, P. G.; Leszczynski, J., Eds.; Linear-Scaling Techniques in Computational Chemistry and Physics.; Springer: New York, 2011. (2) Reimers, J. R., Ed.; Computational methods for large systems: Electronic structure approaches for biotechnology and nanotechnology; Wiley: New York, 2011. (3) Nakata, A.; Bowler, D. R.; Miyazaki, T., Efficient calculations with multisite local orbitals in a large-scale DFT Code CONQUEST. J. Chem. Theory Comput. 2014, 10, 4813-4822. (4) Wilkinson, K. A.; Hine, N. D. M.; Skylaris, C.-K. Hybrid MPI-OpenMP parallelism in the ONETEP linear-scaling electronic structure code: application to the delamination of cellulose nanofibrils. J. Chem. Theory Comput. 2014, 10, 4782-4794. (5) Warshel, A.; Karplus, M. Calculation of ground and excited state potential surfaces of conjugated molecules. I. Formulation and parametrization. J. Am. Chem. Soc. 1972, 94, 5612-5625. (6) Dapprich, S.; Komáromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J. A new ONIOM implementation in Gaussian98. Part I. The calculation of energies, gradients, vibrational frequencies and electric field derivatives. J. Mol. Str.: THEOCHEM 1999, 461, 1-21. (7) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation methods: A route to accurate calculations on large systems. Chem. Rev. 2012, 112, 632–672. (8) 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 CH2 O molecules. Chem. Phys. 1975, 8, 192–200. (9) Gao, J. Toward a molecular orbital derived empirical potential for liquid simulations. J. Phys. Chem. B 1997, 101, 657–663.

31



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 32 of 51

(10) He, X.; Merz, K. M. Divide and conquer Hartree-Fock calculations on proteins. J. Chem. Theory Comput. 2010, 6, 405–411. (11) 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. (12) Friedrich, J.; Yu, H.; Leverentz, H. R.; Bai, P.; Siepmann, J. I.; Truhlar, D. G. Water 26mers 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, 666-670. (13) Liu, J.; Wang, X.; Zhang, J. Z. H.; He, X. Calculation of protein-ligand binding affinities based on a fragment quantum mechanical method. RSC Adv. 2015, 5, 107020-107030. (14) Andreji´c, M.; Ryde, U.; Mata, R. A.; Söderhjelm, P. Coupled-cluster interaction energies for 200-atom host–guest systems. Chem. Phys. Chem. 2014, 15, 3270-3328. (15) Steinmann, C.; Fedorov, D. G.; Jensen, J. H. The effective fragment molecular orbital method for fragments connected by covalent bonds. PLOS One 2012, 7, e41117. (16) Bates, D. M.; Smith, J. R.; Tschumper, G. S. Efficient and accurate methods for the geometry optimization of water clusters: Application of analytic gradients for the two-body:Many-body QM:QM fragmentation method to (H2 O)n . J. Comput. Theor. Chem. 2011, 7, 2753-2276. (17) Huang, L.; Massa, L. Quantum kernel applications in medicinal chemistry. Future Medicinal Chemistry 2012, 4, 1479–1494. (18) Gurunathan, P. K.; Acharya, A.; Ghosh, D.; Kosenkov, D.; Kaliman, I.; Shao, Y.; Krylov, A. I.; Slipchenko, L. V. Extension of the effective fragment potential method to macromolecules. J. Phys. Chem. B 2016, 120, 6562-6657. 32


Page 33 of 51

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


(19) Kiewisch, K.; Jacob, C. R.; Visscher, J. Quantum-chemical electron densities of proteins and of selected protein sites from subsystem density functional theory. J. Chem. Theory Comput. 2013, 9, 2425-2440. (20) 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. (21) Sahu, N.; Gadre, S. R. Vibrational infrared and Raman spectra of polypeptides: Fragmentsin-fragments within molecular tailoring approach. J. Chem. Phys. 2016, 144, 114113. (22) Collins, M. A. Systematic fragmentation of large molecules by annihilation. Phys. Chem. Chem. Phys. 2012, 14, 7744–7751. (23) Ab initio molecular dynamics with intramolecular noncovalent interactions for unsolvated polypeptides. Zhang, L.; Li, W.; Fang, T.; Li, S. Theor. Chem. Acc. 2016, 135, 34. (24) Merz, K. M. Jr. Using quantum mechanical approaches to study biological systems. Acc. Chem. Res. 2012, 47, 2804-2811. (25)

Ryde, U.; Söderhjelm, P.

Ligand-binding affinity estimates supported by quantum-

mechanical methods. Chem. Rev. 2016, 116, 5520-5566. (26) Gordon, M. S.; Smith, Q. A.; Xu, P.; Slipchenko, L. V. Accurate first principles model potentials for intermolecular interactions. Ann. Rev. Phys. Chem. 2013, 64, 553-578. (27) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999-3094. (28) Klamt, A.; Schuurmann, G. COSMO: A new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc., Perkin Trans. 1993, 2, 799-805.

33



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 34 of 51

(29) Ten-no, S.; Hirata, F.; Kato, S. A hybrid approach for the solvent effect on the electronic structure of a solute based on the RISM and Hartree-Fock equations. Chem. Phys. Lett. 1993,

214, 391 - 396. (30) Vreven, T.; Mennucci, B.; da Silva, C. O.; Morokuma, K.; Tomasi, J. The ONIOMPCM method: Combining the hybrid molecular orbital method and the polarizable continuum model for solvation. Application to the geometry and properties of a merocyanine in solution. J. Chem. Phys. 2001, 115, 62-72. (31) Cui, Q. Combining implicit solvation models with hybrid quantum mechanical/molecular mechanical methods: A critical test with glycine. J. Chem. Phys. 2002, 117, 4720-4728. (32) Li, H.; Pomelli, S. C.; Jensen, H. J. Continuum solvation of large molecules described by QM/MM: a semi-iterative implementation of the PCM/EFP interface. Theor. Chem. Acc. 2003, 109, 71–84. (33) Mei, Y.; Ji, C.; Zhang, J. Z. H. A new quantum method for electrostatic solvation energy of protein. J. Chem. Phys. 2006, 125, 094906. (34) Söderhjelm, P.; Kongsted, J.; Ryde, U. Ligand affinities estimated by quantum chemical calculations. J. Chem. Theory Comput. 2010, 6, 1726–1737. (35) 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. (36) Fedorov, D. G.; Kitaura, K., Eds.; The fragment molecular orbital method: Practical applications to large molecular systems; CRC press: Boca Raton, FL, 2009. (37) 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.

34


Page 35 of 51

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


(38) Fedorov, D. G.; Nagata, T.; Kitaura, K. Exploring chemistry with the fragment molecular orbital method. Phys. Chem. Chem. Phys. 2012, 14, 7562–7577. (39) 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. (40) Anzaki, S.; Watanabe, C.; Fukuzawa, K.; Mochizuki, Y.; Tanaka, S. Interaction energy analysis on specific binding of influenza virus hemagglutinin to avian and human sialosaccharide receptors: Importance of mutation-induced structural change. J. Mol. Graph. Model. 2014, 53, 48-58. (41) Öberg, H.; Brinck, T. Fragment molecular orbital study of the cAMP-dependent protein kinase catalyzed phosphoryl transfer: a comparison with the differential transition state stabilization method. Phys. Chem. Chem. Phys. 2016, 18, 15153-15161. (42) 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. (43) Otsuka, T.; Okimoto, N.; Taiji, M. Assessment and acceleration of binding energy calculations for protein-ligand complexes by the fragment molecular orbital method. J. Comput. Chem. 2016, 36, 2209-2218. (44) Prato, G.; Silvent, S.; Saka, S.; Lamberto, M.; Kosenkov, D. Thermodynamics of binding of di- and tetrasubstituted naphthalene diimide ligands to DNA G-quadruplex. J. Phys. Chem. B 2015, 119, 3335-3347. (45) 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.

35



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 36 of 51

(46) Nagata, T.; Fedorov, D. G.; Li, H.; Kitaura, K. Analytic gradient for second order MøllerPlesset perturbation theory with the polarizable continuum model based on the fragment molecular orbital method. J. Chem. Phys. 2012, 136, 204112. (47) Nakata, H.; Fedorov, D. G.; Kitaura, K.; Nakamura, S. Extension of the fragment molecular orbital method to treat large open-shell systems in solution . Chem. Phys. Lett. 2015, 635, 86 - 92. (48) 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. (49) Nagata, T.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. A combined effective fragment potential–fragment molecular orbital method. I. The energy expression and initial applications. J. Chem. Phys. 2009, 131, 024101. (50) Watanabe, H.; Okiyama, Y.; Nakano, T.; Tanaka, S. An interpretation of positional displacement of the helix12 in nuclear receptors: Preexistent swing-up motion triggered by ligand binding. Chem. Phys. Lett. 2010, 500, 116-119. (51) Yoshida, N. Efficient implementation of the three-dimensional reference interaction site model method in the fragment molecular orbital method. J. Chem. Phys. 2014, 140, 214118. (52) Fedorov, D. G.; Ishida, T.; Uebayasi, M.; Kitaura, K. The fragment molecular orbital method for geometry optimizations of polypeptides and proteins. J. Phys. Chem. A 2007,

111, 2722–2732. (53) Okamoto, T.; Ishikawa, T.; Koyano, Y.; Yamamoto, N.; Kuwata, K.; Nagaoka, M. Minimal implementation of the AMBER-PAICS interface for ab initio FMO-QM/MM-MD simulation. Bull. Chem. Soc. Japan 2013, 86, 210-222.

36


Page 37 of 51

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


(54) 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. (55) Nakata, H.; Fedorov, D. G.; Yokojima, S.; Kitaura, K.; Nakamura, S. Simulations of Raman spectra using the fragment molecular orbital method. J. Chem. Theory Comput. 2014,

10, 3689-3698. (56) Avramov, P. V.; Fedorov, D. G.; Sorokin, P. B.; Sakai, S.; Entani, S.; Ohtomo, M.; Matsumoto, Y.; Naramoto, H. Intrinsic edge asymmetry in narrow aigzag hexagonal heteroatomic nanoribbons causes their subtle uniform curvature. J. Phys. Chem. Lett. 2012, 3, 2003–2008. (57) Ishikawa, T.; Yamamoto, N.; Kuwata, K. Partial energy gradient based on the fragmet molecular orbital method: Application to geometry optimization. Chem. Phys. Lett. 2010,

500, 149-154. (58) Tsukamoto, T.; Mochizuki, Y.; Watanabe, N.; Fukuzawa, K.; Nakano, T. Partial geometry optimization with FMO-MP2 gradient: Application to TrpCage. Chem. Phys. Lett. 2012,

535, 157-162. (59) 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. (60) Steinmann, C.; Fedorov, D. G.; Jensen, J. H. Mapping enzymatic catalysis using the effective fragment molecular orbital method: Towards all ab initio biochemistry . PLOS One 2013, 8, e60602. (61) Christensen, A. S.; Steinmann, C.; Fedorov, D. G.; Jensen, J. H. Hybrid RHF/MP2 geometry optimizations with the effective fragment molecular orbital method. PLOS One 2014, 9, e88800.

37



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 38 of 51

(62) Nakata, H.; Fedorov, D. G.; Nagata, T.; Kitaura, K.; Nakamura, S. Simulations of chemical reactions with the frozen domain formulation of the fragment molecular orbital method. J. Chem. Theory Comput. 2015, 11, 3053-3064. (63) Fedorov, D. G.; Ishida, T.; Kitaura, K. Multilayer formulation of the fragment molecular orbital method (FMO). J. Phys. Chem. A 2005, 109, 2638-2646. (64) Sun, Q.; Phan, J.; Friberg, A. R.; Camper, D. V.; Olejniczak, E. T.; Fesik, S. W. A method for the second-site screening of K-Ras in the presence of a covalently attached firstsite ligand. J. Biomol. NMR 2014, 60, 11–14. (65) 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–480. (66) Nakata, H.; Fedorov, D. G.; Nagata, T.; Kitaura, K.; Nakamura, S. Simulations of chemical reactions with the frozen domain formulation of the fragment molecular orbital method. J. Chem. Theory Comput. 2015, 11, 3053-3064. (67) Nagata, T.; Brorsen, K.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. Fully analytic energy gradient in the fragment molecular orbital method. J. Chem. Phys. 2011, 134, 124115. (68) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. J.; 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. (69) Gordon, M. S.; Schmidt, M. W. Advances in electronic structure theory: GAMESS a decade later. In Theory and applications of computational chemistry, the first forty years; Dykstra, C. E.; Frenking, G.; Kim, K. S.; Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005. (70) Fedorov, D. G.; Kitaura, K. The importance of three-body terms in the fragment molecular orbital method. J. Chem. Phys. 2004, 120, 6832–6840. 38


Page 39 of 51

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


(71) 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. (72) Goerigk, L; Collyer, C. A.; Reimers, J. R. Recommending Hartree-Fock theory with London-dispersion and basis-set-superposition corrections for the optimization or quantum refinement of protein structures. J. Phys. Chem. B 2014, 118, 14612-14626. (73) Conrad, J. A.; Gordon, M. S. Modeling systems with π − π interactions using the HartreeFock method with an empirical dispersion correction. J. Phys. Chem. A 2015, 119, 53775385. (74) Ono, T.; Ohta, M.; Iseda, K.; Sada, K. Counter anion dependent swelling behaviour of poly(octadecyl acrylate)-based lipophilic polyelectrolyte gels as superabsorbent polymers for organic solvents. Soft Matter 2012, 8, 3700-3704. (75) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements HPu. J. Chem. Phys. 2010, 132, 154104. (76) Wang, J.; Cieplak, P.; Kollman, P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 2000, 21, 1049–1074. (77) Case, D. A.;

Betz, R. M.;

Botello-Smith, W.;

Cerutti, D. S.;

Cheatham, III, T. E.;

Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W. et al. AMBER 2015, University of California: San Francisco., 2015. (78) J-OCTA, version 2.0 http://www.jocta.com/index.html. . (79) Aoyagi, T.; Sawa, F.; Shoji, T.; Fukunaga, H.; Takimoto, J.; Doi, M. A general-purpose coarse-grained molecular dynamics program. Comp. Phys. Comm. 2002, 145, 267 - 279. 39



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 40 of 51

(80) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general Amber force field. J. Comput. Chem. 2004, 25, 1157–1174. (81) Steinmann, C.; Ibsen, M. W.; Hansen, A. S.; Jensen, J. H. FragIt: A tool to prepare input files for fragment based quantum chemical calculations. PLOS ONE 2012, 7, e44480. (82) Suenaga, M. Development of GUI for GAMESS/FMO calculation. J. Comput. Chem. Jpn. 2008, 7, 33–53. (83) Su, P. F.; Li, H. Energy decomposition analysis of covalent bonds and intermolecular interactions. J. Chem. Phys. 2009, 130, 074109. (84) 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, 22182231. (85) 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. (86) Fedorov, D. G.; Kitaura, K. Pair interaction energy decomposition analysis. J. Comput. Chem. 2007, 28, 222–237. (87) Alexeev, Y.; Mazanetz, M. P.; Ichihara, O.; Fedorov, D. G. GAMESS as a free quantummechanical platform for drug research. Curr. Top. Med. Chem. 2012, 12, 2013–2033. (88) Mazanetz, M. P.; Chudyk, E.; Fedorov, D. G.; Alexeev, Y. Computer aided drug discovery; Zhang, W. Ed.; Springer: New York, 2016.

40


Page 41 of 51

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


TOC Graphic

41



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

TOC


Page 42 of 51

Page 43 of 51

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


Figure1



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

Figure2


Page 44 of 51

Page 45 of 51

Figure3

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




Figure4

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 46 of 51

Page 47 of 51

Figure 5

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




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

Figure6


Page 48 of 51

Page 49 of 51

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


Figure7



Figure 8

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 50 of 51

Page 51 of 51

Figure 9

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



Efficient Geometry Optimization of Large Molecular Systems in

Recommend Documents