Fragmentation: Toward Accurate Calculations on Complex Molecular Systems

Fragmentation: Toward Accurate Calculations on Complex Molecular Systems introduces the reader to the broad array of fragmentation and embedding methods that are currently available or under development to facilitate accurate calculations on large, complex systems such as proteins, polymers, liquids and nanoparticles. These methods work by subdividing a system into subunits, called fragments or subsystems or domains. Calculations are performed on each fragment and then the results are combined to predict properties for the whole system.

Topics covered include:

• Fragmentation methods
• Embedding methods
• Explicitly correlated local electron correlation methods
• Fragment molecular orbital method
• Methods for treating large molecules

This book is aimed at academic researchers who are interested in computational chemistry, computational biology, computational materials science and related fields, as well as graduate students in these fields.

• Language: English
• Publisher: Wiley
• Release date: Aug 4, 2017
• ISBN: 9781119129264

Book preview

List of Contributors

Emily A. Carter

School of Engineering and Applied Science, Princeton University, USA

Garnet K.L. Chan

Frick Chemistry Laboratory, Department of Chemistry, Princeton University, USA

Ajitha Devarajan

Office of University Development, University of Michigan, USA

Johannes M. Dieterich

Department of Mechanical and Aerospace Engineering, Princeton University, USA

Dmitri G. Fedorov

Research Center for Computational Design of Advanced Functional Materials (CD-FMat), National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba, Japan

Alexander Gaenko

Advanced Research Computing, University of Michigan, USA

Kandis Gilliard

Department of Chemistry, University of Illinois at Urbana–Champaign, USA

Mark S. Gordon

Ames Laboratory of United States Department of Energy, USA

Department of Chemistry, Iowa State University, USA

Pradeep K. Gurunathan

Department of Chemistry, Purdue University, USA

Xiao He

School of Chemistry and Molecular Engineering, East China Normal University, China

NYU-ECNU Center for Computational Chemistry, NYU Shanghai, China

So Hirata

Department of Chemistry, University of Illinois at Urbana–Champaign, USA

Jan H. Jensen

Department of Chemistry, University of Copenhagen, Denmark

Carlos A. Jiménez-Hoyos

Frick Chemistry Laboratory, Department of Chemistry, Princeton University, USA

K. V. Jovan Jose*

Department of Chemistry, Indiana University, USA

Murat Keçeli

Department of Chemistry, University of Illinois at Urbana–Champaign, USA

Argonne National Laboratory, USA

Kazuo Kitaura

Graduate School of System Informatics, Kobe University, Japan

Christoph Köppl

Institute for Theoretical Chemistry, University of Stuttgart, Germany

Caroline M. Krauter

Department of Mechanical and Aerospace Engineering, Princeton University, USA

Jinjin Li

Department of Chemistry, University of Illinois at Urbana–Champaign, USA

National Key Laboratory of Science and Technology on Micro/Nano Fabrication, Department of Micro/Nano Electronics, Shanghai Jiao Tong University, China

Jinfeng Liu

School of Chemistry and Molecular Engineering, East China Normal University, China

Qianli Ma

Institute for Theoretical Chemistry, University of Stuttgart, Germany

Hiromi Nakai

Department of Chemistry and Biochemistry, School of Advanced Science and Engineering, Waseda University, Japan

Krishnan Raghavachari

Department of Chemistry, Indiana University, USA

Michael A. Salim

Department of Chemistry, University of Illinois at Urbana–Champaign, USA

Max Schwilk

Institute for Theoretical Chemistry, University of Stuttgart, Germany

Lyudmila V. Slipchenko

Department of Chemistry, Purdue University, USA

Olaseni Sode

Department of Chemistry, University of Illinois at Urbana–Champaign, USA

Department of Chemistry, Biochemistry, and Physics, The University of Tampa, USA

Casper Steinmann

Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, Denmark

Hans-Joachim Werner

Institute for Theoretical Chemistry, University of Stuttgart, Germany

Theresa L. Windus

Ames Laboratory of United States Department of Energy, USA

Department of Chemistry, Iowa State University, USA

Sebastian Wouters

Center for Molecular Modelling, Ghent University, Belgium

Frick Chemistry Laboratory, Department of Chemistry, Princeton University, USA

Kiyoshi Yagi

Department of Chemistry, University of Illinois at Urbana–Champaign, USA

Theoretical Molecular Science Laboratory, RIKEN, Japan

Takeshi Yoshikawa

Department of Chemistry and Biochemistry, School of Advanced Science and Engineering, Waseda University, Japan

Kuang Yu

Department of Mechanical and Aerospace Engineering, Princeton University, USA

John Z. H. Zhang

School of Chemistry and Molecular Engineering, East China Normal University, China

NYU-ECNU Center for Computational Chemistry, NYU Shanghai, China

Department of Chemistry, New York University, USA

