Benchmarking Electron Densities and Electrostatic Potentials of

Article pubs.acs.org/JCTC

Benchmarking Electron Densities and Electrostatic Potentials of Proteins from the Three-Partition Frozen Density Embedding Method Albrecht Goez and Johannes Neugebauer* Theoretische Organische Chemie, Organisch-Chemisches Institut and Center for Multiscale Theory and Computation, Westfälische Wilhelms-Universität Münster, Corrensstraße 40, 48149 Münster, Germany S Supporting Information *

ABSTRACT: The fragment-based Three-Partition Frozen Density Embedding (3-FDE) approach [Jacob, C. R.; Visscher, L. J. Chem. Phys. 2008, 128, 155102] is used to generate protein densities and electrostatic potentials, which are critically assessed in comparison to supermolecular Kohn−Sham Density Functional Theory (DFT) results obtained with sophisticated exchange−correlation functionals. The influence of several parameters and user choices is explored with respect to accuracy and reliability. In addition, a recently implemented combination of the 3-FDE scheme with hybrid functionals is applied in production calculations for the first time. We demonstrate that the 3-FDE method not only closely reproduces results from corresponding supermolecular calculations for routine situations (peptides/proteins in solution) but can even surpass conventional Kohn−Sham DFT in accuracy for difficult cases, such as zwitterionic structures in vacuo. This is due to the fact that the fragmentation inherently limits the overdelocalization caused by the self-interaction error in common DFT approximations. The method is thus not only able to reduce the computational effort for the description of large biological entities but also can strongly reduce the artifacts brought about by the SIE. yields unphysical results for extended systems.8 This error can be reduced by using hybrid functionals, which include a certain amount of exact exchange from Hartree−Fock (HF) theory. Of course, full elimination of the SIE is only guaranteed in case of 100% HF exchange, and the best-suited functional in practice clearly depends on the problem at hand.9 Especially sophisticated range-separated hybrids have gained much attention in this context.10,11 However, the evaluation of exact exchange again increases the computational scaling and is thus limited to smaller molecular systems. One possibility to alleviate the SIE without resorting to large hybrid calculations is to employ a fragment-based method.12−15 Well-known approaches include the Molecular Tailoring Approach,16,17 the Fragment Molecular Orbital Method,18,19 and the Effective Fragment Potential scheme20 as well as multilevel embedding approaches such as the QM/MMpol method21 or different flavors of the ONIOM scheme.22 A comprehensive overview can be found in two recent reviews.2,23 With all these methods, the total problem is broken down into smaller parts, which can be solved independently. The determining factor regarding the accuracy then becomes the approximate description of the interactions between the

1. INTRODUCTION The theoretical modeling of biomacromolecules such as proteins or DNA structures has become an indispensable part of modern biochemistry and biophysics. Traditionally, mainly force field approaches have been employed to describe the energetics and properties of large biomolecular structures.1 However, in addition to the errors introduced by the inherently parameter-dependent nature of these methods, they are fundamentally incapable of describing certain phenomena, such as electronic transitions. Due to the ever-increasing computational power becoming available in the last decades, electronic-structure theories such as wave function methods or Density Functional Theory (DFT) have become more and more applicable to protein modeling (see, e.g., the reviews in refs 2 and 3). The unfavorable scaling of ab initio methods unfortunately limits their application to rather small systems, even though local wave function approaches (for instance the DL-PNO coupled-cluster method4,5) or massively parallel GPU implementations6,7 have made impressive progress in the past years. Furthermore, although Kohn−Sham (KS) DFT approaches using pure density functionals combined with density fitting or the Resolution of the Identity approximation formally scale only as 6(N3), the well-known self-interaction error (SIE) © XXXX American Chemical Society

Received: June 9, 2016

A

DOI: 10.1021/acs.jctc.6b00590 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

simulating a continuum around the system of interest, the overdelocalization reported in refs 11 and 15 can be tremendously reduced, leading to a good description both with supermolecular DFT and the 3-FDE scheme. However, in contrast to regular DFT methods, 3-FDE is shown to also yield excellent results if applied in vacuo, which can be attributed to the fragmentation procedure. In both cases, the efficiency gain due to the linear-scaling nature of the underlying FDE approach allows the application to proteins of much larger size than accessible with regular DFT. The method is therefore proven to be an all-round tool for the investigation of very different biological problems. The remainder of this article is organized as follows: In Section 2, a short summary of the 3-FDE method is given, followed by computational details in Section 3. Section 4 addresses the determination of a supermolecular reference method in order to judge the quality of the fragment-based calculations. A comprehensive assessment of the densities and ESPs generated with the 3-FDE method in solution is presented in Section 5. Furthermore, some tests in vacuo are presented in Section 6. The article is concluded with a summary and outlook in Section 7.

subsystems. Fragment methods have the advantage that the overdelocalization of charge, which is commonly caused by the SIE, is automatically restricted, because the electrons are constrained to their respective fragments. For instance, Olsen et al. could show in a recent article15 that MP2 reference data for the electrostatic potential (ESP) of insulin can be reproduced quite accurately using a fragment approach based on small hybrid DFT calculations, whereas supermolecular B3LYP calculations perform much worse. Another possibility to alleviate the SIE effects would be to employ constrained DFT24 to assign a fixed number of electrons to a certain spatial area. One particularly useful fragment method is the Frozen Density Embedding (FDE) formalism,25 which is based on Subsystem DFT26,27 (for recent reviews, see refs 28 and 29). In FDE, the interaction between the subsystems is described by an effective embedding potential, which is constructed from the electron densities of the individual fragments. In principle, this can be done in an exact fashion; however, the contribution of the kinetic energy of the electrons presents a problem, which in practice necessitates the approximation of an additional functional covering the nonadditive part of the kinetic energy. Unfortunately, an exact, density-dependent expression for this term is not available, which effectively makes Subsystem DFT an approximation to supermolecular KS-DFT. For structures consisting of a single covalent unit (such as most proteins), an additional problem arises. This is due to the fact that the currently available approximations to the nonadditive kinetic energy functional only work well in case of noncovalently bound fragments. However, proteins have to be fragmented in order to gain any advantage over a regular calculation. A pragmatic solution to this problem is the ThreePartition FDE approach (3-FDE) by Jacob and co-workers.30 Here, the problem of covalent cuts is avoided by employing capping groups in the spirit of the Molecular Fractionation with Conjugate Caps (MFCC) approach.31 The 3-FDE method has been tested for several polypeptides and proteins,30,32−34 but the results have only been compared to supermolecular KS-DFT calculations with a GGA exchange− correlation (XC) functional so far. Because of the SIE, this is a questionable reference. In the present article, we take a closer look at several aspects regarding the 3-FDE approach and thoroughly assess the quality of the generated protein densities and electrostatic potentials in comparison with suitable reference methods. In addition, the recently implemented combination of 3-FDE with hybrid XC functionals35 is applied. It must be noted that an efficient as well as consistent treatment of both additive and nonadditive XC contributions is impossible in this context, because such an approach would require supersystem orbitals. A common strategy is thus to use a hybrid functional for the individual fragment contributions, while evaluating the nonadditive part with a purely density-dependent GGA form.36 The implications of this approach for the 3-FDE method are also investigated in the present work. An important consideration concerns the representation of the environment around the protein. Because of the presence of charged side chains, a description in vacuo is far from the physical reality of biological systems. Furthermore, it has already been shown that the SIE can have particularly grave consequences for calculations on charged systems in vacuo, which can be strongly reduced or even eliminated by employing a stabilizing continuum.10 We first demonstrate that by

