Further analysis and comparative study of intermolecular interactions using dimers from the S22 database

Competition Between Halogen Atom and Ring of Halobenzenes as Hydrogen Bond Electron Donor Sites

A halobenzene molecule contains several sites that are capable of acting in an electron-donating capacity within a H-bond. One set of such sites comprise the lone electron pairs of the halogen (X) atoms on the periphery of the ring. The π-electron system above the ring plane can also fulfill this function in many cases. DFT calculations are applied to compare and contrast the propensity of these two site types to engage in such a H-bond within the context of mono, di, tri, tetra, and hexasubstituted halobenzenes. The X atoms chosen for study comprise the full set: F, Cl, Br, and I. It is found that even when the electrostatic potential of the X lone pair is more negative than that above the ring, it is the latter position which is the preferred binding site of HCl in most cases. This preference switches over to the X lone pair only for higher order of substitution, with n=4 or 6. This pattern is explained in large measure by the higher contribution of dispersion when the proton donor is located above the ring.

Ability of the Spectroscopic Properties of the P═Se Bond of a Base to Assess Noncovalent Bond Strength

The viability of the P═Se bond to serve as a monitor of the strength of a noncovalent bond was tested in the context of the (CH3)3PSe molecule. Density functional theory (DFT) computations paired this base with a collection of Lewis acids that spanned hydrogen, halogen, chalcogen, pnicogen, and tetrel bonding interactions and covered a wide range of bond strengths. A very strong linear correlation was observed between the interaction energy and the nuclear magnetic resonance (NMR) 1J(PSe) coupling constant, which could serve as an accurate indicator of bond strength. Also correlating very well with the interaction energy is the stretch of the P═Se bond caused by complexation and the red shift of its stretching frequency. Moderate correlations arise in the chemical shifts of the P and Se nuclei. The σ-hole depth on the Lewis acid is poorly correlated with the energetics, and the same is true for the full electrostatic contribution to the bond energy. Of the various components, it is the polarization energy that correlates most closely with the interaction energy.

Modulating the Competition between Different Atoms to Form Halogen Bonds

I and Br atoms are placed on opposite ends of a n-butyl group, with each allowed to form a halogen bond (XB) with NH3. DFT calculations show that the intrinsic preference of the nucleophile for the heavier I over Br can be reversed by the proper placement of substituents on the alkyl chain. A similar reversal occurs for NH2 and OH groups on the alkyl chain, where substituents make the O a better electron donor than N in an XB to an electrophilic ICCH. The highly mobile π-electron cloud of an aromatic ring makes such reversals much more difficult when the pair of competing atoms are placed on, or within, such a ring.

Participation of transition metal atoms in noncovalent bonds

The existence of halogen, chalcogen, pnicogen, and tetrel bonds as variants of noncovalent σ and π-hole bonds is now widely accepted, and many of their properties have been elucidated. The ability of the d-block transition metals to potentially act as Lewis acids in a similar capacity is examined systematically by DFT calculations. Metals examined span the entire range of the d-block from Group 3 to 12, and are selected from several rows of the periodic table. These atoms are placed in a variety of neutral MXn molecules, with X = Cl and O, and paired with a NH3 nucleophile. The resulting M⋯N bonds tend to be stronger than their p-block analogues, many of them with a substantial degree of covalency. The way in which the properties of these bonds is affected by the row and column of the periodic table from which the M atom is drawn, and the number and nature of ligands, is elucidated.

Spodium Bonding to Dicoordinated Group 12 Atoms

DFT calculations consider the interactions between linear MR2 and a series of N-bases, where M is Hg or Zn and its R substituents are CCH, CN, or NO2. NCH, NH3, and NMe3 were considered as three different N-bases. Zn forms stronger bonds with the N bases than does Hg, and they strengthen along with the electron-withdrawing power of the R substituent, varying over a wide range from 3.4 to 43.9 kcal/mol. Another factor contributing to the bond strength is the nucleophilicity of the base: NCH < NH3 < NMe3. All MR2 Lewis acids can bind at least two bases, which are situated along the R-M-R bisecting plane, fairly close to one another, with θ(N-M-N) angles between 67° and 117°. The presence of a more electron-withdrawing substituent R and more powerful nucleophile allows up to 4 bases to bind to M. The properties of these bonds place them along a continuum, some clearly noncovalent, while other contain a good deal of covalent character.

Halogen Bonding to the π-Systems of Polycyclic Aromatics