Tong Zhu

School of Chemistry and Molecular Engineering, East China Normal University, China

NYU-ECNU Center for Computational Chemistry, NYU Shanghai, China

Note

*Current address: School of Chemistry, University of Hyderabad, India

Preface

Electronic structure theory, also referred to as ab initio quantum chemistry (QC), has attained a high level of maturity and reliability for gas-phase molecules of modest size. Unfortunately, the formal scaling of these methods such as Hartree–Fock (HF), density functional theory (DFT), second-order perturbation theory (MP2), coupled cluster theory (CC), and multi-reference (MR) methods hinder their application to large molecules, to condensed phase systems or to excited electronic state potential energy surfaces. These limitations are especially severe for methods that account for electron correlation, such as MP2, CC, and MR methods, since their scaling with system size is steeper than for the simpler HF and DFT methods. There is therefore a need for computational strategies that nearly retain the accuracy of the most reliable methods while greatly reducing the scaling of these methods as a function of system size. While researchers who are interested in simulations of large molecular systems have often turned to classical molecular mechanics (MM) force fields, MM methods are limited in their applicability. While there are a few exceptions, classical MM cannot realistically treat bond making/bond breaking (the essence of chemistry) or excited state phenomena.

One effective QC approach that has become increasingly popular is referred to as fragmentation (broadly defined) or embedding theory. Fragmentation commonly refers to the physical subdivision of a large molecule into fragments, each of whose energy can be computed on a different compute node, thereby making the overall computation highly parallel. Fragmentation methods of this type scale nearly linearly with system size and can take advantage of massively parallel computers. Fragmentation methods of this type are discussed in Chapters 3, 5, 6, 7, 10, and 11. An alternative approach to physical fragmentation of a molecule is to fragment the wave function, by employing localized molecular orbitals to separate the wave function into domains that can be separately correlated. This approach is based on the fact that electron correlation is short-range. Chapter 1 provides an excellent discussion of local electron correlation methods by one of the leaders in the field.

Embedding methods are similar to fragmentation methods in that a total system is partitioned into multiple subsystems, in a manner that allows the incorporation of interactions among the subsystems. Like fragmentation and local orbital approaches, embedding methods reduce the steep scaling of traditional electronic structure methods. Embedding methods frequently involve multiple levels of theory. Approaches to embedding methods are discussed in Chapters 2, 4, 8, and 9.

The methods that are discussed in this book provide an exciting path forward to the accurate study of large molecules and condensed phase phenomena.

1

Explicitly Correlated Local Electron Correlation Methods

Hans-JoachimWerner, Christoph Köppl, Qianli Ma, and Max Schwilk

Institute for Theoretical Chemistry, University of Stuttgart, Germany

1.1 Introduction

Accurate wave function methods for treating the electron correlation problem are indispensable in quantum chemistry. A well-defined hierarchy of such methods exists, and in principle, these methods allow to approach the exact solution of the non-relativistic electronic Schrödinger equation to any desired accuracy. A much simpler alternative is density functional theory (DFT), which is probably most often used in computational chemistry. However, its failures and uncertainties are well known, and there is no way for systematically improving or checking the results other than comparing with experiment or with the results of accurate wave function methods.

Due to the steep scaling of the computational resources (CPU-time, memory, disk space) with the molecular size, conventional wave function methods such as CCSD(T) (coupled-cluster with single and double excitations and a perturbative treatment of triple excitations) can only be applied to rather small molecular systems. For example, the CPU-time of CCSD(T) scales as , where is a measure of the molecular size (e.g., the number of correlated electrons) and even the simplest electron correlation method, MP2 (second-order Møller-Plesset perturbation theory) scales as . This causes a scaling wall that cannot be overcome. Even with massive parallelization and using the largest supercomputers, this wall can only be slightly shifted to larger systems. However, it is well known that electron correlation in insulators is a short-range effect. The pair correlation energies decay at long-range with R− 6, where R is the distance between two localized spin orbitals. Therefore, the steep scaling is unphysical. It results mainly from the use of canonical molecular orbitals, which are usually delocalized over larger parts of the molecule.

The scaling problem can be much alleviated by exploiting the short-range character of electron correlation using local orbitals and by introducing local approximations. This was first proposed in the pioneering work of Pulay et al. [1–6], and in the last 20 years enormous progress has been made in developing accurate local correlation methods. There are two different approaches, both of which are based on the use of local orbitals. The traditional one is to treat the whole molecule in one calculation and to apply various approximations that are based on the fast decay of the correlation energy. We will denote such methods local correlation methods. A large variety of such approaches has been published in the past [7–59].