2. THEORETICAL BACKGROUND The most important features of the 3-FDE method will be briefly discussed in the following. For more details, the reader is referred to ref 30. The first step toward density-based embedding calculations is the fragmentation of the total system into a number of subsystems. As mentioned in the Introduction, cutting through covalent bonds normally leads to situations where the FDE approach fails due to the currently available approximations for the nonadditive kinetic energy term.37,38 Therefore, some kind of capping scheme is mandatory. In the simplest case, a hydrogen atom could be added as a link atom, as is often done in QM/MM approaches39,40 (for a comprehensive overview, see ref 41). In order to preserve the nature of the backbone as much as possible, the 3-FDE approach however uses a partitioning according to the MFCC scheme.31 Although originally intended for the calculation of interaction energies, the method was quickly extended to allow for the determination of total densities as well.42 In the MFCC approach, the protein is split into subsystems by cutting through amide bonds at specific places in the backbone. The resulting fragments are subsequently capped with suitable terminal groups, usually methylated amino acid termini. This leads to a number of (partly overlapping) so-called capped f ragments, which do not feature any dangling bonds. The electron densities of the latter are determined in isolatedmolecule calculations. In addition, pairs of neighboring capping groups are conjoined to form cap molecules. The electronic structure of these small entities is determined in isolatedmolecule calculations as well. A first approximation to the supermolecular density ρtot can be obtained by simply summing over the individual densities of the capped fragments, ρcapped , A and subtracting those of the cap molecules, ρcap B NF tot

ρ ( r ⃗) =

∑ A

F

NC

ρAcapped ( r ⃗)

−

∑ ρBcap ( r ⃗) B

(1)

C

where N and N are the number of capped fragments and cap molecules, respectively. The resulting density will be referred to as the MFCC density in the following. B


Article

Journal of Chemical Theory and Computation λ( r ⃗) = 0 for r ⃗ ∉ W Kcap

The MFCC density is a rather crude approximation, which especially suffers from errors close to the severed bonds.30,32,33 Furthermore, the subtraction procedure can lead to unphysical areas of negative electron density. The central idea of the 3FDE approach is to improve upon the MFCC density by applying the FDE scheme. The functional principle of the latter is to relax the density of a single (“active”) fragment in an embedding potential generated by all other (“frozen”) fragments. This is accomplished by setting up Kohn−Sham equations with constrained electron density (KSCED)25 for the active system K

is introduced. The densities in the cap region (but not outside of it) are thus guaranteed to cancel exactly in the subtraction procedure upon solution of eq 3. The KSCED can then be rewritten to read ⎡ ∇2 ⎤ + v 3 ‐ FDE[ρK , ρ tot ]( r ⃗)⎥ψK , i( r ⃗) = ϵK , iψK , i( r ⃗) ⎢− ⎣ 2 ⎦

v 3 ‐ FDE[ρK , ρ tot ]( r ⃗) eff emb ⎧ [ρK , ρ tot ]( r ⃗) for r ⃗ ∉ W Kcap ⎪ v [ρ ]( r ⃗) + v K =⎨ ⎪ cap for r ⃗ ∈ W Kcap ⎩ v ( r ⃗)

(2)

Here, ψK,i and ϵK,i represent Kohn−Sham-like orbitals and their corresponding energies, respectively, for fragment K. veff K has the same form as the effective potential in a regular KS-DFT calculation and only depends on the active fragment density ρK. The additional embedding potential is represented by vemb K and contains Coulomb interactions as well as kinetic and XC contributions stemming from all other fragments. By solving the KSCED for the active system, an improved density for the latter is obtained, whereas the inactive density remains unchanged. However, by interchanging the role of active and inactive subsystem(s), the total density can be optimized in an iterative fashion. This procedure has been called “Freeze-andThaw” (FT) in the literature.43 In addition to the regular XC functional, which is required for any DFT calculation, two central quantities have to be approximated in practical FDE calculations. These are the nonadditive contributions (as opposed to the fragment contributions) to the kinetic and XC energy. The nonadditive XC term is usually described by the same XC functional as the fragment contributions. An exception are hybrid functionals, where the lack of supersystem orbitals prevents such an approach. In this case, either a purely density-dependent form must be chosen, or a potential reconstruction scheme44,45 can be employed. The nonadditive term pertaining to the kinetic energy is routinely approximated through a Thomas−Fermilike form or a gradient-corrected expression such as the PW91k functional,46,47 which has been shown to work well for noncovalently bound subsystems47−51 as well as for properly capped MFCC fragments.30,32,33 In contrast to the MFCC approach, the 3-FDE scheme elegantly circumvents the problem of negative densities by applying a special constraint, which forces the density in the cap regions of a capped fragment and its corresponding cap molecules to be equal. This is achieved by adding an additional constraint to the energy minimization condition30

(6)

Inside the cap region Wcap K , the total potential is replaced by a cap potential vcap, which is responsible for enforcing the equality constraint, whereas the regular sum of effective and embedding potential is applied outside this region. In the current implementation,32 the cap region consists of all points of the numerical integration grid which are within 3 a0 of a cap atom and not closer to any noncap atom. Although the size of the cap radius clearly has an effect on the results, the latter was found to be minor as long as a sufficiently large value is chosen.30 We have verified that this finding holds also for the present study (see the Supporting Information for full details). As a first approximation for the cap potential, the one obtained from the corresponding cap molecule calculation is used. Upon self-consistent solution of the KSCED for the capped fragment, this choice will clearly not lead to exactly the same density as in the cap molecule. Therefore, an iterative update procedure is employed. As soon as the SCF convergence measure drops below a certain threshold, the cap potential is updated by means of the potential reconstruction scheme proposed by van Leeuwen and Baerends52 to steer toward the desired cap density. The error in the cap density is measured by calculating the integrated root-mean-square deviation (RMSD) of the electron density in that region. When the latter falls below a user-defined cap density convergence threshold, the cap density is fixed, and the SCF cycle proceeds as normal until total convergence is achieved. In a system with three or more fragments, there are two distinct strategies to carry out FT cycles. After an updated density for the first subsystem has been determined, the embedding potential acting on the next fragment can either directly use the updated density of the first one or employ only the initial densities until a full FT cycle has been completed. The first approach has been called “sequential FT”,30 since each calculation depends on the previous one and the individual fragments cannot be updated in parallel. Accordingly, the other strategy is referred to as “parallel FT”, because different processes can take care of the individual fragments and communication is only required at the end of each cycle. It has been argued that sequential FT cycles might exhibit faster overall convergence,34 which is confirmed in this work (see Section 5).

