Review pubs.acs.org/CR
Large-Scale Computations in Chemistry: A Bird’s Eye View of a Vibrant Field Alexey V. Akimov and Oleg V. Prezhdo* Department of Chemistry, University of South California, Los Angeles, California 90089, United States 3.1.9. Timeline of Semiempirical Methods 3.2. Density-Based Methods 3.2.1. Empirical Density Embedding Schemes: EAM and MEAM 3.2.2. Orbital-Free DFT (OF-DFT) 3.2.3. Timeline of Density-Based Methods 3.3. Bond Order Methods 3.3.1. Bond Order Concept 3.3.2. Bond Order Conservation 3.3.3. Bond Order Potentials and Reactive Force Fields 3.3.4. Construction of Reactive Bond-Order Potentials as a Phenomenological Variational Principle 3.3.5. Timeline of Bond Order Based Methods 4. Computationally Motivated Approximations 4.1. Classification of Computationally Motivated Approaches 4.2. Fragmentation-Based Approaches 4.2.1. Density-Based Fragmentation 4.2.2. Energy-Based Fragmentation 4.2.3. Dynamical Growth with Localization 4.2.4. Diabatic Approaches 4.3. Direct Optimization Methods 4.3.1. Density Matrix 4.3.2. Wave Function or Density Matrix as Dynamical Variables 4.3.3. Orbital Minimization (OM) Methods 4.3.4. Density Matrix Minimization (DMM) Methods 4.3.5. Fermi Operator Expansion (FOE) Methods 4.4. Quantal Force Fields 4.4.1. Basics of QM/MM 4.4.2. Polarizable Force Fields 4.4.3. Effective Fragment Potential (EFP) Method 4.5. Embedding Schemes 5. Conclusions and Outlook 5.1. Methods Diversity: New and Old 5.2. Force Field Development: New Dimensions: Energy → Excited States → Spin? 5.3. Are Semiempirical/Force Field Methods Needed? 5.4. Linear-Scaling Methods from a Different Perspective, Dynamic Programming
CONTENTS 1. Introduction 1.1. The Meaning of “Large Scale” 1.1.1. Large Size 1.1.2. Long Time Scale 1.1.3. Combinatorial Complexity 1.1.4. Use of Many Compute Nodes 1.1.5. Combined Strategies 1.2. When Large-Scale Computations Are Needed 1.3. What Is Large? A Survey of the Records 1.4. Scope and Philosophy of the Review 2. Basic Grounds 2.1. Wave Function Theory (WFT) 2.1.1. Wave Functions and Transformations 2.1.2. Hamiltonian and Energy 2.1.3. Variational Principle and the Eigenvalue Problem 2.1.4. Representation of Operators in Different Bases 2.1.5. Useful Properties 2.2. Density Functional Theory (DFT) 2.3. Limitations of WFT and DFT, and Approaches To Overcome Them 2.3.1. Scaling Law 2.3.2. Classification of Approximations 2.3.3. Sources of Performance Bottlenecks 3. Physically Motivated Approximations 3.1. Semiempirical MO Methods 3.1.1. CNDO and CNDO/2 Methods 3.1.2. INDO and NDDO Methods 3.1.3. MINDO, MNDO, AM1, and PMn Methods 3.1.4. SINDO, SINDO1, and MSINDO Methods 3.1.5. ZINDO Method 3.1.6. Sparkle Model for f-Elements 3.1.7. DFTB and Derived Methods 3.1.8. Extended Hü ckel Theory © 2015 American Chemical Society
5798 5798 5798 5798 5798 5799 5799 5799 5799 5800 5802 5802 5802 5803 5805 5806 5806 5806 5807 5807 5808 5808 5809 5809 5809 5811 5812 5814 5815 5816 5816 5818
5821 5821 5821 5824 5825 5825 5825 5825 5828
5832 5833 5833 5834 5835 5836 5843 5849 5852 5855 5856 5857 5858 5859 5860 5861 5861 5864 5865 5866 5868 5868 5869 5870 5870
Special Issue: Calculations on Large Systems Received: September 18, 2014 Published: April 8, 2015 5797
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews 5.5. Operator Representations, Projectors, Huzinaga−Cantu Equation, and Separability Theory 5.6. Importance of Diabatic Approaches 5.7. Importance of Software Author Information Corresponding Author Notes Biographies Acknowledgments Abbreviations References
Review
the size of a system is size-independent, or constant, although this scaling is never achieved if the atomistic details are accounted for. Thus, the best practical scaling method through which most of the developments are realized is the linear scaling of computational expenses versus the system size. Practically, though, this is also not always possible because the scaling follows a power law or linearity, which becomes evident only for very large systems in which the prefactor of the relationship itself is large. The scaling quickly follows high-power and exponential complexity for more accurate approaches. We will outline the methodology developments on different levels of theory that aim at a linear-scaling relationship in the discussion of the large-size simulations. 1.1.2. Long Time Scale. The group of methods in the longtime scale branch of the large-scale simulation methods mainly addresses the time scale of the processes that can be modeled within a reasonable time frame. This is in contrast to the previous group of methods, which concerns the size of the system that can be handled efficiently. The intrinsic scaling relationship in this branch of methods is linear: the computation of an N-timeslonger process will require an N-times-longer simulation time. Thus, unlike the size-scaling methods, the problem may be perceived as easily overcome. However, the main challenge here lies in the extremely large range of time scales of the processes, varying over many orders of magnitude. If the simulations on the femtosecond (10−15) or picosecond (10−12) time scales are easily achievable nowadays, the simulations on the microsecond (10−6) to second (100) time scales often become prohibitively expensive, if the level of theory remains the same. Of course, with alternative techniques to describe dynamics one can easily overcome this limit. The challenge here then is to encode as much atomistic (or subatomic) details as possible into such methodologies. We want to emphasize that the linear scaling described above does not address the complications arising from statistical samplings and ergodicity issues. For instance, if the thermodynamic or kinetic properties are needed for complex systems, such as proteins, that could exist in a variety of conformations and can be trapped in many local minima, a large number of trajectories are required to attain statistical relevance. A similar situation is encountered in semiclassical approaches to quantum dynamics, which represent nuclear distributions in terms of swarms of trajectories. A large number of such trajectories are typically needed for accurate results, which slows down calculations by several orders of magnitude. An even more problematic situation is observed in the coupled/entangled trajectories methods or in exact quantum mechanical propagation schemes. The computational resources may scale as a power of the number of coupled trajectories, basis states, or grid points, restricting the time scale limits to tens of picoseconds at most. 1.1.3. Combinatorial Complexity. Yet another aspect of large-scale simulations is a situation in which both the size of the system and the time scale of the process (if present) are suitable for existing standard methodologies, but the number of such systems or processes is very large, or multiple variations need to be considered. This can occur with a typical “screening” process, in which the number of possible candidates (e.g., drug molecules, photovoltaic or photocatalytic materials, etc.) is combinatorially large.1−3 Calculations of this type can be organized as highthroughput calculations on a variety of distributed computing systems that need not be homogeneous. Performing the screening of a large database of candidates for the desired target property or application can be prohibitively expensive. This can also be considered a variant of large-size computations.
5872 5872 5872 5873 5873 5873 5873 5873 5873 5874
1. INTRODUCTION Progress in computational chemistry has primarily been motivated by research and development in technological or biomedical applications, as work in these areas can lead to solutions for practical problems that have a direct impact on society. This work is driven by fundamental questions, which the researcher inevitably has to address when explaining the mechanisms of the processes observed experimentally. Although the particular questions and types of systems studied change over time, with currently “fashionable” topics supplanted by newer problems, the general trend in the field of theoretical and computational chemistry remains constant. Namely, there is great interest in performing increasingly larger simulations, including treatment of more complex systems, simulation of longer dynamics, resolution of a larger number of microscopic details, and increase in the diversity of studied systems. The growing capabilities of computational techniques create interesting opportunities for studying a greater range of smallerscale systems, which can be investigated in a brute-force fashion through the use of large-scale computational facilities and resources. The possibility of performing such types of simulations is also motivated by steady progress in computer technologies, encompassing new techniques for mathematical operations in processors, and the use of parallel and highthroughput computing architectures. On the other hand, there is great interest in studying complex systems as a whole, processes that occur in large spatial regions, and processes that occur very slowly. These systems and processes are not amenable to the existing standard techniques, and the need for truly large-scale simulation methodologies cannot be overemphasized. As a result, the computational chemistry community continually faces the stimulating challenge of attempting to reach beyond its current limits. 1.1. The Meaning of “Large Scale”
In general, the term “large scale” can have one of the following five meanings in computational chemistry: (1) large size: power-law and exponential scaling (depending on level of theory) (2) long time scale: linear scaling (yet, in practice, many orders of magnitude would have to be overcome) (3) combinatorial complexity: many variations (distributed computations) (4) use of many compute nodes (5) any combinations of the above 1.1.1. Large Size. As the size of a system increases, the number of interactions increases as well, leading to a proportional growth of the CPU time required for calculations. In an ideal case, the relationship between computational resources and 5798
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
1.1.4. Use of Many Compute Nodes. By this we mean the straightforward utilization of high-performance supercomputers and similar distributed systems, as well as new technologies such as graphical processing units (GPUs). The approach might be thought of as “brute force”, not necessarily employing advanced modeling techniques, but rather attacking the problem’s complexity by efficient parallelization strategies. Of course, advanced methodologies can be and are often augmented with these types of large-scale calculations to make them even more productive and insightful; however, we recognize this type of simulation in our classification to highlight the computational/ parallelization strategies rather than sophisticated physically motivated methodologies. 1.1.5. Combined Strategies. It is often difficult to classify an approach as belonging solely to one of the categories above. Highly productive and insightful simulations do, in fact, combine all or several of the above listed strategies, and are often nonseparable from each other. In the following discussion, we will adopt the term “large-scale” as a generic term for describing one of the five approaches listed above. The first two groups of large-scale computations imply mostly algorithmic and methodological developments, while groups 3 and 4 rely mostly on utilization of large-scale computational facilities, although for group 3 some algorithm/methodology aspects may be important.
using relatively small systems. Other properties that rely on longdistance correlation and complexity of possible reactive channels require large-scale atomistic or coarse-grained models. Examples include reactive dynamics at interfaces,17−23 chemical transformations in energy materials,24−27 combustion28 and catalysis,29,30 and mechanical properties of materials.31−36 Large atomistic models are often needed for studies of transport properties, such as directed motion of biological37 and artificial molecular machines,38−40 surface diffusion via a longrange hopping mechanism,41 and mass and charge transfer in nanoscale structures.42−47 At the quantum level, when quantum dynamics or static ab initio electronic structures are concerned, modeling of reactive processes or large adsorbates requires extended atomistic represention of the substrate, which may be difficult to handle with standard tools. Long-range electrostatic effects may play an important role in charge carrier transport and separation. Similar to avoidance of self-interaction effects in interaction potentials in systems subjected to periodic boundary conditions, artificial quantum confinement should also be minimized, unless the realistic system is very small. Representation of a bulk material should be larger than the exciton size. Calculation of charge carrier transport properties would require inclusion of at least several localized exciton sites, which would further increase the size of the minimal setup.48−50
1.2. When Large-Scale Computations Are Needed
The meaning of the word “large” differs across subfields of computational chemistry and is constantly redefined. In one extreme of the example of a large system is the simulation of large biomolecular complexes. Steady progress has been made over the past few decades in this area. The size and time scales of simulations have increased from approximately 500 atoms simulated over t = 10 ps in 70 s to greater than 2.6 × 106 atoms simulated over t > 1.5 ns in 2000 s. The latter comprise simulations of the following macromolecules: complete satellite tobacco mosaic virus51 (N = 1.07 × 106, t ∼ 10 ns), bacterial flagellar filament52 (N = 2.38 × 106, t ∼ 1.6 ns), and 70S ribosome with mRNA53 (N = 2.64 × 106, t ∼ 4.0 ns). Recently, an allatomic simulation of the mature HIV-1 capsid (N = 60 × 106, t ∼ 100 ns)54 was reported. It constitutes the largest all-atomic system modeled to date, to the best of our knowledge. Type-specific records of systems include the modeling of lipid bilayers and related structures. In 2003 a micelle simulation was reported,55 which included 23 775 particles simulated over a ∼5 ns time scale. A study on ion dynamics and binding in membranes was conducted in 2004, which included about 20 000 atoms simulated over a 200 ns time interval. More recently, in 2011, an even larger all-atomic emulsion simulation was performed.56 This simulation included 63 816 atoms confined in a 60 × 60 × 160 Å box. The use of coarse-grain (CG) approaches allows for the modeling of interesting processes of even larger sizes and time scales, which are comparable to realistic values achievable through experiments. For example, in 2006 a rotation of bacterial flagellum was modeled.57 The simulation included greater than 3400 CG particles, which are equivalent to greater than 15 × 106 atoms in all-atomic model, with a time scale spanning over t > 30 μs. More recently, the liposome has been studied as a potential drug-delivery shuttle for biomedical applications.58 The model included approximately 2500 lipid CG particles and 160 000 water CG particles and the simulation spanned 10 μs. Finally, advances in memory utilization and parallelization strategies allowed simulation of the dynamics in huge molecular systems
1.3. What Is Large? A Survey of the Records
The physics and chemistry of many processes in large-scale systems can be understood in terms of a hierarchy of approximations. This is the typical philosophy of multiscale models, in which effects stemming from lower levels of organization can be incorporated into higher levels via effective mapping. For example, when it is possible to introduce electronic and quantum effects at the atomic level with effective charges4−6 and other interatomic parameters, the atomistic resolution and high-frequency modes can be hidden in a coarse-grain description by constructing effective potentials derived from statistical considerations (e.g., potentials of mean force).7−15 Although multiscale modeling is a well-defined strategy for handling large systems, it is not always a straightforward approach because of the theoretical, computational, and technical difficulties of defining the mapping between different representative levels. In addition, it may not always be suitable in cases when a detailed description must be retained for the entire system of interest at all times. In these situations, the utilization of large systems with explicit representation of all details is necessary by the very nature of the object or the process. For example, biological macromolecules have to be considered explicitly or with a minimal degree of coarse-graining when one is interested in nonlocal effects, such as electron−nuclear correlation16 or reactivity. In other cases smaller atomistic models can be adopted, but at the cost of deteriorated accuracy. For instance, representation of QDs by small clusters can significantly affect their electronic structure, as well as introduce quantum confinement effects not present in a realistic system. Although this rough representation may be the only tractable model, using it can lead to qualitative misrepresentation of the system and its properties. This approach may require careful interpretation of the obtained results or introduction of empirical corrections. The scale of the problem at hand may drive the need for large simulation setups. Some properties, such as heat capacity, mass density, or bulk radial distribution function, can be well described 5799
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
containing 125 nm bilayer vesicles with a 2 μm × 2 μm planar bilayer.59 The system contained 58.5 × 106 CG particles, which is likely the largest simulated system with such a level of theory to date. It should be noted that such a simulation was only possible with the utilization of large supercomputer resources. Extensive reviews of the advances in simulations of biomolecular systems can be found elsewhere,53,60 which include that of coarse-grain simulations.61 The most notable breakthroughs that significantly contribute to this steady progress are the efficient linear-scaling techniques for electrostatic interactions, as well as efficient parallelization methodologies. Generally, computational methods for molecular-mechanicsbased calculations are well-established for most of the central algorithms for the all-atomic interactions already available. The methodological challenges reside mostly in development of accurate potentials, techniques for long-time simulations and efficient sampling, parallelization strategies with favorable scaling, and in flexible and transferable implementation in the community-supported software. Although the challenge of dimensionality is largely resolved, undoubtedly new interesting developments will emerge yet. On the other side of the spectrum are the ab initio and quantum dynamics calculations. Atoms numbering in multiples of tens can very quickly become a really “large” system. This becomes the case when correlated and many-body methods or fully quantal propagation of wave functions is used. On the other hand, in the more popular wave function (WFT) and density functional theory (DFT) methods, the “large” size limit varied through the decades from about 10−20 atoms in the 1960s to currently up to hundreds of atoms nowadays. This transition was led mainly by advances in hardware and computing technologies. Further extensions of this limit will require the utilization of more sophisticated approaches and approximations. One of the first developments in this area was the utilization of an approximation within the WFT or DFT frameworks, which leads to efficient semiempirical and tight-binding approaches.62−64 These methods allow the modeling of systems as large as a few thousand atoms. The use of particle-based quantum dynamics allowed for efficient modeling of dynamics of excited states in large nanoscale structures.5,6 The orbital-free DFT formulation allowed a record simulation of systems with more than 106 atoms.65 Using the effective Hamiltonian approach combined with the classical molecular dynamics (MD) and semiempirical quantum-chemistry methods, the quantum dynamics in soft matter systems have been studied in the Schulten group. The exemplary works include studies of the role of the environment on coherence and charge transfer dynamics in the Fenna−Matthews−Olson complex in the glycerol−water mixture,66 and the record-breaking study of energy transfer in the lamellar chromatophores in Rhodospirillum photometricum.67 The former featured more than 1.4 × 106 QM/ MM computations for each solvent and temperature, and extensive ensemble averaging in the large biomolecular complex. The latter study is remarkable for its size: it comprised 20 × 106 atoms, and the quantum dynamics simulations were performed over 50 ps. The advent and further elaboration of linear-scaling methods68−76 has made it feasible to apply the “standard” ab initio and DFT methods, including some correlated approaches, to systems having over 10 000 atoms. Using one of these methods, a record-breaking simulation of a fullerite system containing over 106 atoms has become possible.77 A single point calculation takes just slightly more than an hour of wall time on a
128 core Xeon cluster. A number of groups are actively developing these methodologies, although the broader scientific community does not yet routinely use them. Likewise, some of these approaches are implemented in a number of available programs, but are not yet common parts of most packages. The development of linear-scaling techniques and their further applications are expected to significantly contribute to advances in computational chemistry in the near future. Finally, another important group of methods is comprised of reactive28,78 and quantal force fields.79−84 These approaches have the combined advantage of being extremely fast through the use of molecular mechanics (classical force field), and being able to describe chemical transformations and excited electronic states (only for quantal FF) via quantum methods. Because of these combined effects, long reactive dynamics of thousands of atoms can be efficiently simulated. Some of the applications for the use of reactive force fields include studies of oxidation and catalytic processes, fracture of nanoscale objects, and interfacial processes. Quantal force fields are more commonly used for biophysical applications, for instance, studies of enzyme reactivity, charge transfer and excitation relaxation in photosynthetic complexes, and reactivity in condensed phases. 1.4. Scope and Philosophy of the Review
It becomes a difficult task to track novel developments and approaches as the number of methods of theoretical and computational chemistry rapidly grows. As such, many groups throughout the world may be focusing on specific methodologies of colleagues in similar disciplines, while overlooking the developments in areas outside their expertise. This narrow approach is further reinforced through the fundamental grounds that the methodologies are based upon, as well as the restricted types of systems and processes under investigation. For example, membrane transport, protein folding, and biomolecule aggregation are often based on coarse-grained techniques with occasional overlap with all-atomic force-field models. Studies of protein−ligand interactions and docking are most commonly based on all-atomic classical molecular mechanics, with an occasional utilization of combined QM/MM schemes. Studies of reactive dynamics in the catalytic centers of enzymes often combine QM/MM and high-level ab initio techniques, while the reactive dynamics in condensed matter systems is typically found through reactive force field approaches. Photoinduced processes require wave function based techniques or perturbative DFT methods. The utilization of the methodologies outside the above spectrum may not always be considered. We believe that an interdisciplinary approach through the use of hybrid techniques analogous to those used in other subfields can be beneficial for scientists having different areas of expertise. There are many excellent reviews and feature papers61,62,76,85−102 on specific techniques targeted to large-scale simulations coming from different perspectives. The recent issue of the Accounts of Chemical Research on fragmentation methods for linear-scaling computations103−111 warrants particular attention. The papers cited above provide essential details, examples of applications, and critical evaluation of the corresponding methodologies. The goal of this review will be to summarize the existing approaches to large-scale calculations from a rather general perspective. We aim to introduce the basic ideas in these approaches and to provide a minimal yet sufficient review on the latest developments across various subdisciplines. We also present a critical evaluation of existing methods and applications, and outline possible future directions in these fields. 5800
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
recent review of Ochsenfeld et al.139 for a more detailed discussion of the FMM-based methods. The details of the FMM implementation in the GAMESS program can be found in ref 140. We also point the reader to the work on matrix reformulation of the FMM method.141 The approaches for accelerated time evolution (group d) are often developed in conjunction with the classical MD methods for systems in the ground electronic state. The techniques include temperature-accelerated dynamics (TAD),142−145 metadynamics,146 hyperdynamics147 and the related bond-boost method, 1 4 8 , 1 4 9 replica-exchange molecular dynamics (REMD),150 and its combination with the temperatureenhanced sampling.151 The discrete MD method developed in a number of groups152−156 obtains notable acceleration of the standard MD calculations. We refer the reader to the existing reviews94,157 for a broader outlook on the enhanced sampling methods. The above techniques achieve long time scales through thouroghly constructed sampling schemes that close the gap between realistic rare events and the atomistic dynamics achievable in computer simulation. It is important to emphasize that the acceleration in this group of methods is not a result of reduction or simplification of the interaction potential. The methods are limited by system size, especially if quantum calculations are to be done. Progress toward joining the two worldsenhanced sampling techniques and linear-scaling methods for quantal calculationshas started only recently,158 and developments in this direction are expected to come out in the near future. Finally, group (e) is represented by the methods for handling larger sizes and time scales via utilization of distributed computations and efficient parallelization stratagies. A notable progress has been achieved within the last subcategory by the rise and widespread utilization of graphical processing units (GPUs). The leading groups developing computational codes and ́ methods with GPUs in mind include the Martinez group (quantum chemisty and dynamics),159−162 the Travesset163 and Schulten164 groups (classical MD), and the Aspuru-Guzik165 and Yasuda groups166 (accelerated computation of molecular integrals in quantum calculations), to name a few. We aim to overview the efficient techniques and strategies of such calculations for researchers studying dynamics, particularly those who are interested in quantum (charge, energy, spin, coherence) and reactive dynamics. Therefore, we pay special attention to presenting methods that have applications to the above dynamics problems. The subjects of extensive timedomain (“long time”) simulations are somewhat beyond the scope of the present review, although they are mentioned where appropriate. While compiling this review, we found the historic method for analysis of various methodologies very fruitful. To understand the diversity of the currently available methodologies, it is often helpful to go back to their origins and study the evolution of the original ideas. The common roots of seemingly different methodologies can be found this way, and other related techniques can be identified. It is common that some of the original ideas are developed quite extensively, while others attract somewhat less interest. However, the less popular ideas still possess potential for further developments, especially in light of the changed paradigm to focus on large-scale systems. We will outline some of those ideas, which might be useful to revitalize. When possible, we provide sufficient details of formulas and derivations. Although in some cases such derivations are rather
In this paper, we aim to avoid repetition of material already presented in specialized reviews. Our strategy is to give an overview of work from seemingly different areas, and present a comparative analysis to discover their common features and possible interplay, as well as their differences. Our goal is to elucidate the hierarchy of the approximations that connects the various techniques. Often, this hierarchy is rather straightforward and follows the natural historic evolution of the methods. Also frequently, different methods may appear at different times and are derived from varying perspectives. Still, they share a common ground and can be related to each other. Sometimes, application of a certain technique developed for specific problems to questions outside its intended scope may yield fruitful results. We outline several such ideas, which we suggest as possible directions for future studies. Although our intent is to present a wide range of methodologies that have applications for systems of differing sizes, we restrict our discussion to methods used for calculation of reactive interactions and excited-state properties in large-scale systems. Other related techniques are mentioned, but are not the main focus. Particularly, we do not discuss any of the following areas: (a) acceleration of computations with correlated electronic structure methods; (b) improved mathematical algorithms for molecular integral calculations and advanced convergence/ summation methods; (c) acceleration of computations with classical force fields; (d) methods for faster time evolution and sampling; (e) parallelization and hardware-related strategies. For completeness, we provide below relevant literature leads for the interested reader. Significant computational advances can be achieved within groups (a) and (b) by utilizing the resolution of identity (RI) technique, also known as density fitting.112−115 The method expands the two-center products entering the four-center molecular integrals either using explicit auxiliary functions, or implicitly. The singular value116 or Cholesky117 decompositions are the most frequent choices in the latter case. An extensive review of the Cholesky method and its utility in the molecular integrals calculations and reduction of the molecular system dimensionality is given by Aquilante et al.118 Beebe et al.117 present a fine description of the decomposition procedure and its application to accelerated computation of the integrals. One of the important aspects of such computations is the fact that numerically insignificant information can be discarded by neglecting the tensor elements smaller than a specified threshold. This simplification can drastically reduce the rank of the problem and lead to improved scalability of the computational costs. The numerically insignificant integrals can be discarded based on the Schwartz inequality119 or more sophisticated criteria, such as mulipole-based integral estimates (MBIE).120 The mathematical grounds of the RI can be related to the inner projection procedure, discussed in the early works of Löwdin.121 A thorough discussion of the RI method is given by Whitten122 and Dunlap et al.123 A number of improvements and extensions of the RI method have been reported, including extension to correlated methods124 and combinations with the divide-andconquer scheme to achieve linear scaling.125 The advanced algorithmic approaches within groups (b) and (c) primarily aim to overcome quadratic scaling of the electrostatic potential computations. The most popular methodology is the one derived from fast multipole moments (FMM).126−128 It includes various methods for classical charged particles,126,129−136 as well as extensions and generalizations for quantum distributions.127,137,138 We refer the reader to the 5801
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
straightforward and can be found in numerous sources, including textbooks and other reviews or original papers, we believe it is crucial to present a consistent and self-contained formulation of the discussed methods. Ideally, such details should be sufficient for the implementation of the methods in computer programs. Due to size restrictions, we will focus on the most central and important details, provide references to more elaborate discussions, and omit some intermediate derivations. One of the motivations behind this approach is to help the interested reader to bridge the gap between theory and practical implementation. We believe that a rigorous and consistent mathematical formulation and intuitively transparent use of notation and illustrations are crucial for proper description of various theories. Finally, we minimize discussion of specific applications, and focus primarily on the theoretical formulations of the methods. Our work puts emphasis on comparative and historic analysis of different methods in one place, providing the reader with a broader outlook and new ideas for future developments of the large-scale computation methodologies.
letters are used to denote AOs and the upper case letters are used to denote atoms. 2.1. Wave Function Theory (WFT)
2.1.1. Wave Functions and Transformations. The simplest ab initio wave function theory (WFT) is the Hartree− Fock method. According to this approach, the wave function of a system with N electrons, Ψ(r1⃗ , ..., rN⃗ ,σ1,σ2, ..., σN), is approximated by a single Slater determinant (SD), Φ0(r1⃗ , ..., rN⃗ ,σ1,σ2, ..., σN): Ψ( r1⃗ , ..., rN⃗ , σ1 , σ2 , ..., σN ) = Φ0( r1⃗ , ..., rN⃗ , σ1, σ2 , ..., σN ) (2.1.1)
Φ0( r1⃗ , ..., rN⃗ , σ1 , σ2 , ..., σN ) 1 N!
=
ψ1( r1⃗ , σ1)
ψ2( r1⃗ , σ1)
... ψN ( r1⃗ , σ1)
ψ1( r2⃗ , σ2) ...
ψ2( r2⃗ , σ2) ...
... ψN ( r2⃗ , σ2) ... ...
ψ1( rN⃗ , σN ) ψ2( rN⃗ , σN ) ... ψN ( rN⃗ , σN ) 1 Â [ψ1( r1⃗ , σ1) ψ2( r2⃗ , σ2) ... ψN ( rN⃗ , σN )] N!
=
2. BASIC GROUNDS In this section, we present the mathematical grounds of the basic theory of molecular interactions. Although the theory is wellknown to most theoretical/computational chemists, it is worth repeating here for the following reasons: (a) Universal notation: We present a common notation for the rest of the paper and facilitate the correlation of quantities presented with those known from other resources. This is motivated by the fact that there is a wide diversity of notations used by different authors, which may be confusing. (b) Clarity: We wish to put all readers, possibly with different levels of familiarity with the subject, on common ground. This may especially be helpful to students and junior scientists, and it might serve to refresh the knowledge of more experienced readers. (c) Foundational knowledge: We present the basic foundations of the molecular interactions methods in order to foresee the possible bottlenecks that could limit the applicability of the methods. This section will help to guide the reader along the hierarchies of further developments in that area of research. (d) Terminology, concepts, and mathematics: Presenting a convenient and systematic notation and basic concepts of quantum chemistry in one place facilitates discussion of the methods to follow, without the need to reintroduce variables, notation, and formulas multiple times. Because the present review covers a variety of advanced linear-scaling methods, the theory recalled in the present section has relevance to the following sections of the review. We adopt the following labeling convention in the present section and the following chapters. The lower case Latin letters such as i, j, k, etc. correspond to general running indices, and also to molecular orbitals. The lower case Latin letters such as a, b, c, etc. correspond to localized (atomic) orbitals. The Greek indices are associated with spin components or Cartesian projections. The upper case Latin letters denote atomic species. The above reservations vary slightly depending on the context of the discussion. In the sections discussing fragment-based methods we use the lower case Latin letters a, b, c, etc. to refer to atoms, and the upper case letters to refer to fragments. In the sections describing atomic orbital (AO) and molecular orbital (MO) transformations and semiempirical methods, the lower case
(2.1.2)
where
 =