The second approach is the so-called fragmentation methods [60–87], in which the system is split into smaller pieces. These pieces are treated independently, mostly using conventional methods (although the use of local correlation methods is also possible). The total correlation energy of the system is then assembled from the results of the fragment calculations. Various methods differ in the way in which the fragments are chosen and the energy is assembled. A special way of assembling the energy using a many-body expansion is used in the so-called incremental methods [88–96], but these also belong to the group of fragmentation methods. Fragmentation methods will be described in other chapters of this volume and are therefore not the subject of this chapter. However, in Section 1.13, we will comment on the relation of local correlation and fragmentation methods.

Another problem of the CCSD(T) method is the slow convergence of the correlation energy with the basis set size. Very large basis sets are needed to obtain converged results, and this makes conventional high-accuracy electron correlation calculations extremely expensive. This problem is due to the fact that the wave function has a cusp for r12 → 0, where r12 is the distance between two electrons. The cusp is due to the singularity of the Coulomb operator , and cannot be represented by expanding the wave function in antisymmetrized products of molecular orbitals (Slater determinants). This leads to the very slow convergence of the correlation energy with the size of the basis set, and in particular with the highest angular momentum of the basis functions. This problem can be solved by including terms into the wave function that depend explicitly on the distance r12, and these methods are known as explicit correlation methods [97–155].

The combination of explicit correlation methods with local approximations has been particularly successful [140–153]. As will be explained and demonstrated later in this chapter, this does not only drastically reduce the basis set incompleteness errors, but also strongly reduces the errors caused by local approximations. Local correlation methods employ two basic approximations. The first is based on writing the total correlation energy as a sum of pair energies. Each pair describes the correlation of an electron pair (in a spin-orbital formulation), or, more generally, the correlation of the electrons in a pair of occupied local molecular orbitals (LMOs). Depending on the magnitude of the pair energies, it is possible to introduce a hierarchy of strong, close, weak, or distant pairs [7,18,31,32]. Different approximations can be introduced for each class, ranging from a full local coupled-cluster (LCCSD) treatment for strong pairs to a non-iterative perturbation correction for distant pairs, which can be evaluated very efficiently using multipole approximations [12, 13]. We will denote such approximations as pair approximations. The second type of local approximations is the domain approximation, which is applied to each individual pair. A domain is a subset of local virtual orbitals which is spatially close to the LMO pair under consideration. Asymptotically, the number of orbitals in each pair domain (the domain sizes) become independent of the molecular size. Also the number of pairs in each class (except for the distant pairs) becomes independent of the molecular size. This leads to linear scaling of the computational effort as a function of molecular size, as has already been demonstrated for LMP2 and up to the LCCSD(T) level of theory more than 25 years ago [12–18].

The critical question is, of course, how quickly the correlation energy as well as relative energies (e.g., reaction energies, activation energies, intermolecular interaction energies, and electronic excitation energies) converge with the domain sizes and how they depend on the pair approximations. The domain sizes which are necessary to reach a certain accuracy (e.g., 99.9% of the canonical correlation energy) depends sensitively on the choice of the virtual orbitals. As is known since the 1960s, fastest convergence is obtained with pair natural orbitals (PNOs) [156], and this has first been fully exploited in the seminal PNO-CI and PNO-CEPA methods of Meyer [157, 158], and somewhat later also by others [159–163]. The problem with this approach is that the PNOs are different for each pair and non-orthogonal between different pairs. This leads to complicated integral transformations and prevented the application of PNO methods to large molecules for a long time. The method was revived by Neese and coworkers in 2009 and taken up also by others (including us) later on [32,33,48–57,146–150]. The problem of evaluating the integrals was overcome by using local density-fitting approximations [22]. Furthermore, the integrals are first computed in a basis of projected atomic orbitals (PAOs), which are common to all pairs, and subsequently transformed to the pair-specific PNO domains [54,146,147,150]. Also, hybrid methods, in which so-called orbital-specific virtuals (OSVs) [164–167] are used at an intermediate stage, have been proposed [53, 146, 147]. Later sections of this chapter will explain these approaches in some detail.

Local approximations have also been developed for multi-reference wave functions [168–176]. The description of these methods is beyond the scope of the current article, but we mention that recently very efficient and accurate PNO-NEVPT2 [175] (N-electron valence state perturbation-theory) and PNO-CASPT2 [176] (complete active space second-order perturbation theory) methods have been described.

In the current article, we will focus on new developments of well-parallelized PNO-LMP2-F12 and PNO-LCCSD-F12 methods recently developed in our laboratory. These methods also have a close relation to the methods of Neese et al. [54–57, 153]. After introducing some benchmark systems, which will be used later on, we will first outline the principles of local correlation and describe the choice of the local occupied and virtual orbitals as well as of the domains. The convergence of the correlation energy as a function of the domain sizes will be demonstrated for various types of virtual orbitals for LMP2. Subsequently, based on these foundations, we will discuss more advanced approximations for distant pairs and close/weak pair approximations used in local coupled cluster methods. Next, we will present an introduction to local explicit correlation methods, and demonstrate the improvements achieved by the F12 approach both for LMP2-F12 and LCCSD-F12. Finally, we will describe some important technical details, such as local density fitting and parallelization. A summary concludes the chapter.