(∫ ρ (r ⃗)dr ⃗ − N )

+

el K

K

⎤

∫ λ( r ⃗)[ρK ( r ⃗) − ρcap( r ⃗)]d r ⃗⎦⎥ =! 0

(5)

where the 3-FDE potential v3‑FDE is defined as

⎡ ∇2 ⎤ + vKeff [ρK ]( r ⃗) + vKemb[ρK , ρtot ]( r ⃗)⎥ψK , i( r ⃗) = ϵK , iψK , i( r ⃗) ⎢− ⎣ 2 ⎦

δ ⎡ E[ρ , ρ tot , ρcap ] + μ δρK ⎢⎣ K

(4)

(3)

Here, μ is the regular Lagrangian multiplier which ensures normalization of the density to the NelK electrons of the respective fragment. The additional multiplier λ(r)⃗ constrains the fragment density to equal the cap molecule density ρcap. In order to confine this constraint to a certain cap region Wcap K , a position dependence of λ in the form C


Article


3. COMPUTATIONAL DETAILS 3.1. Quantum-Chemical Calculations. Calculations with Slater-Type Orbital (STO) basis sets have been performed with a locally modified version of the ADF code.53,54 In all cases, a double-ζ basis set with a single set of polarization functions (DZP) was employed. In order to make sure that our conclusions are independent of the basis set size, one representative example was however also investigated with a triple-ζ basis set with two sets of polarization functions (TZ2P). These results are presented in the Supporting Information. The 3-FDE code in ADF,30 which is based on the regular FDE implementation,55 was used for all fragmentbased calculations. A Voronoi quadrature grid was employed for numerical integration.56,57 The nonadditive kinetic energy functional for FDE was always chosen to be PW91k.46,47 In case of GGA calculations, the nonadditive XC contributions were evaluated with the same functional as the fragment terms. As mentioned above, this consistent approach is not possible if hybrid expressions are employed. A purely density-dependent form was used in these cases (usually a GGA functional) and is explicitly stated in the text where applicable. The PYADF framework58 was used to set up all ADF calculations. PYADF uses OPENBABEL59,60 routines to create the capped fragments. Cuts always occur through a peptide bond between two amino acid residues, and the capping groups were chosen to be either C(O)CH3 or NHCH3 groups for Nterminal and C-terminal residues, respectively. All cap atoms which have a direct equivalent in the full structure (i. e., all backbone atoms as well as those hydrogen atoms connected to the α carbon atom of the adjacent residue) are assigned identical positions as in the original structure. Cap hydrogen atoms replacing either an amino acid side chain or the remainder of the backbone are placed along the bonds to the original connecting atom but at a default distance of 1.07 Å. Even though the cap density cancels out exactly upon full convergence of eq 5, the nature of the capping groups might have an influence on the obtained results due to the iterative numerical solution process. However, in a series of calculations with different varieties of these groups, only mild influences were found. The corresponding tests are discussed in the Supporting Information. When constructing the MFCC density, which takes the role of the initial density for 3-FDE calculations, it is recommended to employ a stabilizing environment already for the fragment calculations, for instance in the form of individual COSMO cavities.33 In contrast, it was recently established that for reasons of robustness the electronic structure of the cap molecules should better be determined in vacuo instead.35 Following these recommendations, all calculations for proteins in solution were set up in this way. An exception is the example presented in Section 6 (insulin in vacuo); here, we prepared all fragments and caps without a solvent environment, which forced us to employ rather small fragments to avoid convergence problems. Some reference calculations were carried out with GaussType Orbital (GTO) basis sets. In this case, HF and DFT runs were performed with ORCA 3.0.3,61 whereas MP2 calculations were facilitated by the RICC2 module62−64 of TURBOMOLE 7.0.65−67 The GTO basis set was always chosen to be ccpVDZ.68,69

Calculations featuring an implicit solvent environment were made possible through the COSMO algorithm.70 The permittivity was always chosen to correspond to water, and the solvent cavity was created with the default settings of the respective program. To guarantee that the slightly different settings between codes do not introduce any artifacts, a corresponding test was carried out, which is presented in the Supporting Information. XC potentials employed in the scope of this work included the BLYP,71,72 BP86,71,73 PBE,74,75 B3LYP,72,76−78 BHandHLYP,71,72,79 PBE0,74,75,80 M06-2X,81 and SAOP82−84 formulations. A long-range corrected version85 of the BP86 functional is denoted by LR-BP86. In addition, the range-separated CAMB3LYP functional86 was used for some tests. In its standard formulation, the error function is used to distinguish between the short- and long-range regimes. Since this function does not offer advantages with STO basis sets, an analogous functional featuring the Yukawa potential87 has been implemented in ADF under the name CAMY-B3LYP.85,88 Apart from this switching function, the two functionals are identical, since all parameters have been assigned the same values. 3.2. Test Systems. 3.2.1. Linear Polyglycines. In analogy to ref 11, a set of linear polyglycines in zwitterionic form (charged termini) was created. A completely linear configuration (ϕ = ψ = ω = 180°) was enforced, and the number of glycine units was varied between 2 and 20. The PROTEIN module of TINKER 6.3.289 was used to obtain initial geometries, which were further optimized with MP2/cc-pVDZ in a COSMO cavity. All structures were constrained to the point group Cs. An example is shown in Figure 1(a). 3.2.2. Miniproteins. In addition to polyglycines, three different miniproteins were investigated. In all cases, a solution-state NMR structure was obtained from the Protein Data Bank, which was subsequently optimized with BP86/DZP and COSMO. The heavy-atom positions were kept fixed throughout the optimization.

Figure 1. Graphical representation of the peptide/protein test set employed in this work. Protein secondary structure elements are shown in a transparent cartoon representation, and hydrogens have been omitted for clarity. The figure was created with VMD99 (secondary structure assignment based on STRIDE100). D


Article