The propensity of the π-electron system lying above a polycyclic aromatic system to engage in a halogen bond is examined by DFT calculations. Prototype Lewis acid CF3I is placed above the planes of benzene, naphthalene, anthracene, phenanthrene, naphthacene, chrysene, triphenyl, pyrene, and coronene. The I atom positions itself some 3.3-3.4 Å above the polycyclic plane, and the associated interaction energy is about 4 kcal/mol. This quantity is a little smaller for benzene, but is roughly equal for the larger polycyclics. The energy only oscillates a little as the Lewis acid slides across the face of the polycyclic, preferring regions of higher π-electron density over minima of the electrostatic potential. The binding is dominated by dispersion which contributes half of the total interaction energy.

Anions as Lewis Acids in Noncovalent Bonds

The ability of an anion to serve as electron-accepting Lewis acid in a noncovalent bond is assessed via DFT calculations. NH3 is taken as the common base, and is paired with a host of ACln - anions, with central atom A=Ca, Sr, Mg, Te, Sb, Hg, Zn, Ag, Ga, Ti, Sn, I, and B. Each anion reacts through its σ or π-hole although the electrostatic potential of this hole is quite negative in most cases. Despite the contact between this negative hole and the negative region of the approaching nucleophile, the electrostatic component of the interaction energy of each bond is highly favorable, and accounts for more than half of the total attractive energy. The double negative charge of dianions precludes a stable complex with NH3.

Tetrel Bonding of the Carbenium Ion Forms a Pentacoordinate Carbon Atom

As a flat trigonal species, the CR3 + carbenium ion contains a pair of deep π-holes above and below its molecular plane. In the case of CH3 + a first base will form a covalent bond with the central C, making the combined species tetrahedral. Approach of a second base to the opposite side results in a longer but rather strong noncovalent tetrel bond (TB). While CMe3 + can also form a similar asymmetric complex with a pair of bases, it also has the capacity to form a pair of nearly equivalent TBs, such that the resulting symmetric trigonal bipyramid configuration is only slightly higher in energy. When the three substituents on the central C are phenyl rings, the symmetric configuration with two TBs predominates. These tetrel bonds are quite strong, reaching up to 20 kcal/mol. Adding OPH2 or OCH substituents to the phenyl rings permits the formation of intramolecular C⋅⋅O TBs to the central C, very similar in many respects to the case where these TBs are intermolecular.

Transition from covalent to noncovalent bonding between tetrel atoms

The strength and nature of the bonding between tetrel (T) atoms in R2T⋯TR2 is examined by quantum calculations. T atoms cover the range of Group 14 atoms from C to Pb, and substituents R include Cl, F, and NH2. Systems vary from electrically neutral to both positive and negative overall charged radicals. There is a steady weakening progression in T-T bond strength as the tetrel atom grows larger, transitioning smoothly from a strong covalent to a much weaker noncovalent bond for the larger T atoms. The latter have some of the characteristics of a ditetrel bond, but there are also significant deviations from a classic bond of this type. The T2Cl4- anions are more strongly bonded than the corresponding cations, which are in turn stronger than the neutrals.

Wolfium bonds in homodimers of MX4Y (M = Mo, W; X = F, Cl, Br; Y = O, S, Se)

The term "wolfium bond" has been recently introduced to describe the noncovalent attraction between an atom of group 6 and a nucleophile via a σ-hole binding site. Crystal structures commonly contain a motif wherein two MX4Y units are arranged in close proximity, where M represents either Mo or W, and X and Y refer to halogen and chalcogen atoms respectively. DFT calculations were thus applied to a wide range of homodimers of these molecules so as to assess their preferred arrangements, and to characterize the types of bonding that are present in each in a systematic manner. The most stable Dual-X configuration is symmetric and contains a pair of equivalent M⋯X bonds. The interaction energies range from -8 to -29 kcal mol-1, and are largest for X = F, Y = O, and M = W. The X electron donor is replaced by Y, and the two wolfium bonds are reduced to one, in the less stable Mono-Y structure, with interaction energies between -2 and -10 kcal mol-1. There is some question as to whether the weaker bonds of this type constitute true wolfium bonds.

Influence of Internal Angular Arrangement on Pnicogen Bond Strength

The three Z-X covalent bonds of a ZX3 unit (Z = P, As, Sb, Bi) are normally arranged in a pyramidal structure. Quantum chemical calculations show that pnicogen bonds (ZBs) to the central Z are weakened if ZX3 is flattened, as in the opening of an umbrella. The partial closing of the umbrella has the opposite effect of substantially strengthening these ZBs, even amounting to a 2- or 3-fold magnification in certain cases. The strongest such bonds, wherein Sb and Bi are in a strained configuration within a ZO3CH model system, have interaction energies of 20 kcal/mol with an NH3 base. Most of these systems, whether flattened or more pyramidal, are capable of engaging in three ZBs simultaneously, despite a certain amount of negative cooperativity.