1.2 Benchmark Systems

Some large molecules and reactions, which we have used extensively to benchmark our methods [32, 145, 147, 148], are shown in Figures 1.1 and 1.2. For easy reference, we have given short names to some of the molecules, which are shown in the figure and will be used throughout this article. Reaction I is the last step in the synthesis of androstendione. In reaction II, testosterone is esterified to make it more lipophilic for a longer retention time in body tissues. Reaction III is the dissociation of a gold(I)-aminonitrene complex (AuC41H45N4P, for simplicity denoted Auamin, see Figure 1.1). This reaction is taken from Ref. [177] and plays an important role in catalytic aziridination and insertion reactions. The Auamin molecule has three phenyl and three mesityl groups and therefore strong long-range dispersion interactions are expected.

Diagrams show chemical reactions and molecules with benchmark with three reactions and plots for precursor, androstendione, testosterone, ester, auamin, auamin2s, mes equals mesityl, and P(Ph)3.

Figure 1.1 Benchmark molecules and reactions.

Image described by caption and surrounding text.

Figure 1.2 Visualization of a selection of large molecules mentioned in Tables 1.3 and 1.4.

In Figure 1.2, some representative medium-sized organic molecules with about 100–150 atoms are shown. The selection of molecules contains a modified crown-ether (C52O6H48) as well as the biologically active molecules nonactin and elaiophylin. We use these and some other molecules to demonstrate the convergence of our local correlation methods toward the canonical results and to perform benchmark calculations.

Most calculations have been carried out using the cc-pVDZ-F12 (VDZ-F12) and cc-pVTZ-F12 (VTZ-F12) basis sets [154]. For the gold atom, the ECP60MDF effective core potential [178] for the inner 60 electrons along with the aug-cc-pVDZ-PP and aug-cc-pVTZ-PP basis sets [179] were used. Density-fitting (DF) approximations were employed throughout this work for the evaluation of 2-electron integrals, and the aug-cc-pVTZ/MP2FIT basis set [180] was used as the DF auxiliary basis. We also used the aug-cc-pVTZ/JKFIT basis, which was derived from the cc-pVTZ/JKFIT basis set [181] by adding for each angular momentum another shell of diffuse functions, as the auxiliary basis for the resolution of the identity (RI) in F12 calculations. Unless otherwise stated, the parameters for local approximations are our program defaults, as described in the following sections.

For reaction III, the computed dissociation energy can be directly compared to an experimental gas-phase value [177] of 196.5 ± 11.2 kJ mol−1, which has been obtained by subtracting the PW91/cc-pVTZ-pp zero-point correction [177] of −8.2 kJ mol−1 from the measured value. The Hartree–Fock value is computed to be 92 kJ mol−1. Thus, the correlation contribution to the dissociation energy is very large, and, as will be demonstrated later, very sensitive to local approximations. The reaction therefore provides an excellent difficult benchmark system.

All methods described and used in this article are implemented in the development version of the MOLPRO quantum chemistry package [182,183], and will be made available for general use in the next release.

1.3 Orbital-Invariant MP2 Theory

The basic approximations in local correlation methods can be well understood by considering the closed-shell MP2 method. We will therefore discuss local MP2 (LMP2) methods first; more advanced coupled cluster methods will be presented in Section 1.9.

Since we intend to use localized orbitals, the theory must be formulated in a way that is invariant with respect to unitary transformations within the occupied valence and virtual orbital spaces. This means that the Fock matrix is non-diagonal in each of these subspaces. The choice of local occupied and virtual orbitals will be discussed in the subsequent sections.

We will introduce quite a number of different orbital spaces and, for the sake of clarity, the index notation is summarized in Table 1.1.

Table 1.1 Index notation for various orbital spaces.

The first-order MP2 wave function is written as

(1.1)

numbered Display Equation

where i, j are occupied and a, b are virtual orbitals. The summation limit val over the first summation indicates that the sum runs only over the correlated valence orbitals (excluding uncorrelated core orbitals). are doubly excited configurations, which are obtained by applying spin-summed excitation operators to the Hartree–Fock reference wave function |0⟩ ≡ |ΦHF⟩:

(1.2) numbered Display Equation

The MP2 amplitudes Tijab can be determined by minimizing the Hylleraas functional

(1.3)

numbered Display Equation

which leads to the first-order amplitude equations

(1.4)

numbered Display Equation

The quantities Rijab are residuals, which must vanish for the optimized amplitudes. Since the configurations |Φabij⟩ are pairwise non-orthogonal for ij and ji, it is convenient to use in equation (1.4) contravariant configurations