Table 1. Normalized, Integrated RMSDs in the Electron Density (Unitless) for Different Test Systems in Solution: (a) Differences between MP2 and Various Methods (First Column) Obtained with the Same Basis Set (cc-pVDZ) and (b) Differences between DZP and cc-pVDZ Basis Sets Obtained with the Same Method (First Column) (a) Method/cc-pVDZ − MP2/cc-pVDZ

a

method

(Gly)2

(Gly)5

(Gly)10

BLYP B3LYP M06-2X CAM-B3LYP HF

0.011 0.008 0.005 0.008 0.016

0.011 0.009 0.005 0.008 0.016

method

(Gly)2

(Gly)5

(Gly)10

BLYP B3LYP M06-2X CAM-B3LYPa HF

0.027 0.027 0.025 0.026 0.025

0.027 0.026 0.025 0.026 0.024

0.026 0.026 0.025 0.026 0.024

(Gly)20

TC5b

IGF-F1-1

RCB

0.012 0.010 0.005 0.009 0.016

0.012 0.009 0.005 0.008 0.015

0.012 0.009 0.005 0.009 0.015

(Gly)20

TC5b

IGF-F1-1

RCB

0.026 0.026 0.024 0.025 0.024

0.025 0.024 0.023 0.024 0.022

0.024 0.023 0.022 0.023 0.021

0.024 0.023 0.022 0.023 0.021

0.011 0.012 0.009 0.009 0.005 0.005 0.008 0.008 0.017 0.017 (b) Method/DZP − Method/cc-pVDZ

Comparison between CAM-B3LYP/cc-pVDZ and CAMY-B3LYP/DZP.

TC5b. TC5b is an artificial protein, which was first synthesized by Neidigh et al. (PDB ID:1L2Y).90 It features the Trp-cage motif and is one of the smallest proteins stable in solution (20 amino acid residues, 304 atoms). In spite of its small size, it already contains a complete α-helix. The structure used in this work carries a total charge of +1 (five charged side chains) and is shown in Figure 1(b). IGF-F1-1. IGF-F1-1 is a protein fragment which is thought to inhibit the interaction between insulin-like growth factor 1 (IGF-1) and IGF-binding proteins. Its structure was determined by Deshayes et al. (PDB ID:1LB7)91 and also contains a prominent α-helix. The system comprises 16 amino acid residues (255 atoms) and possesses a total charge of +2 (five charged side chains). The structure employed in this article is shown in Figure 1(c). RCB-1. RCB-1 contains 19 amino acid residues (283 atoms) and acts as a biomarker in Ricinus communis seeds. In ref 92, a synthesis strategy was presented along with a determination of its structure (PDB ID:2MTM). An overall charge of +2 (four charged side chains) was obtained for the structure employed in the present study, which is graphically represented in Figure 1(d). 3.2.3. Insulin. The monomeric form of insulin was already the subject of similar studies,11,15 although it should be pointed out that these were investigations in vacuo (instead of within a solvent environment). In order to be consistent with these previous studies, the monomer (chains A and B) was extracted from the corresponding crystal structure (PDB ID:1MSO),93 which already contained hydrogen positions. Where applicable, alternative location A was selected. It should be noted that one histidine residue was protonated at both potential sites, which was not altered for the sake of consistency with refs 11 and 15. The resulting structure contains 787 atoms and bears an overall charge of −1, resulting from 11 charged residues. Due to its size, no further geometry optimization was carried out. The structure is shown in Figure 1(e). 3.3. Observed Properties. The electron density was directly used as a measure of accuracy in the form of an RMSD over all points of an atom-centered grid, which was further normalized to the number of electrons Nel in the respective system

ΔRMS(ρ) =

1 N el

∫ (ρref ( r ⃗) − ρapprox ( r ⃗))2 d r ⃗

(7)

The grid was chosen to be the Voronoi numerical integration grid employed by ADF. Furthermore, the ESP was determined on a density isosurface according to ρ = 10−4 e·a−3 o , which roughly corresponds to the distance of van der Waals contacts.11 The set of points was obtained by calculating the electron density on an equidistant grid with a spacing of 0.2 a0 and subsequently selecting all points where 8 × 10−5 e·a−3 0 ≤ ρ ≤ 1 × 10−4 e·a−3 . This assignment was based on the MP2/cc0 pVDZ density for all systems but insulin, where the CAMYB3LYP/DZP density was used instead. Differences in the ESP are given in the form of Mean Absolute Deviations (MADs), averaged over the resulting set of Npnt points MAD

Δ

1 (ϕ) = pnt N

N pnt

∑ i=1

(ϕref ( ri ⃗) − ϕapprox ( ri ⃗))2

(8)

4. DETERMINATION OF A SUITABLE REFERENCE The choice of reference densities and ESPs is crucial when assessing the quality of FDE results. In the past, mainly comparisons to supermolecular DFT calculations with GGA functionals have been presented.30,32,33 However, this might not be the best choice, since the latter have been shown to be affected by the SIE, which might be at least partially cured by fragment-based methods. MP2 has been suggested as an alternative,11 since the underlying HF wave function is inherently SIE-free, and the deficiencies of the HF approach are alleviated by the post-SCF procedure. However, a direct comparison of 3-FDE densities and ESPs to MP2 results is not straightforward. This is due to the fact that 3-FDE is currently only implemented in ADF, which employs STOs, whereas common MP2 codes only work with GTOs. Thus, any comparison would include a deviation due to the basis set difference. The magnitude of this error can be estimated from Table 1. Here, integrated RMSDs of the densities obtained for some test systems in solution are presented. While the upper part of the table lists differences due to the employed method (both calculations using the same basis set), the lower part gives the differences caused by the basis set inequality (both calculations using the same method). E


Article


2 have been employed for the following calculations. In addition, the initial tests were confined to the TC5b protein.

It is evident that the errors caused by the basis set are often more than twice as large, which completely obscures comparisons between MP2 with a GTO basis set and HF/ DFT calculations with an STO set. Similar results were obtained for comparisons between other STO and GTO sets (not shown). Interestingly, the basis set effect is almost equal for all tested methods and systems. The behavior is similar for the ESP (see Table S1 in the Supporting Information), although a slightly larger spread in the magnitude of the basis set effect is observed. The ESP errors due to method and basis set difference are visualized in Figure 2 for two molecular systems. It is clearly

Table 2. Default Settings for 3-FDE Calculations Presented in This Work parameter

default setting

XC functional fragment size FT mode cap convergence threshold nonadditive XC functional

PBE n=2 sequential 5 × 10−4 PBE