Transition between the Noncovalency and Covalency of σ-Hole Bonds

The properties of the bond between a N-ligand and a Lewis acid containing a σ-hole are studied by quantum chemical methods. Interactions considered include pnicogen bonds involving SbX5, PX5, and PX3, where X represents any of the halogen atoms F, Cl, Br, or I. Also studied are the tetrel bonds of PbX4 and SiX4, as well as the chalcogen bond involving TeOX4. Both NH3 and NCH are applied as two possible bases of differing potency. Some of the bonds are very strong with interaction energies easily exceeding 25 kcal/mol and with AIM bond critical point densities much higher than 0.04 au, suggesting their classification as coordinate covalent bonds. The pentavalent SbX5 and PX5 fall into this category when combined with NH3, as does TeOX4. Although the tetrel bonds involving PbX4 are only slightly weaker, they are probably better viewed as a strong noncovalent bond on the cusp of covalency. Changing the internal bonding of hypervalent SbX5 to the more conventional SbX3 weakens the interaction to a classical noncovalent pnicogen bond. Reducing the base nucleophilicity from NH3 to NCH weakens the bonds so that they are clearly noncovalent.

Heavy pnicogen atoms as electron donors in sigma-hole bonds

DFT calculations evaluate the strength of σ-hole bonds formed by ZH3 and ZMe3 (Z = N, P, As, Sb) acting as electron donor. Bond types considered include H-bond, halogen, chalcogen, pnicogen, and tetrel bond to perfluorinated Lewis acids FH, FBr, F2Se F3As, F4Ge, respectively, as well as their monofluorinated analogues. All of the Z atoms can engage in bonds of at least moderate strength, varying from 3 to more than 40 kcal mol-1. In most cases, N forms the strongest bonds, but the falloff from P to Sb is quite mild. However, this pattern is not characteristic of all cases, as for example in the halogen bonds, where the heavier Z atoms are comparable to, or even stronger than N. Most of the bonds are strengthened by replacing the three H atoms of ZH3 by methyl groups, better simulating the situation that would be generally encountered. Structural and NMR shielding data ought to facilitate the identification of these bonds within crystals or in solution.

Competition between Binding to Various Sites of Substituted Imidazoliums

The imidazolium cation has a number of different sites that can interact with a nucleophile. Adding a halogen atom (X) or a chalcogen (YH) group introduces the possibility of an NX···nuc halogen or NY···nuc chalcogen bond, which competes against the various H-bonds (NH and CH donors) as well as the lone pair···π interaction wherein the nucleophile lies above the plane of the cation. Substituted imidazoliums are paired with the NH₃ base, and the various different complexes are evaluated by density functional theory (DFT) calculations. The strength of XB and YB increases quickly along with the size and polarizability of the X/Y atom, and this sort of bond is the strongest for the heavier Br, I, Se, and Te atoms, followed by the NH···N H-bond, but this order reverses for Cl and S. The various CH···N H-bonds are comparable to one another and to the lone pair···π bond, all with interaction energies of 10-13 kcal/mol, values which show very little dependence upon the substituent placed on the imidazolium.

What Next for Quantum Mechanics in Structure-Based Drug Discovery?

There is significant potential for electronic structure methods to improve the quality of the predictions furnished by the tools of computer-aided drug design, which typically rely on empirically derived functions. In this perspective, we consider some recent examples of how quantum mechanics has been applied in predicting protein-ligand geometries, protein-ligand binding affinities and ligand strain on binding. We then outline several significant developments in quantum mechanics methodology likely to influence these approaches: in particular, we note the advent of more computationally expedient ab initio quantum mechanical methods that can provide chemical accuracy for larger molecular systems than hitherto possible. We highlight the emergence of increasingly accurate semiempirical quantum mechanical methods and the associated role of machine learning and molecular databases in their development. Indeed, the convergence of improved algorithms for solving and analyzing electronic structure, modern machine learning methods, and increasingly comprehensive benchmark data sets of molecular geometries and energies provides a context in which the potential of quantum mechanics will be increasingly realized in driving future developments and applications in structure-based drug discovery.

Accurate prediction of energetic properties of ionic liquid clusters using a fragment-based quantum mechanical method

