Molecular Electrostatic Potential and Electron Density of Large

Article Cite This: J. Phys. Chem. A 2019, 123, 6281−6290

pubs.acs.org/JPCA

Molecular Electrostatic Potential and Electron Density of Large Systems in Solution Computed with the Fragment Molecular Orbital Method Dmitri G. Fedorov,*,† Anton Brekhov,‡ Vladimir Mironov,‡ and Yuri Alexeev§

Downloaded via UNIV OF SOUTHERN INDIANA on July 25, 2019 at 08:35:54 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

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 ‡ Department of Chemistry, Lomonosov Moscow State University, Moscow, 119991, Russian Federation § Argonne Leadership Computing Facility and Computational Science Division, Argonne National Laboratory, Argonne, Illinois, 60439, United States S Supporting Information *

ABSTRACT: A solvent screening model for the molecular electrostatic potential is developed for the fragment molecular orbital method combined with the polarizable continuum model at the Hartree−Fock and density functional theory levels. The accuracy of the generated potentials is established in comparison to calculations without fragmentation. Solvent effects upon the molecular electrostatic potential and electron density of solute are discussed. The method is applied to two proteins: chignolin (PDB: 1UAO) and ovine prostaglandin H(2) synthase-1 (1EQG).

1. INTRODUCTION Quantum-mechanical (QM) methods can be used for computing the energy and a variety of other properties derived from the wave function, such as electron density. Driven by the progress in the development of linear-scaling methods1,2 and computational hardware, large-scale simulations calculations including solvent are now possible. A very efficient way of performing large system calculations3 is accomplished by fragmentation, and many fragment-based methods have been developed.4−19 The fragment molecular orbital (FMO) method20−24 is an efficiently parallelized25 fragmentation approach applied to study protein−ligand binding26,27 and inorganic systems.28,29 Interaction analyses in solution have been developed30−33 for FMO and applied to study binding in solution.34,35 Although fragment-based methods typically focus on the energy and its derivatives, a number of methods are explicitly formulated for electron density.10,36−41 Molecular electrostatic potential (MEP) calculations have been done for some fragment-based methods.38,42 Both electron density24,43,44 and MEP44−46 have been evaluated using FMO in the gas phase. There are methods to add electrostatic screening in solution, many of which are based on solving the Poisson−Boltzmann (PB) equation47,48 on a grid. In PB models, the protein and solvent are treated as uniform dielectric media, and the potential is obtained from the charge density. Several programs such as DelPhi,49−51 UHBDP,52 and APBS53 can be used to solve PB equations. There is also a variation of the polarizable © 2019 American Chemical Society

continuum model (PCM), where a part of the protein is treated as a dielectric continuum with a dielectric constant.54 Accurate electron density calculations are needed for a number of applications. For example, they are critical in X-ray crystallography,55−57 for quantum refinement of structures with QM methods,58−60 and in ion mobility spectroscopy.61 MEP is also commonly used in drug discovery and ligand docking.62 The purpose of this work is to add solvent effects to MEP calculations in FMO and to analyze solvent effects on electron density and MEP.

2. METHODOLOGY A concise description of FMO combined63,64 with PCM65 is as follows. The solute is divided into fragments (monomers), which are calculated in the presence of the electrostatic potential (embedding) from the whole system. Consequently, pairs of fragments (dimers) are calculated also in the embedding potential. The total properties are obtained in the many-body expansion.66−69 The solvent in PCM is represented by a continuum with a cavity, the surface of which is divided into small elements called tesserae, with a point charge placed in the center of each tessera. The cavity is defined as a union of atom-centered spheres with overlapping parts removed. The solute and solvent charges mutually polarize each other, and the electronic state of the solute and Received: May 24, 2019 Revised: June 28, 2019 Published: June 28, 2019 6281

DOI: 10.1021/acs.jpca.9b04936 J. Phys. Chem. A 2019, 123, 6281−6290

Article

The Journal of Physical Chemistry A ΔDIJ = DIJ − (DI ⊕ D J )

the values of solvent charges are determined self-consistently. Other solvent models have also been used with FMO.70−72 The MEP of solute monomers (X = I) or dimers (X = IJ) φsolute (r) at some point r can be computed as follows X φXsolute(r)

=

∑ A∈X

ρ X (r) =

ZA − |r − RA|

∫

ρ X (r′) dr ′ |r − r′|

Then, the dimer term in eq 8 can be written as (the nuclear contributions cancel out) ΔφIJsolute(r) = −

(1)

where ZA and RA are the nuclear charge and coordinates of atom A, respectively; r and r′ are electron coordinates and χμ are atomic orbitals (AO). The first term in eq 1 is due to protons (ZA protons in atom A). The electron density ρX(r) of fragment X in the AO basis is DXμν. The solvent screening of the solute electrostatic potential is due to induced solvent charges. In this work, the following solvent contribution is evaluated for MEP of fragment X qi φXsolvent(r) = ∑ ∑ |r − R i| (3) A∈X i∈A

W IJμν =

φ

(r) =

N

V=

=

∑ i=1

qi |r − R i|

(13)

(14)

N

∑ V I + ∑ (V IJ − V I − VJ ) I>J

(15)

X

where V is the contribution of fragment X to the total potential. The two-body terms in eq 15 describe the CT effect in the solute on the potential exerted on the solvent at tessera i ViIJ − ViI − ViJ = −Tr(ΔDIJ w i)

(16)

1 |ν ⟩ |r − R i|

(17)

i wμν = ⟨μ|

In PCM[1] and PCM⟨1⟩, the second sum in eq 15 is neglected, and the induced charges q1 (the superscript 1 signifies the one-body truncation of eq 15) are obtained as

∑ [φIJ (r) − φI (r) − φJ (r)]