Tests concerning the other miniproteins can be found in the Supporting Information. In all cases, the MFCC density/ESP as well as the results from several FT cycles will be assessed. 5.1. XC Functional. The electron density and ESP for TC5b in solution were determined in 3-FDE calculations with the default settings and three different functionals (PBE, B3LYP, CAMY-B3LYP). These can be regarded as representative of the most important classes of density functionals (GGA, hybrid, range-separated hybrid). The resulting deviations from the reference (supermolecular M06-2X/CAMY-B3LYP calculations) are shown in Figure 3. Clearly convergent behavior is observed for all three functionals. After two FT cycles, no changes can be spotted in the deviations. For the density, the order of the functionals with respect to accuracy is preserved through all cycles. In contrast, PBE actually produces a better MFCC ESP than the more sophisticated formulations but is again surpassed by the latter during the FT cycles. For both quantities, CAMY-B3LYP performs best, even though the difference to the next best candidate is small. It should be noted that all tested functionals approach the supermolecular result obtained with the same formulation upon performing FT cycles (see also Section 5.2 and Figure S1 in the Supporting Information). However, a remaining deviation is always present, which can be attributed to the nonadditive kinetic energy approximation. 5.2. Fragment Size. Four different fragment sizes (n = 1,2,5,10) were tested. In the present case, larger fragments are irrelevant, since TC5b is comprised of only 20 amino acid residues. The results are shown in Figure 4. As expected, the deviations can be systematically decreased by increasing the fragment size. Still, a difference to the supermolecular results always remains. An interesting point concerns the converged results. In case of pure density functionals (such as the PBE functional in Figure 4), the density errors become independent of the fragmentation size. This is not true for hybrid functionals (see Figure S2 in the Supporting Information), which is due to the nonadditive XC component. As mentioned above, an inconsistent formulation has to be chosen for the latter due to the orbital dependence of hybrid functionals (in this case the PBE functional). Better approximations might be found by potential reconstruction approaches, but the latter cannot yet be considered suitable for routine calculations.37 The ESP also converges upon increasing the fragment size and/or number of FT cycles, but a larger difference remains even in case of GGA functionals. 5.3. FT Mode. It has been argued that the sequential FT scheme should lead to faster convergence of the total density compared to the parallel approach.34 However, this has not yet been systematically investigated. Figure 5 confirms this assumption for different fragmentation sizes, both for the

Figure 2. Difference ESP for (Gly)5 and TC5b from supermolecular calculations in solution with different methods and basis sets, represented as color-coded density isosurfaces (ρ = 10−4 e·a−3 0 ). Red corresponds to areas of overly negative charge, whereas blue signifies an excess of positive charge. Left column: BLYP/cc-pVDZ − MP2/ccpVDZ (method error). Right column: BLYP/DZP − BLYP/cc-pVDZ (basis set error). Color scale: ± 0.005 au. The figure was created with VMD.99

visible that the basis set error not only is much larger but also affects different parts of the molecule, rendering quantitative comparisons very difficult. These results call for a reference method within ADF (which excludes MP2), in order to guarantee that the same basis set can be chosen. The smallest density deviations from MP2 were found for M06-2X (see Table 1), while the ESP is best described by CAM-B3LYP (see Table S1 in the Supporting Information). Therefore, the reference in all tests presented in the following sections is constituted by supermolecular calculations with M06-2X/DZP (density) or CAMY-B3LYP/ DZP (ESP).

5. RESULTS IN SOLUTION Several important user choices have to be made when preparing 3-FDE calculations. The following five parameters were chosen for closer evaluation: XC functional, fragment size n (number of amino acid residues per fragment), FT mode (parallel or sequential), cap density convergence threshold, and nonadditive XC functional (in case of hybrids). Since the parameter space is rather large, not all possible combinations will be discussed. Instead, one specific setting will be varied at a time starting from a standard configuration. Unless explicitly specified otherwise, the defaults listed in Table F


Article


Figure 3. Influence of the XC functional on density and ESP errors in 3-FDE calculations for TC5b in solution over the course of five FT cycles (cycle 0 = MFCC). Errors are given with respect to supermolecular results obtained with M06-2X (density) or CAMY-B3LYP (ESP).

fragments at somewhat reduced cost compared to the default. The effect of the precise value on the obtained densities and ESPs was found to be very weak (see Figure S3 in the Supporting Information). 5.5. Nonadditive XC Potential. So far, the 3-FDE approach has not been systematically tested with hybrid functionals for the description of the fragment XC contributions, although some studies used regular FDE with hybrid functionals in the framework of so-called disconnectedenvironment models.94−96 To evaluate the nonadditive XC contributions, either a purely density-dependent formulation or a potential reconstruction approach44,45 must be chosen, since no supermolecular orbitals are available. In this work, only standard LDA and GGA formulations have been tested. The results are presented for different combinations of fragment and nonadditive XC potentials in Figure 6. No clear trends can be deduced from the results, but the effect of this choice is in general very weak. Either all three nonadditive formulations (LDA, BLYP, PBE) lead to the same result (B3LYP density/CAMY-B3LYP ESP), or LDA performs slightly different from the GGA formulations (CAMY-B3LYP density/B3LYP ESP). Furthermore, the dependence com-

density and the ESP. In each case, the errors drop slightly faster when the sequential FT scheme is employed. However, the differences are small, and after two FT cycles essentially the same results are obtained. It can be inferred that the choice of FT mode can be made mostly based on the available computer infrastructure, because it does not affect the final results significantly. This is especially true if comparatively large fragments are used or several FT cycles are carried out. 5.4. Cap Density Convergence Threshold. The cap density convergence threshold is used to test for equal electron density between the cap regions of the active capped fragment and those of the corresponding cap molecules. It is evaluated in the form of an integrated RMSD over all grid points belonging to the cap region. The current default in ADF is 1 × 10−4, which proved to be already rather strict. Choosing a tighter threshold (5 × 10−5 and 1 × 10−5 were also tested) led to extremely slow SCF convergence (a limit of 300 cycles was not sufficient for most fragments). On the other hand, convergence problems from the second FT cycle on were encountered if a weaker threshold of 1 × 10−3 was applied. This can be traced back to remaining cap density located inside the active fragment. A threshold of 5 × 10−4 was found to reliably converge all G


Article


Figure 4. Influence of the fragment size on density and ESP errors in 3-FDE calculations for TC5b in solution over the course of five FT cycles (cycle 0 = MFCC). Errors are given with respect to supermolecular results obtained with M06-2X (density) or CAMY-B3LYP (ESP).