Accurate prediction of physicochemical properties of ionic liquids (ILs) is of great significance to understand and design novel ILs with unique properties.

A look at the density functional theory zoo with the advanced GMTKN55 database for general main group thermochemistry, kinetics and noncovalent interactions

We present the updated and extended GMTKN55 benchmark database for more accurate and extensive energetic evaluation of density functionals and other electronic structure methods with detailed guidelines for method users.

Highly Selective Halide Receptors Based on Chalcogen, Pnicogen, and Tetrel Bonds

The interactions of halides with a number of bipodal receptors were examined by quantum chemical methods. The receptors were based on a dithieno thiophene framework in which two S atoms can engage in a pair of chalcogen bonds with a halide. These two S atoms were replaced by P and As atoms to compare chalcogen with pnicogen bonding, and by Ge which engages in tetrel bonds with the receptor. Zero, one, and two O atoms were added to the thiophene S atom which is not directly involved in the interaction with the halides. Fluoride bound the most strongly, followed by Cl^- , Br^- , and I^- , respectively. Replacing S by the pnicogen bonds of P strengthened the binding, as did moving down to As in the third row of the periodic table. A further large increment is associated with the switch to the tetrel bonds of Ge. Even though the thiophene S atom is remote from the binding site, each additional O atom added to it raises the binding energy, which can be quite large, as much as 63 kcal mol^-1 for the Ge⋅⋅⋅F^- interaction. The receptors have a pronounced selectivity for F^- over the other halides, as high as 27 orders of magnitude. The data suggest that incorporation of tetrel atoms may lead to new and more powerful halide receptors.

Crystal structure prediction of rigid molecules

A non-polarizable force field based on atomic multipoles fit to reproduce experimental crystal properties andab initiogas-phase dimers is described. The Ewald method is used to calculate both long-range electrostatic and 1/r6dispersion energies of crystals. The dispersion energy of a crystal calculated by a cutoff method is shown to converge slowly to the exact Ewald result. A method for constraining space-group symmetry during unit-cell optimization is derived. Results for locally optimizing 4427 unit cells including volume, cell parameters, unit-cell r.m.s.d. and CPU timings are given for both flexible and rigid molecule optimization. An algorithm for randomly generating rigid molecule crystals is described. Using the correct experimentally determined space group, the average and maximum number of random crystals needed to find the correct experimental structure is given for 2440 rigid single component crystals. The force field energy rank of the correct experimental structure is presented for the same set of 2440 rigid single component crystals assuming the correct space group. A complete crystal prediction is performed for two rigid molecules by searching over the 32 most probable space groups.

Semiempirical Quantum Mechanical Methods for Noncovalent Interactions for Chemical and Biochemical Applications

Semiempirical (SE) methods can be derived from either Hartree-Fock or density functional theory by applying systematic approximations, leading to efficient computational schemes that are several orders of magnitude faster than ab initio calculations. Such numerical efficiency, in combination with modern computational facilities and linear scaling algorithms, allows application of SE methods to very large molecular systems with extensive conformational sampling. To reliably model the structure, dynamics, and reactivity of biological and other soft matter systems, however, good accuracy for the description of noncovalent interactions is required. In this review, we analyze popular SE approaches in terms of their ability to model noncovalent interactions, especially in the context of describing biomolecules, water solution, and organic materials. We discuss the most significant errors and proposed correction schemes, and we review their performance using standard test sets of molecular systems for quantum chemical methods and several recent applications. The general goal is to highlight both the value and limitations of SE methods and stimulate further developments that allow them to effectively complement ab initio methods in the analysis of complex molecular systems.

A computational study of the interaction between dopamine and DNA/RNA nucleosides

The interaction between protonated dopamine and neutral RNA and DNA nucleosides was studied by means of density functional theory calculations in vacuum and in implicit water. On the most stable complexes formed with each of the nucleosides, the vertical absorption excitation energies were evaluated and compared with the values of separated dopamine and corresponding nucleoside. The most stable complex was formed with guanosine and the spectral changes in this complex resulted in a significant reduction of the oscillator strength of the first dopamine's transition. In the first guanosine's transition, a redshift of 0.2 eV was found combined with a reduction of the oscillator strength.

On the origin of the substantial stabilisation of the electron-donor 1,3-dithiole-2-thione-4-carboxyclic acid···I2 and DABCO···I2 complexes