N ji zy q1 = −C−1jjjj∑ V I zzzz j z k I=1 {

(7)

Using eqs 4, 7 can be rewritten as N

∑ ΔφIJsolute(r)

(18)

The difference between PCM[2] and PCM[1(2)] approaches is that in the former, the solvent charges (q in eq 14) are obtained fully self-consistently by iterating monomer and dimer calculations until convergence. It is accurate, but requires multiple calculations of the whole system. In contrast, the simpler PCM[1(2)] method requires only two calculations of the whole system. The difference between PCM[1] and PCM⟨1⟩ is in the definition of the solvent contribution to the energy,75 and both methods have the same induced solvent charges, electron density, and MEP. The accuracy of all FMO/ PCM methods is compared in Table 1.

(8)

where Δφsolute (r) = φsolute (r) − φsolute (r) − φsolute (r). There is IJ IJ I J no direct solvent contribution to MEP at the two-body level in eq 8. The electron density transfer between fragments I and J is Δρ IJ (r) ≡ ρ IJ (r) − ρ I (r) − ρ J (r) IJ ΔDμν χμ (r)χν (r)

μν ∈ IJ

∑

where Vi is the solute potential on tessera i. C is a square matrix related to the inverse intertessera distances, with the linear size of NTS. There are four main variations63 of PCM in FMO, denoted by PCM[1], PCM⟨1⟩, PCM[1(2)], and PCM[2]. The most accurate method is the latter, for which

I>J

∑

NTS

φIsolvent(r)

q = −C−1V

N

=

(12)

where the sum is over all tesserae (N ). The implementation follows eq 13 rather than eq 3, which is more efficient. Solvent charges qi are obtained by solving the PCM equations64

(6)

I>J

μν ∈ IJ

|ν ⟩ ,

TS

N

φ FMO2(r) = φ FMO1(r) +

|r − R i|

I=1

The potentials φX(r) are added in FMOn using the manybody expansion for n = 1, 2 for a molecular system divided into N fragments

φ FMO2(r) = φ FMO1(r) +

solvent

I=1

I=1

qi

N

(5)

∑ φI (r)

(11)

The total solvent contribution in eq 8 comes from monomers (see eqs 5 and 6),

The total potential due to both solute and solvent is

φ FMO1(r) =

NTS

∑ ⟨μ| i=1

where induced solvent charges qi of tessera i on atom A. qi are obtained from a continuum solvent model such as PCM or solvent model density.73,74 The charges qi in eq 3 are independent of X. Taking eq 3 for X = IJ and splitting the sum over atoms A in dimer IJ into monomers I and J, one obtains qi qi φIJsolvent(r) = ∑ ∑ = ∑∑ A ∈ IJ i ∈ A |r − R i| A ∈ I i ∈ A |r − R i| qi + ∑∑ = φIsolvent(r) + φJsolvent(r) A ∈ J i ∈ A |r − R i| (4) φX (r) = φXsolute + φXsolvent

′

Thus, the meaning of is very simple: it is the effect of CT between solute fragments I and J upon MEP. This CT also contributes to the solvent screening35 of the interaction between I and J as ΔECT·es = −1/2Tr(ΔDIJWIJ), IJ IJ where Wμν are the electrostatic integrals due to solvent charges

(2)

μν ∈ X

IJ

∫ Δ|rρ−(rr′|) dr′ Δφsolute (r) IJ

X Dμν χμ (r)χν (r)

∑

(10)

(9)

where ΔDXμν is the density transfer matrix 6282


Article

The Journal of Physical Chemistry A

Table 1. Accuracy of FMO/PCM Methods (a.u.), in Terms of Maximum (max) Deviations and rmsd to Full Unfragmented Calculations methoda

dimer calculations

density (max)

density (rmsd)

MEP (max)

MEP (rmsd)

FMO1/PCM[1] FMO2/PCM[1] FMO2/PCM[1(2)] FMO2/PCM[2]

none once (in the presence of q1 obtained for monomers) twice (first with q1, then with q) multipleb (first, with q1, after that, with q)

0.1053 0.0273 0.0246 0.0242

0.0000013 0.0000002 0.0000002 0.0000002

0.1122 0.0501 0.0405 0.0389

0.0000066 0.0000043 0.0000022 0.0000019

a

For properties in this table, PCM⟨1⟩ gives the same results as PCM[1]. bTypically, about 10 iterations are needed for convergence.

The electron density of the solute in FMO can be calculated

It has been suggested to use the Gaussian charge smearing in PCM.80 In this work, the smearing is implemented for computing the solvent contribution (eq 13). It is reasonable because the Dirac delta function, representing a point-like distribution, is the weak limit of a sequence of normalized Gaussians

as N

ρ(r) =

N

∑ ρ I (r) + ∑ Δρ IJ (r) I=1

I>J

(19)

ij α yz jj zz kπ {

An analysis of the many-body contributions to density in FMO can be done24 defining polarization and charge transfer contributions. The solvent effects on MEP are (a) the polarization of the electron density DXμν of the solute by the solvent (via the Fock matrix contribution63 due to the solvent), affecting φsolute (r) in eq 1 and (b) the explicit contribution of X the solvent charges to MEP, φsolvent(r) in eq 13. The latter is typically the dominant contribution. In the accuracy evaluation, the following very simple rootmean-square deviation (rmsd) criterion for a property B (B is electron density or MEP) was used for a set of Ngrid grid points 1 E B = grid N

1/2

iαy qiδ(r − R i) → qijjj zzz kπ {

(22)

3/2

2

e−α | r − R i|

(23)

The electrostatic potential of this Gaussian (Gauss) distribution is the Coulomb potential scaled by the error function (erf) factor81 qi VGauss(r) = erf( α |r − R i|) |r − R i| (24)

∑ [BFMO(ri) − Bfull (ri)]2 (20)

where “full” refers to the property obtained without fragmentation. The value of the Coulomb potential of a point charge nucleus (eq 1) is infinite at r = RA for any atom A. Normally, this is not a problem because one is not interested in MEP near nuclei when analyzing chemical properties. To obtain MEP near nuclei, protons in all atoms can be described quantummechanically76,77 like electrons, yielding a proton density ρ̅(r), which can be integrated like the electron density (the second term in eq 1) to obtain the proton contribution. One can approximate the proton density by a set of Gaussians centered on atoms, as in the finite nucleus model in relativistic theories.78 Likewise, the value of MEP is infinite at tessera centers, r = Ri. In PCM, the atomic radii are multiplied with a factor of α = 1.2 pushing the tesserae away from the nuclei. Very large MEP values may appear if the surface used for plotting MEP is very close to a tessera. If |r − RA| or |r − Ri| is smaller than a threshold (taken to be 0.001 au), then the contribution of that A or i is ignored and a warning is printed. In order to calculate MEP for grid points r outside the solvent cavity in PCM (i.e., for a point r inside the continuum solvent), one could use79 φX (r) = φXsolute(r)/ε

as α → ∞

The three-dimensional Gaussian distribution of charges for tessera i is

N grid i=1

2 w e−α(x − x0) → δ(x − x0)

The erf factor approaches 1 as the distance between a tessera center Ri and grid point r increases and the potential becomes nearly identical to the Coulomb potential of a point charge. On the other hand, in the vicinity of Ri, where the Coulomb potential has a singularity, the potential in eq 24 reduces to a value which depends only on α and qi qi lim VGauss = VCoulomb = |r − R i|→∞ |r − R i| (25) lim

|r − R i|→ 0

VGauss = 2qi

a π

(26)

If α is large, then the potential in eq 24 is very similar to the point-charge model, failing to smoothen out the discreet nature of charges. If α is small, then the charge smearing is too broad. Numerical tests suggested that it is reasonable to use 0.5 ≤ α ≤ 5. The asymptotic behavior of VGauss facilitates replacing the computationally expensive error function by a constant of 1 (compare eqs 24 and 25) for almost all grid points. Only in the vicinity of each tessera Ri one has to evaluate erf explicitly. By introducing the approximation error threshold ϵ, one can obtain the following

(21)

where ε is the dielectric constant of the solvent. The potential in eq 21 may have a discontinuity at the cavity surface if the values in eqs 5 and 21 differ. In this work, MEP is calculated according to eq 5 for all grid points (that are typically inside the cavity anyway). Also, for points near the surface, the problem of defining MEP is exacerbated by the charge escape.65

VCoulomb − VGauss =ϵ VCoulomb

(27)

1 − erf( α |r − R i|) = erfc( α |r − R i|) = ϵ

(28)

|r − R i| = 6283

erf c−1ϵ = rthr α

(29) DOI: 10.1021/acs.jpca.9b04936 J. Phys. Chem. A 2019, 123, 6281−6290

Article


Figure 1. MEPs drawn on isosurfaces of the electron density with the cutoff T (a.u.): the upper row in the gas phase and the lower row in solution (FMO2/PCM[1]). The potential range is from −0.2 to 0.2 a.u.

where erfc−1(x) is the inverse complementary error function. The relative accuracy cutoff of ϵ = 10−8 was chosen. In the MEP calculation, if |r − Ri| > rthr (see eq 29), then the point charge potential contribution due to tessera i is computed; otherwise, the Gaussian model in eq 24 is used. The physical picture of the MEP calculation in FMO/PCM is as follows. The polarization of fragments in FMO is reflected in the fragment densities. The electron density is affected by CT as shown in eq 9. The solute contribution to MEP is affected by CT as shown in eq 11. The solvent contribution to MEP is affected by CT via the solute potential contribution in eq 16 that affects the solvent charges (eq 14) and thus MEP according to eq 13.

values of FMO thresholds were used. The 6-31(++)G** basis set was used for chignolin calculations (diffuse functions are added only on negatively charged carboxylate groups). An orthogonal grid surrounding the molecule with an extra space of 2 Å and spacing between points of 0.2 Å was used. C-PCM calculations were done with 60 tesserae per atom and van-derWaals atomic radii. Proteins were fragmented as one residue per fragment. Water was used as the solvent, although the method is general and can be applied to any solvent.

4. RESULTS AND DISCUSSION 4.1. Effect of the Density Surface. The isosurface cutoff T defining the density surface was varied to see how much the potential is affected by the cutoff and how far the tesserae are from the density surface. For these tests, the level of FMO2/ PCM[1] was chosen. There is one cationic (Gly-1) and three anionic (Asp-3, Glu-5 and Gly-10) fragments in chignolin (see Figure 1). It is clear that density cutoff does have a strong effect on the value of the potential. The general trend is that for a larger cutoff the surface is closer to the solute atoms, and thus solute has a stronger effect. This is especially pronounced for MEP in solution, where the solvent screening makes the total potential small (a nearly white surface) for T = 0.001, but for T = 0.025 there is a visible red and blue coloring, corresponding to the negative and positive potentials, respectively. The density surface for the cutoff of T = 0.001 nearly coincides with the solvent surface as shown in Figure 2. 4.2. Accuracy of the Density and Electrostatic Potential in FMO. The accuracy results for chignolin are given in Table 1. For electron density, the level of PCM makes a little difference, and the main change comes from adding interfragment charge transfer in FMO2, reducing the maximum deviation 4 times and the rmsd 6 times. The maximum deviations and rmsd for MEP in FMO2 are reduced 1.5−3 times compared to FMO1. Adding the effect of charge transfer to induced charges (PCM[1] → PCM[1(2)] and PCM[2]) decreases rmsd in MEP by a factor of 2; the approximate PCM[1(2)] approach shows a minor difference to the self-consistent PCM[2].

3. COMPUTATIONAL DETAILS MEP calculations in this work were implemented in GAMESS82 and parallelized using the generalized distributed data interface.83 In each run, two cube files were generated, one with density and another with MEP, for the purpose of plotting MEP on a surface defined by density. Density and MEP were stored as separate shared arrays distributed across all nodes. For a medium size system, chignolin (PDB: 1UAO) was computed with the experimental NMR structure. This is a small protein that consists of 10 residues, of which 4 are charged. Three kinds of calculations were done on chignolin: (a) verifying the effect of threshold T used to generate density surfaces, (b) comparing MEP of FMO-RHF/PCM at various levels in Table 1 to full RHF/PCM without fragmentation, and (c) comparing MEP in gas phase and solution, using FMO. A comparison of electron density in gas phase and solution was also performed. MEP visualization was done with Pymol. Secondly, MEP was computed for a larger system, ovine prostaglandin H(2) synthase-1 (1EQG). The structure was taken from the previous study in the complex with ibuprofen,84 and for the MEP calculation the ligand was removed. Calculations were mostly performed at the Hartree−Fock level except for the 1EQG computations performed with density functional theory (DFT). In DFT, SG-1 grid was used with the CAM-B3LYP functional. Hybrid orbital projection operators were employed for fragment boundaries; the default 6284


Article


Figure 2. Separation of density (blue) and solvent tesserae (shown as gray balls) for chignolin. The density surfaces are plotted for three cutoff values.

Figure 3. ϕFMO(r) − ϕfull(r) plotted on the ρfull(r) surface with the cutoff T = 0.025 au for chignolin at different levels of theory.

The conclusions are that for electron density any PCM level can be used, whereas for the MEP, PCM[1(2)] is preferred. It should be noted, however, that the accuracy estimates are done for the whole cube file and thus do not reflect the effect of choosing a particular surface. They also include the areas near nuclei and tesserae, where MEP may have substantial deviations. However, for practical purposes, these areas are not important, as the potential is usually plotted sufficiently far from the poles. For a better understanding of the accuracy of FMO versus full unfragmented calculations, the values of ϕFMO(r) − ϕfull(r) are plotted on the ρfull(r) surface with the cutoff T equal to 0.025 and 0.005 au shown in Figures 3 and 4, respectively. The worst accuracy is for the lowest level of FMO1/PCM[1] as expected, followed by FMO2/PCM[1]. The charge transfer between fragments is important for the grid points closer to atoms. If the surface is sufficiently separated from atoms as for T = 0.005, then neglecting charge transfer between fragments in FMO1 is a good approximation. For all cutoffs, FMO2/PCM[1(2)] has a comparable accuracy to FMO2/PCM[2], showing insignificance of the full solvent charge relaxation at the two-body level for describing MEP. To get a better understanding of FMO errors, ρFMO2/PCM[2](r) − ρFMO1/PCM[1](r) is plotted in Figure 5a for the cutoff of 0.005. In FMO, proteins are fragmented between Cα and C′ atoms, where C′ is the carboxyl carbon (in C′O). Cα is the bond detached atom (from which the C−C bond is detached), which is most severely affected by fragmentation. Indeed, one can see in Figure 5a that the largest two-body effects are for Cα atoms. The bond fragmentation errors are reduced in FMO2 by considering fragment dimers, in which the bond between Cα and C′ atoms is intact. This correction is

missing in FMO1 leading to the mentioned deviations of electron density. Another source of error is lone pairs on O (in CO and OH groups). Lone pairs are prone to be affected by interfragment charge transfer accounted for in FMO2. However, the density error of FMO1 is mostly localized near atoms. For an electron density isosurface used to plot MEP, errors of this kind are insignificant. It is illustrated in Figure 5b, where ρFMO2/PCM[2](r) − ρFMO1/PCM[1](r) is plotted on the ρfull(r) surface (the cutoff of 0.025). The density errors are noticeable around the fragmentation points, but they are relatively small (on the order of 0.005 au or 0.5% of one electron, as seen by the value range in the legend). FMO1 delivers reasonable values of density and MEP when plotted on surfaces not too near the nuclei. Next, the error behavior of the FMO2 density is analyzed in Figure 5c. The values of ρfull(r) − ρFMO2/PCM[2](r) are plotted with the cutoff of 0.001. This plot reflects the total impact of nbody charge transfer and exchange-repulsion effects24 for n > 2. They are most significant at fragmentation points (the backbone) and near lone pairs of oxygen atoms, because of charge transfer. Moreover, the peptide-bond between Asp-3 and Pro-4 is also affected, mainly due to the lone pairs on nitrogen. Likewise, there is a substantial many-body effect on the lone pair of nitrogen for the peptide bond between Tyr-2 and Asp-3. These effects may be enhanced by the negative charge on Asp-3. However, on the absolute scale, the manybody effects are small, on the order of 0.1% of one electron. 4.3. Solvent Effects on the Electrostatic Potential and Electron Density. The electronic density of the solute is affected by the solvent because of the polarization effects. The 6285


Article


Figure 4. ϕFMO(r) − ϕfull(r) plotted on the ρfull(r) surface with the cutoff T = 0.005 au for chignolin at different levels of theory.

electrostatic potential coming from the solvent contributes to the Fock matrix of solute fragments, thereby polarizing their densities. There is no charge transfer between solute and solvent in PCM. The results are shown in Table 2. It can be seen that rmsd of density change is nearly zero, which is related to the solute charge conservation; the effect of the polarization shows up in the slightly larger value of the maximum deviation. It can be expected that the density near charged residues is strongly affected by the solvent. The solvent (water) makes a large contribution to MEP, as seen in Figure 1 and Table 2. At each density threshold level, solvent makes a large change to MEP. The cationic fragment Gly-1 (N-terminus) has a positive (denoted by blue color) potential around it, which is screened by the solvent unless the surface is very tight, as in the case of T = 0.025 in the gas phase. For T = 0.005, the negatively charged fragments Asp-3, Glu-5, and Gly-10 (C-terminus) have a negative (red) potential around them. The charged residues Asp-3 and Glu5 induce substantial charges on the solvent cavity surface of nearby fragments. This is the reason why not only anionic fragments but also nearby neutral fragments, such as Thr-6, have a positive potential around them. It is interesting to observe that for the high density threshold level of 0.025 the potential has a more pronounced atomistic structure. This reflects the increased importance of individual atomic charges, rather than the overall fragment charge (as seen for T = 0.005 and 0.001) for the areas closer to nuclei and further from tesserae. The solvent effects for this tight surface are weaker as the surface is further away from the solvent. One can clearly see individual functional groups. For example, Gly-

Figure 5. Chignolin results: (a) ρFMO2/PCM[2](r) − ρFMO1/PCM[1](r) plotted with the cutoff of 0.005 (positive and negative values are plotted in blue and red, respectively), (b) ρFMO2/PCM[2](r) − ρFMO1/PCM[1](r) plotted on the ρfull(r) surface (the cutoff of 0.025), the plotted range is from −0.005 to 0.005 a.u., and (c) ρfull(r) − ρFMO2/PCM[2](r) plotted with the cutoff of 0.001.

Table 2. Solvents Effects (a.u.) Upon the Density and MEP of Chignolin at the Level of PCM[1(2)]: Maximum (max) Deviations and rmsd between the Gas Phase and Solution density MEP

max

rmsd

0.057810 0.596261

0.000001 0.000155

10 has an anionic carboxyl group, whose negative (red) potential is seen as red spheres (negative oxygen atoms, the Mulliken charges of −0.84 each) connected to the blue area for carbon (the charge of +0.76). The N-terminus (Gly-1) has a positive potential because of −NH3+, whereas Thr residues have −OH groups that have a negatively charged oxygen. To gain more insight, the values of ρFMO2/PCM[2](r) − FMO2 ρ (r) (the solvated density minus the gas phase density) are plotted in Figure 6a with the cutoff of 0.005. This plot 6286


Article


a small set of values) for correlation studies in drug discovery.85 4.4. Effect of Using Gaussians for Solvent Potentials. Point charge representation of the solvent electrostatics creates a few noticeable patches on the surface with the smallest cutoff (0.001) in solution because the surface gets close to the PCM solvent cavity (Figure 1). Higher density levels correspond to isosurfaces which are closer to the atomic nuclei. In this case, the point charge model for solvent potential works well. The effect of replacing point charges with Gaussian distributions using a uniform exponent α is shown in Figure 7 for the cutoff of 0.001. A zoomed-in view of the surface was

Figure 6. Solvents effects on the electron density in chignolin: (a) ρFMO2/PCM[2](r) − ρFMO2(r) plotted with the cutoff of 0.005 (positive and negative values are plotted in blue and red, respectively) and (b) values of ρFMO2/PCM[2](r) − ρFMO2(r) plotted on the ρfull(r) surface with the cutoff of 0.025.

shows the electron density polarization due to solvent. The magnitude of the density change due to the solvent polarization is on the order of 0.5% of one electron. Aromatic hydrophobic sidechains of Tyr-2 and Trp-9 are not affected by solvent. The largest polarization is for the COO− group of Asp3. All other charged residues (Gly-1 and Gly-10) are affected, as well as neutral residues adjacent to them (Tyr-2 and Trp-9). However, all these changes are on the backbone and local to the nuclei. They have a small effect at the area where density or MEP is usually plotted. The values of ρFMO2/PCM[2](r) − ρFMO2(r) are plotted on the ρfull(r) surface in Figure 6b with the cutoff of 0.025. However, the magnitude of these values is only 0.1% of one electron, and all atoms in this small protein are affected. As a general trend, oxygen atoms on the protein surface tend to gain electron density due to solvent polarization, while hydrogen atoms lose some density. What is the “right” value of T? As a ligand approaches the protein during binding, at first the ligand feels the potential of the protein at a large separation corresponding to small values of T, whereby the combined potential of atoms blurs into residue regions. However, when the ligand locks in and docks into the binding pocket, it sees individual atoms as for larger values of T. Clearly, there is a desolvation accompanying the binding, which changes the potential as well. Thus, a range of T values are important to protein−ligand binding, when considered as a dynamic process with a varying protein−ligand distance. Choosing a single value of T is not enough to describe the full complexity of the protein−ligand binding. However, it may be possible to find an optimum value of T (or

Figure 7. Comparison of MEP of chignolin obtained with point charge and Gaussian models (with the uniform exponent of α), at the level of FMO2/PCM[1] and the density cutoff of 0.001.

chosen where the patches from point charges are most noticeable. As the exponent increases, the potential becomes more smeared and the patches disappear. Based on these results, the value of α = 1 may be a good compromise between a too strong smearing of the charges and generating a clean surface if one is interested in a density surface (cutoff of 0.001) near the PCM surface. 4.5. MEP for Proteins in Solution. The MEP of prostaglandin H(2) synthase-1 (1EQG) was computed at the level of FMO1-CAM-B3LYP/6-31G*/PCM[1]. In this system, there are 17 734 atoms divided into 1092 fragments. In PCM, there are 83 096 tesserae. The results are shown in Figure 8. The calculation with the grid spacing of 0.5 Å took 16.6 h on a PC cluster (300 CPU cores total). The calculation cost is dominated by the second term in eq 1. 6287



■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Dmitri G. Fedorov: 0000-0003-2530-5989 Vladimir Mironov: 0000-0002-9454-5823 Notes

The authors declare no competing financial interest.

■

Figure 8. Structure and MEP of prostaglandin H(2) synthase-1 (1EQG) computed at the level of FMO1-CAM-B3LYP/6-31G*/ PCM[1] with the cutoff of 0.005 and α = 1. The potential range is from −0.2 to 0.2 a.u.

ACKNOWLEDGMENTS D.G.F. thanks Dr. Garib Murshudov (LMB) for many fruitful discussions about MEP. This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Contract No. DE-AC02-06CH11357. This research used the resources of the Argonne Leadership Computing Facility, which is a U.S. Department of Energy (DOE) Office of Science User Facility supported under Contract DE-AC02-06CH11357. A.B. and V.M. thank the Russian Foundation for Basic Research for funding (project 19-03-00043). This work was supported by JSPS KAKENHI grant 19H02682.

5. CONCLUSIONS In this work, the solvent contribution to the MEP has been derived and implemented for the FMO method. The manybody contributions due to charge transfer between solute fragments have been derived. The accuracy of MEP has been established in comparison to the results without fragmentation. It has been shown that the interfragment charge transfer is described in FMO through both direct solute terms (eq 11) and indirect terms via the induced solvent charges (eq 16). Several levels of FMO combined with PCM can be used, providing options to speed up the calculations at a small penalty of accuracy. The dependence of MEP on the isodensity surface controlled by the cutoff has been shown. The closer the surface is located to the solute (the larger the cutoff is), the more atomistic the electrostatic potential looks and the less pronounced the solvent effects are. The role of the cutoff on the potential representation was discussed in the context of considering a binding process dynamically when a ligand approaches a protein. The solvent contribution to MEP was implemented for the point charges and Gaussian distributions. The latter representation provides a better MEP for the small cutoffs defining the density surface. The same Gaussian model can be used to remove the divergence of the potential near nuclei. The developed mechanism for evaluating MEP in solution is expected to have a large impact on the applications of the FMO in drug discovery studies,86,87 in particular for those using the charge complementarity.85 The earlier applications44−46 of FMO to calculate the MEP did not include the solvent contribution and thus resulted in values that are too large, especially near charged residues. The potential can be used to analyze infrared and Raman spectra in solution.88 The formulation developed here for FMO can be used as a basis for evaluating MEP in other fragment-based methods.89,90

■

■

REFERENCES

(1) Linear-Scaling Techniques in Computational Chemistry and Physics; Zalesny, R., Papadopoulos, M. G., Mezey, P. G., Leszczynski, J., Eds.; Springer: Berlin, 2011. (2) Computational Methods for Large Systems: Electronic Structure Approaches for Biotechnology and Nanotechnology; Reimers, J. R., Ed.; Wiley: New York, 2011. (3) 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. (4) 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. (5) Otto, P.; Ladik, J. Investigation of the interaction between molecules at medium distances. Chem. Phys. 1975, 8, 192−200. (6) Gao, J. Toward a Molecular Orbital Derived Empirical Potential for Liquid Simulations. J. Phys. Chem. B 1997, 101, 657−663. (7) Söderhjelm, P.; Ryde, U. How Accurate Can a Force Field Become? A Polarizable Multipole Model Combined with FragmentWise Quantum-Mechanical Calculations. J. Phys. Chem. A 2009, 113, 617−627. (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) Kobayashi, M.; Nakai, H. How Does it Become Possible to Treat Delocalized and/or Open-Shell Systems in FragmentationBased Linear-Scaling Electronic Structure Calculations? The Case of the Divide-and-Conquer Method. Phys. Chem. Chem. Phys. 2012, 14, 7629−7639. (11) Truhlar, D. G.; Yu, H.; Leverentz, H. R.; Bai, P.; Siepmann, J. I. 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, 666−670. (12) Gao, J.; Truhlar, D. G.; Wang, Y.; Mazack, M. J. M.; Lö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. (13) Jose, K. V. J.; Raghavachari, K. Evaluation of Energy Gradients and Infrared Vibrational Spectra through Molecules-in-Molecules

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.9b04936. Timing comparison of MEP/PCM for various levels of FMO and full unfragmented calculations; a comparison of the point charge and Gaussian-based solvent models (PDF) 6288


Article

The Journal of Physical Chemistry A Fragment-Based Approach. J. Chem. Theory Comput. 2015, 11, 950− 961. (14) 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. (15) 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. (16) 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. (17) 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. (18) Fang, T.; Li, Y.; Li, S. Generalized Energy-Based Fragmentation Approach for Modeling Condensed Phase Systems. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2017, 7, No. e1297. (19) Kobayashi, R.; Amos, R.; Collins, M. A. Microsolvation within the Systematic Molecular Fragmentation by Annihilation Approach. J. Phys. Chem. A 2017, 121, 334−341. (20) 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. (21) 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. (22) Fedorov, D. G.; Nagata, T.; Kitaura, K. Exploring Chemistry with the Fragment Molecular Orbital Method. Phys. Chem. Chem. Phys. 2012, 14, 7562−7577. (23) 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. (24) Fedorov, D. G. The Fragment Molecular Orbital Method: Theoretical Development, Implementation in GAMESS and Applications. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2017, 7, No. e1322. (25) Mironov, V.; Alexeev, Y.; Fedorov, D. G. Multithreaded Parallelization of the Energy and Analytic Gradient in the Fragment Molecular Orbital Method. Int. J. Quantum Chem. 2019, 119, No. e25937. (26) 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. (27) Heifetz, A.; James, T.; Southey, M.; Morao, I.; Aldeghi, M.; Sarrat, L.; Fedorov, D. G.; Bodkin, M. J.; Townsend-Nicholson, A. Characterising GPCR-ligand interactions using a fragment molecular orbital-based approach. Curr. Opin. Struct. Biol. 2019, 55, 85−92. (28) 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. (29) Kato, K.; Fukuzawa, K.; Mochizuki, Y. Modeling of Hydroxyapatite-Peptide Interaction Based on Fragment Molecular Orbital Method. Chem. Phys. Lett. 2015, 629, 58−64. (30) Nagata, T.; Fedorov, D. G.; Sawada, T.; Kitaura, K. Analysis of Solute-Solvent Interactions in the Fragment Molecular Orbital Method Interfaced with Effective Fragment Potentials: Theory and Application to a Solvated Griffithsin-Carbohydrate Complex. J. Phys. Chem. A 2012, 116, 9088−9099. (31) 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. (32) Okiyama, Y.; Nakano, T.; Watanabe, C.; Fukuzawa, K.; Mochizuki, Y.; Tanaka, S. Fragment Molecular Orbital Calculations

with Implicit Solvent Based on the Poisson-Boltzmann Equation: Implementation and DNA Study. J. Phys. Chem. B 2018, 122, 4457− 4471. (33) Okiyama, Y.; Watanabe, C.; Fukuzawa, K.; Mochizuki, Y.; Nakano, T.; Tanaka, S. Fragment Molecular Orbital Calculations with Implicit Solvent Based on the Poisson-Boltzmann Equation: II. Protein and Its Ligand-Binding System Studies. J. Phys. Chem. B 2019, 123, 957−973. (34) Fedorov, D. G.; Kitaura, K. Subsystem Analysis for the Fragment Molecular Orbital Method and its Application to ProteinLigand Binding in Solution. J. Phys. Chem. A 2016, 120, 2218−2231. (35) Fedorov, D. G.; Kitaura, K. Pair Interaction Energy Decomposition Analysis for Density Functional Theory and Density-Functional Tight-Binding with an Evaluation of Energy Fluctuations in Molecular Dynamics. J. Phys. Chem. A 2018, 122, 1781−1795. (36) Stoll, H.; Preuß, H. On The Direct Calculation of Localized HF Orbitals in Molecule Clusters, Layers and Solids. Theor. Chim. Acta 1977, 46, 11−21. (37) Yang, W. Direct Calculation of Electron Density in DensityFunctional Theory. Phys. Rev. Lett. 1991, 66, 1438−1441. (38) Babu, K.; Ganesh, V.; Gadre, S. R.; Ghermani, N. E. Tailoring Approach for Exploring Electron Densities and Electrostatic Potentials of Molecular Crystals. Theor. Chem. Acc. 2004, 111, 255−263. (39) Szekeres, Z.; Exner, T.; Mezey, P. G. Fuzzy Fragment Selection Strategies, Basis Set Dependence and HF-DFT Comparisons in the Applications of the ADMA Method of Macromolecular Quantum Chemistry. Int. J. Quantum Chem. 2005, 104, 847−860. (40) He, X.; Zhang, J. Z. H. A New Method for Direct Calculation of Total Energy of Protein. J. Chem. Phys. 2005, 122, 031103. (41) Cembran, A.; Bao, P.; Wang, Y.; Song, L.; Truhlar, D. G.; Gao, J. On the Interfragment Exchange in the X-Pol Method. J. Chem. Theory Comput. 2010, 6, 2469−2476. (42) Le, H.-A.; Lee, A. M.; Bettens, R. P. A. Accurately Reproducing Ab Initio Electrostatic Potentials with Multipoles and Fragmentation. J. Phys. Chem. A 2009, 113, 10527−10533. (43) Inadomi, Y.; Nakano, T.; Kitaura, K.; Nagashima, U. Definition of Molecular Orbitals in Fragment Molecular Orbital Method. Chem. Phys. Lett. 2002, 364, 139−143. (44) Ishikawa, T. Ab Initio Quantum Chemical Calculation of Electron Density, Electrostatic Potential, and Electric Field of Biomolecule Based on Fragment Molecular Orbital Method. Int. J. Quantum Chem. 2018, 118, No. e25535. (45) 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. (46) Mazanetz, M. P.; Chudyk, E.; Fedorov, D. G.; Alexeev, Y. Applications of the Fragment Molecular Orbital Method to Drug Research. In Computer Aided Drug Discovery; Zhang, W., Ed.; Springer: New York, 2016; pp 217−255. (47) Gilson, M. K.; Rashin, A.; Fine, R.; Honig, B. On the Calculation of Electrostatic Interactions in Proteins. J. Mol. Biol. 1985, 184, 503−516. (48) Gilson, M.; Sharp, K. A.; Honig, B. Calculating the electrostatic potential of molecules in solution: Method and error assessment. J. Comput. Chem. 1988, 9, 327−335. (49) Gilson, M. K.; Honig, B. Calculation of the Total Electrostatic Energy of a Macromolecular System. Solvation Energies, Binding Energies and Conformational Analysis. Proteins: Struct., Funct., Bioinf. 1988, 4, 7−18. (50) Nicholls, A.; Honig, B. A Rapid Finite Algorithm, Utilizing Successive Over-relaxation to Solve the Poisson−Boltzmann Equation. J. Comput. Chem. 1991, 12, 435−445. (51) Rocchia, W.; Alexov, E.; Honig, B. Extending the Applicability of the Nonlinear Poisson−Boltzmann Equation: Multiple Dielectric Constants and Multivalent Ions. J. Phys. Chem. B 2001, 105, 6507− 6514. 6289


Article


(73) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (74) Fedorov, D. G. Analysis of Solute-Solvent Interactions Using the Solvation Model Density Combined with the Fragment Molecular Orbital Method. Chem. Phys. Lett. 2018, 702, 111−116. (75) 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. (76) Ishimoto, T.; Tachikawa, M.; Nagashima, U. A Fragment Molecular-Orbital-Multicomponent Molecular-Orbital Method for Analyzing H/D Isotope Effects in Large Molecules. J. Chem. Phys. 2006, 124, 014112. (77) Auer, B.; Pak, M. V.; Hammes-Schiffer, S. Nuclear-Electronic Orbital Method within the Fragment Molecular Orbital Approach. J. Phys. Chem. C 2010, 114, 5582−5588. (78) Malkin, E.; Repiský, M.; Komorovský, S.; Mach, P.; Malkina, O. L.; Malkin, V. G. Effects of Finite Size Nuclei in Relativistic FourComponent Calculations of Hyperfine Structure. J. Chem. Phys. 2011, 134, 044111. (79) Lange, A. W.; Herbert, J. M. A Simple Polarizable Continuum Solvation Model for Electrolyte Solutions. J. Chem. Phys. 2011, 134, 204110. (80) Scalmani, G.; Frisch, M. J. Continuous Surface Charge Polarizable Continuum Models of Solvation. I. General Formalism. J. Chem. Phys. 2010, 132, 114110. (81) Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular Integral Evaluation. Molecular Electronic-Structure Theory; John Wiley & Sons, Ltd: Chichester, 2014; pp 336−432. (82) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., Jr. General Atomic and Molecular Electronic Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (83) 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. (84) 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. (85) Bauer, M. R.; Mackey, M. D. Electrostatic Complementarity as a Fast and Effective Tool to Optimize Binding and Selectivity of Protein-Ligand Complexes. J. Med. Chem. 2019, 62, 3036−3050. (86) Choi, J.; Kim, H.; Jin, X.; Lim, H.; Kim, S.; Roh, I.; Kang, H.; Tai No, K.; Sohn, H. Application of the Fragment Molecular Orbital Method to Discover Novel Natural Products for Prion Disease. Sci. Rep. 2018, 8, 13063. (87) Heifetz, A.; James, T.; Southey, M.; Morao, I.; Aldeghi, M.; Sarrat, L.; Fedorov, D. G.; Bodkin, M. J.; Townsend-Nicholson, A. Characterising GPCR-ligand interactions using a fragment molecular orbital-based approach. Curr. Opin. Struct. Biol. 2019, 55, 85−92. (88) Nakata, H.; Fedorov, D. G. Simulations of Infrared and Raman Spectra in Solution Using the Fragment Molecular Orbital Method. Phys. Chem. Chem. Phys. 2019, 21, 13641−13652. (89) Mei, Y.; Ji, C.; Zhang, J. Z. H. A new Quantum Method for Electrostatic Solvation Energy of Protein. J. Chem. Phys. 2006, 125, 094906. (90) Söderhjelm, P.; Kongsted, J.; Ryde, U. Ligand Affinities Estimated by Quantum Chemical Calculations. J. Chem. Theory Comput. 2010, 6, 1726−1737.

(52) Madura, J. D.; Briggs, J. M.; Wade, R. C.; Davis, M. E.; Luty, B. A.; Ilin, A.; Antosiewicz, J.; Gilson, M. K.; Bagheri, B.; Scott, L. R.; McCammon, J. A. Electrostatics and Diffusion of Molecules in Solution: Simulations with the University of Houston Brownian Dynamics Program. Comput. Phys. Commun. 1995, 91, 57−95. (53) Holst, M.; Baker, N.; Wang, F. Adaptive Multilevel Finite Element Solution of the Poisson−Boltzmann Equation I. Algorithms and Examples. J. Comput. Chem. 2000, 21, 1319−1342. (54) Si, D.; Li, H. Heterogeneous Conductorlike Solvation Model. J. Chem. Phys. 2009, 131, 044123. (55) Murshudov, G. N.; Skubák, P.; Lebedev, A. A.; Pannu, N. S.; Steiner, R. A.; Nicholls, R. A.; Winn, M. D.; Long, F.; Vagin, A. A. REFMAC5 for the Refinement of Macromolecular Crystal Structures. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2011, 67, 355−367. (56) Northey, T.; Zotev, N.; Kirrander, A. Ab Initio Calculation of Molecular Diffraction. J. Chem. Theory Comput. 2014, 10, 4911−4920. (57) Northey, T.; Kirrander, A. Ab Initio Fragment Method for Calculating Molecular X-ray Diffraction. J. Phys. Chem. A 2019, 123, 3395−3406. (58) Cao, L.; Caldararu, O.; Rosenzweig, A. C.; Ryde, U. Quantum Refinement does not Support Dinuclear Copper Sites in Crystal Structures of Particulate Methane Monooxygenase. Angew. Chem., Int. Ed. 2018, 57, 162−166. (59) Merz, K. M. Using Quantum Mechanical Approaches to Study Biological Systems. Acc. Chem. Res. 2014, 47, 2804−2811. (60) Zheng, M.; Reimers, J. R.; Waller, M. P.; Afonine, P. V. Q|R: Quantum-Based Refinement. Acta Crystallogr., Sect. D: Struct. Biol. 2017, 73, 45−52. (61) Alexeev, Y.; Fedorov, D. G.; Shvartsburg, A. A. Effective Ion Mobility Calculations for Macromolecules by Scattering on Electron Clouds. J. Phys. Chem. A 2014, 118, 6763−6772. (62) Pagadala, N. S.; Syed, K.; Tuszynski, J. Software for Molecular Docking: a Review. Biophys. Rev. 2017, 9, 91−102. (63) 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. (64) Li, H.; Fedorov, D. G.; Nagata, T.; Kitaura, K.; Jensen, J. H.; Gordon, M. S. Energy Gradients in Combined Fragment Molecular Orbital and Polarizable Continuum Model (FMO/PCM) Calculation. J. Comput. Chem. 2010, 31, 778−790. (65) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3094. (66) Fedorov, D. G.; Kitaura, K. The Importance of Three-Body Terms in the Fragment Molecular Orbital Method. J. Chem. Phys. 2004, 120, 6832−6840. (67) Watanabe, C.; Fukuzawa, K.; Okiyama, Y.; Tsukamoto, T.; Kato, A.; Tanaka, S.; Mochizuki, Y.; Nakano, T. Three- and FourBody Corrected Fragment Molecular Orbital Calculations with a Novel Subdividing Fragmentation Method Applicable to StructureBased Drug Design. J. Mol. Graphics Modell. 2013, 41, 31−42. (68) Fedorov, D. G.; Asada, N.; Nakanishi, I.; Kitaura, K. The Use of Many-Body Expansions and Geometry Optimizations in FragmentBased Methods. Acc. Chem. Res. 2014, 47, 2846−2856. (69) Richard, R. M.; Lao, K. U.; Herbert, J. M. Aiming for Benchmark Accuracy with the Many-Body Expansion. Acc. Chem. Res. 2014, 47, 2828−2836. (70) 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. (71) Watanabe, H.; Okiyama, Y.; Nakano, T.; Tanaka, S. Incorporation of Solvation Effects Into the Fragment Molecular Orbital Calculations with the Poisson-Boltzmann Equation. Chem. Phys. Lett. 2010, 500, 116−119. (72) 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. 6290


Molecular Electrostatic Potential and Electron Density of Large

Recommend Documents