charged side chains and/or terminal groups were present in the test systems, no stabilizing environment of any kind was employed. It was already shown by Rudberg8 as well as Antony and Grimme10 that such an approach leads to vanishing HOMO−LUMO gaps and difficult SCF convergence. Furthermore, it should be pointed out that a description in vacuo cannot be considered a realistic physiological environment for proteins. As long as equilibrium properties in a natural environment are of interest, it is therefore mandatory to employ at least a qualitative environmental model such as the COSMO method. Table 1 clearly establishes that the SIE is hardly a problem in this case. This can be inferred from the fact that even if pure density functionals are used and the polyglycine chain length is varied between two and 20 units, the errors are completely preserved. This is in stark contrast to the results without solvent from ref 11, where large errors had been obtained already for five glycine units. It is thus confirmed (as reported in refs 8 and 10) that a solvent model − such as COSMO − counteracts the consequences of the SIE and stabilizes the HOMO−LUMO gap.

pletely vanishes when larger fragments are used (not shown). This allows to routinely employ hybrid functionals for the fragment XC contributions in 3-FDE calculations. However, it should be mentioned that the use of hybrids leads to drastically enhanced requirements in terms of computer time. The construction of the MFCC density and two subsequent FT cycles with the PBE functional took only 2.7 h on a 12-core Intel Xeon machine running at 2.5 GHz, while 5.1 h were required with B3LYP and even 10.9 h with the range-separated CAMY-B3LYP. 5.6. Other Proteins. In order to avoid drawing conclusions based on special properties of TC5b, the IGF-F1-1 and RCB-1 proteins were investigated as well. Only the XC functional and the fragment size were varied in these cases, since the effect of the other variables is expected to be insignificant. Qualitatively similar results were obtained, which can be found in Section 5 of the Supporting Information.

6. RESULTS IN VACUO The performance of different density functionals with respect to the ESP for small peptides and proteins has been investigated by Jakobsen et al. in ref 11. However, even though H


Article


Figure 5. Influence of the FT mode on density and ESP errors in 3-FDE calculations for TC5b in solution over the course of five FT cycles (cycle 0 = MFCC). Errors are given with respect to supermolecular results obtained with M06-2X (density) or CAMY-B3LYP (ESP).

individual fragments and two FT cycles are carried out (also in vacuo), the magnitude of the error can be decreased to almost the same value (0.0019 au). Subsequent cycles again have little effect on the deviations. A graphical representation of the ESP error can be found in Figure 7. The results prove that the overdelocalization, which is caused by the SIE already for very small peptides in vacuo, is efficiently counteracted by the fragmentation. The density deviations exhibit similar trends, although the difference between supermolecular and 3-FDE calculations is smaller than for the ESP (not shown). 6.2. Insulin. A much more challenging test case is presented by the insulin monomer. This protein was also investigated in refs 11 and 15. The density and ESP errors are reported in Figure 8 for a small variety of functionals and fragmentation sizes. All other settings were chosen according to the defaults listed in Table 2. The fragmentation reduces the SIE to such a large extent that similar density errors as for the other systems in solution are obtained, in spite of the challenging electronic situation. Smoothly convergent behavior is again observed upon carrying out FT cycles, with no significant changes after the second one. Astonishingly, almost identical deviations are obtained for PBE

In order to reduce the SIE effects in situations where no environment is desired (e.g., in the context of simulated mass spectrometry experiments97,98), Olsen et al. recently proposed to use a fragment-based method.15 Since the SIE is known to lead to overdelocalization of charge, constraining electrons to individual fragments can be regarded as a pragmatic way to counteract this error. Clearly, FDE offers the same advantage as the purely electrostatic method employed in ref 15, while being free of system-dependent parameters. This has been demonstrated with respect to spin density distributions14 as well as charge-transfer excitations.12,13 In the following, some tests in vacuo are presented. 6.1. Triglycine. In addition to introducing large errors, the vanishing HOMO−LUMO gap caused by the SIE makes SCF convergence very hard, which was already noted by Jakobsen and co-workers.11 In our tests, the largest zwitterionic oligoglycine which could be converged with the PBE functional in a supermolecular vacuum calculation with ADF was triglycine. Despite the short distance between the charged groups, the MAD in the ESP was found to be almost an order of magnitude larger (0.0133 au) than for a corresponding calculation in solution (0.0018 au). When the system is partitioned into three I


Article


Figure 6. Influence of the nonadditive XC potential on density and ESP errors in 3-FDE calculations for TC5b in solution over the course of five FT cycles (cycle 0 = MFCC). The first label represents the intrafragment potential, whereas the second one refers to the nonadditive part. Errors are given with respect to supermolecular results obtained with M06-2X (density) or CAMY-B3LYP (ESP).

Figure 7. Difference ESP for (Gly)3 in vacuo for supermolecular, MFCC and 3-FDE calculations with the PBE functional (compared to a supermolecular CAMY-B3LYP reference), represented as color-coded density isosurfaces (ρ = 10−4 e·a−3 0 ). Red corresponds to areas of overly negative charge, whereas blue signifies an excess of positive charge. Color scale: ± 0.01 au. The figure was created with VMD.99

with n = 1 and CAMY-B3LYP with n = 2 from the second cycle on. Although this result must be partly attributed to error cancellation, it demonstrates again that fragment-based methods can be regarded as an effective and inexpensive remedy for the SIE. Good results are also obtained for the ESP, with PBE even performing slightly better than the two hybrid functionals. In this case, already the first polarization cycle leads to almost converged errors. The RMSD in the ESP (not to be confused

with the MAD) for insulin in a supermolecular CAM-B3LYP calculation was reported to be 7.7 kJ·mol−1 in ref 15 (in comparison to a local MP2 reference). With the present approach (3-FDE/PBE/n = 1), a deviation of 11.6 kJ·mol−1 (compared to a supermolecular CAMY-B3LYP reference) could be obtained already in the first FT cycle. It is difficult to compare these numbers directly due to the distinct reference methods and basis sets, but clearly the same order of magnitude is reached. J


Article


Figure 8. Density and ESP errors in 3-FDE calculations for insulin in vacuo with different XC functionals and fragment sizes over the course of three FT cycles (cycle 0 = MFCC). Errors are given with respect to supermolecular results obtained with M06-2X (density) or CAMY-B3LYP (ESP).

allow for a direct comparison between DFT results produced with the former and MP2 results obtained with the latter. Instead, the sophisticated M06-2X and CAMY-B3LYP functionals (both available in ADF) were established as reference methods for the electron density and the ESP, respectively. Finally, some additional tests in vacuo were carried out. This can be regarded as a very challenging situation for DFT methods, since the SIE leads to vanishing HOMO−LUMO gaps and convergence problems in this case. It could be shown that by constraining the electrons to their respective fragments, the problematic overdelocalization is strongly reduced, and 3FDE densities and ESPs of similar quality as in solution can be obtained. On the one hand, the method thus reduces the computational cost, especially for larger proteins, which are out of reach for supermolecular calculations. On the other hand, the errors resulting from a purely density-based XC treatment can be strongly reduced, eliminating one of the most severe problems of such formulations. It can be concluded that the 3-FDE method is an important tool with respect to the calculation of first-principles electron densities and produces good results even in challenging situations, such as those presented by charged side chains