The stabilisation energies of the crystal structures of 1,3-dithiole-2-thione-4-carboxyclic acid···I2 and DABCO···I2 complexes determined by the CCSD(T)/CBS method are very large and exceed 8 and 15 kcal mol(-1), respectively. The DFT-D method (B97-D3/def2-QZVP) strongly overestimates these stabilisation energies, which support the well-known fact that the DFT-D method is not very applicable to the study of charge-transfer complexes. On the other hand, the M06-2X/def2-QZVP method provides surprisingly reliable energies. A DFT-SAPT analysis has shown that a substantial stabilisation of these complexes arises from the charge-transfer energy included in the induction energy and that the respective induction energy is much larger than that of other non-covalently bound complexes. The total stabilisation energies of the complexes mentioned as well as of those where iodine has been replaced by lighter halogens (Br2 and Cl2) or by hetero systems (IF, ICH3, N2) correlate well with the magnitude of the σ-hole (Vs,max value) as well as with the LUMO energy. The nature of the stabilisation of all complexes between both electron donors and X2 (X = I, Br, Cl, N) systems is explained by the magnitude of the σ-hole but surprisingly also by the values of the electric quadrupole moment of these systems. Evidently, the nature of the stabilisation of halogen-bonded complexes between electron donors and systems where the first non-zero electric multipole moment is the quadrupole moment can be explained not only by the recently introduced concept of the σ-hole but also by the classical concept of electric quadrupole moments.

Different aggregation dynamics of benzene–water mixtures

The differences between the molecular aggregations in benzene–water mixtures are identified using all-atom molecular dynamics simulations.

Conformational analysis and parallel QM/MM X-ray refinement of protein bound anti-Alzheimer drug donepezil

The recognition and association of donepezil with acetylcholinesterase (AChE) has been extensively studied in the past several decades because of the former's use as a palliative treatment for mild Alzheimer disease. Herein we examine the conformational properties of donepezil and we re-examine the donepezil-AChE crystal structure using combined quantum mechanical/molecular mechanical (QM/MM) X-ray refinement tools. Donepezil's conformational energy surface was explored using the M06 suite of density functionals and with the MP2/complete basis set (CBS) method using the aug-cc-pVXZ (X = D and T) basis sets. The donepezil-AChE complex (PDB 1EVE) was also re-refined through a parallel QM/MM X-ray refinement approach based on an in-house ab initio code QUICK, which uses the message passing interface (MPI) in a distributed SCF algorithm to accelerate the calculation via parallelization. In the QM/MM re-refined donepezil structure, coordinate errors that previously existed in the PDB deposited geometry were improved leading to an improvement of the modeling of the interaction between donepezil and the aromatic side chains located in the AChE active site gorge. As a result of the re-refinement there was a 93% reduction in the donepezil conformational strain energy versus the original PDB structure. The results of the present effort offer further detailed structural and biochemical inhibitor-AChE information for the continued development of more effective and palliative treatments of Alzheimer disease.

First steps in growth of a polypeptide toward β-sheet structure

The full conformational energy surface is examined for a molecule in which a dipeptide is attached to the same spacer group as another peptide chain, so as to model the seminal steps of β-sheet formation. This surface is compared with the geometrical preferences of the isolated dipeptide to extract the perturbations induced by interactions with the second peptide strand. These interpeptide interactions remove any tendency of the dipeptide to form a C5 ring structure, one of its two normally stable geometries. A C7 structure, the preferred conformation of the isolated dipeptide, remains as the global minimum in the full molecule. However, the stability of this structure is highly dependent upon interpeptide H-bonds with the second chain. The latter forces include not only the usual NH···O interaction, but also a pair of CH···O H-bonds. The secondary minimum is also of C7 type and likewise depends in part upon CH···O H-bonds for its stability. The latter interactions also play a part in the tertiary minimum. A two-strand β-sheet structure is not yet in evidence for this small model system, requiring additional peptide units to be added to each chain.

Open-shell pair interaction energy decomposition analysis (PIEDA): formulation and application to the hydrogen abstraction in tripeptides

An open-shell extension of the pair interaction energy decomposition analysis (PIEDA) within the framework of the fragment molecular orbital (FMO) method is developed. The open-shell PIEDA method allows the analysis of inter- and intramolecular interactions in terms of electrostatic, exchange-repulsion, charge-transfer, dispersion, and optional polarization energies for molecular systems with a radical or high-spin fragment. Taking into account the low computational cost and scalability of the FMO and PIEDA methods, the new scheme provides a means to characterize the stabilization of radical and open-shell sites in biologically relevant species. The open-shell PIEDA is applied to the characterization of intramolecular interactions in capped trialanine upon hydrogen abstraction (HA) at various sites on the peptide. Hydrogen abstraction reaction is the first step in the oxidative pathway initiated by reactive oxygen or nitrogen species, associated with oxidative stress. It is found that HA results in significant geometrical reorganization of the trialanine peptide. Depending on the HA site, terminal interactions in the radical fold conformers may become weaker or stronger compared to the parent molecule, and often change the character of the non-covalent bonding from amide stacking to hydrogen bonding.