(1.5) numbered Display Equation

which have the properties

(1.6)

numbered Display Equation

(1.7) numbered Display Equation

Expanding the first-order wave function as

(1.8) numbered Display Equation

and equating this to equation (1.1) yields the corresponding contravariant amplitudes

(1.9) numbered Display Equation

It is then straightforward to derive that

(1.10) numbered Display Equation

where

(1.11)

numbered Display Equation

are the electron repulsion integrals (ERIs) in Mulliken notation and ρai(r) = φa(r)φi(r) are 1-electron charge distributions. In the following, we will often consider the integrals for a fixed pair ij as a matrix

(1.12) numbered Display Equation

Similarly, the amplitudes Tijab = [Tij]ab and residuals Rijab = [Rij]ab are matrices, where superscripts denote different matrices, and subscripts their indices. Note that the matrices are all defined in the virtual orbital space, that is, their indices always refer to virtual orbitals. The Hylleraas functional (equation 1.3) can then be written as

(1.13)

numbered Display Equation

For optimized amplitudes the residuals Rijab vanish, and the Hylleraas energy E2 becomes equal to the standard second-order energy expression

(1.14) numbered Display Equation

The orbital-invariant form of the MP2 residuals is

(1.15)

numbered Display Equation

where [F]ab and fkl are the external (virtual) and internal (occupied) blocks of the closed-shell Fock matrix. In general, the linear equation system Rijab = 0 (cf. equation 1.4) must be solved iteratively. In the special case where canonical HF orbitals are used (i.e., fij = δijεi, Fab = εaδab), the solution reduces to

(1.16) numbered Display Equation

If the occupied orbitals are local but the virtual orbitals are still canonical, one can write

(1.17)

numbered Display Equation

Approximate amplitudes can then be obtained by neglecting the coupling terms Gijab

(1.18) numbered Display Equation

This is denoted a semi-canonical amplitude approximation. The energies obtained with semi-canonical amplitudes are usually too inaccurate for practical applications, but these amplitudes will be useful for generating PNOs, cf. Section 1.6.3, or for distant pair approximations, cf. Section 1.8.

Similarly, the linear amplitude equations can be solved iteratively by computing in each iteration the residuals and from these, by applying first-order perturbation theory, the amplitude updates

(1.19) numbered Display Equation

The semi-canonical amplitudes are the first approximation in this iterative algorithm.

1.4 Principles of Local Correlation

Let us now assume that |i⟩ is a local orbital with charge center at position Ri. We also assume that |i⟩ can be expanded in a local subset of atomic orbitals (AOs) which are spatially close to Ri. Then the charge distributions ρai(r) will be local around Ri. This is true even for canonical (delocalized) virtual orbitals |a⟩, since only those parts of them will contribute to the charge distribution ρai which overlap with |i⟩. If the virtual orbitals are also local, the locality of the ρai is further improved. If the occupied and virtual orbitals are both expanded in AOs that are close to Ri and Ra, respectively, it follows from the exponential decay of the AOs that the integral (ai|bj) is a sum of contributions that decay exponentially with the distances between the basis functions of |i⟩ and |a⟩.

Similar considerations hold for b and j, and so the integral (ai|bj) will decay exponentially with the distances Rai = |Ra − Ri| and Rbj = |Rb − Rj| (if these are sufficiently large). This forms the basis for domain approximations, which means that to each LMO |i⟩, a subset of local virtual orbitals a can be assigned. We will denote such orbital domains as a ∈ [i]. Similarly, to each pair ij of LMOs, a pair domain [ij] ≡ [i]∪[j] can be assigned. In the LMP2 wave function, the excitations |Φijab⟩ are restricted to the domain a, b ∈ [ij]. The sizes of these domains that are necessary to reach a certain accuracy depend crucially on the choice and construction of the virtual orbitals. This will be discussed in more detail in Section 1.6.

If LMO |i⟩ is far from |j⟩, the local charge distribution ρai will also be far from ρbj. If the virtual orbitals |a⟩, |b⟩ are local with |a⟩ being close to |i⟩ and |b⟩ close to |j⟩, the exchange-type integrals (bi|aj) will decay exponentially and can be neglected at large distances of |i⟩ and |j⟩. Using the semi-canonical amplitude approximation (cf. equation 1.18) the distant pair energies Edistij can then be approximated as

(1.20)

numbered Display Equation

The integrals (ai|bj) can be efficiently computed using local density-fitting techniques (cf. Section 1.12.1) or approximated by a multipole expansion [12, 13, 54] (cf. Section 1.8). Since the densities ρai(r) carry no charge (∫ρai(r)dr = 0), the lowest-order contribution is the dipole–dipole interaction, which decays with R− 3, where Rij is the distance between the charge centers of the LMOs i and j. Consequently, in a local orbital basis, the amplitudes and pair energies decay with R− 3ij and R− 6ij, respectively. This is the basis of pair approximations, which simplify the treatment of distant pairs ij, cf. Section 1.8.