∑ (−1)[P]P ̂ (2.1.3)
P
is the antisymmetrizer, P̂ is the permutation operator, and [P] denotes the parity of the permutation P. The lower-case vector symbol ri⃗ is reserved for the position of the ith electron and the Greek letters σi ∈ {α,β} denote spin variables of ith electron. The one-particle functions ψi(r,⃗ σ) are known as molecular spin− orbital and are commonly represented as a product of spatial and spin functions:
ψi( r ⃗ , σ ) = ψi( r ⃗)σ
(2.1.4)
The spin function, σ, typically takes one of the two values:
⎧ |α ⟩ σ=⎨ ⎩ |β ⟩ ⎪ ⎪
(2.1.5)
The vectors |α⟩ = 10 and |β⟩ = 10 are orthonormal: ⎛0⎞ ⟨α|β⟩ = (1 0)⎜ ⎟ = 0 ⎝1 ⎠ (2.1.6a)
()
()
⟨α|α⟩ = ⟨β|β⟩ = 1
(2.1.6b)
and span a two-dimensional space of the spin functions. In a more general treatment, the spin−orbitals transform to spinor functions that can be represented using the above basis as ⎛ ψ α( r ⃗)⎞ i ⎟ = ψ α( r ⃗)α + ψ β( r ⃗)β ψi( r ⃗ , σ ) = ⎜ i i ⎜ ψ β( r )⎟ ⃗ ⎝ i ⎠
(2.1.7)
In the following discussion we will assume the most common representation of the spin−orbitals; ψi′(r,⃗ σ) will imply either ψi(r)⃗ α or ψi(r)⃗ β. To distinguish between the indexing of spin− orbitals and their spatial components (orbitals), we will utilize primed indices for the former and unprimed indices for the latter. The primed indices will run over a wider range and can be mapped to the unprimed indices using the pair convention: i′ = (i, σi′). For example, for a system with Nα spin-up and Nβ spin5802
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
|ψσ ⟩ = ( ψσ ,1 ψσ ,2 ... ψσ , NM )
down electrons, the index i will belong to the range [1, Nα] or to the range [1, Nβ], depending on which quantity is indexed, while the index i′ will belong to the range [1, Nα + Nβ]. The indexed quantity will be additionally indexed with the spin component. Note that when not dealing with orbitals/spin−orbitals (e.g., when enumerating electrons, atoms, excitations, etc.), the primed and unprimed indices do not take effect, unless otherwise noted. The set of 2NM orthonormal spin−orbitals form a Hilbert space, Ω = {|ψi′⟩ |i′ = 1, ..., 2NM}, with the scalar product of two vectors in this space f i′, f j′ ∈ Ω defined by ⟨fi ′ |f j ′ ⟩ =
∫ dσ ∫ d r ⃗ f i*′ ( r ⃗ , σ ) f j′ ( r ⃗ , σ )
and the AO basis (spin-independent functions) as |χ ⟩ = ( χ1 χ2 ... χN ) A
|ψσ ⟩ = ( ψσ ,1 ψσ ,2 ... ψσ , NM ) = |χ ⟩Cσ ⎛Cσ ,11 ⎜ ⎜Cσ ,21 χ χ χ ... =( 1 2 NA )⎜ ... ⎜ ⎜C ⎝ σ , NA1
and with the basis vectors being orthonormal: (2.1.9)
The above space, Ω′, can be factorized into two subspacesone for each spin component: Ω′ = Ωα ⊗ Ωβ
⎛C * ⎜ σ ,11 ⎜ * ⎜Cσ ,21 ⎜ ... ⎜ * ⎝Cσ , N 1
(2.1.10b) Ωβ = {|i′⟩ = |ψi ′⟩ = |ψ(i , β)⟩ = |ψβ , iβ⟩ = |ψβ , i⟩|β⟩ |i = 1, ..., NM}
(2.1.10c)
To clarify the notation, the pair indexing in |ψ(i,α)⟩ denotes the i′ = (i,σi′)th (i′ ∈ [1,2NM]) spin−orbital from the entire set |ψi′⟩ ∈ Ω′, while the double index in |ψσ,i⟩ denotes the ith orbital from the subset for a given spin component σ:|ψσ,i⟩ = |(ψσ)i⟩ = |ψσ⟩i ∈ Ωσ. We will keep various indices in the subscript, reserving the superscript for mathematical operations (inverse, transpose, etc.). The spatial components of the spin−orbitalorbitals |ψα,i⟩ and |ψβ,i⟩need not be the same for any given MO index i. All components should belong to the space of square-integrable functions (according to the basic requirements of quantum mechanics). In addition, we assume that each subset {|ψσ,i⟩ |i = 1, ..., NM forms an orthonormal basis:
A
N
Ĥ el(1, ..., N ) =
(2.1.11)
∑ h(̂ i) + i=1
∑ Ca(i ,σ)|χa ⟩ = ∑ Cσ ,ai|χa ⟩, a
... Cσ*,1NM ⎞ ⎟ ⎟ ... Cσ*,2NM ⎟ ... ... ⎟ ⎟ * ... Cσ , N N ⎠ A M
(2.1.15b)
where NA is the number of atomic orbitals. In general, the number of AOs is not less than that of the MOs, while both numbers are larger than the number of electrons, N. In the following discussion we will refer to the number of occupied alpha orbitals, Nocc,α, and the number of occupied beta orbitals, Nocc,β. In a general unrestricted case, when all one-electron functions are understood as spin−orbitals, the number of occupied orbitals with each spin is equal to the number of electrons with that spin, Nocc,σ = Nσ. 2.1.2. Hamiltonian and Energy. The nonrelativistic fieldfree electronic Hamiltonian of the N-electron system is written:
which is implicitly assumed in eq 2.1.9. Any other set of linearly independent but nonorthonormal functions, {|χa⟩ |a = 1, ..., NM, can be transformed to the orthonormal MO basis by forming appropriate linear combinations: |ψσ , i⟩ =
(2.1.15a)
= ( χ1* χ2* ... χN*A )
(2.1.10a)
∫ d r ⃗ ψσ*,i( r ⃗) ψσ ,j( r ⃗) = δij
... Cσ ,1NM ⎞ ⎟ ... Cσ ,2NM ⎟ ⎟ ... ... ⎟ ... Cσ , NANM ⎟⎠
( ψσ*,1 ... ψσ*, NM ) = |χ * ⟩Cσ*
Ωα = {|i′⟩ = |ψi ′⟩ = |ψ(i , α)⟩ = |ψα , iα⟩ = |ψα , i⟩|α⟩ |i = 1, ..., NM}
⟨ψσ , i|ψσ , j⟩ =
(2.1.14)
the MO-LCAO equation can be written as
(2.1.8)
⟨ψi ′|ψj ′⟩ = δijδσi′σj′
(2.1.13)
1 2
N
∑ i,j=1 i≠j
1 rij
1 h(̂ i) = − ∇i 2 + V (i) 2 Nnucl
σ = α, β
V (i ) = − ∑
a
n=1
(2.1.12)
(2.1.16a)
(2.1.16b)
Zn | ri ⃗ − R⃗ n|
(2.1.16c)
In these equations we utilize shorthand notation for electronic degrees of freedom: i ≡ (ri⃗ ,σi) denotes both spin and coordinate variables of the ith electron. Because the nonrelativistic field-free Hamiltonian, eqs 2.1.16, does not depend on spin variables, the above notation will most often refer to spatial variables only: i ≡ (ri⃗ ). In eqs 2.1.16, rij = |ri⃗ − rj⃗ | is the distance between electrons i and j
Similar to the notation defined for the vectors in eqs 2.1.10, the notation for matrices is as follows: Ca(i,σ) refers to an element from the ath row and i′ = (i, σ)th column of a single combined matrix C with dimensions NA × 2NM, while Cσ,ai refers to an element from the ath row and ith column of a distinct submatrix Cσ with dimensions NA × NM. In the following discussion we will make an association between the basis {|χa⟩} and the AO. Equation 2.1.12 represents the central assumption of the molecular orbital as a linear combination of atomic orbitals (MO-LCAO). For practical purposes and mathematical clarity, the above transformation, as well as many others to follow, can conveniently be written in matrix notation. Defining the MO basis for a spin component σ as
⎛ ⎞ ∂2 ∂2 1 2 1 ∂2 ∇i = ⎜⎜ 2 2 + 2 2 + 2 2 ⎟⎟ 2 2 ⎝ ∂ xi ∂ yi ∂ zi ⎠
is the kinetic energy operator of the ith electron, V(i) is the external potential of nuclei acting on the ith electron, Zn and R⃗ n are the charge and position of the nucleus n, respectively, and 5803
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Nnucl is the number of nuclei in the system. In this paper we use atomic units: ℏ = 1, e = 1, 1/(4πε0) = 1, and me = 1. The operators ĥ(i) are the one-particle operators, describing the energy of the ith electron interacting with an external potential, but not with other electrons. The operators 1/rij are two-particle operators describing the interaction between the pair of electrons i and j. According to the rules of quantum mechanics, the electronic energy of the system is given by the expectation value of the electronic Hamiltonian, eq 2.1.16a. In the Hartree−Fock method, the expectation value is computed using the SD wave functions eq 2.1.2. Using the fact that the SD is normalized to unity, the final expression for energy is N
1 Eel = ∑ ⟨i′|h|̂ i′⟩ + 2 i ′= 1 N
=
∑ (i′|h|̂ i′) + i ′= 1 N
=
∑ hi′ i′ i ′= 1
1 + 2
1 2
describes the interaction between the charge densities created by the one-particle orbitals ψi′ and ψi′, as is particularly apparent in the chemists’ notation, eq 2.1.19a: Ji ′ j ′ =
[⟨i′j′|i′j′⟩ − ⟨i′j′|j′i′⟩]
[(i′i′|j′j′) − (i′j′|j′i′)]
i ′ , j ′= 1
N
[Ji ′ j ′ − K i ′ j ′] (2.1.17)
i ′ , j ′= 1
The angle brackets in eq 2.1.17 represent physicists’ notation of molecular integrals:
≡
∫ d r1⃗ ∫ d r2⃗
i′
1 dσ2 d r2⃗ ψi*′ (1) ψ j*′ (2) ψi ′(1) ψj ′(2) r12
(2.1.21a)
i′
12
j′
j′
σi ′ , σj ′ ∈ {α , β}
(2.1.21b)
K i ′ j ′ = δσi′, σj′K σi′, σj′, ij
(2.1.21c)
∫ d r1⃗ ∫ d r2⃗ ψσ*,i(1) ψσ ,j(1) r1 ψσ*,j(2) ψσ ,i(2) i′
j′
j′
12
i′
The one-particle molecular integrals appearing in eq 2.1.17 hi ′ i ′ = ⟨i′|h|̂ i′⟩ ≡ hσi′i
∫ dσ1 d r1⃗ ∫ dσ2 d r2⃗ ψi*′ (1) ψj*′ (2) r1 ψj′(1) ψi′(2)
(2.1.22)
are called core integrals. The energy, eq 2.1.17, can be written in the basis of molecular orbitals (not spin−orbital) in a more common form:
while parentheses are used in chemists’ notation:
Nβ
Nα
Ji ′ j ′ ≡ (i′i′|j′j′)
Eel =
⎛ ⎞ 1 ψj ′ψj ′⎟⎟ ≡ ⎜⎜ψi ′ψi ′ r12 ⎝ ⎠
∑ hα ,i + ∑ hβ ,i + i=1
+
∫ dσ1 d r1⃗ ∫ dσ2 d r2⃗ ψi*′ (1) ψi′(1) r1 ψj*′ (2) ψj′(2) 12
(2.1.19a)
+
K i ′ j ′ ≡ (i′j′|j′i′)
1 2 1 2
i=1 Nβ
1 2
Nα
Nα
∑ ∑ (Jα ,α ,ij − Kα ,α ,ij) i=1 j=1
Nβ
∑ ∑ (Jβ ,β ,ij − Kβ ,β ,ij) + i=1 j=1 Nβ
1 2
Nα
Nβ
∑ ∑ (Jα ,β ,ij ) i=1 j=1
Nα
∑ ∑ (Jβ ,α ,ij ) i=1 j=1
(2.1.23)
In the restricted formulation, the spatial orbitals for a given pair of alpha and beta spins are the same:
⎛ ⎞ 1 ψj ′ψi ′⎟⎟ ≡ ⎜⎜ψi ′ψj ′ r12 ⎝ ⎠
∫ dσ1 d r1⃗ ∫
(2.1.20)
(2.1.21d)
1 |ψ ψ ⟩ r12 j ′ i ′
(2.1.18b)
≡
j′
r12
i′
j′
K σi′, σj′, ij ≡
12
≡
i′
∫ d r1⃗ ∫ d r2⃗ ψσ*,i(1)ψσ ,i(1) r1 ψσ*,j(2) ψσ ,j′(2),
Jσ , σ , ij ≡
K i ′ j ′ ≡ ⟨i′j′|j′i′⟩
≡
ρσ , i ( r1⃗) ρσ ,2 ( r2⃗ )
j′
(2.1.18a)
≡ ⟨ψi ′ψj ′|
j′
r12
j′
i′
1 |ψ ψ ⟩ r12 i ′ j ′
∫ dσ1 d r1⃗ ∫
=
i′
Ji ′ j ′ = Jσ , σ , ij
Ji ′ j ′ ≡ ⟨i′j′|i′j′⟩ ≡ ⟨ψi ′ψj ′|
∫ d r1⃗ ∫ d r2⃗
The exchange integral, Ki′j′, does not have a straightforward classical analogue, but can be interpreted in terms of the interaction of charge density fluxes between a pair of molecular spin−orbitals. The above definitions of the Coulomb and exchange integrals retain factors that depend on the spin components of the corresponding spin−orbitals. When integrated over spin variables, these factors yield the more common Coulomb and exchange integrals defined in spaces of orbitals, Jij and Kij:
N
∑
=
i′
i ′ , j ′= 1
∑
12
|ψσ , i( r1⃗)|2 |ψσ , j( r2⃗ )|2
= Jσ , σ , ij
N
∑
∫ dσ1 d r1⃗ ∫ dσ2 d r2⃗ ψi*′ (1) ψi′(1) r1 ψj*′ (2) ψj′(2)
|ψα , i⟩ = |ψβ , i⟩ = |ψi ⟩
1 dσ2 d r2⃗ ψi*′ (1) ψj ′(1) ψ j*′ (2) ψi ′(2) r12
(2.1.24)
and hence
(2.1.19b)
The molecular integrals Ji′j′ and Ki′j′ are known as the Coulomb and exchange integrals, respectively. The Coulomb integral Ji′j′ 5804
Jα , α , ij = Jβ , β , ij = Jα , β , ij = Jβ , α , ij = Jij
(2.1.25)
Kα , α , ij = Kβ , β , ij = K ij
(2.1.26) DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
the Fock matrix for spin component σ, in the AO basis:
If we assume Nα = Nβ = N/2, the expression eq 2.1.23 simplifies to the more commonly used form: N /2
Fσ , ab = Hab + Gσ , ab ⇔ Fσ = H + Gσ
N /2 N /2
and the density matrix for spin component σ, in the AO basis:
∑ hi + ∑ ∑ (2Jij − K ,ij)
Eel = 2
i=1
Nσ
(2.1.27)
i=1 j=1
Pσ , ab
Using eqs 2.1.19a and 2.1.19b, the energy expression, eq 2.1.17, can be transformed to a more convenient and commonly used orbital (unprimed) notation: N
Eel =
i ′= 1 N
=
∑
=
N
1 Cai*′Cbi ′(a|h|̂ b) + 2 i ′= 1 a , b = 1
∑ ∑
Nα
NA
∑ ∑ Nβ
Ca*(i , α)Cb(i , α)(a|h|̂ b) +
NA
∑∑
Ca*(i , β)Cb(i , β)(a|h|̂ b)
⎛1 ⎜ ⎜0 Oσ = ⎜ 0 ⎜ ⎜0 ⎝0
i=1 a,b=1
i=1 a,b=1 Nα
NA
∑ ∑
[Ca*(i , α)Cb(i , α)Cc*(j , α)Cd(j , α)(ab|cd)(αα|αα)
i,j=1 a,b,c ,d
− Ca*(i , α)Cb(j , α)Cc*(i , α)Cd(j , α)(ab|cd)(αα|αα)] +
Nβ
Nα
1 2
NA
∑∑ ∑
[Ca*(i , α)Cb(i , α)Cc*(j , β)Cd(j , β)(ab|cd)(αα|ββ)
Nβ
1 2
Nα
NA
∑∑ ∑
[Ca*(i , β)Cb(i , β)Cc*(j , α)Cd(j , α)(ab|cd)(ββ|αα)
i=1 j=1 a,b,c ,d
− Ca*(i , β)Cb(j , α)Cc*(i , β)Cd(j , α)(ab|cd)(βα|βα)] +
Nβ
1 2
[Ca*(i , β)Cb(i , β)Cc*(j , β)Cd(j , β)(ab|cd)(ββ|ββ)
⟨ψσ , i|ψσ , j⟩ =
∑
Pα , abHab +
+
+ 1 = 2 =
1 2
=
1 2 1 2
∑
=
∑ Cσ+,iaSabCσ , bj
= (Cσ+SCσ )ab
NA
= δab
Pβ , ab[(Pα , cd + Pβ , cd)(ab|cd) − Pβ , cd(ad|cb)]
a,b,c ,d
NA
∑ Pα , ab[2Hab + Gα , ab] + 1 2 a,b=1 NA
∑
Pα , ab[Hab + Fα , ab] +
a,b=1
1 2
∑
Cσ+SCσ = I
Pβ , ab[2Hab + Gβ , ab]
a,b=1 NA
∑
Sab = ⟨χa |χb ⟩
a,b=1
Hab = (a|h|̂ b)
(2.1.35b)
The matrix S contains the overlaps of AOs:
Pβ , ab[Hab + Fβ , ab]
(2.1.36)
To find the orbitals that minimize the energy subject to the orthonormalization condition, eqs 2.1.35, one needs to minimize the extended Lagrangian:
(2.1.28)
The letters a, b, c and d denote atomic orbitals, while indices i and j represent MOs. Here we introduce the core Hamiltonian in the AO basis, {|a⟩ ≡ a ≡ χa ≡ |χa⟩}:
L = Eel − tr(Eα(Cα+SCα − I )) − tr(Eβ (Cβ+SCβ − I )) (2.1.37)
(2.1.29)
with respect to each of the sets of coefficients Cσ and Cβ independently. The matrices Eσ contain undetermined Lagrange multipliers and, in general, may have an arbitrary structure, but eventually one often chooses a diagonal form (canonical MO). Omitting the details of the derivation, the result is the wellknown secular (generalized eigenvalue, Hartree−Fock−Roothan, Roothan−Hall) equation:
the two-particle integral (Coulomb and exchange interactions) matrix for spin component σ, in the AO basis: NA
∑ [Pcd(ab|cd) − Pσ ,cd(ad|cb)] c ,d
(2.1.35a)
or in matrix form:
NA
1 1 = tr(Pα(H + Fα)) + tr(Pβ(H + Fβ)) 2 2
Gσ , ab =
∑ Cσ*,aiSabCσ ,bj a,b
Pα , ab[(Pα , cd + Pβ , cd)(ab|cd) − Pα , cd(ad|cb)]
a,b,c ,d
∑
(2.1.34)
∑ Cσ*,ai⟨χa |χb ⟩Cσ ,bj
NA
∑
0⎞ ⎟ 0⎟ 0⎟ ⎟ 0⎟ ... ⎠
a,b
Pβ , abHab
a,b=1
a,b=1
0 0 0 0 0
a,b
NA
NA
(2.1.33)
0 0 1 0 0
i,j=1 a,b,c ,d
− Ca*(i , β)Cb(j , β)Cc*(i , β)Cd(j , β)(ab|cd)(ββ|ββ)] =
0 ... 0 0 0
NA
∑ ∑
(2.1.32)
2.1.3. Variational Principle and the Eigenvalue Problem. To find the wave function that describes the system with the Hamiltonian given by eqs 2.1.16a, the energy, eq 2.1.28, can be minimized with respect to either the expansion coefficients, Cσ, or with respect to the corresponding density matrices, Pσ. The variation must be performed with an additional constraintthe orthonormality of MOs must be preserved:
i=1 j=1 a,b,c ,d
− Ca*(i , α)Cb(j , β)Cc*(i , α)Cd(j , β)(ab|cd)(αβ|αβ)] +
⇔ Pσ = CσOσ Cσ+
The matrix Oσ in eq 2.1.32 is the density matrix for spin component σ (dimension Nσ × Nσ, σ ∈ {σ,β}) in the MO basis, also known as the occupation number matrix. The first Nσ diagonal elements, which correspond to occupied orbitals, are set to 1:
[Cai*′Cbi ′Ccj*′Cdj ′(ab|cd)
i ′ , j ′= 1 a , b , c , d
NA
1 + 2
i=1
(CσOσ Cσ+)ab
P = Pα + Pβ
i ′ , j ′= 1
NA
∑∑
∑ Cσ ,biOσ ,iiCσ+,ia = (CσOσ Cσ+)ba
and the total density matrix: [(i′i′|j′j′) − (i′j′|j′i′)]
− Cai*′Cbj ′Ccj*′Cdi ′(ab|cd)] =
NA
= ∑ Cσ*, aiCσ , bi = i=1
N
1 2
∑ (i′|h|̂ i′) +
(2.1.31)
(2.1.30) 5805
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Fσ Cσ = SCσEσ
X̂ =
(2.1.38)
The operator representation, eq 2.1.42, is convenient because it explicitly shows in which basis the operator is represented. A change of basis can be done using mathematical properties of ket and bra vectors (e.g., transformation eqs 2.1.15) and standard matrix transformations of the matrices representing operator elements in the given basis, eq 2.1.43. 2.1.5. Useful Properties. Discussions in the context of MO theory often involve a number of properties that can be computed from MOs as target quantities or as useful auxiliary properties. Below we present a brief summary of the definition of these quantities and describe some of their properties. (a) Charge density matrix in general basis set {|χi⟩}:
(2.1.39)
One can then obtain eq 2.1.38 by projecting both sides of eq 2.1.39 on the AO basis functions: ⟨χa |Fσ̂ |ψσ , i⟩ = ⟨χa ||ψσ , i⟩Eσ , i ⇔
∑ ⟨χa |Fσ̂ |χb ⟩Cσ ,bi b
=
∑ ⟨χa ||χb ⟩Cσ ,biEσ ,i ,
∀ i ∈ [1, NM ]
b
(2.1.40)
2.1.4. Representation of Operators in Different Bases. In some derivations, mathematical operations may be significantly simplified by utilizing an explicit representation of any arbitrary operator, X̂ , in an arbitrary basis {|χi⟩} in terms of their matrix elements of the operator in the basis, Xij = ⟨χi|X̂ |χj⟩, where the ket, |χi⟩, and the bra, ⟨χj|, vectors are introduced. In this subsection we summarize some useful results related to such a representation and show some properties of the constructed operators. This math will become useful in the discussion of projectors in block-localized MO and DFT approaches, and for Huzinaga−Cantu equations. We assume that the basis {|χi⟩} is not orthonormal, in general, so ⟨χi |χj ⟩ = Sij
∑ |χa ⟩[∑ (S a,b
(2.1.43)
being the matrix elements of the operator in the given basis. Indeed, computing matrix elements in a given basis, we obtain
EF , σ :
∑ (SAa(S−1XS−1)ab SbB)
(2.1.49b)
where f is the Fermi distribution function having width ΔE, gi,σ f(Ei,σ; EF,σ) gives the population (number of electrons) of the gi,σ degenerate energy level Ei,σ for the value of Fermi energy EF,σ.
(2.1.44)
In particular, the identity operator is given by
2.2. Density Functional Theory (DFT)
∑ |χa ⟩[∑ (S−1)ai δij(S−1)jb ]⟨χb | i,j
(2.1.49a)
−1 ⎛ ⎛ Ei , σ − E F, σ ⎞⎞ fΔE (Ei , σ ; E F, σ ) = ⎜1 + exp⎜ ⎟⎟ ⎝ ⎠⎠ ΔE ⎝
= (S ·S −1·X ·S −1·S)AB
a,b
∑ gi ,σ f (Ei ,σ ; EF ,σ ) = Nσ i=1
i,j
a,b
Î =
(2.1.48b)
Nσ
⟨χA |X̂ |χB ⟩ = ⟨χA |(∑ |χa ⟩[∑ (S −1)ai Xij(S −1)jb ]⟨χb |)|χB ⟩
= XAB
∑ ni
where A represents a set of orbitals of interest (group). Often, such groups consists of the atomic orbitals centered on a given atom, in which case eq 2.1.48b gives the atomic Mulliken populations. (c) The Fermi energy, or chemical potential, EF, is one of the central quantities characterizing the energetics of electrons and their mobility in a given system. For a system of Nσ electrons of a given spin σ, the Fermi energy (thus, separate for each spin) is given by
with
=
(2.1.48a)
i∈A
(2.1.42)
a,b
(2.1.47b)
nA =
)ai Xij(S )jb ]⟨χb |
Xij = ⟨χi |X̂ |χj ⟩
D = Dα + Dβ
which is the electron density on the orbital |χi⟩ (often AO). The group-resolved Mulliken populations are the sums of the orbitalresolved populations:
−1
i,j
(2.1.47a)
ni = Dii
(2.1.41)
−1
Dσ = PσS
Pσ and P are the density matrices (for each spin component and the total), eqs 2.1.32 and 2.1.33; S is the overlap matrix, eq 2.1.36 or 2.1.41. The charge density matrix represents the density of electronic charge rather than the density of electronic states, given by the density matrices Pσ and P. In the orthogonal basis, density and charge density matrices are equivalent. (b) Mulliken populations can appear in different projections. For instance, the orbital-resolved Mulliken populations are
Then one can show that the representation of the arbitrary operator X̂ in this basis is X̂ =
(2.1.46)
a,b
The matrix Eσ contains eigenvalues of the Fock operator for the spin channel σ, F̂σ, (with Fσ being the matrix representation of the Fock operator in AO basis). These eigenvalues have a commonly accepted interpretation of the energies of the oneparticle states (MO) |ψσ,i⟩. The Fock operator plays the role of an effective one-particle Hamiltonian. The corresponding oneparticle Schrödinger equation (SE) is then Fσ̂ |ψσ , i⟩ = Eσ , i|ψσ , i⟩
∑ |χa ⟩Xab⟨χb |
The DFT was founded in 1964 by Hohenberg and Kohn (HK)167 and consists of four main postulates/assumptions: (1) Charge density, ρ(r)⃗ , is the main variable. (2) Uniqueness theorem (mapping statement): There exists a 1-to-1 correspondence between the ground state charge density, ρ(r)⃗ , and the external potential acting on it, v(r)⃗ .
(2.1.45)
with the matrix elements ⟨χi|I|̂ χj⟩ = δij being the Kronecker delta symbols. If the basis {|χj⟩} is orthonormal (e.g., MO basis), the expression 2.1.44 simplifies: 5806
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
δEṽ [ρ] = δρ
(3) Universal functional: There exists a universal functional, F[ρ(r)⃗ ], of charge density, ρ(r)⃗ , that is independent of the external potential, v(r)⃗ , such that the electronic energy of the system for a given charge density can be computed by Ev[ρ( r ⃗)] =
∫ v( r ⃗) ρ( r ⃗) d r ⃗ + F[ρ( r ⃗)]
(2.2.8a)
(2.2.1)
∫ ρ( r ⃗) d r ⃗ = N
δTnint + veff ([ρ]; r ⃗) = μ δρ
veff ([ρ]; r ⃗) = v( r ⃗) +
(2.2.2)
F[ρ] = G[ρ] + J[ρ]
∫
The effective one-electron Hamiltonian, H , is known as the KS Hamiltonian. It serves the same role in DFT as the Fock operator does in WFT. The meaning of the KS orbitals, ψKS i , is subject to multiple interpretations. Despite having multiple interpretations of their meaning, they can always be thought of as the auxiliary functions that describe the density of the system of interacting particles: N
ρ( r ⃗) =
2.3. Limitations of WFT and DFT, and Approaches To Overcome Them
2.3.1. Scaling Law. In this section we present sources of the computational limitations of the WF and DF theories that prevent their application to large-scale simulations, or make such simulations difficult. We then outline the general approaches to solving these problems and describe the classification of such approaches, which are adopted in further discussions in this paper. The scaling law often characterizes the performance of most computational chemistry methods:
(2.2.5)
R = A·N B
=0
(2.3.1)
reflects the fact that the computational resources, R (time and memory), increase as the Bth power of some measure of the system’s size, N, with the proportionality prefactor A. The scaling law eq 2.3.1 is often written in the form O(NB), to emphasize the power-law dependence, which is critically important for largescale applications. The dependence on the prefactor A is often weak and can be reduced by various approximations, and at some point (so-called break-even points), the resources required will be dominated by the size of the system, N. The choice of the size measure, N, can vary depending on the particular method used. In molecular mechanics it is typically set to the number of atoms, N = Nat, while in WFT and DFT it is often set to the number of basis (molecular) orbitals, N = NMO. If MOs are represented in the basis of localized AOs, the proportionalities NMO ∼ NAO ∼ Nat hold. If the MOs are represented in the basis of plane waves (PW), the number of the PW orbitals, NPW, is determined by the maximal kinetic energy,
(2.2.6)
+
(2.2.12)
in which the summation runs over N lowest energy orbitals.
where μ is a Lagrange multiplier. Using eqs 2.2.3−2.2.5, one obtains the Euler−Lagrange equation: ⎡
∑ |ψiKS|2 i=1
∫ v( r ⃗) ρ( r ⃗) d r ⃗ + F[ρ( r ⃗)] − μ[∫ ρ( r ⃗) d r ⃗ − N ]
∫ δρ( r ⃗)⎢⎣v( r ⃗) + ∫ | rρ⃗ −( r ′⃗ r)′⃗ | d r ′⃗ +
(2.2.11) KS
Variation of eq 2.2.1 subject to condition eq 2.2.2 is equivalent to the unconstrained variation of the functional:
δEṽ [ρ] = δρ
(2.2.9)
(2.2.10)
1 HKS = − ∇2 + veff ( r ⃗) 2
For simplicity here and in further discussions we will omit the dependence of the charge density on the coordinates of all electrons, keeping in mind that the variable ρ represents a function. This is emphasized by the square brackets that denote a functional dependence on the argument. The functional G[ρ] is unique and, in general, unknown. Finding this functional is the main challenge of DFT. As a result, although the theory described above is formally exact, it has relatively little practical power on its own and critically depends on practical approximations of the universal functional, F[ρ], especially its non-Coulombic part, G[ρ]. This part contains all correlation and exchange effects. Elaboration of the G[ρ] functional was completed in 1965 by Kohn and Sham (KS).168 According to the KS formulation, the G[ρ] functional can be split into the functional that describes the kinetic energy of noninteracting electrons, Tnint[ρ], and the remainder that describes the exchange and correlation of interacting electrons, Exc[ρ] (which, in addition, includes effects arising from the kinetic energy of the interacting electrons):
Eṽ [ρ] =
δExc[ρ] δρ
HKSψiKS = EiKSψiKS
(2.2.4)
G[ρ] = Tnint[ρ] + Exc[ρ]
∫ | rρ⃗ −( r ′⃗ r)′⃗ | d r ′⃗ +
The Schrödinger equations for each of the noninteracting single particles moving in the field of effective potential, eq 2.2.9, are known as the KS equations:
(2.2.3)
ρ( r ⃗) ρ( r ′⃗ ) d r ⃗ d r ′⃗ | r ⃗ − r ′⃗ |
(2.2.8b)
)with
where N is the total number of electrons in the system. The energy functional, eq 2.2.1, can further be elaborated by separating the Coulombic interactions of the electrons, Vee[ρ], from the universal functional:
1 2
δρ
⎤ + veff ( r ⃗) − μ⎥ d r ⃗ = 0 ⎦
or
(4) Variational principle: The ground state charge density is obtained by minimization of the energy functional, eq 2.2.1, with respect to the charge density, subject to conservation of the total number of electrons:
J[ρ] =
⎡ δTnint
∫ δρ( r ⃗)⎢⎣
δTnint δρ
⎤ δExc − μ⎥ d r ⃗ ⎦ δρ (2.2.7)
Equation 2.2.7 can be represented in a more convenient form that corresponds to the system of noninteracting electrons moving in the field of effective potential, veff(r)⃗ : 5807
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Ecut, of plane waves and by the volume of the system, V: NMO ∼ NPW ∼ VEcut2, and is formally independent of the number of atoms in the system. In many cases, the volume of the simulation cell is proportional to the number of atoms, making the base linearly dependent on the number of atoms in a system as well: NMO ∼ NPW ∼ NatEcut2. Often, not all MOs are required in computations, but only the occupied and, perhaps, some small fraction of the nearby-unoccupied MOs. In this case the base may be reduced, N = NMO(occ), which can be particularly important as the size of the system increases. 2.3.2. Classification of Approximations. Depending on the specific choice of N, other proportionality factors can be adsorbed into the prefactor A. The efforts to minimize the base parameter, N, or the prefactor A can be classified as intensive, because they typically do not take advantage of the opportunities given by the size of the system. These approaches can speed up calculations, but only up to a certain extent, after which they can be prohibitive with the current state of computational technologies. The main challenge of most computational chemistry methods applied to large-scale systems is in reducing the exponent, B, to produce linear (B = 1) or sublinear (B < 1) scaling methods. The desired techniques would take advantage of the opportunities provided by the system size. Therefore, these methods can be classified into a group of extensive approximations. The simplest hierarchy of extensive approximations is the following: correlated WFT methods (B > 4) − Hartree−Fock (B = 4) − DFT (with pure functionals)/semiempirical methods (B = 3) − force fields with full electrostatics (FF, B = 2) − linear-scaling methods (N). Methods that minimize A, N, and B can be classified according to the nature of the approximations and techniques used. In many cases the approximations are based on the use of computationally efficient alternatives for some of the quantities involved in calculations. For example, certain integrals may be approximated by functions that are not computationally intensive, which behave similarly to the integrals of interest. The methods can use various distance cutoffs or value thresholds to neglect some of the terms. Purely mathematical or programming tricks can be used to accelerate some calculations. For example, efficient schemes can be used to accelerate the convergence of self-consistent field iterations. Numerous partitioning methods (such as divide-and-conquer, subsystem DFT, or fragment MO) can also be used. Unlike the previous methods, which minimize the parameters A and N, the latter are used to minimize parameter B. Despite the differences in the ways in which these methods affect the scaling law, eq 2.3.1, none of them alter the physics of or the level of description of these interactions. In other words, they introduce minimal approximation regarding the physics of the interactions. We classify these techniques as computationally motivated. The complementary group of approximations can be classified as physically motivated ones. The approximations in this group are mostly driven by more fundamental assumptions, which are intended to maintain an accurate account of the physics of the process, but describe it in a different way. Often, this alternative way of representing the physics uses more computationally efficient formulations. The best example of physically motivated approximations is DFT, which appeared as an inexpensive alternative to correlated WFTs. Another example of physically motivated models is a transformation from the ab initio Hartree− Fock theory to semiempirical methods. Finally, one can think of the model Hamiltonian approaches as yet another level of
physically motivated models that attempt to mimic more sophisticated all-atomic orbital-based models. The two groups of approximationscomputationally and physically motivatedcan both lead to a substantial acceleration of computations, but through different means. Often, the two types of approximations are present in an interconnected way, making it difficult to distinguish and classify them into one particular group. For example, the formulations of semiempirical methods make use of both computationally motivated approximations, such as the Nishimoto−Mataga or Ohno formulas to represent two-electron integrals. They may also neglect orbital overlaps and physically motivated approximations such as introduction of suitable empirical functions to describe nuclear−nuclear or electron−nuclear interactions or ionization potentials and electronic affinities to approximate certain core matrix elements. 2.3.3. Sources of Performance Bottlenecks. In the simplest quantum mechanical treatment of the electronic structurethe Hartree−Fock techniquethe scaling is quartic, with the inclusion of the number of orbitals, O(Norb4). This scaling arises at the Fock matrix formation step because it is necessary to calculate all of the four-center, two-electron (Coulomb and exchange) integrals, eqs 2.1.18. The scaling of the higher order correlated methods is even more unfavorable, restricting the size of the systems to a few tens of atoms at most; examples include the Møller−Plesset perturbation theory (MP2/MP4), configuration interactions (CI), coupled cluster (CC), and complete active space SCF (CASSCF). In the present review we do not discuss approaches to improving the scalability of those methods. Rather, we will only discuss the methodologies that scale better than the Hartree−Fock O(Norb4) method. The next level of scaling, O(Norb3), is achieved in the semiempirical WFT and DFT (with pure functionals). Both techniques avoid the computation of all or part of the exchange and Coulomb integrals that enter the Fock matrix. In the semiempirical branch this is achieved by neglecting the terms involving the products of the atomic orbitals localized on different atoms, as in the complete neglect of differential overlap (CNDO) technique. Alternatively, in the neglect of diatomic differential overlap (NDDO) approximation, only those exchange integrals that involve products of the orbitals centered only on two atoms (or one atom) are retained. In the KS-DFT the exchange functional depending only on the electronic density effectively represents exchange integrals. In addition, the correlation absent at the HF level can also be introduced in a similar manner, via the correlation functional. The cubic scaling arises because of two reasons: (a) the necessity of computing the density matrix, eqs 2.1.32 and 2.1.33, which involves dense matrix multiplications; (b) the diagonalization step. Although finding only the eigenvalues of the N × N matrix may be achieved in linear time, the eigenvectors (MOs) are typically required, which makes the algorithm scale as O(N3). The next level of scaling, O(N2), can typically be achieved in classical force field/molecular mechanics methods. This scaling arises from pairwise long-range interactions (Coulombic), while the so-called bonded (bonds, angles, torsions, dihedral) and short-range nonbonded (vdW) interactions can be computed in linear time, O(N). The quadratic scaling can be reduced to linear, O(N), or quasi-linear, O(N log N), by utilizing efficient summation techniques (Ewald, particle-mesh Ewald, etc.) for long-range electrostatic interactions, especially in periodic systems. Therefore, the standard (nonpolarizable) force field methods show the promising quasi-linear scaling O(N) (often 5808
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
O(N log N) making it possible to study huge atomistic models comprising over several million atoms with only modest requirements for memory and CPU time. The MM methodologies become more complicated when reactive many-body (such as Sutton−Chen, MEAM, Tersoff− Brenner, ReaxFF, etc.) potentials are used or polarization is included. For example, potentials such as MEAM involve nonseparable nonpairwise terms and are to some extent similar to density functional methods, although they do not require diagonalization/electron density optimization. One of the simplest ways of treating the polarization of classical atoms can be accomplished with charge equilibration (QEq) schemes. Such schemes require the solution of the system of linear equations that introduces cubic scaling O(N3). Although N in this case is notably smaller than Norb, the complexity is still cubic, comparable to that obtained from DFT and semiempirical methods. Finally, we arrive at a group of methods that scale sublinearly with the number of atoms. A variety of coarse-graining (CG) force fields (CG-FF) has been developed, primarily in the biophysical community. These methods treat groups of atoms as effective particles, parametrized to thermodynamic and structural data that are obtained either from all-atomic methods or experimentally. CG-FFs relate to all-atomic FFs (AA-FF), analogous to the relationship of the fragment MO (FMO) or subsystem DFT methods and standard MO and KS-DFT methods, respectively. Alternatively, in the mapping approaches, smaller subsets of effective collective coordinates can be found that describe processes of interest in large-scale systems. In the field of classical MM simulations this is known, for example, as the diffusion map MD. This approach effectively reduces the dimensionality of the problem at hand, leading again to sublinear scaling. In the field of quantum dynamics this approach has widely been used in the framework of model Hamiltonians. Many of the topics mentioned above have been extensively reviewed and discussed in detail previously. As we stated in section 1.4, this review will discuss only some of the promising techniques for computing electronic structure and reactive interactions in large-scale systems. Our review is motivated by the desire to understand the potential of these methods for dynamical problems, including the reactivity on the ground state PES and the nonadiabatic quantum dynamics on the manifold of excited state PESs.
grounds, choosing one of them depending on a desired degree of acceleration. 3.1. Semiempirical MO Methods
This group of methods appeared first, in order to circumvent the quartic scaling of the Hartree−Fock method. The main assumption that makes the semiempirical methods efficient is the neglect of a large portion of molecular integrals from eq 2.1.30 or its approximation with more computationally efficient formulas. Typically, the calculation of the Fock matrix scales quadratically, due to the necessity to account for the pairwise electrostatic interactions of the Coulomb type. In the tightbinding methods without electrostatic effects, such as the extended Hückel theory (EHT) or the density functional tightbinding (DFTB), formation of the Fock (or effective oneelectron Hamiltonian) matrix can be linear, due to the shortrange nature of interactions. The purpose of this section is to summarize some of the central steps in the development of the semiempirical methods, minimizing discussion of the details. Such discussion can be found in the dedicated reviews on the semiempirical methods.62 Rather, we demonstrate evolution of the semiempirical methods and discuss their current status. In certain cases, when reviews are lacking (such as on the EHT), a slightly more extended discussion will be presented. 3.1.1. CNDO and CNDO/2 Methods. 3.1.1.1. Basic Formulation of the CNDO Method. The complete neglect of the differential overlap (CNDO) method, developed by Pople and co-workers,169,170 was one of the first approaches derived directly from the HF theory. The method aimed to reduce the number of two-electron integrals in the Fock matrix, and to preserve the rotational invariance of the resulting approximate Fock operator. The main assumptions of the method can be summarized as follows: (1) Neglect of differential overlap is the general (main) approximation: +
Ca⃗ C⃗b = δab ⇔ S = I
(3.1.1)
Under this assumption, some Hamiltonian matrix elements become zero (see next), and the generalized eigenvalue problem, eq 2.1.38, becomes the standard eigenvalue problem: (3.1.2)
FC = CE
(2) Approximations regarding the two-electron integrals are the following. Because of eq 3.1.1, all two-electron integrals (ab| cd) are zero, unless a = b and c = d:
3. PHYSICALLY MOTIVATED APPROXIMATIONS This group appeared relatively early in the development of quantum chemical methods. Some methods were invented very early, as the most natural and logical ways of reducing computational costs of higher-level approaches. Other methods were developed relatively recently, although the basis for their rise has been forged in earlier developments. Some methods in this group were developed as nonquantum methods, relatively independently of the evolution of the mainstream quantum chemistry methods. We list them in this section as well, because they can be derived from the high-level quantum approaches, by a sequence of appropriate approximations, and because they show big promise for large-scale computations. In contrast to the computationally motivated approaches, which utilize mathematical advantages of method reformulation (although often guided by chemical and physical concepts), the group of methods discussed in this section starts from the physical interpretation first and often relies on less rigorous mathematical or physical
⎧(aa|cc) ≠ 0, a = b , c = d (ab|cd) = δabδcd(aa|cc) = ⎨ ⎩ 0, a ≠ b or c ≠ d ⎪
⎪
(3.1.3)
The nonzero two-electron integrals do not depend on orbital type, to preserve rotational invariance of the Hamiltonian: (aa|bb) = γAB ,
a ∈ A, b ∈ B
(3.1.4)
where A can be the same as B. Here, and in the following discussion, capital letters denote atom labels, and the notation a ∈ Aimplies that the AO a is centered on the atom A. For practical reasons, the integral is evaluated in STO (or in STO-nGTO) using s-orbitals of the valence shell. (3) Approximations regarding the one-electron integrals involve diagonal elements of the core Hamiltonian. Because of 5809
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
⎧ ⎪(a − 1 ∇2 − ∑ V ( r ⃗ , R⃗) b) C ⎪ 2 C ⎪ Hab = ⎨ = Uab − ∑ VAB , a ∈ A , b ∈ A ⎪ B≠A ⎪ ⎪ β 0 S , a ∈ A , b ∈ B, A ≠ B ⎩ AB ab
eq 3.1.1, the one-electron integrals (a|VB|c) will be nonzero only if a = c: ⎧(a|VB|a) ≡ VAB , (a|VB|c) = ⎨ ⎩ 0, a ≠ c ⎪
a ∈ A, A ≠ B
⎪
(3.1.5)
Note that the element(a|VA|a), a ∈ A, is implicitly included in the parametrization of the core Hamiltonian. (4) Off-diagonal matrix elements of the core Hamiltonian are proportional to the AO overlaps, similar to the EHT prescription: 0 Hab = βAB Sab ,
a ≠ b, a ∈ A , b ∈ B
where ⎧Uaa , parametrized 1 Uab ≡ (a − ∇2 − VA( r ⃗ , R⃗) b) = ⎨ 2 ⎩ 0, a ≠ b ⎪
⎪
(3.1.12)
⎧ ⎪(a|V ( r ⃗ , R⃗ )|a), B VAB ≡ (a|VB( r ⃗ , R⃗)|b) = ⎨ ⎪ ⎩ 0, a ≠ b
(3.1.6)
Note that the off-diagonal elements of the overlap matrix are computed only for this purpose, but they are assumed zero to simplify the eigenvalue problem. 3.1.1.2. Explicit Expression of the Fock Matrix Elements. Using the approximations above, the Fock matrix elements can be written as
(3.1.13)
γAB = (nsAnsA|nsBnsB) =
1
To preserve rotational invariance, the integrals are computed using s-type orbitals of the valence shell given by STOs with parametrized exponents:
∑ [Pcdδabδcd(ab|cd) − Pcd , σδadδbc(ad|bc)]
ns(r ) ∼ exp( −ξr )
= Hab + (δab ∑ Pcc(aa|cc)) − Pba , σ(aa|bb) c
= Hab + (δab ∑ ∑ Pcc(aa|cc)) − Pab , σ(aa|bb) C
ξ=
c∈C
c∈C
= Hab + δab ∑ γACPC − Pab , σγAB C
(3.1.7)
or, with a change of indexing: Fab , σ = Hab + δab ∑ γABPB − Pab , σγAB B
∑ Pbb
(3.1.16)
Z
∫ nsA 2(1) | r ⃗ −BR⃗ | d r ⃗ B
(3.1.8)
(3.1.17)
(3) The core integrals Uaaare parametrized from the IP and EA energies and the ERIs, eq 3.1.14. For example, one obtains for the 2sm2pn configuration:
(3.1.9)
is the total electron density on atom B. Equation 3.1.8 can be elaborated further:
Fab , σ
Z′ n
(a|VB( r ⃗ , R⃗)|a) =
where
b∈B
(3.1.15)
where n is the principal quantum number of the valence shell and Z′ is a parameter (with the physical meaning of the screened core charge, but numerically not necessarily equal to it). (2) The nuclear attraction integrals (NAI) of type (a|VB)(r,⃗ R⃗ )| a), a ∈ A, are evaluated explicitly, using valence s-orbitals, as in eq 3.1.15:
= Hab + δab ∑ (aa|cc) ∑ Pcc − Pab , σ(aa|bb) C
2
(3.1.14)
c ,d
PB =
∬ nsA 2(1) | r ⃗ −1 r ⃗ | nsB2(1) d r1⃗ d r2⃗
∑ [Pcd(ab|cd) − Pcd , σ(ad|bc)] c ,d
= Hab +
a∈A
3.1.1.4. Some Technicalities. This section includes some technicalities. (1) The electron repulsion integrals (ERI) of type γAB are computed explicitly:
Fab , σ = Hab + Gab , σ = Hab +
(3.1.11)
U2s,2s = −Is(X , 2sm2pn) − (Zeff, X − 1)γXX
(3.1.18a)
U2p,2p = −Ip(X , 2sm2pn) − (Zeff, X − 1)γXX
(3.1.18b)
The parameters Is(X,2sm2pn) and Ip(X,2sm2pn) can be taken from atomic data (experimental or computational) or can be adjustable. (4) The quantities β0AB are averages of the single-atom parameters β0X:
⎧ H + (P − P )γ + ∑ γABPB , a ∈ A B aa , σ AA ⎪ aa B ≠A =⎨ ⎪ ⎩ Hab − Pab , σγAB , a ≠ b , a ∈ A , b ∈ B
0 βAB =
(3.1.10)
3.1.1.3. Explicit Expression for the Core Hamiltonian Matrix Elements. The explicit expression for the core Hamiltonian matrix elements is
1 0 (β + βB0) 2 A
(3.1.19)
(5) Note that so far we have introduced three types of Z: (a) Z is the nuclear charge. It is used in computation of the NAI, eq 3.1.17, and nuclear−nuclear interaction energy. 5810
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
(b) Zeff = m + n is the number of valence electrons in the 2sm2pn configuration. This quantity plays the role of the effective core charge (but in a crude way, for example, not obtained from the Slater rules). It is used in the core parameter evaluation, eqs 3.1.18. (c) Z′ is the effective core charge, more closely related to the Slater rules and nuclear charge screening. Solving the hydrogenlike problem with this effective charge gives the exponent of the resulting STOs. It is used in the NAI and ERI evaluations as the parameter entering the AO exponents. 3.1.1.5. CNDO/2 Method. The original CNDO approach was valuable for the calculation of energies and electronic structure, but it failed notably in predicting good nuclear geometries. The method has been significantly improved to yield its second version, abbreviated CNDO/2.171 In addition to the approximations made by CNDO, CNDO/2 relies on the following two approximations: (1) Neglect of the so-called “penetration” term: To derive this correction, one can start with the CNDO Fock matrix and regroup the terms in the following way: Faa , σ = Haa + (PB − Paa , σ )γAA +
(ab|cd) = δ[a],[b]δ[c],[d](ab|cd)
where [a] denotes the index of the atom on which the AO χais centered. The Fock matrix can be written with this approximation as Fab , σ = Hab + Gab , σ = Hab + = Hab +
= Uaa + (PB − Paa , σ )γAA +
− Pcd , σδ[a],[d]δ[b],[c](ad|bc)] = Hab + δ[a],[b] ∑ −
∑
Fab , σ = Hab +
∑∑
Pcd(ab|cd) −
∑
Pcd , σ(ad|bc),
c ,d∈A
a, b ∈ A
B≠A
Fab , σ = Hab −
(3.1.26a)
∑ Pcd ,σ(ad|bc),
a ∈ A , b ∈ B, A ≠ B
c∈B d∈A
(3.1.26b) (3.1.20)
with the core Hamiltonian, Hab, given by eq 3.1.11. In contrast to CNDO and CNDO/2, not only the diagonal elements, but all elements Uaband VABin eqs 3.1.12 and 3.1.13 are computed explicitly (or parametrized). This is because the overlap χ*a χb is not neglected. The INDO approximation is essentially the same as the CNDO/2, with the exception that all one-center exchange-type integrals are retained. The approximation 3.1.3 is modified:
(3.1.21)
(2) Utilization of a different strategy for parametrization of the core Hamiltonian, using both IP and EA: −Ia = Uaa + (Zeff, A − 1)γAA
(3.1.22a)
−Aa = Uaa + Zeff, AγAA
(3.1.22b)
(ab|cd) = δ[a],[c]δ[a],[b]δ[c],[d](ab|cd)
(3.1.27)
The Fock matrix can be written as Faa , σ = Uaa +
∑ (PBB − ZB)γAB + ∑ B≠A
so that ⎛ 1 1⎞ Uaa = − (Ia + Aa) − ⎜Zeff, A − ⎟γAA ⎝ 2 2⎠
(3.1.25)
B c ,d∈B
The last sum ∑B≠A(VAB − Zeff,BγAB) has the meaning of “penetration” energy, although there is no strict justification. According to the approximation, this term is neglected. Effectively, the approximation is equivalent to the following:
VAB = Zeff, BγAB
Pcd , σ(ad|bc)
or in the more explicit form:
∑ γAB(PB − Zeff,B)
∑ (VAB − Zeff, BγAB)
Pcd(ab|cd)
c ∈ [b] d ∈ [a]
∑ γABPB − ∑ VAB
B≠A
∑
B c ,d∈B
B≠A
−
∑ [Pcdδ[a],[b]δ[c],[d](ab|cd) c ,d
∑ γABPB B≠A
∑ [Pcd(ab|cd) − Pcd ,σ(ad|bc)] c ,d
B≠A
= Uaa + (PB − Paa , σ )γAA +
(3.1.24)
− Pcd , σ(ac|ad)], Fab , σ = Uab +
(3.1.23)
∑
[Pcd(aa|cd)
c ,d∈A
a∈A
(3.1.28a)
[Pcd(ab|cd) − Pcd , σ(ac|bd)],
c ,d∈A
3.1.2. INDO and NDDO Methods. 3.1.2.1. Basic Formulation. Limited applicability to spin-polarized calculations has been the main difficulty of the CNDO and CNDO/2 models. Because the exchange integrals are completely neglected by the ZDO approximation, the energies of the singlet and triplet configurations are the same. The difficulty can be avoided in the improved approximation schemes: intermediate neglect of differential overlap (INDO)172 and neglect of diatomic differential overlap (NDDO).169 The approaches utilize most of the approximations made in the CNDO/2 method, but they retain a certain fraction of exchange-type integrals. In the NDDO approximation, the four-center integrals (ab|cd) are neglected except when the orbitals a and b are centered on the same atom A, and the orbitals c and d are centered on the atom B (which can be A):
a, b ∈ A , a ≠ b 0 Fab , σ = βAB Sab − Pab , σγAB ,
(3.1.28b)
a ∈ A , b ∈ B, A ≠ B (3.1.28c)
Because of the altered form of the Fock matrix, a slightly different way of parametrization of the core matrix elements is utilized. It takes into account the exchange integrals in the determination of IPs and EA. For instance, the INDO parametrization utilizes the Slater−Condon factors. In the NDDO approximation, the changed parametrization is related to an increased number of two-electron integrals. The integrals can be either computed directly, or treated as adjustable parameters, or both. The increased number of degrees of freedom in the semiempirical NDDO parametrization allows 5811
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
more flexibility and, hence, accuracy of the resulting models. At the same time, the transferability of the obtained parameters (especially if a small training set is used) is relatively low, and the parametrization procedure is more difficult. Historically, the initial version of the NDDO approximation showed somewhat poor performance. In later developments, both the NDDO and INDO Hamiltonians were used as the starting grounds for more accurate and transferable methodologies that will be outlined below. 3.1.3. MINDO, MNDO, AM1, and PMn Methods. The INDO approximation has been extended in a series of modifications, known as the modif ied INDO (MINDO) method.173−175 The following modifications were used in the first version:173 (a) The parameters were fit to experimental data, rather than to exact Hartree−Fock calculations. (b) The core matrix elements were fit to experimental atomic spectra, using transitions among high-spin electronic configurations. (c) An approximate functional form was utilized for estimation of the two-center integrals. The total energy is expressed as the sum of HF-like (electronic) energy, Eel, eq 2.1.28, and terms that represent the core−core interaction, Ecore,AB: Etot = Eel +
∑ Ecore, AB
JAB (RAB) =
JAA = IPA − EAA =
(3.1.30)
RAB2 + (ρA , l + ρB , l )2 2
f2 (RAB) = exp( −αRAB) (3.1.31b)
2
(3.1.34a)
2f JA , l + JB, l 1
(3.1.31c)
2
with an adjustable parameter f = 1.2. (3) The Ohno184 formula: 1
(
2 JAA + JBB
2
)
(3.1.34b)
The resulting MINDO/2 approach was found to be much more accurate and consistent with respect to its predecessor, MINDO (also known as MINDO/1). It gives good estimates for bond lengths, heats of formation, and force constants simultaneously and for a larger set of hydrocarbons. Although the accuracy of the estimated heats of formation improved significantly, it was judged to be “a bit high for chemical purposes”. The method underestimates the strain energy in small rings and overestimates attraction of nonbonded hydrogens, leading to low rotational barriers. Improper description of lone pairs constituted another deficiency. As a consequence, the method had limited applicability to compounds containing the N and O atoms. Finally, the MINDO/3 method was developed.175 The basic methodology was similar to the previous strategies, with an
f
RAB2 +
(3.1.33b)
with the optimal choice for the function f 2(RAB) given by
1 JA , l + JB, l
JAB (RAB) =
(3.1.33a)
⎛Z Z ⎞ Ecore, AB = Eel, AB + ⎜ A B − Eel, AB⎟f2 (RAB) ⎝ RAB ⎠
(2) The Nishimoto−Mataga−Weiss183 formula:
RAB +
(3.1.32)
Balancing of the core−core repulsion terms was the second major modification developed in the MINDO/2. In particular, the core−core repulsion of two atoms, Ecore,AB, should be equal at dissociation to the electron−electron repulsion of the neutral atoms, Eel,AB. There should be no long-range electrostatic interactions between the neutral atoms. To obey this asymptotic property, the following expression has been adapted:
(3.1.31a)
1 1
JAB (RAB) =
∂ 2E ∂qA 2
f1 (R ij) = 1
where ρA,l1 = 1/2JAA and ρB,l2 = 1/2JBB are adjustable in general. An expression similar to eq 3.1.30 was also used to describe the electron−core attraction, VAB. The choice of the function f1(R) is motivated by the proper asymptotic behavior at R → 0 and R → +∞, as well as by computational efficiency. Several other approximations of the kind have been discussed in the literature,181 namely the following three formulas. (1) The Nishimoto−Mataga182 formula: RAB +
(3.1.31e)
with the distance-dependent function f1(Rij). Later, the unit function was found to be the optimal choice:
1 1
JAB (RAB) =
1
+ 2 JBB ek BRAB
Hij = B(Ii + Ij)Sijf1 (R ij)
where ZX, X = A, B, is the core charge, and γAB represents the screened Coulomb potential. The latter was approximated via the Dewar−Sabelli176−178−Klopman−Ohno179,180 expression: γAB = JAB (RAB) =
1 1 J ekARAB 2 AA
The methodology was applied primarily to hydrocarbon molecules and radicals, and was successful in describing their heats of formation. However, the method could not describe bond lengths, and it failed to predict correctly the relative energies of rotational isomers. These shortcomings were addressed in the second version of the MINDO approximation (MINDO/2)174 by developing parametric functions for the core resonance integrals, Hij, and the core−core repulsion terms. Both types of functions depend on the interatomic distance and are chosen to satisfy suitable physical requirements. In particular, the resonance integrals are taken proportional to the AO overlap integrals, Sij, and to the valence-state ionization potentials of orbitals, Ii, similarly to the EHT and it extensions:
The MINDO expression for the core−core repulsion utilized an approximate functional form: Ecore, AB = ZAZBγAB
RAB +
with the adjustable parameter 0.4 ≤ k ≤ 0.8. The parameters JAA describe repulsion of two electrons on the same atom (one-center Coulomb integrals). They are called idempotential or self-Coulomb integrals.186 Having the meaning of the second derivative of energy with respect to atomic charge, they can be evaluated using atomic (or valence state, for orbitalresolved quantities) IP and EA:
(3.1.29)
A 0, γ > 0 (3.3.33)
where β and γ are empirical parameters, and S2(rij) is yet another switching function, typically with a shorter range than S1(rik). S2(rij) is designed to describe charge back-donation from nearby C atoms, to model a charge alternation pattern typical for aromatic systems. Because of the specifics of the interactions in the C60/ Au(surface) system, the explicit covalent bonding potentials similar to those used in reactive force fields were not constructed. Instead, the nearly covalent bonding in this system was explained on the basis of combination of the dispersion and dynamic electrostatic contributions:
(3.3.29)
The split charge qij can be interpreted as a steady-state flux (charge current) between the pair of atoms i and j. Similar to bond orders, it contains information about a pair of atoms rather
E = EvdW + Eelec 5831
(3.3.34a) DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews Eelec =
∑ i ∈ real
⎛ 1 2⎟⎞ ⎜χ q + Jq + ⎝i i 2 ii i ⎠
Review
∑ i , j ∈ real j≠i
qiqj rij
qiqj̃
∑
+
i ∈ real j ∈ image
interfacial systems. It introduced a number of corrections, including new repulsion and bond-bending terms, and overcoordination corrections. It also replaced the point-charge models by the Coulomb integrals over 1s Slater functions. The electrostatic cutoff functions were abandoned in favor of the charge-neutralized real-space direct summation method.136 The second version (COMB2B)21,22 was parametrized against the training set that contained both experimental and high-level ab initio data on small molecules, solid-state systems, and interfacial systems. The data included cohesive energies, formation enthalpies, elastic properties, surface energies, and bond dissociation energies of molecular and anionic species. In addition to the previous formulations, the model included the point dipole terms to improve description of the atomic polarization effects. The third generation of the COMB potential (COMB3) was developed very recently.23 The potential was improved by modifying the expressions for bond order and self-energy. Flexibility of the functional form and a judicious parametrization strategy allowed application of the model to carbon-based materials, hydrocarbons, organometallic compounds, and combinations of these systems. A large-scale (∼5000 atoms) molecular dynamics simulation was reported. 3.3.4. Construction of Reactive Bond-Order Potentials as a Phenomenological Variational Principle. Construction of a suitable functional form that obeys the basic physical considerations and limiting properties, and has good transferability, is certainly the art of computational chemistry. Some aspects of the analytic potential design were extensively discussed by Brenner.86 Although a specific choice is rather arbitrary, it is always motivated by the problems in hand and by the basic chemical concepts, such as directional bonding and bond order.
rij̃
(3.3.34b)
EvdW =
∑ i∈C j ∈ Au
12 ⎡⎛ ⎛σ ⎞6 ⎤ σC − Au ⎞ ⎢ C − Au ⎟ ⎥ ⎜ ⎟ ⎜ − 2⎜ DC − Au⎢⎜ ⎟⎥ r ⎟ ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠
(3.3.34c)
where the meanings of the χi and Jii parameters are similar to those in the Qeq theory. DC−Au and σC−Au are the standard vdW parameters, qi is the partial charge of the atom i computed according to eq 3.3.33, q̃i = −qi is the image charge induced by the charge of the atom i, rij is the distance between two adsorbate atoms, and r̃ij is the distance between the adsorbate atom i and the image of the atom j. The different types of interactions are illustrated in Figure 4. The parameters are designed to reproduce charge transfer between C60 and the Au surface, and the chemisorption energies at different C60/Au geometries, simultaneously. This is done by fitting to extensive semiempirical calculations and to known experimental data. Because the model uses simple parametrized energy/charge functionals, it is efficient computationally, allowing for nanosecond MD simulations with more than 1000 atoms. The ability to reproduce coadsorption effects is an important qualitative advantage of the above model. The construction given by eqs 3.3.32−3.3.34 can be regarded from the mathematical point of view as a reactive bond-order potential, since the total energy is ultimately expressed in terms of bond-order variables. However, the physical interpretation is differentthe model intrinsically couples electrostatic and covalent bonding contributions. The charges are computed via bond orders and can be used in calculations employing an external electric field.39 The response to an external electric field is absent in the early REBO, BOC-MP/UBI-QEP, and EAM/ MEAM potentials. Note that a simple addition of the QEq/EEM scheme in the ReaxFF to the electrostatic energy does not directly couple charge transfer and bond orders. This coupling is present in the qAIREBO force field, although the strategy is different. 3.3.3.9. COMBx (x = 1, 2A, 2B, 3). A series of charge optimized many-body (COMB) potentials was developed recently by the Sinnott group.20−23,439,440 The method originates from the extended Tersoff potential developed by Yasukawa.441 The latter generalized the Tersoff potential to make the repulsion and attraction terms charge dependent. The charges are determined with a charge equilibration principle, although the functional form of the charge-dependent energy is different from that of Goddard. The dynamic-charge scheme was essential for proper description of the charge redistribution in the interfacial Si/SiO2 system. The initial version of the COMB potential (COMB1)20 was mostly a reparametrization of the Yasukawa potential. It also included additional repulsion and bond-bending terms. The extended Lagrangian approach was adapted to solve for coordinate-dependent charges during molecular dynamics evolution. The atomic charges are treated as dynamical variables in this method. The approach increases the number of equations of motion, but avoids solving systems of linear equations, and leads to linear scaling of the required CPU time. The second generation of the COMB potentials included different parametrization sets and focused on new types of systems. The COMB2A439,440 was optimized for bulk and
Figure 4. Interactions in the coupled charge-bond-order scheme. (a) Auxiliary charge density at the central atom C is given by summation of the contributions from all Au atoms (red arrows). The charge on this atom is determined by linear combination of the auxiliary densities of all adsorbate atoms. (Blue arrows indicate possible contributions to/from the nearby partial charges.) (b) Adsorbate point charges (purple circles) induce image charges (blue circles). The energy is determined by summing electrostatic interactions between adsorbate−adsorbate and adsorbate−image atom pairs. q1 and q2 are partial charges of the adsorbate atoms 1 and 2; d1 and d2 are distances from the atoms 1 and 2 to the surface plane.
Over the course of development, the functional form of the reactive potentials gradually became more and more complicated, and involved a large number of adjustable parameters. Some methods evolved as a hierarchy of improvements built up step-by-step to add new physical effects or to include more general bonding patterns. For example, one can observe clear 5832
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
evolutionary lines, such as Tersoff → REBO → AIREBO → qAIREBO or Tersoff → Yasukawa → COMB → COMB2(A,B) → COMB3. These lines can be thought of as minimization paths in the functional space spanned by all functions approximating the true universal energy functional of an arbitrary system. In other words, the process of developing reactive potentials is essentially an execution of variational principle at the meta level. Therefore, the potential development can be called the empirical (phenomenological) variational principle. The central challenge of the DFT is to determine the universal density functional, which can be used to compute the exact energy of an arbitrary many-body system as a function of its electronic density. This problem can be reduced to finding a suitable exchange−correlation functional, and perhaps a kineticenergy functional (as in the OF-DFT). In all cases, the functional is constructed to represent the fundamental electronic interaction principles. The energy can then be found by optimizing the electron density with the proper functional. The construction of the BO potentials starts from the total energy expression and includes functional variation. The electronic density is not directly involved and can be arbitrary, fixedthis can be rationalized by an appropriate functional transformation. The variation of the functional is aimed to arrive at a proper model representing the reality. The relationship between the DFT and the phenomenological variational principle is illustrated in Figure 5. The efforts in the development of reactive force fields are often directed toward description of specific elements and binding environments. They are designed to model particular processes or systems. A generalization to a large database of properties and to a larger scope of elements and systems is often very involved and complex, and is frequently absent. The empirical reactive force field approaches are somewhat less attractive from this point of view than the semiempirical methods, which are based on less drastic approximations and retain more information, including the electronic structure properties. This deficiency of the reactive potentials is compensated by a much higher computational efficiency, which makes simulation of large-scale systems possible when the semiempirical methods become too time-consuming. The steadily growing complexity of the REBO potentials make them appear somewhat “black box”. Although all the terms are physically motivated, their particular choice is flexible and is rather arbitrary. The large number of these terms and their inhomogeneity make further systematic improvements and extensions to more compounds and elements rather difficult. In this regard, it may be rewarding to reformulate the reactive potentials, or at least some of their parts using artificial neural networks (ANN). A properly constructed and trained multilayer ANN is capable of approximating an arbitrary function of desired complexity. An ANN can also “memorize” large amounts of information on different types and use it for further predictions. Both input and output of ANNs can include multiple channels, which may be used systematically for training on inhomogeneous data sources, e.g., elastic constants, heats of formation, optimized geometrical parameters, spectroscopic data, etc. Although an ANN carries little physical interpretation, it can be used in conjunction with physically motivated terms to hide complex correlation effects. For example, one of ANN outputs may return the ratio between the repulsion and attraction terms with a predefined asymptotic, the exponential parameters entering the bond-order calculations, the values of the damping and switching functions, the factors
encoding coordination (e.g., hybridization, over/undercoordination), etc. In the past, ANNs were utilized in relatively few, but very interesting chemically relevant applications, to approximate density functionals,442 potential energy surfaces and molecular dynamics,443−449 transition moment functions,450 dynamical polarizabilities,451 and molecular electrostatic potentials.452 The use and development of the approach is inhibited somewhat by a lack of clear interpretation of the parameters and functional forms (the so-called excitation function of a neuron). We envision that the growing complexity of the functional form of the reactive force fields can facilitate a conceptually simple, systematically improvable, and computationally efficient ANN formulation. In a similar fashion, the ANN formalism may be very useful in the semiempirical methods, especially for incorporating the correlation and complex many-body effects, yet leaving the computational advantages of the basic methods unchanged. 3.3.5. Timeline of Bond Order Based Methods. Table 3 includes the timeline, from 1984 to 2012, of the development of bond order based methods.
4. COMPUTATIONALLY MOTIVATED APPROXIMATIONS Development of this group of methods started when the physically motivated approaches reached practical limits. Although some physically motivated approximations are very efficient, they often lack the quantum-mechanical level of description. When quantum calculations must be performed on large-scale systems, an alternative approach is needed. The computationally motivated approximations make use of the advantages of alternative mathematical reformulations of the existing problem. Fragmentation is a typical strategy. It involves an initial subdivision of the whole system into smaller subsystems, solving the smaller problems, and subsequent reconstruction of the properties and variables of the combined system. These methods are currently known as the linear-scaling methods of quantum chemistry. Another related group includes iterative (direct optimization) methods, which also show linear scaling but avoid diagonalization of large matrices by using iterative procedures, optimization algorithms, and sparse-matrix techniques to solve the Schrödinger equation.
Figure 5. Relationship between the DFT and the phenomenological variational principle. f [q,x] is an effective functional for the bond-order potential, mapping atomic charges, q, and coordinates, x, onto the total energy. This function is a BOP counterpart of the universal functional, F[ρ], that maps the charge density ρ onto the total energy in DFT. 5833
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
The computationally motivated approaches follow the ideas of parallelismall smaller subproblems can be treated simultaneously, in parallel. On the contrary, the philosophy of the physically motivated approximations is essentially serialthe entire problem is treated as the whole, only using more efficient formulas. (This does not preclude parallelization of realistic calculations based on these methods, of course.) The difference in the philosophies of the two groups of methods also reflects the state of computational technologies available at times of their development. Parallel computations were practically not available in the 1960s−1980s, when most of the basic physically motivated approaches were developed. Parallel computations emerged starting in the 1990s, when the computationally motivated (linear-scaling) methods reached their active development phase. Parallelization facilitated orientation of novel computational methods of quantum chemistry toward a fragmentation-based, divide-and-conquer philosophy.
Table 3. Timeline of Bond Order Based Methods Development year
event
authors
1984
universal scaling of cohesive energies BOC-MP proposed simple many-body reactive potential (Stillinger−Weber potential) three-body potential for Si general reactive potential based on the many-body bond-ordered terms reactive potentials for atomic scattering in BO space EEM bond-order term for covalently bonded materials (Tersoff or Abell−Tersoff potential) UBI-QEP BO potential for hydrocarbons (Brenner, Tersoff−Brenner, or REBO potential) Qeq NBI-MD classical fluctuating charge model adaptive intermolecular REBO (AIREBO)accounting for LJ interactions ReaxFF for hydrocarbons maximal entropy valence bond (MEVB) method REBO2improved accuracy for hydrocarbons ReaxFF extended to transition metals split charge equilibration (SQE) charge transfer with polarization current equalization (QTPIE) COMB1 coupling of the split-charge equilibration scheme with with BOPa prototype of qAIREBO COMB2A COMB2B qAIREBO COMB3
Rose et al.336 Shustorovich380 Stillinger and Weber400
1985
1986 1988 1990
4.1. Classification of Computationally Motivated Approaches
1991 1994
In this section, we discuss the linear-scaling techniques of quantum methods. We use the computationally motivated approaches synonymously with the linear-scaling method in this context. To start our review of this topic, it is helpful to give a broader look at the existing approaches and to classify the methods. A detailed description of the methods and their acronyms will be presented after that. The wide variety of existing linear-scaling methods can be classified according to different principles. The following are the most popular among existing classifications: (1) Li et al.454 classified the methods into two main groups, depending on the principle that guides fragmentation: (a) density matrix based methods (D&C, ELG, MFCC)
2000 2001
2002 2005 2006 2007
(b) energy-based approaches (MFCC, MTA, SFM) In the first group of methods, the main goal is to approximate the density matrix of the whole system by the density matrices of the fragments. The former is not explicitly constructed, but approximations are used to compute energy and other molecular properties from the densities (density matrices) of the subsystems. The methods from the second group access the energy of the whole system directly, by representing it with a sum of energies of the fragments and additional nonadditivity terms. (2) Fedorov and Kitaura455 suggested three main categories: (a) divide-and-conquer methods (D&C)
2009 2010 2011 2012
Biswas and Hamann401 Abell394 Garcia and Laganà396 Mortier et al.418 Tersoff453 Shustorovich382 Brenner408 Rappe and Goddard186 Sellers391 Rick et al.4 Stuart et al.412 van Duin et al.28 ́ 426 Morales and Martinez Brenner et al.409 Nielson et al.29 Nistor et al.437 ́ 427 Chen and Martinez Yu et al.20 Mikulski et al.428 Shan et al.439,440 Devine et al.21,22 Knippenberg et al.429 Liang et al.23
• one step
(b) many-body molecular interaction methods (FMO, MFCC, etc.) (c) transferable approaches (SFM) The principle of classification is quite similar to that of Li et al.454 Additional complexity arises from the differentiation of the categories b and c. The difference lies in the way the interactions between fragments are added. (3) Suárez et al.456 classified the many-body expansion (MBE) methodologies into two groups, according to the topology of fragmentation: (a) overlapping fragments (MTA, MFCC, FEM, SFM)
• two step The one-step group essentially substitutes the energy-basedmethods group. It corresponds to the methods in which a property of the composed system (e.g., energy) is directly computed from the properties of the fragments. The two-step group substitutes the density-matrix-based methods group of the Li et al. classification. In these methods, an intermediate property of the whole system is first computed from the corresponding pieces of the fragments (e.g., density matrix). The intermediate property is then used to compute the required properties (e.g., energy) of the entire system. At the second level of the classification hierarchy, each of the first-level groups is subdivided into three groups: • one body
(b) disjoint fragments (FMO, KEM) (4) Gordon et al.76 recognized the difficulties of the former classification schemes and the need for a more complex one. They generalized the scheme of Li et al.454 by a two-level classification scheme. At the first level, a method can be classified as
• many body • conglomerate As the number of new methods for computations on largescale systems grows, the distinctions become more indefinite, fuzzy, and a straight linear classification becomes problematic. This problem is already seen in the classification of Li et al., in 5834
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
which a single method (e.g., MFCC) can be classified into more than one group. The classification scheme of Gordon et al. partially solves this ambiguity; however, different versions of the same family of methods (e.g., FMO) can eventually be classified into different groups, which may become distinct at higher levels of hierarchy. In our opinion, a modern successful classification scheme cannot be linear. Instead, it is best represented by an interconnected tree (Figure 6). The tree mimics the history of the methods development, and their branching and combination. If in the beginning of the era of large-scale calculations the distinction was very clear, as progress continues, new methods appear by combination of ideas from different groups of the earlier methods. The history process acts as a variational, genetic optimization procedure, pushing researchers to combine the most successful ideas and approaches found in different aspects of the large-scale computations. Less successful techniques eventually lose popularity. Therefore, it can be expected in the near future that the broad variety of existing approaches may merge into a few most optimal methods, which will combine the traits of many developed techniques. For the purpose of presentation and organization of the material in this section, we will adopt the following classification scheme. All methods for quantum-mechanical calculations of large systems can be grouped into three large categories: • fragmentation-based methods • direct optimization methods • quantal force fields It is possible that techniques major for a given group of methods can also be used in another group. For example, direct optimization can be applied within the fragmentation-based schemes, such as D&C. Quantal force fields are derivatives of the energy-based or other fragmentation methods and, in principle, can be classified into that category. However, due their distinct features and potential diversity in the future, we prefer to separate them into an independent group, forgetting for a moment their true origin. This separation will also help to structure the presentation and provide the capability to start a new classification within the group. Since seemingly distinct methods may have large similarities in their foundation, or formulation, or both, it is often difficult to classify them into a single group. Instead, we choose to tag the fragmentation-based methods according to the following three types of properties: Property 1: • density partitioning and derived schemes (approximation of density matrix) • Hamiltonian partitioning derived schemes (approximation of Hamiltonian) Here we essentially follow the original classification of Li et al., but generalizing energy-based partitioning to the Hamiltonianbased partitioning. Property 2: • static partitioning (working with the whole pool of fragments at once) • dynamic partitioning (working only with a subset of fragments at given time) The methods of the first type utilize a predefined set of fragments, all at once. Most methods can be classified into this category. The dynamic partitioning schemes utilize only a subset
Figure 6. Tree structure of the existing linear-scaling methods. Offspring methods may be derivatives of some parent method. Additionally, an offspring method may be a derivative of more than one parent method. The direction of the arrow indicates gradual complication and specialization of the methods. Separation of the arrows along the xscale characterizes the similarity of the approaches. Initially distant methods spread out by acquiring new traits and can eventually overlap, because the new traits may be common across the newly developed methods. See text for abbreviations.
of all available fragments at each iteration. Some property is dynamically adjusted and truncated to maintain constant complexity. The ELG and Stewart’s LMO methods are the clearest examples of the methods in this category. Property 3: • adiabatic partitioning (fragments are close to what they are in the combined system) • diabatic partitioning (fragments are diabatic states and can be rather arbitrary) Most of the discussed methods are classified as adiabatic partitioning. As an example, conjugate caps are chosen in the MFCC method to represent the environment of a given fragment in a realistic system as accurately as possible. In the D&C scheme, one chooses a buffer or double-buffer region to minimize spurious effects at the boundary. In this sense, a judicious choice of fragmentation is very important for obtaining desirable accuracy. The methods of the second group are less popular, but are actively developed too. In the diabatic partitioning methods, spurious boundary effects are relatively unimportant. The same applies to an unphysical charge distribution over the fragments into which the whole system is partitioned. One seeks to construct a complete basis of diabatic configurations and to represent the whole system as their linear combination. 4.2. Fragmentation-Based Approaches
We start by presenting yet another practical classification. It is specific to fragment-based approaches and is designed for structuring the current presentation. We adopt the following classification: (1) density-based fragmentation (D&C, MEDLA, etc.) (2) energy-based fragmentation (FMO, MFCC, SMF, MTA, E-EDC) (3) dynamical growth with localization (LMO, ELG) (4) diabatic approaches (FMO-MS-RMD, SFS) Note that the three types of properties discussed in section 4.1 may be attributed to any group of methods in this classification. 5835
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Figure 7. Illustration of the charge density definition that avoids using KS orbitals (a) and the partitioning of the total charge density into two fragments Aand B, using the smooth switching functions defined by eq 4.2.4 (b). The switching functions in (b) are shifted along the Y coordinate for clarity. The positioning along the X coordinate corresponds to that in the combined system, outlined in the middle drawing. {pα,σ|α = A,B} is the density partitioning scheme for the spin channel σ, H is the Hamiltonian, EF,σ is the Fermi energy for the spin channel σ, and Nσis the number of electrons with a particular spin z-projection σ.
4.2.1. Density-Based Fragmentation. 4.2.1.1. Divideand-Conquer (D&C). One of the first linear-scaling approaches is the divide-and-conquer (D&C) technique. It was proposed originally by Yang68,69 within the DFT framework. Later it was extended to semiempirical70,72 and ab initio Hartree−Fock75 molecular orbital methods by Dixon, He, Merz, and many other authors. The MO-based versions of the D&C algorithm were extended to MP2 and CCSD levels of theory. The development of the D&C algorithm was preceded by a range of related methodologies, not necessarily based on self-consistent formulation. Below we summarize the history of the method development and describe the basics of the major techniques. 4.2.1.1.1. Charge Density Partitioning. The D&C method68,69 represents one of the first approaches to achieve linear scaling in quantum mechanical calculations. The central D&C idea is similar to that of any other density partition schemeto “abuse” the main assumption of the DFT, which utilizes the charge density as the main variable describing the state of the system, and to make a good use of its additivity: ρA ∪ B ( r ⃗) = ρA ( r ⃗) + ρB ( r ⃗)
KS expression for the total density, eq 2.2.11, for all electrons with the spin component, σ, using the Heaviside function: ρσ ( r ⃗) = ⟨r |η(E F, σ − Ĥ )|r ⟩
(4.2.2)
where ⎧1, η (x ) = ⎨ ⎩ 0,
x>0 x≤0
(4.2.3)
The definition eqs 4.2.2 and 4.2.3 can be interpreted as the sum of squares of the amplitudes of all eigenvectors of the Hamiltonian Ĥ with eigenvalues lower than the Fermi energy EF,σ, eqs 2.1.49 (Figure 7a). The electron density partitioning scheme is defined by a set of smooth switching (weighting, partitioning) functions: {pA , σ ( r ⃗):
∑ pA ,σ ( r ⃗) = 1} A
(4.2.4)
Each function pA,σ(r)⃗ = pA,σ(r;⃗ {R⃗ }) defines the region that encloses the fragment A in the space of Cartesian coordinates of all atoms, {R⃗ }. For example, the function can be 1 within the interior of the fragment, 0 outside of the fragment, and an intermediate value within some region around the fragment boundary (Figure 7b). The only requirement is that the sum of contributions from all fragments should be 1 at each coordinate vector r.⃗ Apart from this requirement, the partitioning functions, pA,σ(r)⃗ , can be chosen quite arbitrarily. A variant suggested by Yang69 is based on spherical atomic densities pa,σ(|r ⃗ − R⃗ a|):
(4.2.1)
where A ∪ B represents the union of the sets (fragments) A and B. The assumption eq 4.2.1 makes the partitioning of a complex system into fragments easy, appealing, and natural, in contrast to the wave function based composition/partitioning schemes. In contrast to charge densities, wave functions are nonadditive, in general. Thus, one needs to account for renormalization and boundary effects in the MO-based approaches, when reconstructing the wave function of the entire system from those of the fragment MOs. For these reasons, the MO-based D&C methodologies were developed later and utilize an object that is similar to the charge density in the DFT, namely, the density matrix, eqs 2.1.32 and 2.1.33, to perform partitioning. The density matrix can be considered the WFT counterpart of the charge density in the DFT. Therefore, it can be used naturally for constructing WFT-based D&C partitioning schemes. These methods will be described in the following sections. The main purpose of Yang’s D&C formulation of the standard DFT was to avoid using KS orbitals, and to work only in terms of electron density. The starting point is to reformulate the standard
pA , σ ( r ⃗ ; {R⃗}) =
gA , σ ( r ⃗ ; {R⃗}) ∑A gA , σ ( r ⃗ ; {R⃗})
(4.2.5)
with gA , σ ( r ⃗ ; {R⃗}) =
∑ ρa ,σ (| r ⃗ − RA⃗ |) a∈A
(4.2.6)
Here and in the following discussion of fragmentation/partition/ subsystem methods, we will use the capital letters A, B, C, etc. to denote fragments/sets/subsystems and small letters a, b, c, etc. to denote atomic species. 5836
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Using the switching functions eq 4.2.4, one can divide the total density, eq 4.2.2, into a sum of the fragment densities: ρσ ( r ⃗) = ⟨r |η(E F, σ
− Ĥ )|r ⟩ =
∑ pA ,σ ( r ⃗)⟨r|η(EF,σ
HÂ =
∑ ρA ,σ ( r ⃗)
The FMOs, {|ψA,σ,i⟩}, are determined by solving the set of coupled equations, eq 4.2.11. Because this procedure accounts for electronic density deformations of one fragment induced by all other fragments, we classify these FMOs as adiabatic, in analogy with the terminology utilized in the field of quantum dynamics. The opposite limit is to utilize the basis of FMOs obtained from calculation of noninteracting (isolated) fragments. Such approaches are utilized by a number of authors and will be discussed later. The FMOs can be classified as diabatic in this situation. Finally, many intermediate types of coupling between fragments can be included into the procedure used to determine FMOs. Common to the QM/MM methods, they will be discussed later. 4.2.1.1.3. Subsystem Equilibrium and D&C Algorithm. Combining eqs 4.2.7−4.2.9 and keeping in mind the localization approximations, eqs 4.2.10−4.2.11, the total density of the whole system can be computed by
(4.2.7)
A
with the charge density on the fragment A defined by ρA , σ ( r ⃗) = pA , σ ( r ⃗)⟨r |η(E F, σ − Ĥ )|r ⟩
(4.2.8)
4.2.1.1.2. Hamiltonian Localization and Fragment MOs (FMOs). The evaluation of ⟨r|η(EF,σ − Ĥ )r⟩ requires knowledge of the Hamiltonian of the whole system, which depends on the densities of both spin-up and spin-down electrons. In order for the D&C approach to be useful, the second main assumption is made. Namely, the global Hamiltonian, Ĥ , is approximated by a set of local Hamiltonians for subsystems, Ĥ A. In addition, for practical purposes, a smooth equivalent of the Heaviside function, eq 4.2.3, is introduced. With these approximations, eq 4.2.8 turns into a practically more appealing expression: ρA , σ ( r ⃗) = 2pA , σ (r )⟨r |fΔE (HÂ ; E F, σ )|r ⟩
ρσ ( r ⃗) =
∑
≈
∑
The Fermi energy, EF,σ, is assumed to be the same for all subsystems (fragments), which is consistent with its physical meaning as the chemical potential. The self-consistent solution corresponds to the equilibrium state, namely to the energy minimum, which follows from the variational principle. Therefore, all subsystems must be in equilibrium with each other. In other words, the self-consistent solution is obtained when the Fermi energies of all subsystems are the same. The numerical value is obtained from the normalization:
i,j
Nσ =
|χA , a ⟩[∑ (SA−1)ai HA , ij(SA−1)jb ]⟨χA , b |
a,b∈A
(4.2.10a)
i,j
=
∫ ρσ ( r ⃗) d r ⃗ ∑ ∑ fΔ (EF,σ − EA ,σ ,i) ∫ pA ,σ ( r ⃗)|ψA ,σ ,i( r ⃗)|2 d r ⃗ A
and (4.2.10b)
The overall D&C scheme is iterative, since the total density is computed by summing up densities of all subsystems, and the density of each subsystem depends on the Hamiltonian of the whole system. The principal scheme of the D&C algorithm is illustrated in Figure 8. 4.2.1.1.4. Source of Computational Savings. Unlike the conventional DFT, in which the KS Hamiltonian is represented by a single NM × NM matrix, the dimensionality of the fragment KS Hamiltonian is reduced to NMA × NMA, where NMA is the number of basis states (and, equivalently, MOs) associated with the fragment A. The number of the SCF equations one needs to solve increases and is equal to the number of subsystems, NS. However, this increase is only linear with NS. The reduced dimensionality of the KS Hamiltonian that must be diagonalized has a much more drastic effect on performance, reducing the CPU time cubically. In total, the scaling of computational s resources reduces to O(∏NA=1 NMA3). Assuming that sizes of all subsystems are comparable to each other (even partitioning), NMA = NM/NS, ∀ A = 1, ..., NS, one obtains O(NS(NM/NS)3). Choosing NS proportional to the system size, NS = NM/α, one obtains a linear-scaling performance: O(α2NM) = O(NM). 4.2.1.1.5. Anatomy of Fragmentation and Terminology. In his original DFT-based D&C scheme, Yang68,69 introduced the
The first double summation in eq 4.2.10a runs over all AOs localized on the fragment A. It is important to emphasize that eq 4.2.10b defines matrix elements of the f ull Hamiltonian. It depends on properties of all fragments/subsystems, therefore leading to a set of coupled SCF (KS) equations of eq 2.1.38 type, but with additional indexing due to fragments: HA , σCA , σ = SACA , σEA , σ
(4.2.11)
where HA,σ is the KS Hamiltonian, eq 2.2.11, of the fragment A, and SA is the overlap matrix of the AOs localized on the fragment A: ⟨χA , i |χA , j ⟩ = SA , ij
(4.2.12)
The matrix CA,σ contains coefficients that transform the set of nonorthogonal fragment AOs, {|χA,i⟩}, to the set of orthogonal fragment-localized MOs (also known as fragment MOs, FMOs), {|ψA,σ,i⟩}:
∑ CA ,σ ,ai|χa ⟩, a
i
(4.2.16)
HA , ab ≡ ⟨χA , i |Ĥ |χA , j ⟩
|ψA , σ , i⟩ =
i
(4.2.15)
|χA , a ⟩[∑ (S −1)ai Hij(S −1)jb ]⟨χA , b |
a,b∈A
∑ pA ,σ ( r ⃗) ∑ fΔ (EF,σ − EA ,σ ,i)|ψA ,σ ,i( r ⃗)|2 A
(4.2.9)
where fΔE(Ĥ A; EF,σ) is the Fermi function, eq 2.1.49b, of the operator argument. A definition of the subsystem Hamiltonians, Ĥ A, is the last ingredient needed to complete the D&C scheme. The subsystem Hamiltonian Ĥ A is defined as the projection of the full KS Hamiltonian on the subspace spanned by the basis orbitals (AOs) that are localized on the corresponding subsystem A, {|χA,i⟩}. Using eq 2.1.42, this can be written formally as HÂ =
(4.2.14)
i,σ
− Ĥ )|r ⟩
A
=
∑ |ψA , σ ,i⟩EA , σ ,i⟨ψA , σ ,i|
σ = α, β (4.2.13)
The projected Hamiltonian, eqs 4.2.10, simplifies significantly in the FMO basis: 5837
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Figure 8. Algorithm of the D&C method based on charge density partitioning in DFT. The arrow that closes the loop represents the step in which verification for self-consistency among all subsystems is performed, as judged on the basis of computed Fermi levels of all subsystems. See text for detailed discussion of the terms in the scheme on the right.
Figure 9. Illustration of partitioning in the D&C scheme. (a) Subdivision of the system on the fragments (subsystems) and their componentscore and buffer regions. (b) Weights of different types of matrix elements in the density matrix of a given fragment. {φi} are basis functions localized in different spatial regions.
concept of buffer atoms. According to Yang’s definition, these atoms are not part of the subsystem (in our further terminology called “core region”), but they are in its neighborhood and contribute to the subsystem orbitals. The anatomy of the partitioning of a molecular system into fragments is illustrated in Figure 9. The molecular system is composed of generally overlapping, localized regions, also called fragments or subsystems, FA:
(4.2.18a)
CA ∩ BA = ϕ
(4.2.18b)
Division into the core and buffer sets must satisfy the condition that the core regions of different fragments do not overlap, for any pair of distinct fragments: CA ∩ CB = ϕ
NF
M = ∪ FA
(4.2.17a)
FA ∩ FB ≠ ϕ
(4.2.17b)
A=1
FA = CA ∪ BA
(4.2.19)
The overlap may be nonzero in some definitions.70 This deviation from the definition of the core region, eqs 4.2.18−4.2.19, can be resolved by introducing additional partitioning to separate nonoverlapping core sets that obey eq 4.2.19.72 The sets FI, CI, and BI in the partitions above can be understood as the sets of atoms, in the spirit of the original Yang68,69 terminology. It is more general and convenient for further discussion to think about them as sets of localized
Here, NF is the number of fragments into which the system is subdivided, M indicates the entire molecular system, and ϕ denotes the empty (null) set. Each localized region, FA, is composed of the nonoverlapping core, CA, and buffer, BA, regions: 5838
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Figure 10. Illustration of the partitioning scheme in the space of orbital indices with a three-orbital example: (a) original Yang scheme; (b) improved scheme by Merz. The construction of the total density matrix from the density matrices of the fragments is also illustrated. Panel a shows partitioning of the orbital space M = {1, 2, 3} into two fragments, F1 = {1, 2} and F2 = {1, 2, 3}. (This example only illustrates the principle, and does not provide much computational saving.) The first fragment is composed of the core, C1 = {1}, and buffer, B1 = {2}, regions. The second fragment is composed of the core, C2 = {2, 3}, and buffer, B2 = {1}, regions. The core regions of the two fragments do not overlap. The partitioning in (b) is F1 = {1, 2, 3} and F2 = {1, 2, 3} with C1 = {1, 2}, B1 = {3}, and C2 = {2, 3}, B2 = {1}. The overlap of the cores of the two subsystems in this case is not empty: C1 ∩ C2 = {1, 2} ∩ {2, 3} = {2} ≠ ϕ. Pij, Pij′, and Pij″ are the elements of the density matrix. Prime and double prime indicate fragment matrices; unprimed quantities correspond to the reconstructed density of the entire system.
regions of a given fragment, eq 4.2.18b, there is no additional restriction on overlaps of the partitions. In particular, the core of one fragment can overlap with the buffer of another fragment:
orbitals. Most often, the bijection between the two sets can be established: Satomic, A = {I : i ∈ Sorbital, A , i ∈ I }
(4.2.20a)
CA ∩ BA ′ ≠ ϕ
and inversely Sorbital, A = {i: I ∈ Satomic, A , i ∈ I }
(4.2.21a)
Moreover, the core of one fragment can be the buffer of another, as illustrated in Figure 9a:
(4.2.20b)
CA = BA ′
where Satomic,A and Sorbital,A are the sets of atomic and orbital indices of the fragment A, respectively, I is the atom index, i is the orbital index, and the notation i ∈ I implies that the orbital i is localized on the atom I. A unique mapping from the atomic set to the orbital set may not be possible in a more general case. This is why it is more general to think of the sets in the above partitioning in terms of orbital indices (essentially Hilbert subspaces), not just the sets of atomic fragments. It is because we think about the partitions in terms of sets of orbitals rather than atoms that we call the core what Yang originally called a subsystem. Our definition goes along the lines used later by other authors.72 We hope this definition, together with the above partitioning anatomy, can help to reduce confusion and possible ambiguity and lead to more accurate terminology. The main purpose of the buffer atoms is to supplement AOs at the subsystem boundary. This requirement stems from the approximation in eq 4.2.10a. Namely, the global Hamiltonian is defined on the Hilbert space of all AOs of the entire system. However, because of the localization approximation, eq 4.2.10a, which essentially disregards overlaps of AOs belonging to different fragments, ⟨χA,i|χB,j⟩ = δABSij, the fragment-localized Hamiltonian defined on the space of only those AOs that belong to a given fragment is not coupled to Hamiltonians of all other fragments. In other words, it is defined on one of the mutually orthogonal Hilbert subspaces of the original Hilbert space of the entire system. In order for the subsystems to be coupled, their orbital spaces should overlap, which is represented by eq 4.2.17b. Apart from the assumptions of nonoverlapping cores, eq 4.2.19, and the absence of overlaps between core and buffer
(4.2.21b)
4.2.1.1.6. Density Matrix Partitioning. The original formulation of the D&C was based on division of the system in the Euclidian space. The method operated directly on the charge density. In other words, the charge density was the partitioned quantity. The need for computing three-dimensional integrals of the type ∫ pA,σ(r)⃗ |ψA,σ,i(r)⃗ |2 dr ⃗ in eq 4.2.16 was one of the drawbacks of such an approach. In addition, the formulation was restricted to the DFT. In the refined formulation, Yang and Lee457 considered partition in the space of basis orbitals, making the density matrix the main partitioned object, eqs 2.1.32 and 2.1.33. The partition scheme, eq 4.2.4, becomes independent of the coordinate system and is substituted by a discrete analogue: {pA , σ , ij :
∑ pA ,σ ,ij
= 1, ∀ i , j ; σ ∈ {α , β}} (4.2.22)
α
The set of the piecewise functions depending on atomic orbital indices provides the simplest choice of partition functions:
pA , σ , ij
⎧1 i ∈ CA , j ∈ CA ⎪ = ⎨1/2 i ∈ CA , j ∈ BA or j ∈ CA , i ∈ BA ⎪ i ∉ CA , j ∉ CA ⎩0 (4.2.23a)
where CA denotes the core region of the fragment A, and i ∈ CA indicates that the orbital i belongs to the set CA. The orbital-based partitioning given by eq 4.2.23a is schematically illustrated in Figure 10a. 5839
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
In a later work, Dixon and Merz70 proposed an improved D&C density matrix partitioning scheme: pA , σ , ij
⎧ i ∈ BA j ∈ BA ⎪0 =⎨ ⎪ ⎩1/nij otherwise
Nσ
PA , σ , ij = pA , σ , ij
∑ fΔ (FF ,σ − EA ,σ ,k)CA*,σ ,ikCA ,σ ,jk k=1
(4.2.24)
(4.2.23b)
where BA is the buffer region of fragment A, and nij is the number of subsystems in which both basis functions i and j appear as nonbuffer functions. The scheme eq 4.2.23b was improved further72 by splitting the overlapping “core” parts into the true nonoverlapping core set (see our definition according to eqs 4.2.18−4.2.19) and the inner buffer (first buffer layer, buffer 1). The original buffer became the outer buffer (second buffer layer, buffer 2), Figure 10a:
pA , σ , ij
⎧1/n i , j ∈ C or i ∈ C , j ∈ B or j ∈ C , i A A 1, A A ⎪ ij =⎨ ∈ B1, A ⎪ ⎩0 otherwise
Figure 11. (a) Double buffer partitioning scheme by Merz72 (b) Illustration of the success (blue double arrow) and failure (red double arrow) of the scheme eq 4.2.23a to describe partitioning of the weights across the subsystems of different sizes. See text for discussion of system partitioning into different regions.
(4.2.23c)
Treatment of the interface between the subsystems is the main distinguishing feature of the new formulation. The partitioning eq 4.2.23b requires an overlap of the subsystems in their nonbuffer atoms (Figure 10b). The partitioning eq 4.2.23c is more consistent with the previous definition of the core regions, eqs 4.2.18−4.2.19. The core region contains only a high-quality density matrix. The inner buffer serves the purpose of coupling different subsystems, and in general, contains a less accurate density matrix than the core region. The outer region serves only the purpose of insulating the fragment from all other subsystems. The orbitals of that region are used in the SCF equations, but are not used for the density matrix refinement. The outer buffer in the partitioning scheme eq 4.2.23c corresponds to the standard buffer of the scheme eq 4.2.23b. The union B1,A ∪ CA of eq 4.2.23c corresponds to the core region of eq 4.2.23b. Finally, the union of the two buffer regions of eq 4.2.23c, B1,A ∪ B2,A, corresponds to the buffer of the scheme eq 4.2.23a. Utilization of the orbital-dependent normalization factor, nij, instead of the fixed factor 1/2 used in the original Yang scheme, eq 4.2.23a, presents another important modification. This approach is important if the system is split into fragments of notably different size, say fragments A and B, such that the core of the fragment A overlaps with the buffer region of the fragment B differently from how the core of the fragment B overlaps with the buffer of the fragment A (Figure 11a). According to the partitioning eq 4.2.23a, the partition factor for the density matrix element corresponding to the pair of orbitals i ∈ CA and i ∈ CB will be the same for each of the two fragments, pA,ij = pB,ij = 1/2, no matter how the orbitals relate to the sets BA and BB. In reality, equal weights of 1/2 are obtained only if the orbitals are symmetrically distributed (Figure 11a, blue double arrow), that is, i ∈ CA, i ∈ BB, and j ∈ CB, j ∈ BA. Both i and j belong to the two subsystems in this case. The partitioning factor is given by the normalization coefficient 1/nij = 1/2, which is the same as in the original formula, eq 4.2.23a. In a more complex situation, the orbital i ∈ CA may be outside the buffer of the fragment B, in which the partition of the small fragment gets zero weight, pB,ij = 0. From the point of view of the subsystem A, both orbitals still belong to a single fragment, and eq 4.2.23c yields pA,ij = 1. Thus, although the sum of the partitioning factors remains constant, the distribution of the weight across the subsystems is changed.
where {(EA,σ, CA,σ)} are the eigenvalues and corresponding eigenvectors of the interacting fragment, determined from the solution of eq 4.2.11. The total density matrix is given by the sum of the fragmental density matrices, as illustrated in Figure 10:
Pσ =
∑ PA,σ
(4.2.25)
A
The Fermi energy is obtained in a manner similar to the charge density partitioning scheme, using the following normalization condition: Nσ =
∑ Pσ ,ijSij = ∑ PA ,σ ,ijSij i,j
A ,i,j
(4.2.26)
Note that eq 4.2.26 does not involve any numerical threedimensional integration, leading to efficiency improvement. The density matrix based partition scheme can be applied to both the DFT and WFT methods, ab initio or semiempirical. The orbital-based partitioning leads to the charge density less localized in the physical space than that derived from the direct space based partitioning. This drawback is emphasized as the orbital size increases (e.g., for diffuse orbitals). The difference between the two partitioning schemes is easy to understand in terms of the atomic limit. In the case of atomic subdivision, the partitioning schemes eqs 4.2.4 and 4.2.23 become the Hirshfeld458 and Mulliken327,459 population analyses, respectively. 4.2.1.1.7. Elaborations and Extensions of the D&C Method. The basic D&C methodology described in detail above was extended to a variety of other electronic structure methods. In particular, the density matrix formulation of the D&C method was applied to accelerate calculations based on semiempirical methods, to study biomolecules (the groups of W. Yang, D. York, K. Merz, G. Scuseria, etc.)70,72,79,80,460 and nanoscale materials, such as fullerene clusters.461 Other variations include the implementations within DFT (the St-Amant group),462−464 HF and HF/DFT hybrids (the groups of H. Nakai, K. Merz, etc.),75,465,466 MP2 (the Nakai group),467−469 and CC (the Nakai group).470,471 The approach was adapted to unrestricted formulations, including HF/DFT472 and MP2.473 Error sources of the method were analyzed on a model polypeptide system, 5840
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Figure 12. Illustration of the MEDLA. P is the density matrix. The charge density distributions are reprinted from ref 476. Copyright 1993 American Chemical Society.
the “parent” molecules is essentially equivalent to that of the buffer regions in the D&C scheme. Moreover, the “parent” molecules can be constructed automatically for each given system of interest, as outlined in the works of Exner and Mezey.479−481 Once the charge density (or density matrix) for the parent molecule is obtained, the fragment densities (density matrices) can be computed, using the inversion of the D&C density matrix partitioning scheme of Yang et al., eq 4.2.23a. Each fragment is equipped with an internal coordinate system for further transformation of charge densities. By performing computations on a set of parent molecules, a fragment density database can be constructed. Finally, the precomputed fragment charge densities are used to reconstruct a charge density of large molecules. This is done by translations and rotations of the building blocks to match their internal coordinate system with the position and orientation of the blocks in the corresponding fragments of the target molecule. Summation of the charge densities follows, again using eq 4.2.23a, to construct the density of the whole system from the fragment densities, as in the D&C scheme. Utilization of fuzzy charge density boundaries is one of the central principles of the MEDLA approach. This is in contrast to defining a strict atoms-in-molecule partitioning done by Bader. In fact, the MEDLA partitioning is a close analogue of the Mulliken population analysis, extended to fragments. We want to emphasize that, although the partitioning is fuzzy in terms of the charge density, it is discrete in terms of the density matrix partitioning. We have already encountered a similar discussion of the charge density versus density matrix partitioning and their similarities to the Bader versus Mulliken population analyses in the context of the D&C approach. The same opportunities are present in the MEDLA as well. Although the approach is conceptually appealing and is useful for quick reconstruction of the charge densities of large systems, and even for express estimation of the energy of the total system using the D&C DFT formulations, there are some notable limitations. First of all, the method lacks a self-consistent scheme to adjust fragment-in-a-molecule densities. This approach is justified in many cases by the short range of quantum interactions and the localized nature of chemical bonds. In this regard, the methodology is somewhat reminiscent of tight-binding methods, which often disregard long-range interactions. Therefore, one
indicating that a large buffer space cutoff leads to errors due to the loss of variational flexibility of MOs, while a small buffer cutoff restricts the capability of the subsystem MOs to be orthogonal to the true whole-system MOs.462 The D&C scheme was also implemented in parallel computing environments.464 As an example of the remarkable capabilities of the method, a semiempirical calculation on a truly large-scale protein system containing over 9000 atoms was reported.460 The calculations could be performed on a typical workstation, indicating that the method has a large potential for application to even larger systems. 4.2.1.1.8. Implementation. The D&C SCF and correlation methods have been implemented into the GAMESS package474 by the Nakai group. Alternative implementation is available in the DivCon 99 program and its subsequent versions, distributed by the Merz group.475 4.2.1.2. Additive Fuzzy Density Fragmentation (AFDF: MEDLA, ADMA). An alternative approach to partitioning of the charge density (or density matrix), similar in spirit to the D&C method described above, was developed independently by Walker and Mezey, to constitute their molecular electron density Lego approach (MEDLA).476−478 Unlike the self-consistent schemes of Yang et al., the method utilized a significantly simplified, but computationally more efficient scheme. The method was originally targeted on the reconstruction of the electron density for the whole system. Later extensions of the original MEDLA methodology allowed calculation of electrostatic potentials,479,480 energies,480 and dipole moments480 of large systems. The MEDLA abbreviation is often used interchangeably with the abbreviation ADMA, which stands for a more general term: adjustable density matrix approach. Both methods are classified into the group of the additive fuzzy density fragmentation (AFDF) techniques, first introduced by Mezey and co-workers. Several related approaches developed by other authors can also be classified into this group of methods. The MEDLA technique consists of several steps, outlined in Figure 12. The central idea of the method is to compute a set of charge densities of the rigid molecular fragments. This is done by calculations involving small molecules that contain the fragment of interest (so-called parent molecules). Because the patent molecules are typically small, the computations can be done at the ab initio level. It is also worth mentioning that the concept of 5841
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Figure 13. Illustration of the FA-ADMA principle. Molecular structures are adapted from ref 481. Copyright 2004 American Chemical Society.
can expect that systems exhibiting charge transfer and delocalized electronic densities may lead to deteriorated accuracy. Charge transfer and delocalization effects are of special importance in excited states, in particular. Thus, the MEDLA may be difficult to apply to excited states, and is mostly targeted to the ground electronic state. The difficulties may be partially accounted for by choosing sufficiently large fragments, which, however, may lead to loss of computation speed. Second, the construction of a fragment library assumes the utilization of fixed geometries for constituent fragments. Thus, the calculation of molecular dynamics or reactive processes is limited or may require a very large database of the charge densities for different geometries of considered fragments. On the contrary, when the intrafragment dynamics is expected to be negligible, the MEDLA may be very advantageous in terms of computational costs and accuracy. The method can be successful with macromolecular nonreactive systems, such as biopolymers and soft-matter systems, and with a range of materials science and biomedical applications. Finally, one more limitation is associated with the information storage difficulties. As the authors indicate, the molecular charge density is stored in cubic grid files. The memory requirement rises very quickly with the system size, rendering the method inapplicable for sufficiently large systems. The latter difficulty can be mitigated, if one stores density matrix information for fragments only. Then, one can utilize that information to compute properties of the large molecules on-the-fly, similar to the procedure used in the D&C schemes. Below we summarize the set of formulas for computation of some molecular properties within the ADMA.480 These are obtained by representing the total density matrix as the sum of the fragment densities, as in eq 4.2.25. We present these formulas explicitly, for completeness of the present account.
Hartree−Fock energy: Eel =
NAO
Fσ , ab = Hab +
∑ ∑ [PA ,cd(ab|cd) − PA ,σ ,cd(ad|cb)] A
c ,d
(4.2.28b)
electrostatic potential: ϕ( r ⃗) = ϕnucl( r ⃗) − ϕel( r ⃗) Nnucl
ϕnucl( r ⃗) =
∑ a=1
(4.2.29a)
Za | r ⃗ − Ra⃗ |
(4.2.29b)
NAO
ϕel( r ⃗) =
∑ ∑ PA ,ij⟨χi ( r ′⃗ )| ⇀
1
| r − r ′⃗ |
i,j=1 A
|χj ( r ′⃗ )⟩ (4.2.29c)
dipole moment: Nnucl
ϕ( r ⃗) =
∑ a=1
NAO
ZaRa⃗ −
∑ ∑ ⟨χi | r |⃗ χj ⟩ i,j=1 A
(4.2.30)
The lack of long-range interactions is one of the main shortcomings of the ADMA. It has been addressed recently. As argued in most of the early works on the MEDLA and ADMA, the electrostatic interactions can be accounted for by increasing the size of “parent” molecules used to calculate the fragment densities from the first principles. It has been found that high accuracy calculations require surroundings of each fragment in a radius of more than 10 Å.480 This size deteriorates computations significantly, and more efficient approaches were necessary. To address this efficiency issue, Exner and Mezey developed the
NAO A
A
with
∑ χi ( r ⃗) χj ( r ⃗) ∑ PA ,σ ,ij i,j=1
∑ [tr(PA ,α(H + Fα)) + tr(PA ,β(H + Fβ))] (4.2.28a)
electron density: ρσ ( r ⃗) =
1 2
(4.2.27) 5842
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
field-adapted ADMA (FA-ADMA).481 It essentially follows the QM/MM ideas of electrostatic embedding. The schematic procedure of the FA-ADMA is presented in Figure 13. First, the standard ADMA procedure is performed, to initialize the process. Once the charge density of the entire molecule is obtained, the Mulliken charges on all atoms can be easily computed. In the following steps the calculations on all parent molecules are performed, to obtain the density matrices for the corresponding fragments. Unlike the standard procedure, the Hamiltonian of the parent molecule includes the electrostatic field created by the point charge of environment that is not included into the parent molecule itself. The computed fragment densities reflect the existing long-range electrostatic potential created by the environment. The resulting fragment densities are utilized to construct the charge density in the entire molecule and to update all point charges. The cycle is repeated until selfconsistency of charges and fragment densities is achieved. The FA-ADMA has become much more similar to the original D&C approach than to the original ADMA. The ADMA is essentially the tight-binding approximation to the D&C procedure. The FA-ADMA can be considered yet another approximation of the D&C schemethe one that utilized only a conglomerated quantity derived from the density matrix and approximating its effects (Mulliken charges). One can envision further approximations that mimic the effect of the entire density matrix to a higher accuracy. One such approach uses multipole expansions of the charge density, extensively utilized in the quantal force fields. Thus, we have outlined direct connections between the FA-ADMA, the quantal force fields, and the QM/ MM approaches. We have also drawn the hierarchy of approximations: D&C → FA-ADMA → ADMA/MEDLA. Finally, one can find many similarities with the FMO method to be described in section 4.2.2.1. In particular, the embedding part for the incorporation of long-range effects of the environment is the same between the two formulations, as well as the QM/MM approaches. 4.2.2. Energy-Based Fragmentation. 4.2.2.1. Fragment Molecular Orbital (FMO) Method. One of the natural and chemically intuitive ways of accelerating quantum-chemical calculations, in particular for large systems, is to divide the entire molecule into smaller fragments, perform calculations on the smaller fragments, and then obtain the properties (typically energy) of the combined system using the properties of the smaller building blocks. An approach of this type, known as the fragment molecular orbital (FMO) method, was first proposed by Kitaura,482 and then extensively developed by Kitaura and Fedorov.90,483−485 To be fair, this philosophy arose a long time ago in the diatomics-in-molecules (DIM)486 and the pair interaction molecular orbital (PIMO) methods,487 but only the FMO has become a fully functional technique for sufficiently large systems. Further development of the FMO method included its generalization to the ONIOM-type multilayer FMO method (MLFMO),488 which led to acceleration of the molecular integral calculations.483 The multilayer approach allows utilization of alternative levels of theory and different basis sets, to describe complementary parts of the system with varying accuracy and, possibly, to accelerate calculations. The FMO was extended to include correlation at the MP2,489 CIS,490 CC,491 and MCSCF492 levels. The correlation methods were included into the multilayer schemes, allowing further acceleration. A typical MCSCF calculation involved about 1300 basis functions and could be completed in less than 2 h on a standard PC with 1 GB
of memory. A combination of the FMO method with the TDDFT was developed, allowing one to access excited state properties of large systems, including excited state energy gradients.493,494 The method was extended to open-shell systems.494,495 An efficient parallelized implementation of the FMO based on the MP2 level of theory was tested. The calculation involved over 3000 atoms and 44 000 basis functions and was completed in ca. 7 min on ca. 132 000 cores.496 Very recently, an unprecedented record-breaking simulation of the fullerite system containing over 106 atoms has been reported.77 The simulation relied on both FMO for linear scaling of computations and on the DFTB Hamiltonian for a smaller prefactor in the computing time. A single point calculation of such a system takes only slightly more than 1 h of the wall time on a 128-core Xeon cluster. The extensive reviews on applications of the FMO method to large-scale systems and computation of different properties can be found elsewhere.76,90 4.2.2.1.1. Basic Formulation. To illustrate the basic idea of the FMO method, we start with the Hamiltonian of the entire system under consideration, eq 2.1.16a, which can be summarized as N
Ĥ =
⎧ 1 ∇i 2 − 2 ⎩ ⎪
∑ ⎨− ⎪
i=1
Nnucl
∑ a=1
⎫ ⎬+ | ri ⃗ − Ra⃗ | ⎭ Za
⎪
⎪
N
∑ i,j i>j
1 | ri ⃗ − rj⃗| (4.2.31)
One can partition the Hamiltonian of the full system into the Hamiltonians of the subsystems, Ĥ I, and the Hamiltonians describing interactions between the subsystems. In particular, we consider the fragment I, immersed (embedded) into the environment of all other fragments (Figure 14a). The sum N
∑ i,j i>j
1 | ri ⃗ − rj⃗|
can be split into electronic repulsion within the fragment (yellow double arrows in Figure 14a) and interaction with the electrons of all other fragments (green arrows): N
∑ i,j i>j
1 = | ri ⃗ − rj⃗|
NI
∑ i,j∈I i>j
1 + | ri ⃗ − rj⃗| NI
+
NJ
∑∑∑ J≠I i∈I j∈J i>j
NJ
∑∑ J≠I i,j∈J i>j
1 | ri ⃗ − rj⃗|
1 | ri ⃗ − rj⃗|
(4.2.32)
where NI is the number of electrons in the group I. The last two terms in eq 4.2.32 can be grouped to produce
∑ ∫ d r ′⃗ J≠I
ρJ ( r ′⃗ ) | ri ⃗ − r ′⃗ |
where ρJ(r′⃗ ) represents the electronic density due to the fragment J. Finally, the Hamiltonian of the group I can be written as NI
HÎ =
∑ i=1
5843
{
1 − ∇i 2 + VI( ri ⃗) 2
}
NI
+
∑ i,j i>j
1 | ri ⃗ − rj⃗|
(4.2.33)
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review NI + NJ
HIJ =
∑ i=1
⎧ ⎪ 1 ⎨− ∇i 2 − ⎪ ⎩ 2
∑ a=1
Za | ri ⃗ − Ra⃗ | ⎫ ⎪
∑
+
Nnucl
J≠I ,K
NI + NJ
ρ ( r ′⃗ ) ∫ d r ′⃗ | r ⃗K− r ′⃗ | ⎬ + ∑ i,j ⎪
⎭
i
i>j
1 | ri ⃗ − rj⃗|
HIJ ΨIJ = EIJ̃ ΨIJ
Nnucl a=1
Za | ri ⃗ − Ra⃗ |
+
∑ ∫ d r ′⃗ J≠I
ρJ ( r ′⃗ ) | ri ⃗ − r ′⃗ |
(4.2.34)
where VI(ri⃗ ) plays the role of an effective potential that acts on the electrons of the fragment I. It is often called the embedding potential of the fragment I. Note that the motion of the electrons that “belong” to the fragment is affected by the field of all available nuclei, not only by those belonging to the subsystem. It is also affected by the electrostatic potential created by the electrons from other fragments. Thus, the coupling of a given fragment with all other fragments is mostly due to electrostatic interactions. Solving the stationary electronic problem for all fragments {I}: HI ΨI = EĨ ΨI
(4.2.37)
leading to electronic energies of the dimers, Ẽ IJ, and the corresponding dimer-localized fragment MOs, ΨIJ. The total energy of a dimer is obtained by adding the corresponding nuclear repulsion terms, similar to how this is done for monomers. The frozen monomer densities, ρK(r′⃗ ), are typically used for simplicity and computational efficiency. The approximation implies that there is no need for an iterative solution of the fragment dimer stationary Schrödinger equation, eq 4.2.37. It should be noted that although the Hamiltonians of the monomers and dimers include electrostatic interactions with the electronic densities of surrounding fragments via the
Figure 14. Definition of partitioning of the system into monomers (a) and dimers (b) embedded into the environment of all other fragments. Small blue circles represent atoms. Regions represent fragments.
VI( ri ⃗) = − ∑
∑ ∫ d r ′⃗ J≠I ,K
ρK ( r ′⃗ ) | ri ⃗ − r ′⃗ |
terms, the energy of the dimer is not equal to the sum of the monomer energies: EIJ ≠ EI + EJ. The difference originates from the interfragment electronic interaction shown by the white arrows in Figure 14b. It can be expressed as NI + NJ
(4.2.35)
HIJ − (HI + HJ ) =
EI = EĨ +
∑ a,b∈I aj
subject to the constraint of the total number of electrons localized on each fragment, NI, and to spatial localization of the corresponding wave functions on the fragments, one obtains the localized fragment MOs, ΨI, and the electronic energies of monomers, E′I. The total energy of the monomer, EI, is then obtained by adding nuclear repulsion terms to the electronic energy, Nat
(4.2.36)
⎛ ⎜ NI 1 − ⎜∑ + r | − ⎜ i , j ∈ I i ⃗ rj⃗| ⎝ i>j NI
=
ZaZb
NJ
∑∑ i∈I j∈J
|Ra⃗ − R⃗b|
1 | ri ⃗ − rj⃗|
NJ
∑ i,j∈J i>j
⎞ 1 ⎟ ⎟ | ri ⃗ − rj⃗| ⎟ ⎠
(4.2.38)
The types of interactions included in prototypical monomer and dimer fragment Hamiltonians are represented pictorially in Figure 14. The difference between the dimer energy and the sum of monomer energies:
The electron densities of the monomers are obtained in an iterative way, similar to how the electron density matrix is iteratively obtained in a standard SCF scheme in an atomic/ molecular orbital basis. In this way, changes of the electron density of one fragment induce changes of the electron density of the all other fragments in a self-consistent manner. The idea of performing self-consistent calculations on sets of fragments and adjusting their embedding potentials was proposed back in 1975 as the mutually consistent field method.497 The monomer partition scheme presented above helps in reducing the computational costs drastically. Instead of solving the eigenvalue problem for N electrons that scales as O(N3), one solves only smaller problems for Nfr subsystems with {NI|I ∈ 1, ..., Nfr} electrons, respectively. The computation time in the latter N scheme scales as O(∏I=1fr NI3), which is significantly more favorable that that for the complete system, especially when the size of the system increases. In addition to monomers, all possible dimers are computed, using the dimer Hamiltonian defined similarly to eq 4.2.33, but on a bigger set of atomic and orbital centers:
ΔEIJ = EIJ − (EI + EJ )
(4.2.39)
is crucial for an accurate description of the energy of the combined system. It reflects the complex nonadditive effects. The total energy of the system is expressed within the FMO as Nfr
E=
∑
Nfr
EIJ − (Nfr − 2) ∑ EI
I ,J=1 I>J
I=1
(4.2.40)
This expression can be rewritten in a more systematic way: Nfr
E=
Nfr
∑ EI +
∑
I=1
I ,J=1 I>J
ΔEIJ (4.2.41)
Equation 4.2.41 is a special case of the generalized many-body expansion (GMBE) up to m terms: 5844
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews Nfr
Em =
Review Nfr
Nfr
Nfr
Nfr
Nfr
Here, E′X is the intrinsic energy of the fragment (X = I for monomer, X = IJ for dimer) without the embedding potential, VX is the embedding (environment) potential, and PX is the density matrix of the fragment X. Equations 4.2.47 are equivalent to eq 4.2.41, but are much more convenient for introducing approximations to reduce computational costs associated with the embedding potential. For sufficiently separated fragments, the energy of the dimer can be approximated by the sum of the monomer energies plus electrostatic interactions between them. In addition, Mulliken charges can be used instead of the full density matrix, if the interfragmental distance is large enough. For details and discussion of these approximations, we refer the reader to the original work483 and to the recent review.498 We would like to emphasize two points. First, the energy decomposition in the FMO method in eqs 4.2.47 has a useful interpretation. The term
(3) + ... ∑ EI + ∑ ∑ ΔEIJ(2) + ∑ ∑ ∑ ΔEIJK I ,J=1 J=1 J>I
I=1
=1 I ,J=1 J=1 K JJ
E′X = EX − Tr(PXVX ),
ΔPIJ = PIJ − (PI ⊕ PJ )
(4.2.47a)
X = I , IJ
+
∑
ΔE′IJ
I ,J=1 I>J
accounts for implicit many-body effects. This is because the energies of the monomers and dimers are affected by the embedding potential that incorporates the many-body effects. Note that the charge transfer energy is largely composed of the explicit contrbution, but also contains a fraction of the implicit one. The second point concerns the utility of semiempirical schemes. Computation of molecular integrals constitutes a major bottleneck of the FMO method. Therefore, approximate but significantly faster schemes,483 similar to those utilized in the semiempirical methods, are favored over the more accurate, but significantly slower HF-type descriptions. Combination of the FMO method with another semiempirical approachthe DFTB methodallowed reaching the unprecedented scale of quantum calculations employing over 106 atoms.77 By systematically including the higher-order corrections, one can achieve the desired accuracy for large systems. Depending on the perturbation (correction) level, one can identify FMO2, FMO3, etc. methods. The three-body terms may be required sometimes to achieve acceptable accuracy.484,499 In particular, explicit water solvation of biopolymers or accurate calculation of dipole moments may require the use of three-body terms. The role the four-body correction was assessed by Nakano et al.500 They found that such a correction is important, if a nonconventional partitioning of peptides is used (segmentation of residues of side and main chain). The fourth-order FMO was also found essential for accurate modeling of systems having threedimensional structures, such as adamantane-shaped clusters. The higher-order terms are negligible most of the time. As the size of the system increases, the number of fragments into which the system should be partitioned increases as well. This naturally leads to an increased number of possible n-mers that need to be computed. Specifically, the number of n-mers that can be formed from Nfr monomers is given by the combinatorial number CnNfr = Nfr/[n!(Nfr −n)!]. Because the number grows very
(4.2.45)
Nfr
Nfr
Nfr
∑ E′I I=1
m = 3 I1> I2 > ... > Im − 2
aIJ , I1...Im−2 =
[Tr(ΔPIJVIJ )]
accounts for explicit many-body effects due to charge transfer between fragments I and J under the influence of the electrostatic potential of all other fragments. The term
Nfr
I ,J=1 J=1 J>I
∑ I ,J=1 I>J
∑ AI + ∑ ∑ ΔAIJ′ I=1
Nfr
(4.2.47b) (4.2.47c) 5845
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
rapidly, it is common to include only dimers in the calculations on large systems. The approximation provides sufficient accuracy for most applications and, therefore, is well justified. A nearly linear scaling can be achieved using efficient approximations for computation of the energies of far separated dimers.483 Additional savings may be obtained if many dimers are similar to each other or are composed of similar monomers. 4.2.2.1.2. Elaborations and Derivatives. In the original FMO method, as well as in the MBE approach, the fragments into which the system is partitioned are assumed nonoverlapping. Formulations of the fragment-based methods for overlapping fragments were also proposed. For example, Fedorov proposed an extension of the FMO method suitable to partitions in solids.501 The improved scheme was applied to zeolites. 4.2.2.1.3. Implementations. The FMO and its various improvements are implemented in GAMESS and AbInitMP.490 The Web interface, input preparation, and analysis tools were developed by Fedorov and co-workers. 4.2.2.2. Molecular Fractionation with Conjugate Caps (MFCC). The MFCC method was developed by Zhang and coworkers,502 as an approach to the efficient calculation of the interaction of a protein, P, with an arbitrary molecule, M (e.g., solvent). The protein molecule, P, can be represented schematically as a sequence of connected amino acid fragments, Ai: P = nA1−A 2−A3−...−AN c
Figure 15. Schematic representation of the fractionation scheme of a protein, used in the MFCC method. + ⎧ i=0 ⎪ NH (NH 2 3 ), Ci* = ⎨ ⎪ ⎩ NH 2 , i = 1, ..., N − 1
⎧ R i + 1−Cα H 2 , i = 0, ..., N − 1 ⎪ Ci = ⎨ A C = R −C H−COO− (A C N i i α ⎪ i i R C H COOH) , i=N = − − ⎩ N α
Schematically, the choice of the capping fragment for the system shown in Figure 15 is shown in Figure 16. The method was combined with the conductor-like polarizable continuum model (CPCM), to describe solvation of proteins.503 4.2.2.3. Systematic Molecular Fragmentation (SMF) Approach. The systematic molecular fragmentation approach504 was developed by Deev and Collins on the basis of the Chen’s MFCC method.505 Further elaborations were performed by the Collins group and included extensions to treat rings,506 longrange506,507 and nonbonded508 interactions, and to compute molecular potential energy surfaces507,509 and electrostatic potentials.510 The method was extended to periodic systems.507,511 It allows utilization of high-level correlated ab initio approaches, such as MP2 and CCSD(T). The method was applied to model lattice dynamics and to compute phonon spectra and neutron scattering intensities.507 Very recently, an improved version of the algorithm (so-called, SMF by annihilationSMFA) has been applied to compute the relative energies of a series of proteins.512 The idea of the SMF method is very simplethe molecule is hierarchically divided into a set of overlapping fragments:
(4.2.48)
where nA1 and ANc denote N- and C-terminals of the protein: nA1 = NH3+R1−CαH−CONH
(4.2.49a)
AN c = R N −Cα H−COO−
(4.2.49b)
and Ai is the amino acid residue: A i = R i−Cα H−CONH
(4.2.50)
The fractionation scheme for proteins is presented in Figure 15. To compute the protein−molecule interaction energy, E(P ∪ M), the protein is subdivided into fragments, chosen in the MFCC to be the amino acid residues. Each fragment is capped with a pair of conjugate cap groups: C*i and Ci. The interaction energy of the capped fragment with the molecule M, E(Ci−1 * AiCi ∪ M), is computed for each fragment i, along with the energy of interaction of the pair of conjugate caps with the molecule, E(C*i Ci ∪ M). The protein−molecule interaction energy is then computed as N
E (P ∪ M ) =
M 0, n → F1 + F2
(4.2.53)
F1 ∩ F2 = ϕ
(4.2.54)
Each fragment is capped with hydrogens:
Fi → Fi ∪ H
N−1
∑ E(Ci*− 1AiCi ∪ M) − ∑ E(Ci*Ci ∪ M) i=1
(4.2.52)
(4.2.55)
and the resulting smaller molecules are used to repeat the process recursively:
i=1
(4.2.51)
Mi , n + 1 = Fi ∪ H
The role of the conjugate caps is 2-fold: • to better represent the electronic structure of the system around the covalent bond being cut; in this regard, the conjugate caps play the role very similar to the buffer regions in the D&C methods
(4.2.56)
In the above equations, the index n indicates the level of fragmentation hierarchy. The energy of the molecule is approximated by a simple expression: E(Mn) ≈ E(M1, n + 1) + E(M 2, n + 1) − E((F1 ∩ F2) ∪ H ) (4.2.57)
• to better approximate the environment of the correspond-
which can be applied up to a specified level of fragmentation. The two main restrictions are (a) molecules are broken only through single bonds, not multiple bonds and (b) no explicit charge separation across the system occurs. The derivation of the
ing amino acid residue To satisfy the second requirement, the most natural choice for the capping groups is based on the nearby residue: 5846
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
set of overlapping fragments.513 Therefore, the method follows the general scheme of the D&C algorithm. However, unlike the D&C procedure, the electronic densities of the fragments are not affected by those of all other fragments, and no iterative procedure is required. Instead of the partitioning schemes, such as eqs 4.2.23a and 4.2.23b, coupled with the iterative SCF procedure for all fragments, the method simply composes the densities of individual fragments in a proper way, to approximate the density matrix and energy of the whole system. The density matrix composed of the subsystem density matrices stitched together may not produce the population equal to the total number of electrons. An appropriate scaling factor is computed in this case, and the density matrices of the subsystems are rescaled. Because the MTA does not adjust the density matrices of the subsystems self-consistently, the choice of the partitioning and the overlap scheme is especially critical for this approach. In particular, double and triple bonds, as well as aromatic and small rings, should not be broken into different fragments. The fragments should contain a connected graph and should be sufficiently large to simulate accurately the properties of the composed molecule. To facilitate fragmentation, an automatic partitioning algorithm was developed by the Gadre group.514 The non-self-consistent nature of the MTA approach makes it very similar to other methods, notably the density builder approach (MEDLA or ADMA) of Mezey and co-workers, and the FMO of Fedorov and Kitaura. The similarity with the first method comes from the composition procedure for the density matrix (or charge density). The procedure may or may not be self-consistent in the density builder approaches. The similarity with the FMO method comes from the fact that both methods are noniterative at the total energy levels and can be formulated merely in terms of the total energies. The method was extended to the MP2514 and periodic systems.517 The automatic fragmentation procedure was developed.514,518 Following the initial attempts to perform MTA-based optimization,519 a recent modification of the original MTAthe cardinality-guided MTA (CG-MTA)incorporates energy gradients and Hessians for structural optimization.518 The CG-MTA is based on set inclusion and exclusion principles. The molecular system is initially partitioned into fragments of a specified maximal size and the so-called Rgoodness of each atom in a given fragment. The latter parameter is used to assign each atom to a specific fragment. The initially prepared fragments are capped by hydrogen atoms, so as to saturate dangling bonds. Once the set of fragments, {Fi}, is generated, one can utilize Venn diagrams to compute the prefactor specifying the contribution of the energy of each fragment, {Ei}, to the total energy of the system. The prefactor is computed as (−1)K−1, where K is the number of fragments involved in the intersection. The basic scheme of the CG-MTA is presented in Figure 17. The total energy is computed on the basis of the energies of all fragments and their intersections. In the case of two fragments, the formula has a very simple form:
Figure 16. Schematics showing the fragmentation of a polypeptide in the MFCC method. The conjugate caps are marked by blue (Cterminus) and red (N-terminus) boxes. The fragments are marked by green boxes. The molecules to the right of the capped fragments 0 and 1 represent the total addition to the bare fragment. The energy of this additional molecule must be subtracted from the energy of the capped fragment to get the energy contribution of a fragment.
SMF implies that the broken bonds are sufficiently distant from each other. Therefore, the larger are the resulting fragments, the more accurate are the computed energies, although at a cost of reduced computation speed. The choice of the fragmentation scheme is critical. The fragmentation can be performed systematically, at different levels. It is illustrative to consider a chain of five functional groups (fragments), F1F2F3F4F5. The possible division schemes are FF 1 2F3F4F5 → FF 1 2 + F2F3 + F3F4 + F4F5 − F2 − F3 − F4 (level 1)
(4.2.58a)
FF 1 2F3F4F5 → FF 1 2F3 + F2F3F4 + F3F4F5 − F2F3 − F3F4 (level 2) FF 1 2F3F4F5 → FF 1 2F3F4 + F2F3F4F5 − F2F3F4
(4.2.58b)
(level 3) (4.2.58c)
The higher the level, the better the accuracy. Because of the recursive nature, the size of each fragment in the final level can be kept sufficiently small, to make calculations of each subsystem as fast as possible and to achieve linear scaling of CPU time with the size of the whole system. The energy gradients and higher-order derivatives can be easily obtained, making the method applicable to geometry optimization and molecular dynamics problems. On the other hand, the method is nonvariational and is based solely on energies of the system. Therefore, the electronic structure of the original system cannot be obtained, restricting the method to the ground electronic state. In principle, a time-dependent formulation may be envisioned, to access properties of excited states, similarly to the extensions of the FMO method.493 However, excitations are often associated with notable charge delocalization over entire system, which may violate the basic assumptions of the method. 4.2.2.4. Molecular Tailoring Approach (MTA). Initially, the MTA approach was extensively developed by the Gadre group.513−516 They developed a method for theoretical synthesis of molecular electrostatic potentials of large molecules using the information obtained by computing the electronic structure of a
E = E(F1) + E(F2) − E(F1 ∩ F2)
(4.2.59)
For the system shown in Figure 17b the expression also includes a tertiary term: E = E(F1) + E(F2) + E(F3) − E(F1 ∩ F2) − E(F1 ∩ F3) − E(F2 ∩ F3) + E(F1 ∩ F2 ∩ F3) 5847
(4.2.60) DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Figure 17. (a) Fragmentation in the CG-MTA. (b) Venn diagram showing set inclusion and exclusion principle, as applied to overlapping fragments. The regions indicated by the black letters produce positive prefactors in the energy expression, while the regions indicated by the blue letters produce negative energy prefactors. Fi denotes fragment i.
Nonetheless, the point charges provide the most natural interpretation of the quantities ρi in the E-EDC context. The charge transfer, Δρ, is assumed to be relatively small, such that the Taylor expansion eq 4.2.63 can be truncated after the quadratic terms. One considers charge transfer between the sites i and j, and associated with it change of the total energy:
Finally, the expression for the general case is E=
∑ E(Fi) − ∑ E(Fi ∩ Fj) + ∑ E(Fi ∩ Fj ∩ Fk) + ... i
i,j
i ,j,k
(4.2.61)
The expressions for the energy gradients and Hessians are straightforward. They are obtained by differentiation of the energy expression, eq 4.2.61, with respect to the coordinates of interest: ∂E = ∂Xa
∑ i
+
∂E(Fi ) − ∂Xa(Fi )
∑ i ,j,k
∑ i,j
∂E(Fi ∩ Fj)
∂Xa(Fi ∩ Fj ∩ Fk)
(4.2.64a)
ΔEj = aj( −Δρ) + bj( −Δρ)2
(4.2.64b)
ΔE = ΔEi + ΔEj = (ai − aj)Δρ + (bi + bj)(Δρ)2
∂Xa(Fi ∩ Fj)
∂E(Fi ∩ Fj ∩ Fk)
ΔEi = aiΔρ + bi(Δρ)2
(4.2.64c)
Using the variational principle and the DFT, stating that any charge density fluctuation with respect to the ground state density will raise the energy, eq 4.2.64c, the authors show that, under the above approximations, the gradient of the energy with respect to density is constant at any point:
+ ... (4.2.62)
where the coordinate Xa(S) of the variable a belongs to the set S. The energy formula eq 4.2.61 generalizes similar expressions obtained by Deev and Collins in their SMF approach,504 and by Chen et al. in their molecular fractionation with conjugate caps (MFCC).505 4.2.2.5. Energy-Based D&C with Charge Conservation (EEDC). Energy-based divide-and-conquer (EDC) with charge conservation (E-EDC) was developed by Song, Li, and Fan.520 They pointed out that the major errors in the energy-based fragmentation schemes originate from artificial charge transfer between fragments. Eventually, it violates conservation of the total charge of the whole system. The E-EDC method extrapolates the energy value computed with a nonconserved charge to the more accurate value that corresponds to the conserved charge. The mathematical and physical foundation of the method can be summarized as follows. Similar to the SCC-DFTB family of methods, one expands the DFT energy in a Taylor series with respect to the electron density fluctuation, Δρ:
ai = aj = ... ≡ a =
∂E ∂ρ
(4.2.65)
Assuming that the second order terms in Δρ can be neglected, and generalizing eqs 4.2.64 to many fragments, a linear relationship between the total energy change, ΔE = ∑i ΔEi, and the total charge change, ΔQ = ∑i Δρi, can be derived: ΔE ≈ aΔQ
(4.2.66)
The subsequent computations are organized to compute the points necessary for the extrapolation based on eq 4.2.66. First, the system is fragmented in different waysmonomers, dimers, trimers, etc.in a manner similar to the FMO or GMBE formulations. All subsystems are capped with H atoms if chemical bonds are broken. Second, the total energy of the entire system at a given order of the many-body expansion, Em, can be computed according to the general formula (see section 4.2.2.1). In the reported scheme, the authors neglect the effect of point charges of the environment on the energy of the embedded fragment. Each of the fragments (monomers, dimers, etc.) is computed independently, in a vacuum. The total charge of the system at a given level of MBE is computed as
ΔEi = E(ρi + Δρ) − E(ρi ) = aiΔρ + bi(Δρ)2 + ... (4.2.63)
where index i refers to the spatial point at which the charge density, ρi, is changed the most. Apparently, the approximation is derived with the intent to relate this point charge density to the classical point charges. In a more rigorous approach, the electron density fluctuation may be spread out over a notable spatial region. The local approximation may not always be accurate.
Qm =
∑ Q Am(A) A
5848
(4.2.67) DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
The quantity QmA (A) is computed using the MBE, similar to the energies: Nfr
Q Am(A)
= Q A(A) +
∑
Nfr (2) ΔQ AI (A )
I=1
+
aperiodic polymers at the semiempirical, ab initio HF, and DFT levels of theory.541−544 Recently, the method has been extended to include correlation at the MP2 level,545 and to describe excited states at the CIS546 and linear-response time-dependent Hartree−Fock547 levels. Geometry optimization548 and molecular dynamics549 schemes based on the ELG method were developed. The improved version for delocalized systems was also described.550 Before providing the mathematical foundations of the method, we outline its main idea. The short-range nature of chemical bonding constitutes the major assumption of the method. (Strictly) localized molecular orbitals ((S)LMOs) are the central tool of the ELG. The system of interest can be subdivided into the initial fragment F0 and a set of monomers of relatively small size, Mn, n = 1, ..., Nmonomers. The initial fragment F0 is further subdivided into two components, A0 and B0, such that the fragment A0 is sufficiently distant from the site to which the monomers are attached, while the fragment B0 contains that site (Figure 18). It is conventional to assume that the subfragment B0 is larger than the subfragment A0. Having defined the initial fragment F0 = A0 ∪ B0, one obtains its molecular orbitals by solving the conventional SCF equation, eq 2.1.38. The resulting canonical MOs are typically delocalized over the entire subsystem F0. That is, the MO expansions in the basis of standard Cartesian AOs, {(χa)}, involve many AOs often localized in different regions across the system and having comparable weights. In order to make use of the short-range chemical bonding principle, the MOs are localized into the regions A0 and B0. This is achieved by using the uniform localization method,551,552 which is similar to the one proposed by Edmiston and Ruedenberg.553 To start the fragment localization procedure, the nonorthogonal AOs, {|χa⟩}, are preorthogonalized to form a set of hybrid orthogonal AOs, {|χ̃a⟩}. There exist two major orthogonalization schemes: the one initially proposed by Imamura530 and the improved scheme based on the Löwdin orthogonalization.533 It was found that the orthogonalization scheme critically affects the accuracy of the calculations with the elongation method.537 Orthogonalization is achieved via the unitary transformation:
Nfr
(3) ∑ ∑ ΔEAIJ (A ) I=1 J=1
(4.2.68)
+ ... (2) ΔQ AI (A) = Q AI(A) − Q A(A)
(4.2.69)
(3) ΔQ AIJ (A) = Q AIJ(A) − Q AI(A) − Q AJ(A) + Q A(A)
(4.2.70)
QA(A) is the total charge of all atoms in the subsystem A, excluding the cap atoms. QAI(A)and QAIJ(A) are the total charges belonging to the fragment A in the two-body AI and three-body AIJ subsystems, respectively. With all quantities needed for the linear extrapolation of the energies defined, eq 4.2.66, the extrapolation can be written: E ≈ E M = Em −
Em − Em − 1 m (Q − Q M ) Qm − Qm−1
(4.2.71)
where M is the highest order of MBE included in the scheme. In other words, the total energy is approximated by the M-body expansion, but with charge conservation. The total charge is approximated at the M-body level, QM, by the total charge of the system, Q, which is known as an input to the calculations: QM = Q. 4.2.2.6. Other Related Methods. There are many more fragmentation-based linear-scaling methods. They follow the ideas of the methods presented above, but adapt a slightly different formulation and introduce corrections for missing effects. Therefore, rather than discussing these methods in detail, we list them here for completeness of the review: (1) molecules-in-molecules (MIM),521 utilizing the ideas of the ONIOM and energy-based fragmentation schemes like the FMO (2) many-overlapping body (MOB) expansion,522 a two-body expansion using overlapping fragments; the method was utilized for computations of large polypeptide and alkane chains (3) dimer of dimers (DOD) method by Saha and Raghavachari523 (4) generalized energy-based fragmentation (GEBF) by Li et al.454,524,525 (5) kernel energy method (KEM) by Karle et al.526 (6) hybrid many-body interaction (HMBI) by Beran and Nanda;527,528 this method is close to quantum force fields and hybrid QM/MM approaches (7) multilevel fragment-based approach (MFBA) by Ř ezáč and Salahub529 4.2.3. Dynamical Growth with Localization. 4.2.3.1. Elongation Method (ELG). The elongation method was originally proposed by Imamura and co-workers as a computational approach to electronic structure calculations that mimic a theoretical synthesis of polymers.530 Hence, the first applications were to one-dimensional systems, but later extensions to twoand three-dimensional systems531,532 were also developed. The method was also actively developed by the Aoki lab. Further elaborations included an improved localization scheme,533 calculation of the local density of states,534 implementation of an efficient cutoff scheme for additional acceleration of ab initio calculations,535−539 and extensive analysis of performance and efficiency.537 The method was applied to study periodic540 and
NM
|χã ⟩ =
∑ Uba|χb ⟩
(4.2.72)
b=1
where U is the corresponding unitary transformation matrix (UU+ = U+U = I). The canonical MOs of each fragment, which are typically expressed in the nonorthogonal AO basis, eq 2.1.12, can then be represented in the basis of the orthogonal hybridized AOs: NM
|ψσ , i⟩ =
NM
∑ Cσ ,ai|χa ⟩ = ∑ a=1
NM + Cσ , aiUab |χb̃ ⟩ =
a,b=1
σ = α, β
∑ Cσ̃ ,bi|χb̃ ⟩, b=1
(4.2.73)
NM
Cσ̃ , bi =
∑ Cσ ,aiUba+ ⇔ Cσ̃ = U +Cσ a=1
(4.2.74)
The SCF equation, eq 2.1.38, becomes in the new basis ̃ σ̃ Eσ Fσ̃ Cσ̃ = SC
(4.2.75)
with Fσ̃ = U +Fσ U 5849
(4.2.76) DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
bij =
∑
Cσ̃ , ajCσ̃ , bjSab̃ +
a,b∈A
cij =
∑
∑
Cσ̃ , aiCσ̃ , biSab̃ (4.2.82b)
a,b∈B
Cσ̃ , aiCσ̃ , bjSab̃ −
a,b∈A
∑ a,b∈B
Cσ̃ , aiCσ̃ , bjSab̃ (4.2.82c)
Extremum of the localization measure is achieved at π ω θext = − (4.2.83a) 4 2 ⎛ bij − aij ⎞ ⎟⎟ ω = a tan⎜⎜ ⎝ 2cij ⎠
The value that maximizes eq 4.2.81 is chosen from the set of values that deliver the extremum. Initially, a set of canonical MOs is divided into two groups: occupied, {ψocc,σ,i}, and virtual {ψvirt,σ,i}. The sequence of 2 × 2 rotations, eq 4.2.78, is performed to achieve localization. The rotations are carried out separately for pairs belonging to each of the two groups, until any pair of the resulting localized orbitals ψ̃ occ,σ,i, ψ̃ occ,σ,j ∈ {ψ̃ occ,σ,i} or ψ̃ virt,σ,i, ψ̃ virt,σ,j ∈ {ψ̃ virt,σ,i} gives no further localization improvement. The ideal localization cannot be achieved typically, and the iterations are stopped after the interfragment coupling,536 d = ∑i∈B ∑a,b∈A |C*ai SabCbi|, drops below a predefined threshold value. Once the localization within each of the two sets of orbitals is completed, the four sets of LMOs are obtained: {ψ̃ A,occ,σ,i}, {ψ̃ A,virt,σ,i}, {ψ̃ B,occ,σ,i}, and {ψ̃ B,occ,σ,i}. The resulting Hamiltonian has the following block structure:
Figure 18. Basic scheme of the elongation method. The scheme shows how the original molecule is gradually reconstructed by sequential addition of monomers. Mn is the monomer attached at time step n. Fn is the n-mer: the structure being computed after n steps. An and Bn are parts of the n-mer computed. The part A is frozen and effectively excluded from computations. The part B is treated fully at the nth step.
S ̃ = U +SU
(4.2.83b)
(4.2.77)
The fragment localization procedure consists of the following steps. The pair of LMOs ψ̃ σ,i and ψ̃ σ,i is obtained by a rotation of the pair of the canonical MOs, ψσ,i and ψσ,j, in the hyperplane i−j: ⎛ ψσ̃ , i ⎞ ⎛ cos θ ⎞⎛ ψσ , i ⎞ ⎜ ⎟ = ⎜ sin θ ⎜ ⎟ ⎜ ψ̃ ⎟ ⎝−cos θ sin θ ⎟⎠⎜ ψσ , j ⎟ ⎝ ⎠ ⎝ σ ,j ⎠
At this point one can utilize localized MOs (LMOs) of the lead fragment B0 and similarly obtained LMOs (or canonical MOs) of the first monomer fragment, M1, to describe coupling of the fragment F0 with the first monomer. The resulting combined Hamiltonian takes the form
(4.2.78)
The rotation angle θ is chosen to deliver the maximum of the sum of populations localized on each of the regions A and B (localization measure): Lij = ⟨ψà , σ , i|ψà , σ , i⟩ + ⟨ψB̃ , σ , j|ψB̃ , σ , j⟩
(4.2.79)
Because of the spatial separation of the localization regions A0 and M1, the LMOs of the fragment A0 are unlikely to be affected by introduction of the monomer M1, provided the size of the region B0 is sufficiently large. Therefore, the Hamiltonian elements that couple fragments A0 and M1 are negligible, which is indicated by ∼0 in eq 4.2.85. Further, it is assumed that the introduction of the monomer M1 does not affect the coupling between the fragments A0 and B0, and that coupling between these fragment will not have effect on optimization of the wave function for the extended system (with monomer M1). Therefore, the elements F̃ occ,A,B, F̃virt,A,B, F̃ occ,B,A, and F̃ virt,B,A, are set to zero. This approximation turns the Hamiltonian of the extended system into the block-diagonal Hamiltonian:
where ψX̃ , σ , i =
∑ [sin(θ)Cσ̃ ,ai + cos(θ)Cσ̃ ,aj]χã a∈X
+
∑ [−cos(θ)Cσ̃ ,ai + sin(θ)Cσ̃ ,aj]χã a∈X
(4.2.80)
Using eq 4.2.78, the localization measure can be expressed as Lij = aij sin 2(θ ) + 2cij sin(θ ) cos(θ ) + bij cos2(θ ) (4.2.81)
with aij =
∑ a,b∈A
Cσ̃ , aiCσ̃ , biSab̃ +
∑ a,b∈B
⎛ F̃ ⎞ A ,A 0 ⎟ FF̃ + M = ⎜⎜ ⎟ ̃ 0 F ⎝ B+M ⎠
Cσ̃ , ajCσ̃ , bjSab̃ (4.2.82a) 5850
(4.2.86) DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
technique does not help to get rid of the unfavorable cubic scaling on its own. In the most straightforward implementation, it only reduces the prefactor by up to 4 times. Using the principle of occupied−virtual block annihilation,554 one is only interested in considering the pairs composed of one occupied, ψi, and one unoccupied, ψa, orbitals. Moreover, it is always possible to exclude distant pairs because of their small overlap and negligible magnitude of the matrix element. For the pair of orbitals whose Fock matrix element
with ⎞ ⎛ F̃ occ, A , A − E · I 0 ⎟ ⎜ ̃ FA , A = ⎜ ⎟ ̃ 0 F E I − · virt, A , A ⎠ ⎝
(4.2.87)
The canonical orbitals for the new fragment B1 = B0 ∪ M1 are obtained by diagonalizing the corresponding Hamiltonian. The LMOs of the fragment A0 are kept frozen and do not enter the new set of SCF equations, but the density is used to compute the total energy. The new iteration only considers the mixing of the regions B0 and M1 (stage 2). Once the canonical MOs for the combined system B1 are obtained, the procedure described above can be repeatedthe fragment B1 is partitioned into the frozen and active parts, the canonical MOs are localized into the fragments, an additional monomer is added, and the frozen part is dropped (stage 3). The latter feature of the algorithm allows one to keep the size of the active part of the system (the one that determines the dimensionality of the eigenvalue problem) constant at all times, independently of the size of the entire system. This leads to linear scaling of the algorithm in the CPU time. 4.2.3.1.1. Implementation. The elongation scheme with cutoff (ELG/C) was implemented in the GAMESS package. 4.2.3.2. Localized Molecular Orbitals (LMO) Method. In order to reduce the scaling to linear or quasilinear order, Stewart utilized the idea very familiar to most practicing computational chemiststhe localized nature of molecular orbitals. Namely, Stewart proposed the localized molecular orbitals (LMO) method.554 The method has many common points with the ELG method just discussed. Most importantly, both methods utilize dynamical MO growth algorithms and make good use of the localized nature of the MOs. For computations of the MOs, one starts with guess orbitals derived from the standard Lewis diagrams of molecular structure. Each atom is assumed to have one s orbital and three p orbitals. These four orbitals are hybridized and oriented along the directions of possible (sigma) bonding with the nearest neighbors. The hybrids are forced to remain orthogonal to each other. If the number of neighbors is small and some of the hybrids remain unused, their orientation is rather arbitrary, and the orthogonality is preserved. Unlike the conventional MO descriptions, in which MOs are extended to the entire system, the LMO method starts with the localized orbitalsatomic and diatomicthat exist in local regions of space. However, in order to converge the SCF calculations, the LMOs are allowed to grow dynamically. As discussed in the early works of Stewart et al.,555 it is not necessary to compute all eigenvalues and eigenvectors by diagonalizing the entire Fock matrix. It is sufficient to annihilate only those elements of the Fock matrix that couple occupied and virtual MOs. The approach is based on the assumption that the Fock matrix has an approximately diagonal form in the specified basis. Also, it is assumed that the entire matrix of occupied and virtual orbitals is known from the previous calculations: C = (Co Cv ), where Co (nbas × nocc) and Cv (nbas × nvirt) are matrices of occupied and virtual orbitals written in a given basis, respectively. In the LMO approach these orbitals are known by construction. They are hybrids of the corresponding AOs of single atoms and diatomic fragments. It should be noted that the
FiaMO = ⟨ψi|F |̂ ψa⟩ =
∑ ∑ ∑ Ci*λ⟨χλ |F |̂ χμ ⟩Cμa A ,B λ∈A μ∈B
=
∑ ∑ ∑ Ci*λFλμAOCμa A ,B λ∈A μ∈B
(4.2.89)
is not negligible in the MO basis, the annihilation is performed by mixing the two orbitals by ⎛ ψ ′i ⎞ ⎛ α β ⎞⎛ ψi ⎞ ⎜⎜ ⎟⎟ = ⎜⎜ ⎟⎟⎜ ⎟ ⎝ ψ ′a ⎠ ⎝−β α ⎠⎝ ψa ⎠
(4.2.90)
with 1 (1 + D/ 4Fia 2 + D2 ) 2
α=
(4.2.91a)
β = ϕ 1 − α2
(4.2.91b)
D = Fii − Faa
(4.2.91c)
⎧+1, ϕ=⎨ ⎩−1, ⎪
⎪
Fij < 0 F≥0
(4.2.91d)
The LMO growth procedure is shown in Figure 19. Although the MO localization region grows with each annihilation level spanning new atoms and bonds, the contributions of AOs differentiate. For example, the middle atoms give large total wave function amplitudes, while the peripheral atoms give only negligible contributions. As the contributions become smaller than a specified threshold, the AOs are pruned from the LMO. The LMO size becomes relatively constant in this way, typically ranging from 100 to 130 atoms for large systems. Most importantly, the LMO size becomes independent of the size for the entire system, ensuring linear scaling of the computational efforts. The effect is not observed for relatively small systems, because the fully grown LMOs become comparable with the fully delocalized MOs obtained directly. The large-scale simulation of systems with up to 1 000 000 atoms on parallel computers is among the most recent application of the LMO method.556 Although the LMO algorithm is the key to success, efficient parallelization and interprocessor communication strategies, and modern computer architectures involving symmetric multiprocessing (SMP), are essential. 4.2.3.3. Iterative Stochastic Subspace SCF (IS3CF). The annihilation of the occupied−virtual block elements, similar to the one discussed by Stewart, was utilized by Loos, Rivail, and Assfeld.557 All available N orbitals are partitioned in their iterative stochastic subspace SCF method (IS3CF) into K sets of size L = N/K (Figure 20). The annihilation procedure between the orbitals belonging to different sets is then performed in the manner similar to that of the LMO. The choice of a particular 5851
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Figure 20. Schematic representation of the IS3CF procedure. K is the number of sets into which the original pool of N orbitals is partitioned; L = N/K is the number of orbitals in a given set. k enumerates all sets of L orbitals.
Figure 19. Schematic representation of the MO growth procedure in the LMO method of Stewart. ψ and ψ′ are the unmixed and mixed MOs, used in the growth procedure. See eq 4.2.90 and the corresponding discussion in the text.
L × Lk
Fk̃ k
= (BkN × Lk )+ FAOBkN × Lk
(compute the transformed Fock matrix)
orbital pair is stochastic and is made based on the norm of the localization measure, similar to one used in the ELG method. The annihilations are organized in a hierarchy depending on the level of subdivision of the entire set of orbitals. Occupied− virtual orbital pairs are formed at the first level of approximation. The occupied orbital, i, is chosen randomly, while the virtual one, a, is chosen as the orbital that maximizes the Frobenius norm of the product: Ci+FAOCa
F
(4.2.92)
Tr(X +X )
(4.2.93)
L × Lk
Fk̃ k
F
≡
L ×L = C̃k k kEkLk × Lk ,
∀ k = 1, ..., K
(solve the eigenvalue problem)
(4.2.94c)
L × Lk
CkN × Lk = BkN × Lk C̃k k
(transform the eigenvectors to the original basis) (4.2.94d)
The superscripts indicate matrix dimensions. Once the orbitals and their energies are computed for all sets of a given subdivision level, the total matrix of coefficients and the matrix of eigenvalues are reconstructed by concatenation of the matrices for each subset:
with X
L × Lk
C̃ k k
(4.2.94b)
At the first level of partitioning, the matrices Ci and Ca are column vectors containing expansion coefficients of the occupied and virtual orbitals in the AO basis, respectively. FAO is the Fock matrix in the AO basis. Also at the first level, the quantities in eq 4.2.92 are simply the absolute values of the elements of the Fock matrix in the MO basis. The occupied−virtual pairs (and other tuples) are combined in a similar way at the second and higher levels. The first pair (tuple) is chosen randomly, and the second one is chosen to maximize the norm, similarly in eq 4.2.92. Each of the matrices C now contains several column vectors corresponding to the orbitals included into a given group. The stochastic partitioning procedure is continued until the desired number of subsets is formed. At the highest possible level of partitioning, all N MOs are included in a single subset, yielding the conventional SCF method with an N × N Fock matrix. The dimension of the Fock matrices is L × L at the lower levels. Although the number of subsets is larger, the total computational efforts scale more favorably: O((N/L) × L3) = O(N × L2). This is linear in N for a sufficiently large N/L ratio. Once K sets of orbital pairs are constructed stochastically, eigenvalues and eigenvectors for each set are computed by the following set of transformations:
C N × N = C1N × L1 ∪ C2N × L2 ∪ ...
(4.2.95a)
E N × N = E1N × L1 ∪ E2N × L2 ∪ ...
(4.2.95b)
The resulting set of orbitals is utilized to start the stochastic partitioning scheme again, until the desired convergence is achieved. The methodology requires an initial guess of N orthonormal MOs, which can be chosen similarly to the LMO method. 4.2.4. Diabatic Approaches. Identification of an efficient system fragmentation is the key challenge of the methods discussed above. Many ways of partitioning a system into specific fragments have been discussed under different frameworks, each based on its own criterion of accuracy and efficiency. The role of partitioning is less important in the D&C methods, in which the number of electrons in each fragment is determined uniquely from the Fermi energy equilibration principle. At the same time, the partitioning scheme in the D&C method is also a matter of judicious preparation (e.g., how to choose buffer regions, overlaps, etc.). The ambiguity associated with the choice of the partitioning scheme is minimized in the group of methods that we classify as “diabatic”. These methods define a partitioning basis by considering a subspace of various possible charge localizations and molecular partitionings. Essentially, any type of basis configurations can be considered, each corresponding to a smaller subsystem. This is the beauty of the diabatic approach. The configurations may be rather artificial in terms of their
BkN × Lk ← CkN × Lk (definition of the orthogonalization matrix for a given set) (4.2.94a) 5852
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
⎧ δδ , A=B ⎪ ij σσ ′ ⟨ψA , σ , i|ψB , σ ′ , j⟩ = ⎨ ⎪ ⎩ Sij , A ≠ B
physical and chemical interpretation or may have a clear meaning. The efficiency comes from the fact that diabatic states can be computed using relatively small noninteracting fragments of a system. 4.2.4.1. Block-Localized WFT and DFT Methods. The blocklocalized WFT (BLWFT) was developed in 1998 by Mo and Peyerimhoff for studies of electronic delocalization in organic molecules.558 The method was applied to study hyperconjugation,559 and to perform energy decomposition and basis set superposition error analysis.560 The system is partitioned into K subsystems {Sk, k = 1, ..., K} with Mk AO basis states and Nk electrons in the kth group, such that the total number of electrons and AO basis functions are constant:
The expansion of the sth diabatic state, Ψs, in the overall AO basis is ⎛C1 ⎜ 0 C=⎜ ⎜ ... ⎜ ⎝0
(4.2.96a)
k=1
0 C2 ... 0
... 0 ⎞ ⎟ ... 0 ⎟ ... ... ⎟ ⎟ ... CK ⎠
Es = ⟨Ψ|s Ĥ |Ψ⟩ s =
K
∑ Mk = M
One constructs many-electron auxiliary states of each group as the product of the MOs in this group:
P = C(C +SC)−1C +
ΦI , k ( r1⃗ , ..., rN⃗ k , σ1 , σ2 , ..., σNk) = ψk ,1( r1⃗ , σ1) ψk ,2( r2⃗ , σ2) ... ψk , N ( rN⃗ k , σNk) k
∑ Ck ,a(i , σ)|χa ⟩ = ∑ Ck , σ ,ai|χa ⟩, a ∈ Sk
(4.2.97)
σ = α, β
a ∈ Sk
(4.2.98)
The requirement of expansion of MOs for the kth group only in terms of the AOs belonging to the given group is the key assumption of the block-localized WFT, providing the main computational savings of the method. We follow the definitions introduced in section 2.1. Unlike eq 2.1.12, each one-particle MO has an addition index corresponding to the subset. The states, eq 4.2.97, are many-particle basis functions localized on the kth subsystem. They need not be antisymmetric. These states are used to construct physically meaningful many-particle diabatic states, obtained as antisymmetric products of the many-particle auxiliary functions of all groups, eq 4.2.97: Ψs = A[̂ ΦI1,1ΦI2 ,2...ΦIN
K
, NK ]
(4.2.102)
(4.2.103)
The one-electron MOs of each group k, ψk,σ,i, are obtained by solving a set of K coupled Schrödinger equationsone for each group. The coupling between the equations is realized via the density-matrix-dependent Fock operator, but the dimensionality of each electronic problem is smaller than that of the entire system, leading to linear scaling. The organization of computations to find each diabatic state, Ψs, is similar to that of the D&C method. The main focus of the initial works on the BLWFT was on an efficient construction of diabatic states. The authors could then consider several meaningful configurations (partitioning schemes) to analyze reactive processes. A more systematic and general approach to reactive systems, especially to dynamics, was developed by Gao et al., who combined the BLWFT approach for the MO description with the valence bond (VB) theory, leading to the MOVB method.561−563 We refer the reader to the review on the VB theory for more details and discussion of the modern developments of the subject.564 In general, the “true” wave function is given in the VB theory by a linear combination of different resonant configurations. The configurations are chosen to be the diabatic states Ψs in the MOVB theory:
where I is the multi-index defining a particular choice of oneelectron spin−orbitals within the kth group. The one-electron spin−orbitals of the kth group are represented in the basis of the AOs belonging to this group: |ψk , σ , i⟩ =
1 Tr(P· (H + F )) 2
where H and F are the standard one-electron core Hamiltonian and the Fock matrices in the AO basis, and P is the generalized density matrix:
(4.2.96b)
k=1
(4.2.101)
where each Ck is the matrix of coefficients for the MOs of the kth group in the basis of the kth group AOs. Each matrix has dimensions Mk × Nk. We omit the state index s in these matrices for simplicity. The total energy of the determinant is given by the expression similar to the Hartree−Fock, eq 2.1.28:
K
∑ Nk = N
(4.2.100)
ΨMOVB =
∑ as Ψs s
(4.2.99)
(4.2.104)
The linear combination of configurations, eq 4.2.104, is similar to that used in the post-HF correlated methods such as the CI, CC, and CASSCF. In a manner similar to those methods, the expansion, eq 4.2.104, provides a way of incorporating nondynamic correlation. Together with a reasonable choice of the number of configurations and with the BLW approach to compute them, the method is applicable to large systems and reactive processes in a broad meaning. The expansion coefficients as are determined by diagonalization of the Hamiltonian matrix written in terms of the manyelectron configurations:
The construction, eq 4.2.99, is a coarse-grained and generalized version of the conventional Slater determinant, eq 2.1.2. The generalization comes via the antisymmetrizer  , which does not have the simple form of a determinant. Similar to the conventional Slater determinant, there are many possible “excitation” of a reference function, Ψ0. Thus, an additional index s enumerates all possible configurations. These configurations come by choosing one of many possible combinations of the multicomponent indexes I1, ..., INK. The MOs are orthonormal within each group, but they are not required to be orthogonal across different groups:
Ha = SaE 5853
(4.2.105) DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
where HIJ = ⟨Ψ|I H |ΨJ ⟩
(4.2.106)
SIJ = ⟨Ψ|Ψ I J⟩
(4.2.107)
HI,II = ⟨ψI|H |ψII⟩
(2N + 1)!
(4.2.110b)
(2N + 1)!
|φ1D(1) ... φND+ 1(N + 1) φ1A (1) ... (4.2.111a)
and vice versa: 1
ψII ≈
(2N + 1)! φÑ A+ 1(N
|φ1̃ D(1) ... φÑ D(N ) φ1̃ A (1) ...
+ 1)|
(4.2.111b)
The KS orbitals of the noninteracting fragments are, in general, nonorthogonal across the subsets belonging to different fragments. Therefore, they are orthogonalized using the Löwdin technique, which is known to preserve orbital localization character. The coupling matrix element is given by HI,II = ⟨ψI |H |ψII⟩ ≈ ⟨ψI |HIIKS|ψII⟩ = ⟨φND+ 1|hIIKS|φNA+ 1⟩ (4.2.112)
in the basis of the approximate donor and acceptor states. The last equality shows that the coupling matrix elements are determined by the singly occupied molecular orbitals (SOMOs) of the charged donor and acceptor species. A similar technique works for hole diffusion with the only exception that one has to consider positively charged species. The accuracy of the approximations was studied by comparison with the more rigorous constrained DFT calculations, validating the computational methodology.48 The methodology was applied to a variety of organic semiconducting materials.50 Although the FO-DFT is very simplistic, it can be applied to large finite and periodic systems. Many nuclear configurations can be sampled, allowing MD and systematic screening. 4.2.4.3. FMO-MS-RMD. The common idea of the diabatic approaches is well illustrated by the recently developed FMObased multistate reactive molecular dynamics (FMO-MS-RMD) of Lange and Voth.572 The method was motivated by the need to address uniqueness of the partitioning scheme. Instead of developing the “optimal” technique of partitioning a system onto fragments, the authors utilize an idea analogous to that of the multistate valence bond (MS-VB) and multiconfigurational reactive MD573 theories used in reactive force fields. The idea is also analogous to the configuration interaction approach used in the wave function theory. A set of different fragments (obtained by distinct partitioning schemes) is constructed to serve as the basis set for performing the FMO calculations. The wave function of the total molecular system, |Ψtot⟩, is approximated by a linear combination of different fragmentation “states”, |ΨS⟩:
(4.2.109)
The coupling between the initial and final diabatic states is to be computed. A straightforward approach is to construct wave functions for the diabatic states as single Slater determinants composed of the molecular orbitals of the DA complex: |φ1DA (1)...φ2DA (2N + 1)| N+1
|φ1̃ DA (1)...φ2̃ DA (2N + 1)| N+1
φNA (N )|
is one of the most important parameters determining the hopping rates. In the FO-DFT, a system of interest is composed of the donor D and acceptor A species, each containing N electrons. One extra electron (the one that is transferred) is added, to describe the electron transfer reaction:
1
1
ψI ≈
(4.2.108)
D− + A → D + A−
(2N + 1)!
This choice is a special case of the determinants made in the MOVB (or BLDFT) method discussed above. In order to avoid explicit calculations of all molecular orbitals of the combined DA system {φDA i , i = 1, ..., 2N + 1}, the diabatic states are approximated using the orbitals of the noninteracting fragments, {φDi , i = 1, ..., N + 1} and {φAi , i = 1, ..., N}:
The Hamiltonian and overlap matrix elements, eqs 4.2.105 and 4.2.107, respectively, are evaluated using the Slater determinants constructed from the nonorthogonal orbitals. As a result, the final expressions are more involved than in the standard (orthogonal MO) case. The details of the derivations can be found elsewhere.561−563 Thus, the energy evaluation does not involve computation of the adiabatic PESs. Only diabatic energies and couplings are needed, and diagonalization of a relatively smalldimensional Hamiltonian matrix, eq 4.2.106. Recently, Cembran et al. extended the MOVB theory to the DFT in the block-localized DFT (BLDFT) method.565 The method combines the best of the two worldsdynamic correlation from the DFT (via correlation functional) and nondynamics correlation from the VB theory (via multiple configurations). Together with the block-localized scheme, the method has a high potential to become a very accurate and efficient tool for studies of reactive and electronic (e.g., excited state) dynamics in large-scale systems. In passing, we would like to comment on the use of diabatic states. Diabatic states have clear physical meaning, and are convenient for interpretation and construction. They may be defined quite arbitrarily and, hence, may be used as a convenient basis for quantum dynamics, and, as a specific case, for stationary problems. For example, most of the nonadiabatic electronic dynamics methods, such as trajectory surface hopping (TSH)566,567 or Ehrenfest568−571 algorithms, can be formulated in a diabatic basis. They often show results comparable to simulations in the adiabatic basis. The fact that the basis is arbitrary means that each state may be chosen as a MO of an independent fragment. However, one may need many different types of fragments. In particular, it is crucial to include charged configurations A, A+, A−, etc. for the description of charge transfer and polarization phenomena. 4.2.4.2. Fragment Orbital DFT (FO-DFT) Method. The philosophy similar to that of the BLWFT and BLDFT methods was adopted by Oberhofer and Blumberger, who developed the fragment orbital DFT (FO-DFT) method.49 Calculation of charge carrier mobilities in organic materials using surface hopping schemes is the main target of their method. The electronic coupling between the diabatic (e.g., noninteracting donor and acceptor) states
ψI =
1
ψII =
|Ψtot⟩ =
∑ cS|ΨS⟩ S
(4.2.113)
where cS is the expansion coefficient. The authors consider proton transfer dynamics in water as an example of the reactive system which needs a multistate description (Figure 21).
(4.2.110a) 5854
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
⟨Φ1|Φ2⟩ = S12
The method can be considered a model Hamiltonian approach. Unlike the BLDFT and BLWFT methods, the matrix elements of the Hamiltonian are computed using the FMO description. The diagonal Hamiltonian matrix elements, HS,S, are defined as the FMO energies for the corresponding schemes with additional empirical corrections. The off-diagonal elements, HS,S′, correspond to a particular chemical reaction that transforms one set into the other. Once the Hamiltonian supermatrix is formed, the coefficients of the expansion, eq 4.2.113, are obtained by the Hamiltonian diagonalization. The resulting energy of the system is assigned to be the solution that delivers the minimum. Finally, the gradients can be computed using the Hellmann−Feynman theorem, to perform MD simulation or geometry optimization. 4.2.4.4. Superposition of Fragment States (SFS). The superposition of fragment states (SFS) method of Valeev et al.574 utilizes a set of fragment ionized states to represent the adiabatic state of a large-scale system. It was proposed initially as a tight-binding, model Hamiltonian approach for studies of charge transfer in organic solids. Later, it was combined with correlated methods by Zhang and Valeev, to describe ionized states of large clusters.575 The basis idea can be illustrated with the ionized state of an AB dimer. The total wave function of the ionized state, |(AB) +⟩, can be represented by a linear combination of the charge-localized diabatic states:
(4.2.116)
The matrix elements of the resulting coarse-grained (fragment basis) Hamiltonian are determined from the Koopmans theorem: ⎧ e, ⎪ i ̂ ⟨Φ|i Ĥ |Φ⟩ j − EAB ≈ −⟨ψiX |FAB|ψ j ⟩ ≡ ⎨ ⎪J , X ⎩ ij
i=j i≠j (4.2.117)
It is important to emphasize that the Fock operator of the entire system, F̂AB, is utilized to compute the above integrals. The integrals are, in general, different from those obtained for the isolated fragments, ⟨ψiX|F̂X|ψjX⟩, X = A, B. The difference may be notable and, in fact, may play a central role.575 At the same time, the computational cost for such an approach will be comparable to that for the entire system, unless an efficient and accurate separability approximation for such matrix elements is developed. This approximation can be facilitated by the preliminary transformation of the nonorthogonal diabatic basis to the orthogonal one. The Löwdin transformation is the standard choice minimizing distortion of the original wave function. Utilization of the diabatic charge-constrained states, as given by the constrained DFT (CDFT), was combined by van Voorhis and co-workers 576 with the linear-scaling code CONQUEST.87,577,578 The method requires prior knowledge of the charge localization states. These states can be generated in a number of ways. The CDFT by van Voorhis,576,579−583 dimer splitting,584,585 and the fragment orbital586−588 methods are among the most popular. In principle, one can consider a sufficiently large basis of different charge states and obtain an accurate description of the system. This approach may particularly be useful for studying charge transfer dynamics. 4.3. Direct Optimization Methods
Most of the “dynamical” methods were extensively discussed by Goedecker.589 According to his classification the methods fall into one of the following six basic classes: (1) Fermi operator expansion (FOE):590 In this approach the density matrix is constructed directly from the Hamiltonian of the system. (2) Fermi operator projection (FOP): Unlike the FOE method, the entire density matrix is not constructed; only the subspace spanned by the occupied states is searched for. (3) Divide-and-conquer method (D&C): The density matrix for the whole system is computed by combining the density matrices for the subsystems. (4) Density matrix minimization (DMM): The density matrix is determined by optimization of a specially developed functional with respect to all components of the density matrix, each treated as an independent variable. (5) Orbital minimization (OM): The density matrix is computed via the standard definition in terms of localized functions (Wannier functions, in general; eigenfunctions, in particular). The localized functions are found by optimization of the specially derived functional with respect to the expansion coefficients in terms of the basis functions. (6) Optimal basis density matrix minimization (OBDMM): It combines the ideas of the previous two methods. In addition to finding an optimal density matrix, as an expansion in terms of a specified basis, the basis itself is optimized by additional procedures.
Figure 21. Basis of four fragmentation (diabatic) “states” of the protonated water tetramer. |ψn⟩ is the fragmentation set n. Reprinted from ref 572. Copyright 2013 American Chemical Society.
|(AB)+ ⟩ = CA|Φ1⟩ + CB|Φ2⟩
(4.2.114)
|Φ1⟩ = |A+B⟩ = |A+⟩|B⟩ = aî A|A⟩|B⟩
(4.2.115a)
|Φ2⟩ = |AB+⟩ = |A⟩|B+⟩ = aî B|A⟩|B⟩
(4.2.115b)
where âiX is the annihilation operator that removes an electron from the ith orbital of the fragment X. The states |A+B⟩ and |AB+⟩ describe localization of the positive charge on the fragments A and B, respectively. Within the simplest approach, the diabatic states are approximated by the product of the Slater determinants for the individual fragments in the specified charge state. The Slater determinants can be obtained by different methods outlined below. The diabatic states are normalized, but nonorthogonal, initially: 5855
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review Norb
The D&C group of methods has already been discussed in section 3, as one of the static fragmentation methods. It is distinct in our classification, because it relies on the solution of the eigenvalue problem. On the contrary, the dynamical methods avoid solving the eigenvalue problem. Instead, they solve for the density matrix or coefficients by an iterative approach, propagating them from the initial guess to the optimized value. This is typically done by the standard minimization algorithms such as conjugate-gradient or steepest descent, although the propagated variables may be considered dynamical variables, evolving in time. In the following discussion, we will briefly overview the basic formulation of the DMM, OM, and FOE groups of methods, along with the most recent advances along those lines. In addition, we consider the group of methods that treat the wave functions (orbitals) or density matrices as dynamical variables. This group is not specifically classified by Goedecker, although it is reminiscent of the DMM and OM formulations. 4.3.1. Density Matrix. Before discussing direct the optimization methods, it is worth refining the details of the density matrix definitions and transformations under different conditions. This subsection extends the definitions presented in sections 2.1.2 and 2.1.4 to more general cases and puts additional emphasis on the density matrix variables. To simplify the notation, we utilize the spin−orbital indexing, such that each index denotes a unique combination of spin and spatial parts of an orbital. The one-particle charge density matrix, ρ, is defined in the orthonormal MO basis, {ψi}, as
ρ( r ⃗ , r ′⃗ ) =
i=1
∑ niψi*( r ⃗) ψi( r ′⃗ ) i=1
= ⟨r |(∑ ni|ψi ⟩⟨ψi|)|r′⟩ = ⟨r |ρ ̂|r′⟩
gives the density matrix in operator representation: Norb
ρ ̂ = (∑ ni|ψi ⟩⟨ψi|)
Because the MO set, {ψi}, is orthonormal, the matrix representation is given by eq 2.1.46, leading to Norb
ρ ̂ = ( ∑ |ψi ⟩Pij⟨ψj|) ⇔ Pij = Oij = niδij
This trivial result is implied by eqs 2.1.33 and 2.1.34. In short, the density matrix in the orthogonal MO basis is diagonal with the elements equal to the populations. The populations can be generalized further to finite temperatures and fractional values, via the Fermi distribution function. Consider the density matrix representation in the basis of nonorthogonal AOs. Using eq 4.3.7, and expansion eq 2.1.12, we get Norb i,j=1
(4.3.1)
i,j=1
Norb
= ( ∑ |χa ⟩Pab⟨χb |) (4.3.8a)
i,j=1
P = COC +
(4.3.2)
(4.3.8b)
which is the standard expression, eq 2.1.32. Having generalized the standard definition, we can now consider the nonorthonormal MO basis, {ψ̃ i}:
(4.3.3)
⟨ψĩ |ψj̃ ⟩ = SMO, ij
(4.3.9)
The definition of the density matrix, eq 4.3.1, remains the same, when expressed in terms of othogonalized orbitals. The common orthogonalization procedure is via the Lö wdin transformation:
(4.3.4a)
normalization:
|ψi ⟩ =
∑ |ψj̃ ⟩SMO,ji−1/2 (4.3.10)
j
(4.3.4b)
so that
idempotency: ⟨ψi|ψj⟩ =
ρ2 ( r ⃗ , r ′⃗ ) ≡
Norb
ρ ̂ = ( ∑ |ψi ⟩Oij⟨ψj|) = ( ∑ |χa ⟩(COC +)ab ⟨χb |)
hermiticity:
∫ ρ( r ⃗ , r ⃗) d r ⃗ = N
(4.3.7)
i,j=1
A proper one-particle density matrix must satisfy the following basic criteria:
ρ( r ⃗ , r ′⃗ ) = ρ( r ′⃗ , r ⃗)
(4.3.6)
i=1
The diagonal element of the density matrix gives the well-known charge density: ρ( r ⃗) = ρ( r ⃗ , r ⃗) = ρ( r ⃗ , r ′⃗ )|r ⃗ ′ = r ⃗
(4.3.5)
i=1
ni is the occupation number of the orbital i. It is either 0 or 1 in most cases (zero-temperature limit). ψi(r)⃗ is the coordinate representation of the MO |ψi⟩: ψi( r ⃗) ≡ ⟨ r |⃗ ψi ⟩
i=1
Norb
Norb
ρ( r ⃗ , r ′⃗ ) =
Norb
∑ niψi*( r ⃗) ψi( r ′⃗ ) = ∑ ni⟨r|ψi⟩⟨ψi|r′⟩
∫ ρ( r ⃗ , r ″⃗ ) ρ( r ″⃗ , r ′⃗ ) d r ″⃗ = ρ( r ′⃗ , r ⃗)
∑ SMO,ia−1/2⟨ψã |ψb̃ ⟩SMO,bj−1/2 a,b
=
(4.3.4c)
∑ SMO,ia−1/2SMO, abSMO, bj−1/2 = δij a,b
It is often convenient to work with the abstract operator or matrix representations of the density matrix. The expression
(4.3.11)
The operator representation gives 5856
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Norb
The normalization is
ρ ̂ = ( ∑ |ψi ⟩Oij⟨ψj|)
Norb
i,j=1 Norb
Norb
= (∑
∑
tr(ρ ̂) = tr( ∑ |ψã ⟩(SMO−1/2OSMO−1/2)ab ⟨ψb̃ |) a,b=1
|ψã ⟩SMO, ai−1/2OijSMO, jb−1/2⟨ψb̃ |)
Norb
i,j=1 a,b=1
∑
=(
Norb
a ′ , b ′= 1 a,b=1
= ( ∑ |ψã ⟩(SMO−1/2OijSMO−1/2)ab ⟨ψb̃ |)
Norb
a,b=1
=(∑ a,b=1
P̃ =
∑
=
Norb
|ψã ⟩Pab̃ ⟨ψb̃ |)
SMO−1/2OSMO−1/2
(SMOSMO−1/2OSMO−1/2SMO)a ′ b ′
a ′ , b ′= 1 Norb
(4.3.12a)
∑
= (4.3.12b)
(SMO1/2OSMO1/2)a ′ b ′
a ′ , b ′= 1
= Nocc
The occupation matrix, O, is diagonal, and contains 0 and 1 in the zero-temperature limit. If the basis set consists of only occupied orbitals, the expression eq 4.3.12b simplifies: P ̃ = SMO−1/2ISMO−1/2 = SMO−1
(4.3.14b)
The hemiticity, eq 4.3.4a, is trivial. 4.3.2. Wave Function or Density Matrix as Dynamical Variables. One of the first approaches to solving the electronic structure problem via a direct minimization of the total energy with respect to orbitals was developed by Car and Parrinello in 1985.591 Their approach is based on a simple classical mechanical idea of an extended Lagrangian and treatment of the MOs, |ψi⟩, and their time derivative, |ψ̇ i⟩ as coupled dynamical variables. The Lagrangian is
(4.3.12c)
The density matrix operator for nonorthonormal MOs can be expressed in the basis of nonorthonormal AOs: Norb
ρ ̂ = ( ∑ |ψã ⟩Pab̃ ⟨ψb̃ |) a,b=1 Norb
L=K−E+
̃ +)ab ⟨χ |) = ( ∑ |χa ⟩(CPC b
∑ Λij(⟨ψi|ψj⟩ − δij)
Norb a,b=1
P = CSMO−1/2OSMO−1/2C +
where the last term imposes the orthonormalization restriction on the orbitals. K and E are the kinetic and potential energies, respectively:
(4.3.13a) (4.3.13b)
K=
The expression eq 4.3.13b can be transformed further in the zerotemperature limit to
C
which is the expression used earlier in the context of the BLWFT and BLDFT methods. One can verify the three main properties of the density matrix, eqs 4.3.4, as applied to the operator representations discussed above. For example, the idempotency, O2 = O, is satisfied in the nonorthogonal MO case, eq 4.3.13a: ρ 2̂ = ( ∑ |ψã ⟩(SMO−1/2OSMO−1/2)ab ⟨ψb̃ |) Norb
|ψã ′⟩(SMO−1/2OSMO−1/2)a ′ b ′⟨ψb̃ ′|)
Norb a,b=1 a ′ , b ′= 1
|ψã ⟩(SMO−1/2OSMO−1/2)ab SMO, ba ′
(SMO−1/2OSMO−1/2)a ′ b ′⟨ψb̃ ′| Norb
=
∑
|ψã ⟩(SMO−1/2O2SMO−1/2)ab ′⟨ψb̃ ′|
a , b ′= 1
= ρ̂
∑ μi ⟨ψi̇ |ψi̇ ⟩ i
∑ ⟨ψi|∇i 2 |ψi⟩ + U[ρ] i
(4.3.16)
(4.3.17)
δE + δψi*
(4.3.18a)
∑ Λijψj j
(4.3.18b)
Both electrons and nuclei are treated on the same footing. Using eqs 4.3.18, one performs the simulated annealing MD calculations that bring the kinetic energy to a finite temperature. In the limit of zero temperature, the time derivative of orbital velocity approaches zero, ψ̈ i = 0, meaning that the stationary state is found and that the variational procedure converged to a solution given by the matrix diagonalization. In addition, the dynamics converges to the solution that corresponds to orthonormalized orbitals. The expensive matrix diagonalization procedure is substituted with the much more efficient simulated annealing dynamics, allowing modeling of large molecular or solid-state systems.
a ′ , b ′= 1
∑
1 2
μi ψï = −
a,b=1
=
n
1 2
MnR̈ n = −∇R n E
Norb
∑
∑ M n Ṙ n 2 +
where μi is the fictitious mass associated with the orbital |ψi⟩, Mn is the mass of the nth nucleus, Ṙ n is the velocity of the nth nucleus, and U[ρ] is the electronic DFT potential energy that includes all contributions other than the kinetic energy. The Lagrangian, eq 4.3.15, generates coupled dynamics of ionic particles (nuclei) and fictitious electronic degrees of freedom, represented by orbitals:
(4.3.13c)
(
1 2
E=−
P = CSMO−1/2OSMO−1/2C + = CSMO−1C + = C(C +SAOC)−1 +
(4.3.15)
i,j
a,b=1
= ( ∑ |χa ⟩Pab⟨χb |)
⟨ψã ′|ψã ⟩(SMO−1/2OSMO−1/2)ab ⟨ψb̃ |ψb̃ ′⟩)
(4.3.14a) 5857
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
In contrast to the CP method, a strict optimization is performed for each single-point calculation. The basic formulations were given by Mauri et al.600,601 and Ordejón et al.602,603 The method directly minimizes the energy by propagating the MO expansion coefficients using the standard techniques, such as the conjugategradient optimization. The derivation of the OM methods can be systematized using the notation introduced in section 4.3.1. The energy of the system is given in the MO basis by
The simulated annealing may require a large number of steps, creating one of the difficulties associated with the Car−Parrinello (CP) dynamics. Because the fictitious mass of the orbital, μ, is rather small, the integration time step may need to be reduced considerably, if stable dynamics is to be obtained. Both these factors may, in principle, become obstacles, leading to the poorer-than-expected performance boost. Yet anther difficulty is the possibility of converging into one of the local minima. This may become more probable as the size of the system increases, and the number of shallow local minima grows significantly. Using sufficiently slow simulated annealing cooling protocols may help in avoiding the problem, but will require increased computational time. Further improvements of the CPMD method were reported soon after the original work. For example, Marx et al.592,593 developed the extensions to add nuclear quantum effects via path integrals. Doltsinis and Marx594 extended the method to the nonadiabatic dynamics, by combining it with the restricted openshell KS (ROKS) method595 and Tully’s fewest switches surface hopping algorithm.566 To allow even further improvement in computational efficiency, Galli and Parrinello reformulated the method in a nonorthogonal basis.596 The use of a nonorthogonal basis is needed for localization of orbitals and separation of distant localization regions, to yield a linear-scaling approach. This generalization has several important points. First, one introduces the auxiliary function, defined similarly to eq 4.3.10, but with a different power of the overlap matrix:
|ψi̅ ⟩ =
Eel = tr(PMOHMO)
where PMO is the density matrix, and HMO is the Hamiltonian matrix, both given in the MO basis. In the general case of nonorthonormal MOs, in zero-temperature limit, the density matrix is given by the inverse of the MO overlap matrix, eq 4.3.12c. Thus, the energy is given by Eel = tr(SMO−1HMO)
SMO−1 = (I − (I − SMO))−1 ≈ I + (I − SMO) + (I − SMO)2 + ...
Ω = tr(IHMO + (I − SMO)HMO)
Ω = tr(HMO − Λ(SMO − I ))
∑ |ψi⟩Vi ,S⟨ψi̅ | i
Vi , S = A[1 − θ(|r − RS| − σ )]
(4.3.28)
The stationary solution
∂Ω =0 ∂Λ
(4.3.29)
corresponds to the orthonormality: SMO = I. It is achieved for Λ = HMO. Thus, the grand potential for the unconstraint optimization is
(4.3.21)
Finally, to force the MO localization in different spatial regions, a nonlocal potential is added to the Hamiltonian. The potential takes the form similar to that of the core pseudopotentials common to DFT:
V=
(4.3.27)
Alternatively, the orthonormality constraint can be imposed via the Lagrange multipliers method, for example as proposed by Ordejón et al.602,603 In this case, the grand potential is
(4.3.20)
The resulting equation does not contain the orthonormalization term ∑j Λijψj, because the orthonormalization is implicitly incorporated in the definition, eq 4.3.19. Second, to avoid the need for the overlap matrix inversion in eq 4.3.19, the iterative scheme that computes the auxiliary functions in linear time is used: |ψi̅ ⟩(n + 1) = (1 − SMO)|ψi̅ ⟩(n) + |ψĩ ⟩
(4.3.26)
With only the first-order terms, eq 4.3.25 becomes the grand potential:
(4.3.19)
With the new function, the equation of motion for electronic variables, eq 4.3.18b, becomes
δE μi ψï = − δψi̅ *
(4.3.25)
The true solution corresponds to a set of orthonormal MOs, in which case the overlap matrix is diagonal. The orthonormalization of the orbitals is essentially equivalent to the matrix diagonalization and should be avoided if a linear-scaling method is to be formulated. To mitigate the need for the matrix inversion, Mauri et al.600,601 approximated the inverse of the overlap matrix with its Taylor expansion terminated at an odd power:
∑ |ψj̃ ⟩SMO,ji−1 j
(4.3.24)
Ω = tr(HMO − HMO(SMO − I ))
(4.3.30)
The result is exactly the same as that obtained by Mauri et al.600,601 using the first-order Taylor expansion, eq 4.3.27. The constraining parameters Λ in the grand potential of the form given in eq 4.3.28 were used explicitly in the formulation of Wang and Teter.604 They argued that a relatively small magnitude of these parameters (larger but comparable to the magnitude of the Hamiltonian itself) can be used to formulate an efficient gradient molecular optimization strategy for linear scaling. The formulations of Mauri et al. and Ordejón et al. utilized only occupied orbitals. The extension to a more general basis, which can include unoccupied orbitals, was done by Kim et al.605 The extension can be derived using the general constructs and transformations introduced in section 4.3.1. The grand potential can be written in the general case simply as
(4.3.22) (4.3.23)
where A is the strength of the potential, θ(·) is the step function, and RS and σ are the center and the radius of the localization region S, respectively. More recently, an ab initio density matrix propagation (ADMP) method was formulated. It is similar to the CPMD, but uses the density matrix elements as the main dynamical variables.597−599 4.3.3. Orbital Minimization (OM) Methods. The OM methods are very closely related to the CP and similar methods. 5858
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews Ω = tr((Ĥ − εFI )̂ ρ ̂) + εFN
Review
Ω = tr(P(̃ H − εFI )) = tr((3P 2 − 2P 2)(H − εFI ))
(4.3.31)
(4.3.36)
where the last term accounts for the normalization condition of the density matrix, eq 4.3.4b. 4.3.4. Density Matrix Minimization (DMM) Methods. The DMM was formulated in 1993 independently by Li, Nunes, Vanderbilt,606 and Daw.607 In the early works, the scheme was applied mostly to one-element systems at the tight-binding (TB) level of description, which already provided some computational advantages for large-scale calculations of the time. The use of the simple tight-binding or semiempirical schemes seems to be a natural choice for the initial stages of developing a linear-scaling algorithm. This common trend appeared in the discussion of the fragmentation-based approaches. Such choices allow the researcher to focus on the algorithm itself rather than on details of the potential. We illustrate the idea of the DMM methods using the early examples with the tight-binding Hamiltonians. The Fock matrix in eq 4.2.28b reduces to the core Hamiltonian under the tightbinding approximation. Consequently, the whole Hamiltonian is given by the core Hamiltonian, H, independent of the density matrix. The total energy is then simply
Finally, the energy expression, eq 4.3.36, can be optimized using any standard minimization algorithm, such as conjugate gradient or steepest descent, using the gradient of the energy with respect to the variational degrees of freedom: δΩ = 3(PH′ + H′P) − 2(P 2H′ + PH′P + H′P 2) δP (4.3.37)
H′ = H − εFI
The method is variational and can be used to compute energy derivatives with respect to arbitrary parameters λ using the Hellmann−Feynman theorem: ⎛ dH ⎞ dΩ ⎟ = tr⎜P ̃ ⎝ dλ ⎠ dλ
∑
The same results are obtained by Daw from slightly different grounds. He started by considering a smeared step function of the Fermi−Dirac type: θ(β , H′) ≡ θβ ≡
PabHab = tr(PH )
a,b=1
(4.3.32)
(4.3.41)
The above expression can be integrated to obtain θ(β,H′) starting from θ(0,H′). The zero-temperature limit of the function eq 4.3.40 is the density matrix:
(4.3.33)
ρ( r ⃗ , r ′⃗ ) = lim θβ
(4.3.42)
β →∞
The result given by eq 4.3.42 populates the energy levels below the chemical potential (Fermi energy). Using eqs 4.3.41 and 4.3.42, one can obtain the equation of motion for the density matrix: ∂ ρ( r ⃗ , r ′⃗ ) = H′ρ( r ⃗ , r ′⃗ )(1 − ρ( r ⃗ , r ′⃗ )) ∂β
(4.3.43)
Symmetrization of eq 4.3.43 gives eq 4.3.37: ∂ρ 1 1 = (H′ρ + ρH′) − (H′ρ2 + ρH′ρ + ρ2 H′) ∂β 2 3 (4.3.44)
A simplified version of the DMM was developed by Challacombe.610 The method uses several direct minimization steps, followed by the McWeeny purification procedure. The latter shows a quadratic convergence and may help in accelerating computations. The purification procedure requires a smaller number of floating point operations and can be more efficient than a direct minimization. Very recently, different purification schemes were benchmarked and also compared to the gradient minimization.611 The authors show that the purification procedure is much more efficient and can lead to at least an order of magnitude acceleration of calculations. Similarly, Daniels and Scuseria612 found, using a different set of molecular systems containing up to 20 000 atoms, that purification of the density matrix gave the most efficient
(4.3.34)
where P̃ is the density matrix purified from P by one iteration step. A few purification iterations are typically needed after each gradient minimization step. The authors argue that the purified density matrix represents a physically meaningful and converged quantity and that it can be used to compute the expectation values of any arbitrary operator, Â :
̃ ̂) A = tr(PA
(4.3.40)
∂ θβ = H′θβ(1 − θβ) ∂β
Satisfying the idempotency, eq 4.3.4c, is the most difficult challenge. The idempotency condition imposes constraints on the populations of the occupied and unoccupied levels, given by the eigenvalues of the density matrix. Because the DMM does not utilize orbitals explicitly, the occupied and unoccupied states are defined as the energy levels below and above the Fermi energy, respectively. Disregarding the idempotency may allow the occupation numbers of the levels below the Fermi energy to approach +∞, and the occupation numbers of the levels above the Fermi energy to approach −∞. A solution is given by the purification procedure described by McWeeny.608,609 Provided that a reasonable initial guess is used, such that the eigenvalues of the density matrix are confined in the interval [−1/2, 3/2], the following purification procedure will converge the eigenvalues to 0 or 1: P ̃ = 3P 2 − 2P 3
1 1 + exp(βH′)
where β is the reciprocal temperature parameter. The derivative of the function θ(β,H′) with respect to this parameter is expressed in terms of the function itself:
One wishes to minimize the energy functional, eq 4.3.32, with respect to the density matrix. The minimization should maintain the basic properties of the density matrix. The hermiticity is satisfied by the definition of the density matrix. To satisfy the normalization constraint, one considers the unconstrained minimization of the extended (grand) energy: Ω = Eel − εFN = tr(P(H − εFI ))
(4.3.39) 607
Norb
Eel =
(4.3.38)
(4.3.35)
The unpurified density matrix is regarded as an auxiliary (trial) variational parameter. With this assumption, the grand energy of the system, eq 4.3.33, is redefined: 5859
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
procedures. Unlike Rudberg and Rubensson,611 they found that the Fermi operator expansion method was somewhat slower than the conjugate gradient DMM. Larsen et al.613 avoided the need for the purification in the DMM by using exponential parametrization of the atomic density matrix, a set of nonorthogonal projection operators, and a multilevel preconditioning during the optimization. A similar approach was used by Matsuoka et al.,614 who formulated a direct density matrix optimization procedure based on iterative unitary transformations. A proper preconditioning is especially important when the basis set size is large, which leads to the overlap matrices to be ill-conditioned. Mostofi et al.615 showed that this problem is related to kinetic energy contribution, and they proposed a correction scheme. Recently, Sałek et al.616 used the preconditioned conjugate gradient optimization of the density matrix to obtain linearscaling HF and KS-DFT implementations. The method was further extended by augmenting with the linear-response theory, to compute excited state properties.617 The frequency-dependent polarizabilities were computed for polyalanine peptides with up to 1400 atoms, using hybrid functionals. One may observe that the method avoids expensive matrix diagonalization procedures. However, matrix multiplication is also an expensive operation that scales cubically with the matrix size. Thus, by itself the DMM method provides no advantage over the standard methods, in terms of computational efficiency. To make the approach linear scaling, one must utilize a distance cutoff for the density matrix elements and sparse matrix techniques. The cutoff sets matrix elements to zero if the corresponding orbitals are separated by a distance larger than a predefined cutoff value. This makes the density and Hamiltonian matrices sparse, especially for large systems. The matrix information is then stored in a compact and convenient form to reduce memory requirements and CPU time needed to multiply matrices. Another complication is associated with density matrix truncation beyond the cutoff distance. The normalization condition, eq 4.3.4b, may be violated as a result of the truncation. This error is usually small for semiconductors with a large band gap, but may be critical to metallic systems. Qiu et al.618 suggested modifying the chemical potential (Fermi energy) to account for the changed number of electrons. They performed a two-step optimization, alternating the conjugate gradient optimization of the density matrix and the Fermi energy adjustment. The method was utilized to perform molecular dynamics studies of a large carbon-containing system. The authors also suggested utilization of the trial density matrix, P, as a dynamical variable within the Car−Parrinello extended Lagrangian method. The DMM method was generalized to nonorthogonal basis in later works.619 The derivations can be put in the framework introduced in section 4.3.1 of the present review. Millam and Scuseria620 and Daniels and Scuseria71 performed the first applications of the DMM method to very large chemical systems. The approaches utilized semiempirical model Hamiltonians. Minimization was performed with the conjugate gradient technique. The approach was applied initially to a polyglycine chain of about 500 atoms and a water cluster of about 900 atoms.620 Even larger examples were demonstrated later, including a 20 000 atom linear polyglycine chain, an RNA molecule containing ca. 6300 atoms, and a 1800 water cluster.71 The authors investigated different cutoff schemes and found that a soft cutoff allows one to obtain more accurate results and to
avoid unphysical oscillations of the system energy with respect to the cutoff distance. However, the soft cutoff scheme leads to quadratic scaling due to accounting for pairwise electrostatic interactions. In further work, Li et al.621 developed a quasiNewton alternative to the conjugate gradient version. It relied on the direct inversion in iterative space (DIIS) method and was found to be 30% faster than the conjugate gradient version, when applied to polyglycine chains containing 10−100 residues. Simulations of a 300 molecule water cluster showed better convergence. Although the DMM by itself is a linear-scaling technique, the computational costs increase if a very large basis is used. One needs to reduce the number of basis states, or the dimensionality of the density matrix. This is achieved in the group of methods that can be referred to as the optimal basis density matrix minimization (OBDMM) group, according to Goedecker’s classification. For instance, Hernández et al.622,623 simultaneously optimized both the expansion coefficients and the basis set for representation of the density matrix. Independently, Hierse and Stechel624 developed a similar approach. The adaptive strategy was used by Berghold et al.,625 who combined the polarized atomic orbital (PAO) method needed to reduce a large basis set into a smaller basis, and the DMM approach to achieve linear scaling. The authors also compared the performance of different linear-scaling techniques and showed that efficiency increases in the sequence DMM < FOE < canonical purification of the density matrix.626 The DMM method was implemented in a number of codes. For example, Rudberg et al.627 implemented it in the Ergo code and showed linear-scaling performance for the systems with up to 300 000 basis functions with hybrid DFT methods. Hine et al.628 used their modified version in the ONETEP code, and reported simulations with over 10 000 atoms. Challacombe et al.610 presented the FreeON (former MondoSCF) code. 4.3.5. Fermi Operator Expansion (FOE) Methods. The FOE and FOP methods developed by Goedecker et al.590,629 constitute yet another group of methods. Similar to the DMM approaches, the FOE and FOP work directly with the density matrix. However, a projection step is used instead of the matrix minimization. This eliminates the DMM and OM problems associated with multiple local minima. The method relies on a representation of the density matrix by a finite-temperature Fermi matrix: FεF , β = f (β(H − εFI ))
(4.3.45)
where f(·) is the Fermi−Dirac distribution function, eq 2.1.49b. The function of the matrix is defined via its series expansion. Because eq 4.3.45 is a finite temperature approximation of the density matrix (e.g., see eq 4.3.42), the electronic energy is given by E = tr(FεF , βH ) =
∑ ⟨ψi|Fε̂ ,β|ψj⟩⟨ψj|Ĥ |ψi⟩ F
(4.3.46)
i,j
If {ψ} is the basis of eigenfunctions of the system Hamiltonian, Ĥ , the matrix H is diagonal with elements:
⟨ψi|Ĥ |ψj⟩ = εiδij
(4.3.47)
As a result, the matrix FεF,β is also diagonal in this basis: ⟨ψi|Fε̂ F , β|ψj⟩ = f (β(εi − εF))δij 5860
(4.3.48) DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
interactions, including polarization. We have already discussed some early schemes to account for such effects, including various charge and electronegativity equalization methods and fluctuating charge models. Those approaches can be considered coarsegrained versions, with no subatomic electrostatics resolved. The quantal force fields extend these methods to add subatomic details of charge distributionvia molecular orbitals. 4.4.1. Basics of QM/MM. 4.4.1.1. Energy Partitioning. The hybrid quantum mechanics/molecular mechanics (QM/MM) scheme is the earliest version of the quantal force fields. The distinction between the QM and MM regions is more definite here than in later force fields of this sort. The method was first introduced by Warshel and Levitt in 1976.632 The literature on this subject is very diverse and broad. We refer the reader to the general methodological review633 for more details on the method. Various specialized reviews discuss applications, including recent progress on the simulation of enzymatic and organic reactions with QM/MM,93,634,635 combinations of the QM/MM method with the DFTB and its variants,267,268,273 charge transfer in disordered organic semiconductors,636 and the role of electrostatics in polarizable force fields and QM/MM.133 Partitioning of a system into two regionsquantum and classicalconstitutes the main idea of the QM/MM method. The total energy of the hybrid system is represented as
The Fermi matrix is approximated in an arbitrary basis by a matrix polynomial, p(H), of the Hamiltonian matrix, H: Npl
FεF , β ≈ p(H ; εF , β) =
∑ ciPi(H ) i=0
(4.3.49)
where Npl is the power of the polynomial expansion, P(·) is a basis polynomial, and {ci} are the expansion coefficients. Pi(H) = Hi is the most trivial choice of the polynomial functions.590 Chebyshev polynomials were found to produce more stable results.629−631 The expansion coefficients are chosen to satisfy the expansion in any basis. From the practical point of view, the simplest choice of the coefficients is based on approximation of the Fermi matrix in the basis of eigenfunctions of the Hamiltonian, which is simply the Fermi−Dirac distribution function, eq 4.3.48. The Fermi−Dirac function is approximated in the interval spanned by the lowest, εmin, and highest, εmax, eigenvalues of the system Hamiltonian. The Fermi−Dirac function must be approximated only in those regions that show non-negligible density of states. In particular, the Fermi− Dirac function may be rather arbitrary in the gap region of large band gap semiconductors. Therefore, simplified functions such as629 f (ε) =
⎛ ⎛ ε − ε ⎞⎞ 1⎜ F ⎟⎟ ⎜⎜ 1 − erf ⎟⎟ 2 ⎜⎝ ⎝ Δεgap ⎠⎠
E = EQM + EMM + EQM/MM (4.3.50)
(4.4.1)
where the terms on the right-hand side are the energies of the quantum mechanical and classical regions, and the energy of interaction between them, respectively. The former two energies are well-defined. The biggest challenge is to define the last term. Warshel and Levitt computed the latter by considering the diagonal elements of the fully quantum Fock matrix. It can be written explicitly in the zero differential overlap (ZDO) approximation:
can be used to derive expansion coefficients. Finally, the above calculations require knowledge of the properties εmin, εmax, and εF. These are computed not directly, because it would require matrix diagonalization, but via auxiliary procedures.630,631 To show a linear scaling of time, all matrix multiplications must utilize sparse matrix representation and algebra. The sparsity of the matrix follows from the decay properties of the off-diagonal elements, which are set to zero beyond some distance cutoff.
Faa = Faa ,QM −
∑
Q BγAB ,
a∈A
B≠A B ∈ MM
4.4. Quantal Force Fields
The methods classified into this group typically combine several strategies for faster computations, including those discussed in the previous sections. For this reason, we consider them separately, under a different level of classification. The main idea of the quantal force fields is to retain a quantum-mechanical description, as least for a part of the system, but to achieve this via a computationally efficient scheme. Often, a scheme can be formulated using classical force field terms, e.g., electrostatics and dispersion. The methods of this group can combine semiempirical, ab initio, or classical force fields elements, more general physically based approximations, as well as fragmentation and other linear-scaling approaches. In principle, any semiempirical method can already be considered a quantum force field, in the sense that it includes a handful of physically or computationally motivated terms. The language of classical force fields, such as partial atomic charges, dipole moments, and dispersion, is also used in some of the methods discussed above. For example, the MNDO utilizes the multipole expansion for more efficient calculation of twoelectron integrals. Various core repulsion terms, accounting for dispersion interactions, were introduced in the EHT and PMn methods. From the alternative viewpoint, the quantal force fields can be considered derivatives of the classical ones. The main focus in such generalizations is put on a description of electrostatic
Faa ,QM
(4.4.2a)
⎡ ⎢ 1 = ⎢Uaa + Paaγaa − 2 ⎢ ⎣
∑ Pbbγab − ∑ b∈A b≠a
B≠A B ∈ QM
⎤ ⎥ Q BγAB⎥ ⎥ ⎦ (4.4.2b)
where Faa,QM describes only the contributions present in the isolated QM region. The second term in eq 4.4.2a describes interaction of orbitals in the QM region with the classical-point charges, Q, and atoms belonging to the MM region. Using eqs 4.4.2 as the starting ground, one develops an asymptotic expression for the Coulomb integrals, γAB ≈ 1/rAB, and adds the induced dipole interaction: Faa ,QM/MM = Faa ,QM −
∑ B≠A B ∈ MM
QB rAB
+ Find, aa ,
a∈A (4.4.3)
The inductive term, Find,aa, describes the self-consistent interaction of the induced dipoles in the MM region. The induced dipole on each atom A is computed in this region as
μA⃗ = αAE ⃗(RA⃗ ) 5861
(4.4.4a) DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews E ⃗(RA⃗ ) = −∑ B
Review
QB RAB
∑
R⃗ − 3 AB
B ∈ MM B≠A
⎛ μ ⃗ R⃗ ⎞ AB ∇⃗⎜⎜ B 3 ⎟⎟ R ⎝ AB ⎠
molecular electrostatic potentials from the bond fragments in the original formulations.641−646 The idea gets the most popularity in the QM/MM context, especially for treating the boundary of the quantum and classical regions. The concept of the fragment self-consistent method (FSCF) is a special version of the SLMO.640,647−649 The FSCF is implemented in the Sybyl program by the Náray-Szabó group using the NDDO Hamiltonian for quantum systems. Applications included enzyme reactions,650,651 amorphous solid systems,652−654 and liquid phases.655 The FSCF method is formulated as follows. First, the system is divided into three parts: core (C), polarizable (P), and nonpolarizable (N) regions. The total energy is expressed as
(4.4.4b)
where α is the atomic polarizability, and E⃗ (R⃗ A) is the local electric field at the location of atom A. The set of eqs 4.4.4 is solved iteratively for induced dipoles in the MM region. The contribution to the Fock matrix due to induced dipoles, Find,aa, can then be computed as ⃗ μB⃗ RAB
∑
Find, aa =
B ∈ MM B≠A
RAB3
a∈A
,
(4.4.5)
Finally, the total energy of interaction between the QM and MM regions is given by
∑
EQM/MM =
Q AQ B RAB
A ∈ QM B ∈ MM
∑
Edisp =
A ∈ QM B ∈ MM
E ind = −
1 2
+ Edisp + E ind
⎛ A B ⎞ ⎜ AB12 − AB6 ⎟ RAB ⎠ ⎝ RAB
∑ A ∈ QM B ∈ MM
QA
Etot = Eel +
(4.4.6b)
ψπ , l =
∑ cA(pz )A ,l A
Q AQ B
A ,B
RAB
(4.4.8)
Eel = EC + EP + EC − P
(4.4.9)
Keeping in mind the definitions eqs 2.1.28−2.1.30, the electronic energy of the core (QM) region is given by NAO
(4.4.6c)
EC =
The dispersion term, eq 4.4.6b, is introduced empirically, similarly to how it is treated with a standard classical force field model. 4.4.1.2. Boundary Treatment and Strictly Localized Orbitals. The treatment of the boundary region between the QM and MM subsystems is one of the most difficult parts of the QM/MM methodology. Although the basic electrostatic coupling scheme presented above is often accurate for distant parts of the system and weak coupling, it may be insufficient for a description of covalent interactions. Certain care must be taken for the orbitals of the QM region in this case. One of the main approaches to handle the difficulties associated with the partition of the QM and MM regions via a covalent bond is to use strictly localized MOs (SLMOs). The concept of SLMO appeared starting from the early 1950s.637−639 In the early 1980s Náray-Szabó proposed to utilize the strictly localized MOs (SLMOs) for calculations of large systems.640 The simplest set of strictly localized MOs can be constructed on the basis of conventional chemical concepts of bond types. Namely, the lone pair, two-center σ-bond, and many-center π-orbitals can be constructed from atomic orbitals as ψLP, l = χA , l (4.4.7a) ψσ , l = cAχA , l + cBχB , l
∑
where the summation in the second term includes classical electrostatic interaction of atomic cores (in the C and P regions) and point charges (in the N region). The first term represents the electronic energy of the quantum part (the C and P regions). It is decomposed as
(4.4.6a)
⃗ μB⃗ RAB RAB3
1 2
∑ a,b=1 a,b∈C NAO
=
∑ a,b=1 a,b∈C
1 PabHab + 2 1 PabHab + 2
NAO
∑
[Pα , abGα , ab + Pβ , abGβ , ab]
a,b=1 a,b∈C NAO
NAO
∑ ∑
[Pα , ab(Pcd(ab|cd)
a,b=1 c ,d=1 a,b∈C c ,d∈C
− Pα , cd(ad|cb)) + Pβ , ab(Pcd(ab|cd) − Pβ , cd(ad|cb))] NAO
=
∑ a,b=1 a,b∈C
1 PabHab + 2
NAO
NAO
∑ ∑
[PabPcd(ab|cd)
a,b=1 c ,d=1 a,b∈C c ,d∈C
− (Pα , abPα , cd + Pβ , abPβ , cd)(ad|cb)]
(4.4.10)
The expression simplifies for closed-shell systems, Pσ = (1/2) P, σ = α, β: NAO
EC =
∑
PabHab +
a,b=1 a,b∈C
1 2
NAO
NAO
∑ ∑
PabPcd(ab||cd)
a,b=1 c ,d=1 a,b∈C c ,d∈C
(4.4.11a)
with ⎡ ⎤ 1 (ab||cd) = ⎢(ab|cd) − (ad|cb)⎥ ⎣ ⎦ 2
(4.4.7b)
(4.4.11b)
The core integrals, Hab, a, b ∈ C, include interactions with all atomic cores and point charges in the system. The energy of the polarizable region is computed by assuming that both the core one-electron Hamiltonian and the density matrix have blockdiagonal structure that corresponds to nonzero elements only for the hybrid AOs belonging to the same SLMO:
(4.4.7c)
The SLMO approach is very similar to the LMO method described in section 4.2.3.2. It is also used implicitly in the ELG method, when fragment localization is performed. We refer the reader to the corresponding sections as well as to the original literature for the details of the method. It is important to note for our purpose that the method serves as the backbone of an efficient computational strategy for calculations on large molecular systems. It was utilized successfully for computing
NAO
EP =
∑ ∑ α∈P a,b=1 a,b∈α
PabHab +
1 2
NAO
NAO
∑ ∑ ∑
PabPcd(ab||cd)
α ,β∈P a,b=1 c ,d=1 a,b∈α c ,d∈β
(4.4.12) 5862
DOI: 10.1021/cr500524c Chem. Rev. 2015, 115, 5797−5890
Chemical Reviews
Review
Finally, the interaction between the two regions is described by EP =
1 2
NAO
F μSv = Hμv +
λ,η
NAO
∑ ∑ ∑
PabPcd(ab||cd)
a,b=1 α∈P c ,d=1 a,b c ,d∈α
Fab = Hab +
1 2
⎡
∑ Pll⎣⎢(μv|ll) −
+ (4.4.13)
l
∑ ∑ ∑
⎤ 1 (μl|lv)⎥ + ⎦ 2
∑ Q A(μv|sAsA) A∈E
The total energy is then given by the sum of the subsystem, environment, and subsystem/environment interface terms:
NAO
α∈P a,b=1 c ,d=1 a,b∈α c ,d∈β
⎤ 1 (μη|λv)⎥ ⎦ 2
(4.4.15)
Variation of the total energy yields the Fock matrix of the form: NAO
⎡
∑ Pλη⎣⎢(μv|λη) −
PabPcd(ab||cd)
E = ES + EE + EES (4.4.14)
1 2
ES =
The method evolved into the local self-consistent (LSCF) approach656 and was later used as a basis to construct a hybrid classical quantum force field (CQFF).657 Both works applied the novel technique to study proton transfer in biological systems. The approach was implemented in the GEOMOP package. The methodology was developed by the group of Jean-Louis Rivail and is currently being pushed forward by the Xavier Assfeld group. A detailed recent account of the intricacies of the methodology is available.658 In the LSCF method, the entire system is divided into two regions: the quantum subsystem, S, to be treated quantum mechanically, and its environment, E, described with the classical force field (Figure 22). The main difficulty arises at the interface of the two regions. The interface is typically chosen via one or multiple single bonds X−Y. The atom X belonging to the quantum system is called a frontier atom. The AOs localized on this atom are assumed to consist of one s and three p-type functions. Hybridization of the four AOs yields four localized hybrid MOs. One of these MOs, |l⟩, is then frozen and is rotated to align with the vector connecting the two atoms, to represent the bonding orbital along the X−Y bond. The other three orbitals, |i⟩, |j⟩, and |k⟩, are used as the basis for variational calculations of the quantum subsystem S. The quantum subsystem has to be orthogonal to the frozen orbital |l⟩. There may be an arbitrary number, f, of distinct contacts between the MM and QM regions, in general. Thus, if the total number of orbitals in the full QM region is N, the variational subspace consists of N − f orbitals. Similarly to the standard QM/MM scheme, the electrostatic contributions due to point charges in the MM system are included in the Fock matrix formation:
(4.4.16a)
∑ Pμv[Hμv + Fμv] + ∑ A , A ′∈ S A