The accuracy of quantum chemical methods for large noncovalent complexes

We evaluate the performance of the most widely used wavefunction, density functional theory, and semiempirical methods for the description of noncovalent interactions in a set of larger, mostly dispersion-stabilized noncovalent complexes (the L7 data set). The methods tested include MP2, MP3, SCS-MP2, SCS(MI)-MP2, MP2.5, MP2.X, MP2C, DFT-D, DFT-D3 (B3-LYP-D3, B-LYP-D3, TPSS-D3, PW6B95-D3, M06-2X-D3) and M06-2X, and semiempirical methods augmented with dispersion and hydrogen bonding corrections: SCC-DFTB-D, PM6-D, PM6-DH2 and PM6-D3H4. The test complexes are the octadecane dimer, the guanine trimer, the circumcoronene…adenine dimer, the coronene dimer, the guanine-cytosine dimer, the circumcoronene…guanine-cytosine dimer, and an amyloid fragment trimer containing phenylalanine residues. The best performing method is MP2.5 with relative root mean square deviation (rRMSD) of 4 %. It can thus be recommended as an alternative to the CCSD(T)/CBS (alternatively QCISD(T)/CBS) benchmark for molecular systems which exceed current computational capacity. The second best non-DFT method is MP2C with rRMSD of 8 %. A method with the most favorable "accuracy/cost" ratio belongs to the DFT family: BLYP-D3, with an rRMSD of 8 %. Semiempirical methods deliver less accurate results (the rRMSD exceeds 25 %). Nevertheless, their absolute errors are close to some much more expensive methods such as M06-2X, MP2 or SCS(MI)-MP2, and thus their price/performance ratio is excellent.

Electrostatically embedded generalized molecular fractionation with conjugate caps method for full quantum mechanical calculation of protein energy

An electrostatically embedded generalized molecular fractionation with conjugate caps (EE-GMFCC) method is developed for efficient linear-scaling quantum mechanical (QM) calculation of protein energy. This approach is based on our previously proposed GMFCC/MM method (He; et al. J. Chem. Phys. 2006, 124, 184703), In this EE-GMFCC scheme, the total energy of protein is calculated by taking a linear combination of the QM energy of the neighboring residues and the two-body QM interaction energy between non-neighboring residues that are spatially in close contact. All the fragment calculations are embedded in a field of point charges representing the remaining protein environment, which is the major improvement over our previous GMFCC/MM approach. Numerical studies are carried out to calculate the total energies of 18 real three-dimensional proteins of up to 1142 atoms using the EE-GMFCC approach at the HF/6-31G* level. The overall mean unsigned error of EE-GMFCC for the 18 proteins is 2.39 kcal/mol with reference to the full system HF/6-31G* energies. The EE-GMFCC approach is also applied for proteins at the levels of the density functional theory (DFT) and second-order many-body perturbation theory (MP2), also showing only a few kcal/mol deviation from the corresponding full system result. The EE-GMFCC method is linear-scaling with a low prefactor, trivially parallel, and can be readily applied to routinely perform structural optimization of proteins and molecular dynamics simulation with high level ab initio electronic structure theories.

Perspective: Advances and challenges in treating van der Waals dispersion forces in density functional theory

Electron dispersion forces play a crucial role in determining the structure and properties of biomolecules, molecular crystals, and many other systems. However, an accurate description of dispersion is highly challenging, with the most widely used electronic structure technique, density functional theory (DFT), failing to describe them with standard approximations. Therefore, applications of DFT to systems where dispersion is important have traditionally been of questionable accuracy. However, the last decade has seen a surge of enthusiasm in the DFT community to tackle this problem and in so-doing to extend the applicability of DFT-based methods. Here we discuss, classify, and evaluate some of the promising schemes to emerge in recent years. A brief perspective on the outstanding issues that remain to be resolved and some directions for future research are also provided.

JACOB: a dynamic database for computational chemistry benchmarking

JACOB (just a collection of benchmarks) is a database that contains four diverse benchmark studies, which in-turn included 72 data sets, with a total of 122,356 individual results. The database is constructed upon a dynamic web framework that allows users to retrieve data from the database via predefined categories. Additional flexibility is made available via user-defined text-based queries. Requested sets of results are then automatically presented as bar graphs, with parameters of the graphs being controllable via the URL. JACOB is currently available at www.wallerlab.org/jacob.