1.5 Orbital Localization

Localized occupied molecular orbitals (LMOs) |i⟩ ≡ |φloci⟩ can be obtained from the canonical molecular orbitals (CMOs) |φcank⟩ by a unitary transformation

(1.21) numbered Display Equation

Uncorrelated core orbitals are excluded from the localization, since mixing of core and valence orbitals would affect the correlation energy even without any local approximations. The localized valence orbitals |i⟩ can be expressed in the atomic orbital (AO) basis |μ⟩ ≡ |χμ⟩ by a rectangular transformation matrix L (NAO × Nval), which is related to the corresponding canonical MO-coefficient matrix C by the unitary transformation Uloc:

(1.22)

numbered Display Equation

In the following, the indices i, j, k, l will always refer to LMOs.

There are various criteria that can be used to determine the transformation matrix Uloc. For a recent review, see Ref. [184]. Perhaps the most well-known ones are the methods of Foster and Boys (FB), Pipek-Mezey (PM) [185], and Edmiston and Ruedenberg (ER) [186, 187]. ER localization minimizes the overall inter-orbital electron repulsion. This requires the 2-electron integrals, and is therefore quite expensive and rarely used. In FB localization, the sum of the orbital variances

(1.23) numbered Display Equation

is minimized. PM localization maximizes the sum of the squared Mulliken partial charges qiA:

(1.24) numbered Display Equation

(1.25) numbered Display Equation

where the first sum in the latter equation runs over all basis functions μ centered at atom A. Effectively, this means that the number of centers at which the orbitals are localized is minimized. Pipek-Mezey (PM) localization often suffers from the problem that poor localization may result in calculations with large basis sets containing diffuse functions.

The latter problem can be avoided by using intrinsic bond orbitals (IBOs) [188]. In this case, a molecule-intrinsic minimal basis of polarized atomic orbitals (IAOs) |ρ⟩ [188–191] is created, in which the occupied Hartree–Fock orbitals can be expanded exactly. The IAOs are generated from a suitable minimal basis of free-atom AOs using projection operators [188] and are subsequently symmetrically orthogonalized. The free atom AOs are stored in a basis set library (we use a subset of the generally contracted cc-pVTZ basis where available). The partial charges are then defined as

(1.26) numbered Display Equation

and used in the PM localization. Alternatively, Knizia [188] has proposed to replace the exponent of two in equation (1.24) by four, which helps to converge the iterative localization procedure in highly symmetric cases where several equivalent solutions exist. We will denote the IBO localizations with exponents 2 and 4 as IBO2 and IBO4, respectively. The resulting IBOs closely resemble PM orbitals, but are very insensitive to basis set variations; in particular, the commonly encountered artifacts of PM orbitals with diffuse basis sets are absent. Additionally, the IBO construction yields stable partial orbital charges which can be used for defining domains (see Section 1.6.2). Closely related localization schemes were also put forward by Lehtola and Jónsson [192] and West et al. [191]. The relation between these schemes has been investigated by Janowski [189].

The construction of IBOs is simple and efficient, and since the iterative localization is carried out in the minimal IAO basis, it is much faster than the standard PM localization. Even though it scales cubically with molecular size, the time to generate IBOs is found to be short as compared to a PNO-LMP2 calculation, even for cases with about 10,000 basis functions.

A fourth localization scheme, which has been used in the past in local correlation methods [28], is based on the natural bond orbital (NBO) analysis of Weinhold and co-workers [193–196] and yields so-called natural localized molecular orbitals (NLMOs) [194]. Somewhat similar to IBOs, the NLMOs are obtained by a sequence of transformations which first yields natural atomic orbitals (NAOs), then NBOs, and finally the NLMOs. However, these transformations are much more complicated and less uniquely defined than in the IBO scheme. Based on these transformations partial charges can also be defined [28], which are quite stable with respect to basis set variations and are rather similar to IBO partial charges.

Image described by caption and surrounding text.

Figure 1.3 Localized C–C bonding orbitals of naphtalene. FB localization yields pairs of equivalent banana bonds, while the other localization schemes keep σ and π orbitals separate. Only the symmetry unique orbitals are shown. The specifications in curly brackets denote the symmetry operations (mirror planes) which generate equivalent orbitals. The IBO4 orbitals are not fully symmetry-adapted or equivalent. Basis set: VTZ.

FB localization yields qualitatively different localized orbitals for double bonds than the other schemes. Usually, PM, IBO, and NLMO localization keeps σ- and π-orbitals separate (even in non-planar molecules). In contrary, Boys localization often generates two equivalent banana-bond-like orbitals. Some examples are shown in Figure 1.3 for naphtalene. The PM, IBO2, and NBO orbitals are qualitatively similar but not identical. The differences are rather small and therefore only the IBO2 orbitals are shown in Figure 1.3.