The 3-FDE scheme can thus be regarded as a versatile tool for the quantum-chemical investigation of protein properties. Even under very different conditions (in vacuo or with a continuum solvent model), deviations of similar magnitude are obtained, because the method considerably reduces the artifacts caused by the SIE with common DFT approximations.

7. SUMMARY AND OUTLOOK The 3-FDE method has been established as a useful approach to reliably calculate protein densities and ESPs. Several parameters were probed with respect to accuracy of the produced results as well as ease of calculation (SCF convergence, CPU time). In particular, the effects of the XC functional, the fragmentation size, the FT mode, and the role of the nonadditive XC potential were investigated. While the former two largely determine the quality of the results, the FT mode was found to be uncritical, allowing to choose the most suitable mode for the available infrastructure. This holds true for the nonadditive potential as well, especially if larger fragments are used. In the course of the assessment, it was found that the differences between STO and GTO basis sets are too large to K


Article


(22) Svensson, M.; Humbel, S.; Froese, R. D. J.; Matsubara, T.; Sieber, S.; Morokuma, K. J. Phys. Chem. 1996, 100, 19357−19363. (23) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Chem. Rev. 2012, 112, 632−672. (24) Wu, Q.; Van Voorhis, T. Phys. Rev. A: At., Mol., Opt. Phys. 2005, 72, 024502. (25) Wesołowski, T. A.; Warshel, A. J. Phys. Chem. 1993, 97, 8050− 8053. (26) Senatore, G.; Subbaswamy, K. R. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 34, 5754−5757. (27) Cortona, P. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 44, 8454−8458. (28) Jacob, C. R.; Neugebauer, J. WIREs Comput. Mol. Sci. 2014, 4, 325−362. (29) Wesołowski, T. A.; Shedge, S.; Zhou, X. Chem. Rev. 2015, 115, 5891−5928. (30) Jacob, C. R.; Visscher, L. J. Chem. Phys. 2008, 128, 155102. (31) Zhang, D. W.; Zhang, J. Z. H. J. Chem. Phys. 2003, 119, 3599− 3605. (32) Kiewisch, K.; Jacob, C. R.; Visscher, L. J. Chem. Theory Comput. 2013, 9, 2425−2440. (33) Goez, A.; Jacob, C. R.; Neugebauer, J. Comput. Theor. Chem. 2014, 1040-1041, 347−359. (34) Goez, A.; Neugebauer, J. J. Chem. Theory Comput. 2015, 11, 5277−5290. (35) Goez, A.; Neugebauer, J. Mol. Phys. 2016, DOI: 10.1080/ 00268976.2016.1199823. (36) Laricchia, S.; Fabiano, E.; Della Sala, F. Chem. Phys. Lett. 2011, 518, 114−118. (37) Fux, S.; Jacob, C. R.; Neugebauer, J.; Visscher, L.; Reiher, M. J. Chem. Phys. 2010, 132, 164101. (38) Jacob, C. R.; Visscher, L. Towards the Description of Covalent Bonds in Subsystem Density-Functional Theory. In Recent Progress in Orbital-free Density Functional Theory; Wesołowski, T. A., Wang, Y. A., Eds.; World Scientific Publishing: Singapore, 2013. (39) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7, 718− 730. (40) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700−733. (41) Senn, H. M.; Thiel, W. Angew. Chem., Int. Ed. 2009, 48, 1198− 1229. (42) Gao, A. M.; Zhang, D. W.; Zhang, J. Z. H.; Zhang, Y. Chem. Phys. Lett. 2004, 394, 293−297. (43) Wesołowski, T. A.; Weber, J. Chem. Phys. Lett. 1996, 248, 71− 76. (44) Roncero, O.; de Lara-Castells, M. P.; Villarreal, P.; Flores, F.; Ortega, J.; Paniagua, M.; Aguado, A. J. Chem. Phys. 2008, 129, 184104. (45) Roncero, O.; Zanchet, A.; Villarreal, P.; Aguado, A. J. Chem. Phys. 2009, 131, 234110. (46) Lembarki, A.; Chermette, H. Phys. Rev. A: At., Mol., Opt. Phys. 1994, 50, 5328−5331. (47) Wesołowski, T. A.; Chermette, H.; Weber, J. J. Chem. Phys. 1996, 105, 9182−9190. (48) Wesołowski, T. A. J. Chem. Phys. 1997, 106, 8516−8526. (49) Fux, S.; Kiewisch, K.; Jacob, C. R.; Neugebauer, J.; Reiher, M. Chem. Phys. Lett. 2008, 461, 353−359. (50) Götz, A. W.; Beyhan, S. M.; Visscher, L. J. Chem. Theory Comput. 2009, 5, 3161−3174. (51) Schlüns, D.; Klahr, K.; Mück-Lichtenfeld, C.; Visscher, L.; Neugebauer, J. Phys. Chem. Chem. Phys. 2015, 17, 14323−14341. (52) van Leeuwen, R.; Baerends, E. J. Phys. Rev. A: At., Mol., Opt. Phys. 1994, 49, 2421−2431. (53) Amsterdam Density Functional program, Theoretical Chemistry, Vrije Universiteit, Amsterdam. URL: http://www.scm.com (accessed Sept 21, 2016). (54) Te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. J. Comput. Chem. 2001, 22, 931−967.