Assessment of the PW86+PBE+XDM density functional on van der Waals complexes at non-equilibrium geometries

The deficiency of conventional density-functional theory (DFT) in properly describing van der Waals (vdW) (especially dispersion-bound) complexes has been extensively addressed in the past decade. There are now several new methods published in the literature that are capable of accurately capturing weak dispersion interactions in complexes at equilibrium geometries. However, the performance of these new methods at non-equilibrium geometries remains to be assessed. We have previously published [F. O. Kannemann and A. D. Becke, J. Chem. Theory Comput. 6, 1081 (2010); A. D. Becke, A. A. Arabi, and F. O. Kannemann, Can. J. Chem. 88, 1057 (2010)] that the functional PW86+PBE+XDM for exchange + correlation + dispersion, respectively, is a highly accurate functional for general thermochemistry and vdW complexes at equilibrium geometries. Here, we show that this nonempirical, except for two parameters in the dispersion damping part, functional also performs well for vdW complexes at compressed and stretched intermonomer separations. The mean absolute relative error (MARE) is 9.4% overall for vdW complexes in the "S22×5" database incorporating compressed and stretched geometries [J. Rezac, K. E. Riley, and P. Hobza, J. Chem. Theory Comput. 7, 2427 (2011)]. Our largest MARE on the S22×5 database is 13.3% on the compressed geometry set.

Benchmarking van der Waals density functionals with experimental data: potential-energy curves for H2 molecules on Cu(111), (100) and (110) surfaces

Detailed physisorption data from experiment for the H(2) molecule on low-index Cu surfaces challenge theory. Recently, density functional theory (DFT) has been developed to account for nonlocal correlation effects, including van der Waals (dispersion) forces. We show that the functional vdW-DF2 gives a potential-energy curve, potential-well energy levels and difference in lateral corrugation promisingly close to the results obtained by resonant elastic backscattering-diffraction experiments. The backscattering barrier is sensitive to the choice of exchange functional approximation. Further, the DFT-D3 and TS-vdW corrections to traditional DFT formulations are also benchmarked, and deviations are analyzed.

An Efficient Method to Evaluate Intermolecular Interaction Energies in Large Systems Using Overlapping Multicenter ONIOM and the Fragment Molecular Orbital Method

We propose an approach based on the overlapping multicenter ONIOM to evaluate intermolecular interaction energies in large systems and demonstrate its accuracy on several representative systems in the complete basis set limit at the MP2 and CCSD(T) level of theory. In the application to the intermolecular interaction energy between insulin dimer and 4'-hydroxyacetanilide at the MP2/CBS level, we use the fragment molecular orbital method for the calculation of the entire complex assigned to the lowest layer in three-layer ONIOM. The developed method is shown to be efficient and accurate in the evaluation of the protein-ligand interaction energies.

Calculations on noncovalent interactions and databases of benchmark interaction energies