Table 1.2 shows a comparison of domain sizes and correlation energies using various localization schemes for a number of aromatic molecules. Despite the qualitative differences of the various localized orbitals, the localization method has only a minor effect on the domain sizes and correlation energies. Exceptions are special cases in which the LMOs are neither symmetry-adapted nor symmetry-equivalent. It may then happen that the domains are neither symmetry-adapted nor symmetry-equivalent, and this can cause spurious dipole moments perpendicular to a molecular plane and other artifacts. An example is the IBO4 localization of naphtalene, cf. Figure 1.3. In longer polycyclic aromatic hydrocarbons such as tetracene or pentacene, none of the localization methods creates symmetry-adapted or equivalent LMOs. In such cases, the non-equivalent orbital domains should be merged so that, overall, the symmetry properties are recovered.

Table 1.2 Domain sizes and PNO-LMP2 correlation energies (in Eh) using various localization schemes.

Basis set: VTZ, TLMO = 0.2, ToccPNO = 10− 8. The parameter IEXT and ToccPNO determine the PAO and PNO domain sizes, respectively (see Section 1.7).

1.6 Local Virtual Orbitals

In this section, we will discuss various types of virtual orbitals and the corresponding domain constructions. The choice of virtual orbitals is most essential for the accuracy and efficiency of local correlation methods. There are schemes which localize virtual orbitals in a similar way as the occupied orbitals, yielding orthonormal localized virtual orbitals (LVOs) [184]. Even though these orbitals have good locality properties and have been successfully used in fragmentation approaches, we have found that, in local correlation methods, much smaller domains are sufficient for a certain accuracy if projected atomic orbitals (PAOs), orbital-specific virtuals (OSVs), or pair natural orbitals (PNOs) are used. PAOs are obtained by projecting atomic orbitals or basis functions to the virtual space, yielding a single set of non-orthogonal functions. In contrast, in OSV and PNO methods, different orthonormal orbital sets are used for each orbital or orbital pair, respectively. These sets are mutually non-orthogonal. In contrast to LVOs, PAOs, OSVs, and PNOs are not obtained by applying a localization criterion, but they are inherently local by construction. The domain sizes that are necessary to obtain a certain fraction of the canonical correlation energy decrease in the order LVO > PAO > OSV > PNO. As will be shown later, the PNO domain sizes are typically one order of magnitude smaller than the PAO ones, which reduces the number of amplitudes to be optimized by two orders of magnitude. LVOs will not be further considered in this article. Before defining PAOs, OSVs, and PNOs in detail, some general considerations and definitions will be presented.

1.6.1 Pseudo-Canonical Pair-Specific Orbitals

In the most general case, one can use completely independent sets of virtual orbitals for each pair ij. We will denote these as pair-specific virtuals (PSVs). PNOs are a special choice of PSVs, but for the time being, we consider the problem in general. The tilde over the index indicates that the orbitals may be non-orthogonal (but they must always be orthogonal to the occupied orbital space). The PSVs are related to the orthonormal canonical virtuals (CMOs) {b} by transformation matrices :

(1.27) numbered Display Equation

These transformations are of course entirely useless unless they produce more local virtual orbitals so that one can restrict the PSV excitation subspaces to domains [ij]]PSV. This means that we approximate the pair functions as

(1.28)

numbered Display Equation

The transformation matrices then become rectangular, with more rows than columns. To keep the formulation compact and to avoid superscripts on subscripts, we will always use the convention that the subscripts of matrices such as , Tij, or Kij correspond to the pair ij given by the superscript of the matrix (here pair ij). The CMOs can be considered as a special case in which the domain [ij]CMO is identical and complete for all pairs.

The PSVs may in general be non-orthogonal with overlap

(1.29)

numbered Display Equation

Note that here the indices and belong to different pairs, as indicated by the superscripts of the overlap matrix. Inserting equation (1.27) into equation (1.28) yields for the transformation of the amplitudes

(1.30) numbered Display Equation

Integral and Fock matrices transform in the opposite way, for example,

(1.31) numbered Display Equation

In matrix form, these transformations can be written as

(1.32) numbered Display Equation

(1.33) numbered Display Equation

Since the number of PSVs per pair is in general smaller than the number of CMOs, the inverse transformations are not possible. Equation (1.32) shows that one can always transform the amplitudes from the PSV domain back to the CMO basis, but not vice versa. In contrary, integrals can only be transformed from the larger basis to the smaller one, cf. equation (1.33).