without a stabilizing environment. The 3-FDE scheme is proven to be less sensitive to the electronic situation in question than regular KS-DFT and proves to be a robust strategy for the investigation of biomolecular properties. Future work will focus on applying the strategies benchmarked in this work to real problems, e.g., the determination of dynamic excitation energy modulation in pigment−protein complexes.33 The 3-FDE method will be of great help in this respect, since it facilitates a parameter-free, fully quantum-chemical simulation of almost arbitrarily large biological entities.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00590. Further results for the systems discussed as well as for two other proteins and several technical tests (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This work was supported by a VIDI grant (700.59.422) of The Netherlands Organisation for Scientific Research (NWO). REFERENCES

(1) Akimov, A. V.; Prezhdo, O. V. Chem. Rev. 2015, 115, 5797−5890. (2) Raghavachari, K.; Saha, A. Chem. Rev. 2015, 115, 5643−5677. (3) Curutchet, C.; Mennucci, B. Chem. Rev. 2016, DOI: 10.1021/ acs.chemrev.5b00700. (4) Neese, F.; Hansen, A.; Liakos, D. G. J. Chem. Phys. 2009, 131, 064103. (5) Riplinger, C.; Neese, F. J. Chem. Phys. 2013, 138, 034106. (6) Ufimtsev, I. S.; Martinez, T. J. J. Chem. Theory Comput. 2009, 5, 2619−2628. (7) Kulik, H. J.; Luehr, N.; Ufimtsev, I. S.; Martinez, T. J. J. Phys. Chem. B 2012, 116, 12501−12509. (8) Rudberg, E. J. Phys.: Condens. Matter 2012, 24, 072202. (9) Perdew, J. P.; Ernzerhof, M.; Burke, K. J. Chem. Phys. 1996, 105, 9982−9985. (10) Antony, J.; Grimme, S. J. Comput. Chem. 2012, 33, 1730−1739. (11) Jakobsen, S.; Kristensen, K.; Jensen, F. J. Chem. Theory Comput. 2013, 9, 3978−3985. (12) Pavanello, M.; Neugebauer, J. J. Chem. Phys. 2011, 135, 234103. (13) Pavanello, M.; Van Voorhis, T.; Visscher, L.; Neugebauer, J. J. Chem. Phys. 2013, 138, 054101. (14) Solovyeva, A.; Pavanello, M.; Neugebauer, J. J. Chem. Phys. 2012, 136, 194104. (15) Olsen, J. M. H.; List, N. H.; Kristensen, K.; Kongsted, J. J. Chem. Theory Comput. 2015, 11, 1832−1842. (16) Gadre, S. R.; Shirsat, R. N.; Limaye, A. C. J. Phys. Chem. 1994, 98, 9165−9169. (17) Babu, K.; Gadre, S. R. J. Comput. Chem. 2003, 24, 484−495. (18) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Chem. Phys. Lett. 1999, 313, 701−706. (19) Fedorov, D. G.; Kitaura, K. J. Phys. Chem. A 2007, 111, 6904− 6914. (20) Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. J. Phys. Chem. A 2001, 105, 293−307. (21) Thompson, M. A. J. Phys. Chem. 1996, 100, 14492−14507. L


Article

Journal of Chemical Theory and Computation (55) Jacob, C. R.; Neugebauer, J.; Visscher, L. J. Comput. Chem. 2008, 29, 1011−1018. (56) Boerrigter, P. M.; Te Velde, G.; Baerends, E. J. Int. J. Quantum Chem. 1988, 33, 87−113. (57) te Velde, G.; Baerends, E. J. J. Comput. Phys. 1992, 99, 84−98. (58) Jacob, C. R.; Beyhan, S. M.; Bulo, R. E.; Gomes, A. S. P.; Götz, A. W.; Kiewisch, K.; Sikkema, J.; Visscher, L. J. Comput. Chem. 2011, 32, 2328−2338. (59) The OPENBABEL package. URL: http://openbabel. sourceforge.net (accessed Sept 21, 2016). (60) Guha, R.; Howard, M. T.; Hutchison, G. R.; Murray-Rust, P.; Rzepa, H.; Steinbeck, C.; Wegner, J.; Willighagen, E. L. J. Chem. Inf. Model. 2006, 46, 991−998. (61) Neese, F. WIREs Comput. Mol. Sci. 2012, 2, 73−78. (62) Weigend, F.; Häser, M. Theor. Chem. Acc. 1997, 97, 331−340. (63) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294, 143−152. (64) Deglmann, P.; May, K.; Furche, F.; Ahlrichs, R. Chem. Phys. Lett. 2004, 384, 103−107. (65) TURBOMOLE V7.0 2015, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since, 2007. Available from http://www. turbomole.com (accessed Sept 21, 2016). (66) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Chem. Phys. Lett. 1989, 162, 165−169. (67) von Arnim, M.; Ahlrichs, R. J. Comput. Chem. 1998, 19, 1746− 1757. (68) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007−1023. (69) Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1993, 98, 1358− 1371. (70) Klamt, A.; Schüürmann, G. J. Chem. Soc., Perkin Trans. 2 1993, 799−805. (71) Becke, A. D. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098− 3100. (72) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (73) Perdew, J. P. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (74) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (75) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1997, 78, 1396. (76) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (77) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200− 1211. (78) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (79) Becke, A. D. J. Chem. Phys. 1993, 98, 1372−1377. (80) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158−6170. (81) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215−241. (82) Gritsenko, O. V.; Schipper, P. R. T.; Baerends, E. J. Chem. Phys. Lett. 1999, 302, 199−207. (83) Schipper, P. R. T.; Gritsenko, O. V.; van Gisbergen, S. J. A.; Baerends, E. J. J. Chem. Phys. 2000, 112, 1344−1352. (84) Gritsenko, O. V.; Schipper, P. R. T.; Baerends, E. J. Int. J. Quantum Chem. 2000, 76, 407−419. (85) Seth, M.; Ziegler, T. J. Chem. Theory Comput. 2012, 8, 901−907. (86) Yanai, T.; Tew, D. P.; Handy, N. C. Chem. Phys. Lett. 2004, 393, 51−57. (87) Yukawa, H. Proc. Phys. Math. Soc. Jpn. 1935, 17, 48−57. (88) Akinaga, Y.; Ten-no, S. Chem. Phys. Lett. 2008, 462, 348−351. (89) Ponder, J. W.; Richards, F. M. J. Comput. Chem. 1987, 8, 1016− 1024. (90) Neidigh, J. W.; Fesinmeyer, R. M.; Andersen, N. H. Nat. Struct. Biol. 2002, 9, 425−430. (91) Deshayes, K.; Schaffer, M. L.; Skelton, N. J.; Nakamura, G. R.; Kadkhodayan, S.; Sidhu, S. S. Chem. Biol. 2002, 9, 495−505. (92) Boldbaatar, D.; Gunasekera, S.; El-Seedi, H. R.; Göransson, U. J. Nat. Prod. 2015, 78, 2545−2551.

(93) Smith, G. D.; Pangborn, W. A.; Blessing, R. H. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2003, 59, 474−482. (94) Daday, C.; König, C.; Valsson, O.; Neugebauer, J.; Filippi, C. J. Chem. Theory Comput. 2013, 9, 2355−2367. (95) Daday, C.; König, C.; Neugebauer, J.; Filippi, C. ChemPhysChem 2014, 15, 3205−3217. (96) Daday, C.; Curutchet, C.; Sinicropi, A.; Mennucci, B.; Filippi, C. J. Chem. Theory Comput. 2015, 11, 4825−4839. (97) Bauer, C. A.; Grimme, S. Org. Biomol. Chem. 2014, 12, 8737− 8744. (98) Grimme, S.; Bauer, C. A. Eur. Mass Spectrom. 2015, 21, 125− 140. (99) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33−38. (100) Frishman, D.; Argos, P. Proteins: Struct., Funct., Genet. 1995, 23, 566−579.

M


Benchmarking Electron Densities and Electrostatic Potentials of

Recommend Documents