Although covalent interactions determine the primary structure of a molecule, the noncovalent interactions are responsible for the tertiary and quaternary structure of a molecule and create the fascinating world of the 3D architectures of biomacromolecules. For example, the double helical structure of DNA is of fundamental importance for the function of DNA: it allows it to store and transfer genetic information. To fulfill this role, the structure is rigid to maintain the double helix with a proper positioning of the complementary base, and floppy to allow for its opening. Very strong covalent interactions cannot fulfill both of these criteria, but noncovalent interactions, which are about 2 orders of magnitude weaker, can. This Account highlights the recent advances in the field of the design of novel wave function theory (WFT) methods applicable to noncovalent complexes ranging in size from less than 100 atoms, for which highly accurate ab initio methods are available, up to extended ones (several thousands atoms), which are the domain of semiempirical QM (SQM) methods. Accurate interaction energies for noncovalent complexes are generated by the coupled-cluster technique, taking single- and double-electron excitations iteratively and triple-electron excitation perturbatively with a complete basis set description (CCSD(T)/CBS). The procedure provides interaction energies with high accuracy (error less than 1 kcal/mol). Because the method is computationally demanding, its application is limited to complexes smaller than 30 atoms. But researchers would also like to use computational methods to determine these interaction energies accurately for larger biological and nanoscale structures. Standard QM methods such as MP2, MP3, CCSD, or DFT fail to describe various types of noncovalent systems (H-bonded, stacked, dispersion-controlled, etc.) with comparable accuracy. Therefore, novel methods are needed that have been parametrized toward noncovalent interactions, and existing benchmark data sets represent an important tool for the development of new methods providing reliable characteristics of noncovalent clusters. Our laboratory developed the first suitable data set of CCSD(T)/CBS interaction energies and geometries of various noncovalent complexes, called S22. Since its publication in 2006, it has frequently been applied in parametrization and/or verification of various wave function and density functional techniques. During the intense use of this data set, several inconsistencies emerged, such as the insufficient accuracy of the CCSD(T) correction term or its unbalanced character, which has triggered the introduction of a new, broader, and more accurate data set called the S66 data set. It contains not only 66 CCSD(T)/CBS interaction energies determined in the equilibrium geometries but also 1056 interaction energies calculated at the same level for nonequilibrium geometries. The S22 and S66 data sets have been used for the verification of various WFT methods, and the lowest RMSE (S66, in kcal/mol) was found for the recently introduced SCS-MI-CCSD/CBS (0.08), MP2.5/CBS (0.16), MP2.X/6-31G* (0.27), and SCS-MI-MP2/CBS (0.38) methods. Because of their computational economy, the MP2.5 and MP2.X/6-31G* methods can be recommended for highly accurate calculations of large complexes with up to 100 atoms. The evaluation of SQM methods was based only on the S22 data set, and because some of these methods have been parametrized toward the same data set, the respective results should be taken with caution. For really extended complexes such as protein-ligand systems, only the SMQ methods are applicable. After adding the corrections to the dispersion energy and H-bonding, several methods exhibit surprisingly low RMSE (even below 0.5 kcal/mol). Among the various SMQ methods, the PM6-DH2 can be recommended because of its computational efficiency and it can be used for optimization (which is not the case for other SQM methods). The PM6-DH2 is the base of our novel scoring function used in in silico drug design.

Force-field and quantum-mechanical binding study of selected SAMPL3 host-guest complexes

A Merck molecular force field classical potential combined with Poisson-Boltzmann electrostatics (MMFF/PB) has been used to estimate the binding free energy of seven guest molecules (six tertiary amines and one primary amine) into a synthetic receptor (acyclic cucurbit[4]uril congener) and two benzimidazoles into cyclic cucurbit[7]uril (CB[7]) and cucurbit[8]uril (CB[8]) hosts. In addition, binding enthalpies for the benzimidazoles were calculated with density functional theory (DFT) using the B3LYP functional and a polarizable continuum model (PCM). Although in most cases the MMFF/PB approach returned reasonable agreements with the experiment (±2 kcal/mol), significant, much larger deviations were reported in the case of three host-guest pairs. All four binding enthalpy predictions with the DFT/PCM method suffered 70% or larger deviations from the calorimetry data. Results are discussed in terms of the molecular models used for guest-host complexation and the quality of the intermolecular potentials.

Conformational Analysis of Free and Bound Retinoic Acid

The conformational profiles of unbound all-trans and 9-cis retinoic acid (RA) have been determined using classical and quantum mechanical calculations. Sixty-six all-trans-RA (ATRA) and forty-eight 9-cis-RA energy minimum conformers were identified via HF/6-31G* geometry optimizations in vacuo. Their relative conformational energies were estimated utilizing the M06, M06-2x and MP2 methods combined with the 6-311+G(d,p), aug-cc-pVDZ and aug-cc-pVTZ basis sets, as well as complete basis set MP2 extrapolations using the latter two basis sets. Single-point energy calculations performed with the M06-2x density functional were found to yield similar results to MP2/CBS for the low-energy retinoic acid conformations. Not unexpectedly, the conformational propensities of retinoic acid were governed by the orientation and arrangement of the torsion angles associated with the polyene tail. We also used previously reported QM/MM X-ray refinement results on four ATRA-protein crystal structures plus one newly refined 9-cis-RA complex (PDB ID 1XDK) in order to investigate the conformational preferences of bound retinoic acid. In the re-refined RA conformers the conjugated double bonds are nearly coplanar, which is consistent with the global minimum identified by the Omega/QM method rather than the corresponding crystallographically determined conformations given in the PDB. Consequently, a 91.3% average reduction of the local strain energy in the gas phase, as well as 92.1% in PCM solvent, was observed using the QM/MM refined structures versus the PDB deposited RA conformations. These results thus demonstrate that our QM/MM X-ray refinement approach can significantly enhance the quality of X-ray crystal structures refined by conventional refinement protocols, thereby providing reliable drug-target structural information for use in structure-based drug discovery applications.