Within each domain [ij]PSV, the energy is still invariant with respect to transformations among the orbitals . It is useful to apply another transformation within each domain which makes the orbitals pseudo-canonical (PC). This means that in the subspace of the domain [ij]PSV, the PSVs become orthonormal and diagonalize the Fock matrix. The pseudo-canonical pair-specific orbitals (PC-PSVs) are denoted |rij⟩ and are related to the original PSVs by the transformation

(1.34) numbered Display Equation

so that

(1.35)

numbered Display Equation

In practice, this is achieved by first orthogonalizing the orbitals (by singular value decomposition if linear dependencies exist), and then diagonalizing the Fock matrix in the basis of the orthogonal PSVs. Note that the PC-PSVs of different pairs are still non-orthogonal. The total transformation from CMOs to PC-PSVs is then

(1.36) numbered Display Equation

(1.37) numbered Display Equation

If the domain spans the whole virtual space, the PC-PSVs become equal to the CMOs.

It is straightforward to transform the amplitude equations to the PC-PSV basis by forming RijPC-PSV = Qij†RijCMOQij and inserting TijCMO = QijTijPC-PSVQij†. Finally, the products [Qij†Qkl]rs are replaced by [Sij, kl]rs (cf. equation 1.29). This always occurs if an amplitude index is not connected to an integral index. For example, the LMP2 amplitude equations in equation (1.17) become

(1.38)

numbered Display Equation

The optimized amplitudes are now defined by the conditions Rijrs = 0, r, s ∈ [ij]PC-PSV. Due to the domain approximation, the number of amplitudes and corresponding equations is much smaller than in the canonical case. The Hylleraas functional in the PC-PSV basis is computed as

(1.39)

numbered Display Equation

In the following subsections, we will address the question how to define the PSVs and which are the most suitable ones.

1.6.2 Projected Atomic Orbitals

In a seminal series of papers on local correlation methods, Pulay has proposed to use projected atomic orbitals (PAOs) to span the virtual space [1], and this has been widely used later in local correlation methods. The PAOs are labeled by indices and are defined as:

(1.40) numbered Display Equation

where are suitable atomic orbitals and projects out the contributions of the occupied space,

(1.41) numbered Display Equation

so that for all . In general, the atomic orbitals can be represented in the AO basis by a block-diagonal matrix CAO

(1.42) numbered Display Equation

Each non-zero block of the coefficient matrix CAO corresponds to one center. In most cases, the coefficient matrix is assumed to be diagonal, that is, the AO basis functions themselves are projected.

In the AO basis, the PAOs are represented by a coefficient matrix P,

(1.43) numbered Display Equation

with

(1.44) numbered Display Equation

The rectangular ( ) coefficient matrix Cv represents the canonical virtual orbitals, and [S]μν = ⟨μ|ν⟩ is the overlap matrix of the basis functions. Thus, the PAOs can be obtained from the canonical virtual orbitals by a transformation

(1.45) numbered Display Equation

(1.46) numbered Display Equation

Note that the matrix is rectangular (Nvirt × NAO), and the NAO PAOs are non-orthogonal and linearly dependent. Their overlap matrix

(1.47) numbered Display Equation

has Nocc zero eigenvalues. The linear dependencies in the PAO set are removed by singular value decomposition for individual domains.

The PAOs are neither orbital- nor pair-specific, and each one can be related to an AO (or a basis function) at one atom. In most cases, PAO domains are therefore center-based, that is, all PAOs that originate from basis functions at a selection of centers comprise the orbital domains [i]PAO. The domains are then uniquely defined by a list of centers, independent of the basis set. For each LMO i, the centers can, for example, be selected based on their partial IBO charges QiA (cf. equation 1.26). PAO pair domains are then defined as the union of the two related orbital domains, [ij]PAO = [i]PAO∪[j]PAO. These are orthonormalized and made pseudo-canonical, as explained in Section 1.6.1. More details about the choice of PAO domains will be given in Section 1.7.

An advantage of the non-orthogonal PAOs is that they are pair-independent, and the overlap or Fock matrices for a domain [ij]PAO can always be extracted from the full NAO × NAO matrices SPAO or FPAO, respectively. Also the integral transformation is simpler than for truly pair-specific orbitals, in particular for local CCSD methods. In PAO-LMP2 or PAO-LCCSD methods, it is therefore advantageous to compute the residuals first in the non-orthogonal PAO basis, and transform them to the pair-specific PC-PAO basis just for the update of the amplitudes. Subsequently, the amplitude update is transformed back to the non-orthogonal PAO basis:

(1.48) numbered Display Equation

(1.49) numbered Display Equation

These transformations are similar to those in equations (1.32) and (1.33), but with Qij replaced by WijPAO.

Unfortunately, rather large PAO domains (400-600 PAOs per pair for triple-ζ basis sets) are necessary to recover about 99.9% of the canonical correlation energy in large molecules. Thus, the memory and CPU requirements of such calculations can be very large, despite the formally linear scaling (if distant pairs are