# Binary Code Similarity Detection

Zian Liu  Department of Computer Science and Software Engineering
Swinburne University of Technology
Data61, CSIRO
Melbourne, Australia
<EMAIL_ADDRESS>
###### Abstract

Binary code similarity detection is to detect the similarity of code at binary
(assembly) level without source code. Existing works have their limitations
when dealing with mutated binary code generated by different compiling
options. In this paper, we propose a novel approach to addressing this
problem. By inspecting the binary code, we found that generally, within a
function, some instructions aim to calculate (prepare) values for other
instructions. The latter instructions are defined by us as _key_ instructions.
Currently, we define four categories of key instructions: calling
subfunctions, comparing instruction, returning instruction, and memory-store
instruction. Thus if we symbolically execute similar binary codes, symbolic
values at these key instructions are expected to be similar. As such, we
implement a prototype tool, which has three steps. _First_ , it symbolically
executes binary code; _Second_ , it extracts symbolic values at defined key
instructions into a graph; _Last_ , it compares the symbolic graph similarity.
In our implementation, we also address some problems, including path explosion
and loop handling.

###### Index Terms:

Binary code, code analysis, symbolic execution

##  I Introduction

My PhD goal aims at analyzing binary code to improve binary code security.
This includes similarity detection, code diffing, bug finding, etc., at binary
level. I plan to achieve this goal by implementing novel approaches in a
binary code similarity detecting tool. The current process is that the tool’s
framework has finished. I am fine-tuning this tool to conducting on a large
scale dataset.

With the rapid development of the software industry, binary code similarity
detection is playing an increasingly critical role. Software producers release
their products as binary code mainly to protect their source code. Binary code
similarity detection can be used in many fields such as bug search, malware
detection, malware clustering, malware linage, patch generation and analysis,
porting information, software theft detection[1], etc. However, binary code
similarity detection is challenging, which is mainly due to, 1) The same
source code can be compiled on different architectures, resulting in different
binary codes using different instruction sets; 2) During compiling process,
different compilers and compiling options may produce significantly different
binary code for the same source code; 3) Producers can use many code
protecting tools to obfuscate the code. There are some existing solutions to
this problem. These works can be divided into three categories with different
standards. In Section II, we mainly introduce these works according to their
comparison type: syntactic, semantic, and structural similarity. A typical
syntactical approach represents the code by statisticizing the occurrence of
specific strings in the code. Alternatively, one can use machine learning
based method to learn the vectorized representation of one or more sentences
in the code. A semantic approach evaluates whether two codes have similar
functionality or impact. Structural approach compares graph features such as
control flow graph and call graph. However, these solutions still have some
limitations. For example, syntactic similarity focus on the ‘appearance’ of
binary code thus sometimes failing to explore semantic similarity. Semantic
similarity, however, is not scalable and may not cover all of the code or
situations. Moreover, Ren et al. pointed out that existing binary code
similarity detection overlooked the impact of compiling options [2]. They
claimed that the testing set used in current work lacked enough compiling
options mutation. They tested several state-of-the-art binary similarity
detection tools with their crafted dataset. The detecting result all turned
out to be less accurate. Structural similarity is vulnerable for cross-
architecture comparison since they assume that basic blocks remain their
feature and relationship with other basic blocks under all circumstances.
However, by utilizing Ren el al’s method, one can tune the compiling
configuration to have less structural similarity. Therefore, our research
question is how can we detect binary code similarity even when these compiling
options are deployed.

By analyzing the binary code, we observed that some instructions calculate
value for some other instructions. The latter instructions are more important
for deciding the similarity while the former instructions should be given less
consideration. This is because the symbolic value of the former can propagate
to latter instructions. Only inspecting the symbolic value of these latter
instructions can reserve the main symbolic meaning of the function. We thus
define these important instructions as Key Instructions. In our tool, they are
translated to Intermediate Representation (IR) and connected based on control
flow to form a Key IR graph. We implemented our idea in a prototype tool. This
tool is able to lift different architectures’ binary code into Key IR graphs
and compare their similarity. This tool consists of three phases: 1) symbolic
execution, 2) Key IR graph construction, and 3) Key IR graph similarity
detection. In summary, our contributions are: 1) We proposed a Key IR graph to
abstract the semantic of binary code. 2) We implemented a prototype tool,
including the symbolic execution engine, IR lifting engine, and comparison
engine. This tool currently supports two architectures, i.e., x86-64 and ARM.
3) We proposed a novel method to compare the similarity of Key IR graphs.
Currently we have finished the implementation process and on our way to
evaluation.

##  II Background

As [1] mentioned in their survey, existing works on binary code similarity can
be classified according to different standards. In terms of the comparison
granularity, these works can be divided into instruction level, basic block
level, function level, and whole program level. In the binary similarity
comparison, one input is used as a query to compare with the target. Thus
according to input-target number, they can be classified into one-to-one, one-
to-many, and many-to-many comparisons. According to supported architectures,
they can be divided into a single architecture and cross-architecture. Single
architecture only supports one kind of assembly instruction set, while cross-
architecture supports more than one. According to analysis type, they can be
categorized into static analysis, static analysis, and hybrid analysis. Static
analysis analyzes the code without running the code, while dynamic analysis
executes the code with some input. The advantage of static analysis is
scalability and code coverage. Dynamic analysis is resource consuming but
helpful to gain semantic information. According to comparison type, they can
be classified into syntax, structural, and semantics similarity comparison.
Syntax similarity captures instruction representation similarity, while
structural similarity focuses on similarity in terms of graph representation
of the binary code. Semantic similarity highlights the similarity of code’s
impact.

We introduce various techniques used according to different comparison types.
For syntax similarity, common strategies include hashing, embedding, and
alignment. [3, 4, 5] all use hashing technique to output various instructions
sequences into fixed length of output hashing value. Equivalent output value
implies syntax similarity. [6, 7, 8, 9] generate an embedding from sequences.
[10, 11, 12, 13] automatically learn the embeddings for each instruction and
use them to produce the basic-level or function-level embedding. [14, 15, 16]
align two sequences and decide their similarity. For structural similarity,
common methods are optimization solution, k-subgraph matching, path
similarity, and graph embedding. [3, 17, 18, 19] transform the problem into
finding the mapping between two CFGs with minimum cost (optimization
solution). [20, 18, 21, 22] divide the graph into k subgraphs, so that each
one has k connected nodes. Matched number of those subgraphs indicates the
extent of the similarity. [23, 24, 16] determine similarity based on paths.
[19, 25, 26] extract features from graphs into feature vectors and determine
the vector similarity. For semantic similarity, general techniques are
instruction classification, input-output pairs, symbolic execution, theorem
prover, and semantic hashes. [20, 22, 24] classified instruction based on
their semantic in terms of arithmetic, logic or data transfer purpose. [27,
28, 29, 5, 30, 31] check whether output are the same to the input. [32, 27, 8,
23, 28, 33] use symbolic formulas to represent the binary code. With the
symbolic formula, [32] uses Theorem Prover to check whether two different
symbolic formula has similar output. [27, 34, 35] use symbolic hashes as an
alternative to the theorem prover. [36, 37] determine the edit distance of the
tree/graph of the symbolic formula.

##  III Methodology

The comparison granularity of our tool is function level. Given two functions
from two different input binary codes, our tool consists of three phases: 1.
symbolic execution, 2. Key IR graph construction, and 3. Key IR graph
similarity comparison. We firstly symbolically execute the binary code several
times to get their possible symbolic values. Then for the most important
instructions for comparison, we lift them and their symbolic values into Key
IRs. Those Key IRs are then combined to produce the Key IR graph. Lastly, we
compare the similarity of two IR graphs.

###  III-A Symbolic execution

Figure 1: 1st run

Figure 2: 2nd run

In this section, we randomly select paths and symbolically execute the
function to infer the results of each operand in each instruction. As shown in
Figure. 1 and Figure. 2, a function may contain many paths and symbolic
execution might have the path explosion problem. Thus we symbolically execute
the function many times. For each run, we randomly select a main path by using
deep-first algorithm to symbolic execute this path. However, if we only
execute one main path each run, we still may miss many instructions after many
times of run (e.g., we may miss instructions 3, 4, 5 in Figure.1, and
instructions 2, 4, 5 in Figure.2). To cover all the instructions in each run,
except the main path, we also select other paths to cover all the
instructions. We do not propagate to the instructions we already encountered
to avoid the path explosion problem in each run.

Since we randomly select paths in each run, some instructions demonstrate
different symbolic values. In fact, not each instruction can have different
states. Only the instructions at the joint of some paths can hold different
symbolic values (e.g., instruction 6 in Figure.1 and 2). Other instructions
can only hold one possible symbolic value (e.g., all other instructions in
Figure.1 and 2). The percentage of these possible symbolic values revealed is
dependent on the times of run. Generally, the more time the run, the more
chances we can reveal all the possible values.

###  III-B Key IR graph construction

(a) calling subfunction

(b) comparing instruction

(c) returning instruction

(d) memory address writing

Figure 3: Four types of Key Instructions

As we observed in the binary code, many instructions are making preparations
for some other instructions, such as calling subfunction, comparing
instructions, as shown in Figure. 3. In Figure.3(a), parameters are firstly
loaded into registers rcx and edx before calling the EVP_CIPHER_CTX_ctrl
subfunction. In Figure.3(b), values in register rbp are loaded from register
r13 and subtracted with value in register r14. Value in register rax is loaded
from some memory address before comparing with rbp. In Figure.3(c), value in
register eax is set to 0FFFFFFFF before stack balancing instructions (i.e.,
pop instructions) and return instruction. In Figure.3(d), value in somewhere
of the memory address is loaded to register rax and added to the value in
register r14. Then this value is written to some other memory location
[rsp+1D8h+var_170]. As we can see in all these examples, some instructions are
preparing values for some later instructions. We define these preparing
instructions as non-Key Instructions while the other instructions as Key
Instructions. From our observation, during execution, the non-Key Instructions
should propagate their values into the Key Instructions. And the Key
Instructions better describe the behavior of the binary code. Similar binary
code should contain similar Key Instructions with similar values (e.g.,
calling the same subfunction with similar parameters or comparing similar
values to similar values). We define the Key Instructions to contain two
types: control-flow impacting and control-flow irrelevant. Specifically, the
control-flow impacting include three categories: 1. calling subfunction, 2.
comparing instruction, and 3. returning instruction. We select them because
Calling subfunctions leads the execution flow into other parts of the binary
code. The result of comparing instructions decides which branches to take
next. Returning instructions leads execution flow back to the caller function.
The control-flow irrelevant instructions refer to the memory writing
instructions since it does not affect the control flow.

With the result of the first module, symbolic execution, we can translate Key
Instructions into Key IRs with symbolic values. Each instruction in the binary
code corresponds to one Key IR node. Then we connect the Key IRs based on
their control flow to form the Key IR graph. It is important to note that some
nodes only contain one possible symbolic value among those Key IR nodes, while
the other nodes might contain multiple symbolic values because they are the
joint of multiple paths.

###  III-C Key IR graph comparison

Figure 4: Key IR graph comparison

Because of the impact of cross-architecture, different compiling options, and
other factors, the Key IR graph might contain reordered nodes, inserted new
nodes, and duplicated nodes. Sometimes Key IR graph might also lose some nodes
and have other mutations because of the compilation options described in [2].
We propose a fuzzy method to match two Key IR graphs. To compare the
similarity of two Key IR graphs, we check how many similar nodes in two
graphs. To compare the similarity of a pair of node, we divide this process
into two phases: 1. single node textual similarity, and 2. context similarity.
In the first phase, we pick the most similar node pair from two Key IR graphs
as potentially similar pairs if these two nodes’ symbolic value has high
textual similarity. The symbolic values have been simplified beforehand by
msynth 111https://github.com/mrphrazer/msynth. Next, we compare the similarity
of this node pairs’ context. The context here refers to the neighbor nodes of
a node within a given boundary. Again, we compare these neighbors’ simplified
symbolic value textual similarity. An example of this graph comparison is
shown in Figure.4. Suppose we find the potentially similar node pair (blue
circle) with high similarity, and the context boundary is 1. Thus we examine
the neighbor nodes within this bound (green circle) and check their
similarity. When detecting the context similarity, the existence of similar
nodes within the context indicates the similarity. The more similar nodes are,
the more similar these binaries will be. We ignore the relations between these
nodes because the complexity of compilation options can significantly mutate
the relations.

##  IV Implementation

Our implementation is an IDA Pro plugin written in C++ and python. Our tool is
based on IDA pro’s control flow graph and its API. We have finished our
implementation and on our way for evaluation.

###  IV-A Symbolic execution

####  IV-A1 Symbolic values

At the beginning of each function, we assign VAR $0$ to VAR $n$ to parameters
$1$ to $N$ of this function. With the execution of each instruction, those
initial tags will be represented in various expressions and are propagated to
different instructions as their operands’ values. The similarity of the
expressions reveals similar semantics. Thus we transform binary code
similarity detection into expression similarity detection.

####  IV-A2 Handling loops

During the execution, we might encounter loops. Since it is often difficult to
decide the time of execution of a loop, it is challenging to propagate the
values of instructions in the loop to some Key Instructions. Also, it is non-
trivial to analyze which instructions’ operands are invariant during loops and
updated in each loop. To address these problems, we symbolically execute the
instructions in a loop twice. The invariant operands keep their constant
symbolic value, while the updated operands add an ‘ITER()’ notation outside
their value to highlight that their values are changed in the loop.

###  IV-B Key IR graph construction

We recover both control-flow impacting and control-flow irrelevant
instructions into Key IRs. We identify those instructions if they match some
patterns. The identification is implemented as a rule in the tool. For
control-flow impacting instructions, we recover three types, including 1.
calling subfunction instructions, 2. comparing instruction, and 3. returning
instructions. For Type 1, we find the corresponding parameters for the
subfunction. The patterns for matching Type 2 instruction are more complicated
than other types since some compiler options can translate comparison
instructions from source code to equivalent codes in an implicit way as in
[2]. For Type 3, we find the last modification of register rax in x86-64 or
register R0 in ARM before the function returns. For control-flow irrelevant
instruction, i.e., memory writing instruction, we identify them if the
mnemonic is ‘mov’ in x86-64 or ‘STR’ in ARM, and the destination operand is a
memory address.

Each Key Instruction is then transformed to a node in the Key IR graph. It is
important to note that each Key Instruction might have several possible
symbolic values. Thus we reserve all of them for each Key Instruction. After
we recovered all the nodes, we connect them based on the control flow provided
by IDA pro. The result Key IR graph has the same control flow as IDA pro’s
control flow. The difference is that all the non-Key Instructions are removed.

###  IV-C Key IR graph comparison

To compare the similarity of two symbolic formulas, we firstly simplify them
using msynth [msynth]. Msynth is originally used as a binary code
deobfuscation tool. It aims to simplify very complex binary code formulas into
simple formulas. We use it here to facilitate our comparison. Moreover,
various compilation options can be mitigated since the resulting formula will
be more similar, even the original formula may be different.

##  V Evaluation

We have collected many open-source benchmarks used widely in existing works,
i.e., OpenSSL, Coreutils, SPEC CPU2006, and SPEC CPU 2017. We will answer two
research questions in our experiment: 1. Can we effectively detect cross-
architecture binary codes? 2. Can we mitigate the influence of compiling
options? To answer the first question, we aim to compile our benchmark dataset
on ARM and x64 architectures with different compilers such as GCC and Clang.
Then we will prepare random pairs of similar and dissimilar function pairs
with labels indicating their similarity. To answer the second question, we
will compile benchmark dataset by using the compilation options described in
[2] and prepare similar and dissimilar function pairs. In both experiments,
the percentage of correctly detected similar and dissimilar function pairs
will be regarded as the accuracy, which is to indicate the performance of our
work.

## References

  * [1] I. U. Haq and J. Caballero, “A survey of binary code similarity,” ACM Comput. Surv., vol. 54, Apr. 2021.
  * [2] X. Ren, M. Ho, J. Ming, Y. Lei, and L. Li, “Unleashing the hidden power of compiler optimization on binary code difference: An empirical study,” in Proceedings of the 42nd ACM SIGPLAN International Conference on Programming Language Design and Implementation, PLDI 2021, (New York, NY, USA), p. 142–157, Association for Computing Machinery, 2021.
  * [3] X. Hu, T.-c. Chiueh, and K. G. Shin, “Large-scale malware indexing using function-call graphs,” in Proceedings of the 16th ACM Conference on Computer and Communications Security, CCS ’09, (New York, NY, USA), p. 611–620, Association for Computing Machinery, 2009.
  * [4] M. R. Farhadi, B. C. Fung, P. Charland, and M. Debbabi, “Binclone: Detecting code clones in malware,” in 2014 Eighth International Conference on Software Security and Reliability (SERE), pp. 78–87, 2014.
  * [5] Z. Xu, B. Chen, M. Chandramohan, Y. Liu, and F. Song, “Spain: Security patch analysis for binaries towards understanding the pain and pills,” in 2017 IEEE/ACM 39th International Conference on Software Engineering (ICSE), pp. 462–472, 2017.
  * [6] I. Santos, F. Brezo, J. Nieves, Y. Penya, B. Sanz, C. Laorden, and P. G. Bringas, “Idea: Opcode-sequence-based malware detection,” in ESSoS, 2010\.
  * [7] B. Kang, T. Kim, H. Kwon, Y. Choi, and E. G. Im, “Malware classification method via binary content comparison,” in Proceedings of the 2012 ACM Research in Applied Computation Symposium, RACS ’12, (New York, NY, USA), p. 316–321, Association for Computing Machinery, 2012.
  * [8] B. H. Ng and A. Prakash, “Expose: Discovering potential binary code re-use,” in 2013 IEEE 37th Annual Computer Software and Applications Conference, pp. 492–501, 2013.
  * [9] S. H. H. Ding, B. Fung, and P. Charland, “Kam1n0: Mapreduce-based assembly clone search for reverse engineering,” Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2016\.
  * [10] B. Liu, W. Huo, C. Zhang, W. Li, F. Li, A. Piao, and W. Zou, “diff: Cross-version binary code similarity detection with dnn,” in Proceedings of the 33rd ACM/IEEE International Conference on Automated Software Engineering, ASE 2018, (New York, NY, USA), p. 667–678, Association for Computing Machinery, 2018.
  * [11] F. Zuo, X. Li, Z. Zhang, P. Young, L. Luo, and Q. Zeng, “Neural machine translation inspired binary code similarity comparison beyond function pairs,” ArXiv, vol. abs/1808.04706, 2019.
  * [12] S. H. H. Ding, B. Fung, and P. Charland, “Asm2vec: Boosting static representation robustness for binary clone search against code obfuscation and compiler optimization,” 2019 IEEE Symposium on Security and Privacy (SP), pp. 472–489, 2019.
  * [13] L. Massarelli, G. A. D. Luna, F. Petroni, L. Querzoni, and R. Baldoni, “Safe: Self-attentive function embeddings for binary similarity,” ArXiv, vol. abs/1811.05296, 2019.
  * [14] B. Baker and R. Muth, “Compressing differences of executable code,” 2012.
  * [15] Y. David and E. Yahav, “Tracelet-based code search in executables,” in Proceedings of the 35th ACM SIGPLAN Conference on Programming Language Design and Implementation, PLDI ’14, (New York, NY, USA), p. 349–360, Association for Computing Machinery, 2014.
  * [16] H. Huang, A. M. Youssef, and M. Debbabi, “Binsequence: Fast, accurate and scalable binary code reuse detection,” in Proceedings of the 2017 ACM on Asia Conference on Computer and Communications Security, ASIA CCS ’17, (New York, NY, USA), p. 155–166, Association for Computing Machinery, 2017.
  * [17] M. Bourquin, A. King, and E. Robbins, “Binslayer: Accurate comparison of binary executables,” in Proceedings of the 2nd ACM SIGPLAN Program Protection and Reverse Engineering Workshop, PPREW ’13, (New York, NY, USA), Association for Computing Machinery, 2013.
  * [18] S. Cesare, Y. Xiang, and W. Zhou, “Control flow-based malware variantdetection,” IEEE Transactions on Dependable and Secure Computing, vol. 11, no. 4, pp. 307–317, 2014.
  * [19] Q. Feng, R. Zhou, C. Xu, Y. Cheng, B. Testa, and H. Yin, “Scalable graph-based bug search for firmware images,” in Proceedings of the 2016 ACM SIGSAC Conference on Computer and Communications Security, CCS ’16, (New York, NY, USA), p. 480–491, Association for Computing Machinery, 2016.
  * [20] M. Lindorfer, A. Di Federico, F. Maggi, P. M. Comparetti, and S. Zanero, “Lines of malicious code: Insights into the malicious software industry,” in Proceedings of the 28th Annual Computer Security Applications Conference, ACSAC ’12, (New York, NY, USA), p. 349–358, Association for Computing Machinery, 2012.
  * [21] W. M. Khoo, A. Mycroft, and R. Anderson, “Rendezvous: A search engine for binary code,” in 2013 10th Working Conference on Mining Software Repositories (MSR), pp. 329–338, 2013.
  * [22] S. Alrabaee, P. Shirani, L. Wang, and M. Debbabi, “Fossil: A resilient and efficient system for identifying foss functions in malware binaries,” ACM Trans. Priv. Secur., vol. 21, Jan. 2018.
  * [23] L. Luo, J. Ming, D. Wu, P. Liu, and S. Zhu, “Semantics-based obfuscation-resilient binary code similarity comparison with applications to software plagiarism detection,” in Proceedings of the 22nd ACM SIGSOFT International Symposium on Foundations of Software Engineering, FSE 2014, (New York, NY, USA), p. 389–400, Association for Computing Machinery, 2014.
  * [24] S. Alrabaee, P. Shirani, L. Wang, and M. Debbabi, “Sigma: A semantic integrated graph matching approach for identifying reused functions in binary code,” Digit. Investig., vol. 12, pp. S61–S71, 2015.
  * [25] J. Gao, X. Yang, Y. Fu, Y. Jiang, and J. Sun, “Vulseeker: A semantic learning based vulnerability seeker for cross-platform binary,” in Proceedings of the 33rd ACM/IEEE International Conference on Automated Software Engineering, ASE 2018, (New York, NY, USA), p. 896–899, Association for Computing Machinery, 2018.
  * [26] X. Xu, C. Liu, Q. Feng, H. Yin, L. Song, and D. Song, “Neural network-based graph embedding for cross-platform binary code similarity detection,” Proceedings of the 2017 ACM SIGSAC Conference on Computer and Communications Security, 2017.
  * [27] W. Jin, S. Chaki, C. Cohen, A. Gurfinkel, J. Havrilla, C. Hines, and P. Narasimhan, “Binary function clustering using semantic hashes,” in Proceedings of the 2012 11th International Conference on Machine Learning and Applications - Volume 01, ICMLA ’12, (USA), p. 386–391, IEEE Computer Society, 2012.
  * [28] J. Pewny, B. Garmany, R. Gawlik, C. Rossow, and T. Holz, “Cross-architecture bug search in binary executables,” it - Information Technology, vol. 59, pp. 83 – 91, 2017.
  * [29] M. Chandramohan, Y. Xue, Z. Xu, Y. Liu, C. Y. Cho, and H. B. K. Tan, “Bingo: cross-architecture cross-os binary search,” Proceedings of the 2016 24th ACM SIGSOFT International Symposium on Foundations of Software Engineering, 2016.
  * [30] U. Kargén and N. Shahmehri, “Towards robust instruction-level trace alignment of binary code,” 2017 32nd IEEE/ACM International Conference on Automated Software Engineering (ASE), pp. 342–352, 2017.
  * [31] S. Wang and D. Wu, “In-memory fuzzing for binary code similarity analysis,” in Proceedings of the 32nd IEEE/ACM International Conference on Automated Software Engineering, ASE 2017, p. 319–330, IEEE Press, 2017.
  * [32] S.-C. Wang, C. Liu, Y. yu Li, and W. Xu, “Semdiff: Finding semtic differences in binary programs based on angr,” 2017.
  * [33] Y. David, N. Partush, and E. Yahav, “Statistical similarity of binaries,” SIGPLAN Not., vol. 51, p. 266–280, June 2016.
  * [34] A. Lakhotia, M. D. Preda, and R. Giacobazzi, “Fast location of similar code fragments using semantic ’juice’,” in Proceedings of the 2nd ACM SIGPLAN Program Protection and Reverse Engineering Workshop, PPREW ’13, (New York, NY, USA), Association for Computing Machinery, 2013.
  * [35] Y. David, N. Partush, and E. Yahav, “Similarity of binaries through re-optimization,” Proceedings of the 38th ACM SIGPLAN Conference on Programming Language Design and Implementation, 2017.
  * [36] Q. Feng, M. Wang, M. Zhang, R. Zhou, A. Henderson, and H. Yin, “Extracting conditional formulas for cross-platform bug search,” in Proceedings of the 2017 ACM on Asia Conference on Computer and Communications Security, ASIA CCS ’17, (New York, NY, USA), p. 346–359, Association for Computing Machinery, 2017.
  * [37] J. Pewny, F. Schuster, L. Bernhard, T. Holz, and C. Rossow, “Leveraging semantic signatures for bug search in binary programs,” in Proceedings of the 30th Annual Computer Security Applications Conference, ACSAC ’14, (New York, NY, USA), p. 406–415, Association for Computing Machinery, 2014.


# Too big, too small or just right? A benchmark assessment of density
functional theory for predicting the spatial extent of the electron density of
small chemical systems.

Diptarka Hait<EMAIL_ADDRESS>These authors contributed equally to this
work.  Kenneth S. Pitzer Center for Theoretical Chemistry, Department of
Chemistry, University of California, Berkeley, California 94720, USA  Chemical
Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California
94720, USA  Yu Hsuan Liang  These authors contributed equally to this work.
Kenneth S. Pitzer Center for Theoretical Chemistry, Department of Chemistry,
University of California, Berkeley, California 94720, USA  Martin Head-Gordon
<EMAIL_ADDRESS>Kenneth S. Pitzer Center for Theoretical Chemistry,
Department of Chemistry, University of California, Berkeley, California 94720,
USA  Chemical Sciences Division, Lawrence Berkeley National Laboratory,
Berkeley, California 94720, USA

###### Abstract

Multipole moments are the first-order responses of the energy to spatial
derivatives of the electric field strength. The quality of density functional
theory (DFT) prediction of molecular multipole moments thus characterizes
errors in modeling the electron density itself, as well as the performance in
describing molecules interacting with external electric fields. However, only
the lowest non-zero moment is translationally invariant, making the higher-
order moments origin-dependent. Therefore, instead of using the $3\times 3$
quadrupole moment matrix, we utilize the translationally invariant $3\times 3$
matrix of second cumulants (or spatial variances) of the electron density as
the quantity of interest (denoted by $\bm{\mathcal{K}}$). The principal
components of $\bm{\mathcal{K}}$ are the square of the spatial extent of the
electron density along each axis. A benchmark dataset of the principal
components of $\bm{\mathcal{K}}$ for 100 small molecules at the coupled
cluster singles and doubles with perturbative triples (CCSD(T)) at the
complete basis set (CBS) limit is developed, resulting in 213 independent
$\bm{\mathcal{K}}$ components. The performance of 47 popular and recent
density functionals is assessed against this Var213 dataset. Several
functionals, especially double hybrids, and also SCAN and SCAN0 yield reliable
second cumulants, although some modern, empirically parameterized functionals
yield more disappointing performance. The H and Be atoms in particular are
challenging for nearly all methods, indicating that future functional
development could benefit from inclusion of their density information in
training or testing protocols.

##  I Introduction

Quantum chemistry methods have seen increasingly widespread use over the last
few decades, and now permit sufficiently accurate computation of energies
ranging in scale from total atomization energiesKarton, Daon, and Martin
(2011) to extremely weak noncovalent interactionsRezac and Hobza (2013). In
particular, density functional theoryBecke (2014); Jones (2015); Mardirossian
and Head-Gordon (2017) (DFT) is ubiquitously employed to study medium to large
systems as it provides an acceptable balance between accuracy and
computational cost. Although exact in theoryHohenberg and Kohn (1964); Levy
(1979), practical usage of DFT almost always entails use of computationally
tractable density functional approximations (DFAs) within the Kohn-Sham (KS)
formalismKohn and Sham (1965). It is often remarked that DFAs are not really
systematically improvable the way the exact coupled cluster (CC) hierarchy is,
but the best-performing DFAs on a given rung of Jacob’s ladderPerdew and
Schmidt (2001) statistically improve upon lower rungsMardirossian and Head-
Gordon (2017); Goerigk  _et al._ (2017); Mardirossian and Head-Gordon (2018);
Najibi and Goerigk (2018), indicating that greater complexity could (but is
not guaranteed to) lead to better accuracy. Caution is still warranted while
applying KS-DFT to systems with substantial levels of delocalization
errorPerdew  _et al._ (1982); Mori-Sánchez, Cohen, and Yang (2006); Hait and
Head-Gordon (2018a); Lonsdale and Goerigk (2020) or multireference
characterCohen, Mori-Sánchez, and Yang (2008a, b); Hait, Rettig, and Head-
Gordon (2019a), but DFAs are quite effective overall at predicting ground
state relative energies associated with main-group chemistryMardirossian and
Head-Gordon (2017); Goerigk  _et al._ (2017), and even organometallic
chemistry.Dohm  _et al._ (2018); Chan, Gill, and Kimura (2019) Indeed, modern
DFAs often surpass the accuracy of single-reference wave function approaches
like second-order Møller-Plesset perturbation theory (MP2) or CC singles and
doubles (CCSD) for such problemsMardirossian and Head-Gordon (2018); Hait and
Head-Gordon (2018b); Rettig  _et al._ (2020).

This is perhaps unsurprising, as DFA development has historically emphasized
improved prediction of such chemically relevant energies, out of the belief
that this leads to better approximations to the exact functional. However, the
exact functional ought to map the exact density to the exact energy, and it is
entirely possible for a DFA to obtain a reasonable energy from a relatively
inaccurate density via hidden cancellation of errors. It has indeed been
suggested that modern DFAs predict worse densities than older, less
parameterized modelsMedvedev  _et al._ (2017a)-leading to considerable
discussionHammes-Schiffer (2017); Kepp (2017); Medvedev  _et al._ (2017b);
Graziano (2017); Korth (2017); Brorsen  _et al._ (2017); Wang  _et al._
(2017a); Ranasinghe, Perera, and Bartlett (2017); Mayer  _et al._ (2017);
Gould (2017); Hait and Head-Gordon (2018b) in the scientific community about
density predictions from DFT. It is certainly plausible that greater
complexity in the functional form creates additional avenues for incorrect
behavior (as can be seen in wave function theory as wellLeininger  _et al._
(2000); Hait  _et al._ (2019)), even without accounting for overfitting to
empirical data.

Separate from the debate about “true paths” for DFA development,Medvedev  _et
al._ (2017a) we note that densities are important for a very practical reason.
The electron density controls the response of the electronic energy to
external electric fields, such as those that might be encountered in condensed
phases or in spectroscopic simulations. Density errors are thus related to
energy errors under a different one-particle potential (that DFT should still
be formally capable of addressingHohenberg and Kohn (1964)). Application of
DFT beyond ground state gas-phase chemical physics is thus reliant upon DFT
yielding reasonable densities, or at least a characterization of when and how
various DFAs fail. This has led to investigations on the performance of DFT
for predicting dipole momentsHickey and Rowley (2014); Verma and Truhlar
(2017); Hait and Head-Gordon (2018b); Johnson  _et al._ (2019); Grotjahn  _et
al._ (2020); Hait and Head-Gordon (2018c) and polarizabilitiesHickey and
Rowley (2014); Thakkar and Wu (2015); Hait and Head-Gordon (2018d); Grotjahn
_et al._ (2020); Withanage  _et al._ (2019); Hait, Rettig, and Head-Gordon
(2019b), as they represent the linear and quadratic responses, respectively,
of the energy to a spatially uniform electric field.

The natural next step would be to assess the linear response of the energy to
a (spatially) constant electric field gradient, which is equivalent to the
molecular quadrupole moment $\mathbf{Q}$. Indeed, there have been a number of
assessments of the quality of various density functionals for evaluating
molecular quadrupole moments in the pastCohen and Tantirungrotechai (1999); De
Proft, Tielens, and Geerlings (2000); Hofinger and Wendland (2002); Matta
(2010), albeit before the development of many modern functionals. At the same
time, efforts have been made to converge high-quality coupled cluster theory
calculations of quadrupole moments with respect to basis set,Halkier  _et al._
(1997); Zhang, Calvin, and Valeev (2016); Bokhan, Trubnikov, and Bartlett
(2016) setting the stage for benchmark assessments. $\mathbf{Q}$ is however
not translationally invariant for polar systems, making it a somewhat
unsuitable metric. We thus decided to investigate the performance of DFT in
predicting the second spatial cumulant of the electron density
($\bm{\mathcal{K}}$) instead, which is a translationally invariant quantity
connected to $\mathbf{Q}$ and the dipole moment $\bm{\mu}$. $\bm{\mathcal{K}}$
thus relates the density to the response of the energy to electric fields that
linearly vary with spatial coordinates. Mathematically, $\bm{\mathcal{K}}$ is
also the variance of the electron density. To the best of our knowledge, there
have been no previous studies that explicitly focus on prediction of
$\bm{\mathcal{K}}$ by DFT or MP/CC methods. However, a related quantity had
previously been studied with Hartree-Fock (HF)Hollett, Kelly, and Poirier
(2006) and found to be useful in predicting steric effects of substituents.

##  II Second spatial cumulant of electron density

###  II.1 Definition

Let us consider a system with electron density $\rho(\mathbf{r})$ and total
number of electrons $N=\displaystyle\int\rho(\mathbf{r})d\mathbf{r}$. The
spatial probability density $p(\mathbf{r})$ for a single electron is:

$\displaystyle p(\mathbf{r})=\dfrac{\rho(\mathbf{r})}{N}$ (1)

An alternative definition of $p(\mathbf{r})$ in terms of the wave function
$\left|\Psi\right\rangle$ is:

$\displaystyle p(\mathbf{r})$
$\displaystyle=\left\langle\mathbf{r}\right|\left(\underset{N-1}{\mathbf{Tr}}\left[\left|\Psi\right\rangle\left\langle\Psi\right|\right]\right)\left|\mathbf{r}\right\rangle=\dfrac{\Gamma_{1}(\mathbf{r},\mathbf{r})}{N}$
(2)

which corresponds to the scaled diagonal elements of the one particle density
operator, $\Gamma_{1}(\mathbf{r},\mathbf{r^{\prime}})$ that results from
tracing out $N-1$ degrees of freedom from the full density operator
$\Gamma_{N}=\left|\Psi\right\rangle\left\langle\Psi\right|$. For a single
Slater determinant, this is equivalent to averaging the square of each
occupied spin orbital.

The first and second spatial moments of the electron density are consequently:

$\displaystyle\langle\mathbf{r}\rangle$
$\displaystyle=\displaystyle\int\mathbf{r}p(\mathbf{r})d\mathbf{r}$ (3)
$\displaystyle\langle\mathbf{r}\mathbf{r^{T}}\rangle$
$\displaystyle=\displaystyle\int\left(\mathbf{r}\mathbf{r^{T}}\right)p(\mathbf{r})d\mathbf{r}$
(4)

The second spatial cumulant of the electron density is thus:

$\displaystyle\bm{\mathcal{K}}$
$\displaystyle=\langle\left(\delta\mathbf{r}\right)\left(\delta\mathbf{r}\right)^{\mathbf{T}}\rangle=\langle\mathbf{r}\mathbf{r^{T}}\rangle-\langle\mathbf{r}\rangle\langle\mathbf{r}\rangle^{\mathbf{T}}$
(5)

which is equivalent to the spatial (co)variance of the probability
distribution $p(\mathbf{r}$). The individual components are thus of the form:

$\displaystyle\mathcal{K}_{xx}$ $\displaystyle=\langle x^{2}\rangle-\langle
x\rangle^{2}=\displaystyle\int
x^{2}p(\mathbf{r})d\mathbf{r}-\left(\displaystyle\int
xp(\mathbf{r})d\mathbf{r}\right)^{2}$ (6) $\displaystyle\mathcal{K}_{xy}$
$\displaystyle=\langle xy\rangle-\langle x\rangle\langle
y\rangle=\displaystyle\int xyp(\mathbf{r})d\mathbf{r}-\left(\displaystyle\int
xp(\mathbf{r})d\mathbf{r}\right)\left(\displaystyle\int
yp(\mathbf{r})d\mathbf{r}\right)$ (7)

etc., in terms of the individual cartesian directions $x,y,z$.
$\bm{\mathcal{K}}$ is thus translationally invariant, and represents the
“width” (or “spread”) of the electron density. We also note that the
eigenvalues of $\bm{\mathcal{K}}$ are rotationally invariant. In addition,
while we have only focused on the second-order cumulant in this work, higher-
order cumulants of $p(r)$ can similarly be readily defined and utilized, for
problems where they might be relevant.

###  II.2 Connection to Multipole Moments

If the system has nuclear charges $\\{Z_{A}\\}$ at positions
$\\{\mathbf{R_{A}}\\}$, the dipole moment $\bm{\mu}$ and (cartesian/non-
traceless) quadrupole moment $\mathbf{Q}$ are:

$\displaystyle\bm{\mu}$
$\displaystyle=\displaystyle\sum\limits_{A}Z_{A}\mathbf{R_{A}}-\displaystyle\int\mathbf{r}\rho(\mathbf{r})d\mathbf{r}=\displaystyle\sum\limits_{A}Z_{A}\mathbf{R_{A}}-N\langle\mathbf{r}\rangle$
(8) $\displaystyle\mathbf{Q}$
$\displaystyle=\displaystyle\sum\limits_{A}Z_{A}\left(\mathbf{R_{A}}\right)\left(\mathbf{R_{A}}\right)^{\mathbf{T}}-\displaystyle\int\left(\mathbf{r}\mathbf{r^{T}}\right)\rho(\mathbf{r})d\mathbf{r}=\displaystyle\sum\limits_{A}Z_{A}\left(\mathbf{R_{A}}\right)\left(\mathbf{R_{A}}\right)^{\mathbf{T}}-N\langle\mathbf{r}\mathbf{r^{T}}\rangle$
(9)

Therefore the second cumulant can be expressed as:

$\displaystyle\bm{\mathcal{K}}$
$\displaystyle=\dfrac{1}{N}\left(\displaystyle\sum\limits_{A}Z_{A}\left(\mathbf{R_{A}}\right)\left(\mathbf{R_{A}}\right)^{\mathbf{T}}-\mathbf{Q}\right)-\dfrac{1}{N^{2}}\left(\displaystyle\sum\limits_{A}Z_{A}\mathbf{R_{A}}-\bm{\mu}\right)\left(\displaystyle\sum\limits_{A}Z_{A}\mathbf{R_{A}}-\bm{\mu}\right)^{\mathbf{T}}$
(10)

###  II.3 Physical interpretation

Let us begin with the case of atoms, and suppose that the single electron
probability distribution for a (spherical) atom at $(0,0,0)$ is
$p_{1}(\mathbf{r})$. Spherical symmetry indicates that:

$\displaystyle\displaystyle\int xp_{1}(\mathbf{r})d\mathbf{r}$
$\displaystyle=\displaystyle\int
yp_{1}(\mathbf{r})d\mathbf{r}=\displaystyle\int
zp_{1}(\mathbf{r})d\mathbf{r}=0$ (11) $\displaystyle\displaystyle\int
xyp_{1}(\mathbf{r})d\mathbf{r}$ $\displaystyle=\displaystyle\int
yzp_{1}(\mathbf{r})d\mathbf{r}=\displaystyle\int
xzp_{1}(\mathbf{r})d\mathbf{r}=0$ (12) $\displaystyle\displaystyle\int
x^{2}p_{1}(\mathbf{r})d\mathbf{r}$ $\displaystyle=\displaystyle\int
y^{2}p_{1}(\mathbf{r})d\mathbf{r}=\displaystyle\int
z^{2}p_{1}(\mathbf{r})d\mathbf{r}=\eta^{2}$ (13)

This implies that
$\mathcal{K}_{xx}=\mathcal{K}_{yy}=\mathcal{K}_{zz}=\eta^{2}$. Therefore
$\sqrt{3}\eta$ is an effective atomic radius associated with the spatial
extent of the electron density, much like a covalent or van der Waals radius.
In particular for the H atom, $\eta$ is the Bohr radius (1 a.u.=0.529 Å).

For molecules without symmetry, $\bm{\mathcal{K}}$ is not generally a diagonal
$3\times 3$ matrix, and likewise
$\mathcal{K}_{xx}\neq\mathcal{K}_{yy}\neq\mathcal{K}_{zz}$. $\bm{\mathcal{K}}$
can be diagonalized, to yield principal axes, and 3 principal components, very
analogous to the molecular inertia tensor associated with the nuclei,
$\mathbf{I}=\sum\limits_{A}m_{A}\left(\mathbf{R_{A}}\right)\left(\mathbf{R_{A}}\right)^{\mathbf{T}}$.
The eigenvalues of $\bm{\mathcal{K}}$ thereby give the square of the spatial
extent of the electron density along each principal axis. Systems where both a
Cn ($n>1$) axis of rotation and $\sigma_{v}$ plane(s) of symmetry are present
(i.e. Cnv or higher symmetry) have the principal axes be defined by molecular
symmetry, and we shall later choose members of our data set based on this
convenient simplification.

###  II.4 Behavior vs system size

The behavior of $\bm{\mathcal{K}}$ vs system size could be useful in
characterizing the utility of this quantity. A simple case to consider is a 1D
lattice of $M$ noninteracting He atoms, placed at
$(0,0,0),(0,0,a),(0,0,2a)\ldots(0,0,(M-1)a)$. Let each atom have spatial
extent $\eta$ as in Eq. (13).

Subsequently, $p(\mathbf{r})$ for the supersystem is:

$\displaystyle p(\mathbf{r})$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}p_{1}(\mathbf{r}-ma\hat{\mathbf{z}})$
(14)

which is the average of the individual single electron probability
distributions. Transverse to the lattice vector $a\hat{\mathbf{z}}$, we thus
have:

$\displaystyle\langle x^{2}\rangle$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\displaystyle\int
x^{2}p_{1}(\mathbf{r}-ma\hat{\mathbf{z}})d\mathbf{r}=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\displaystyle\int
x_{m}^{2}p_{1}(\mathbf{r}_{m})d\mathbf{r}_{m}=\eta^{2}$ (15)
$\displaystyle\langle x\rangle$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}x\displaystyle\int
p_{1}(\mathbf{r}-ma\hat{\mathbf{z}})d\mathbf{r}=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\displaystyle\int
x_{m}p_{1}(\mathbf{r}_{m})d\mathbf{r}_{m}=0$ (16)

(where $\mathbf{r}_{m}=\mathbf{r}-ma\hat{\mathbf{z}}$). Therefore
$\mathcal{K}_{xx}=\eta^{2}=\mathcal{K}_{yy}$ and is invariant vs system size.

Parallel to the lattice vector $a\hat{\mathbf{z}}$ however, we obtain:

$\displaystyle\langle z^{2}\rangle$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\displaystyle\int
z^{2}p_{1}(\mathbf{r}-ma\hat{\mathbf{z}})d\mathbf{r}$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\displaystyle\int\left(z_{m}+ma\right)^{2}p_{1}(\mathbf{r}_{m})d\mathbf{r}_{m}$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\left(\displaystyle\int
z_{m}^{2}p_{1}(\mathbf{r}_{m})d\mathbf{r}_{m}+m^{2}a^{2}\displaystyle\int
p_{1}(\mathbf{r}_{m})d\mathbf{r}_{m}+2am\displaystyle\int
z_{m}p_{1}(\mathbf{r}_{m})d\mathbf{r}_{m}\right)$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\left(\eta^{2}+m^{2}a^{2}\right)=\eta^{2}+\dfrac{(M-1)(2M-1)}{6}a^{2}$
(17) $\displaystyle\langle z\rangle$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\displaystyle\int
zp_{1}(\mathbf{r}-ma\hat{\mathbf{z}})d\mathbf{r}$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\displaystyle\int\left(z_{m}+ma\right)p_{1}(\mathbf{r}_{m})d\mathbf{r}_{m}$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}\left(\displaystyle\int
z_{m}p_{1}(\mathbf{r}_{m})d\mathbf{r}_{m}+ma\displaystyle\int
p_{1}(\mathbf{r}_{m})d\mathbf{r}_{m}\right)$
$\displaystyle=\dfrac{1}{M}\displaystyle\sum\limits_{m=0}^{M-1}ma=\dfrac{M-1}{2}a$
(18) $\displaystyle\therefore\mathcal{K}_{zz}$ $\displaystyle=\langle
z^{2}\rangle-\langle z\rangle^{2}=\eta^{2}+\dfrac{(M^{2}-1)}{12}a^{2}$ (19)

$\mathcal{K}_{zz}$ thus grows as $O(M^{2})$ vs the system size $M$. This
increase is however entirely due to geometric/structural factors (i.e. only
dependent on the lattice spacing $a$ and not on the electronic contribution
from $p_{1}(\mathbf{r})$). Thus differences between $\mathcal{K}_{zz}$
computed by various methods will be independent of $M$, and will solely be a
function of the subsystem densities $\\{p_{1}(\mathbf{r})\\}$. This is not
true in the interacting subsystems limit, but the analysis nonetheless reveals
a significant contribution to $\bm{\mathcal{K}}$ from molecular geometry
alone.

A similar analysis for 2D (square) and 3D (cubic) lattices shows that
$\mathcal{K}_{ii}$ grows as $O(M)$ and $O(M^{\frac{2}{3}})$ respectively vs
number of identical noninteracting subsystems $M$ (where $\hat{\mathbf{i}}$ is
parallel to the lattice vectors), due to geometric factors, while the
electronic contribution intrinsically arising from $p(\mathbf{r})$ remains
constant. Consequently, differences in $\bm{\mathcal{K}}$ should be size-
intensive in the non-interacting limit. However, relative error in
$\bm{\mathcal{K}}$ would shrink with increasing $M$ (as the reference value in
the denominator would increase). We consequently only consider absolute errors
of the form $\bm{\mathcal{K}}-\bm{\mathcal{K}}_{ref}$, vs a reference value
$\bm{\mathcal{K}}_{ref}$. This stands in contrast to properties like dipole
moments and static polarizabilities, which are properly size-extensive and
thus suitable for relative/percentage error based metrics. It is also possible
to look at standard deviations (i.e. $\sqrt{\mathcal{K}_{ii}}$) instead of
variances, but the geometric factors would prevent the errors from being size-
intensive in that case.

NSP | SP
---|---
AlF | Cl2 | Mg | AlH2 | NH
Ar | ClCN | Mg2 | BH2 | NH2
BF | ClF | N2 | BO | NO2
BH2Cl | FCN | NH3 | BS | NP
BH2F | H2 | NH3O | Be | Na
BH3 | H2O | NaCl | BeH | Na2
BHF2 | HBO | NaH | C2H | NaLi
BeH2 | HBS | Ne | CF2 | O2
C2H2 | HCCCl | OCl2 |  3CH2 | O3
C2H4 | HCCF | PH3 | CH3 | OF2
CH2BH | HCHO | PH3O | CN | P
CH2BO | HCN | SCl2 | F2 | P2
CH3Cl | HCl | SF2 | H | PCl
CH3F | HF | SH2 | H2CN | PF
CH3Li | HNC | SO2 | HCHS | PH
CH4 | He | SiH3Cl | HCP | PH2
CO | LiBH4 | SiH3F | Li | PO2
CO2 | LiCl | SiH4 | N | S2
CS | LiF | SiO | NCl | SO-trip
CSO | LiH |  | NF | SiH3
|  |  | NF2 |
Table 1: The 100 species in the dataset, sorted by whether they are not spin-
polarized (NSP) or spin-polarized (SP).

##  III Dataset

We have investigated $\bm{\mathcal{K}}$ for 100 small main-group systems
(listed in Table 1), for which it was possible to get highly accurate
benchmark values with CC singles and doubles with perturbative triples
(CCSD(T)Raghavachari  _et al._ (1989)) at the complete basis set (CBS) limit.
The systems were chosen from the combined datasets considered in Refs 22 and
43 such that $\bm{\mathcal{K}}$ was diagonal for a fixed coordinate system for
all (spatial symmetry preserving) electronic structure methods. Asymmetric
systems with nondiagonal $\bm{\mathcal{K}}$ were not included for simplicity.
Linear molecules with a single $\pi$ electron/hole (such as OH) were also not
considered, as real valued orbitals would be incapable of predicting a
cylindrically symmetric orbital of $\pi$ symmetry Lee  _et al._ (2019);
Brakestad  _et al._ (2020), and would thus spuriously lead to
$\mathcal{K}_{xx}\neq\mathcal{K}_{yy}$. Complex valued orbitals would thus be
necessary to describe such species with proper symmetry within a DFT
frameworkLee  _et al._ (2019).

The complete dataset consists of 59 not spin-polarized (NSP) and 41 spin-
polarized (SP) systems. The NSP vs SP classification was done on the basis of
whether the stable HF solution has $\langle S^{2}\rangle=0$ or not. NSP
species are thus unambiguously closed-shell, and are thus more likely to be
‘easier’ for single-reference quantum chemistry methods like KS-DFT or MP2.
Indeed, the size of the (T) correction to $\bm{\mathcal{K}}$ was quite small
for nearly all species (NSP or SP, as discussed later in Sec IV.1), indicating
that the multireference character of the chosen systems (if any) did not
strongly influence $\bm{\mathcal{K}}$ predictions and that CCSD(T) is likely
to be an adequately accurate benchmark. Furthermore, none of these systems
have pathological, delocalization driven qualitative failuresDutoi and Head-
Gordon (2006); Ruzsinszky  _et al._ (2006); Hait and Head-Gordon (2018b, c);
Vydrov, Scuseria, and Perdew (2007), making them reasonable choices for
understanding the behavior of DFAs in the regime where they are expected to
work well.

The error $\epsilon_{i,m}$ in an individual component $\mathcal{K}^{m}_{ii}$
for a given molecule $m$ (vs a reference value $\mathcal{K}^{m}_{ref,ii}$) is:

$\displaystyle\epsilon_{i,m}$
$\displaystyle=\mathcal{K}^{m}_{ii}-\mathcal{K}^{m}_{ref,ii}$ (20)

The cumulative errors over all molecules and directions are thus:

  1. 1.

Root mean square error (RMSE):
$\sqrt{\dfrac{1}{3N}\displaystyle\sum\limits_{m=1}^{N}\left(\epsilon_{x,m}^{2}+\epsilon_{y,m}^{2}+\epsilon_{z,m}^{2}\right)}$.

  2. 2.

Mean error (ME):
$\dfrac{1}{3N}\displaystyle\sum\limits_{m=1}^{N}\left(\epsilon_{x,m}+\epsilon_{y,m}+\epsilon_{z,m}\right)$.

  3. 3.

Maximum absolute error (MAX):
$\max\left(\left|\epsilon_{i,m}\right|\right)\forall i\in\\{x,y,z\\}$ and
$\forall m\in\mathbb{Z}^{+},m\leq N$.

All quantities are reported in atomic units (a.u.) unless specified otherwise.

###  III.1 Computational Details

All calculations were performed with the Q-Chem 5 package Shao  _et al._
(2015), using fixed geometries obtained from Refs 22 and 43 (also provided in
the supporting information). We examined the performance of 47 DFAs spanning
all five rungs of Jacob’s ladder (as can be seen from Table 2), ensuring
reasonable representation at each level. Individual DFAs were selected based
on (perceived) popularity, recency and performance over other benchmark
datasetsMardirossian and Head-Gordon (2017); Goerigk  _et al._ (2017); Hait
and Head-Gordon (2018b, d); Mardirossian and Head-Gordon (2018).

$\bm{\mathcal{K}}$ was obtained from $\mathbf{Q}$ computed for this work and
$\bm{\mu}$ obtained from Ref 22. $\mathbf{Q}$ for self-consistent field (SCF)
approaches like HF and non double hybrid DFAs were found analytically via
integrating over $\rho(\mathbf{r})$. $\mathbf{Q}$ for other methods (MP2, CC
or double hybrid DFT) was found via a central, two point finite difference
approach using a constant electric field gradient of $2\times 10^{-4}$ a.u.
(similar to Ref 22). Finite difference errors thus introduced appear to be
quite small, as the RMS deviation between analytic and finite-difference
CCSD/aug-cc-pCVTZDunning Jr (1989); Woon and Dunning Jr (1995); Peterson and
Dunning Jr (2002); Prascher  _et al._ (2011) $\mathcal{K}_{ii}$ is $2.2\times
10^{-4}$ a.u.

Comparison between wave function theory and various DFAs are only meaningful
at the CBS limit, due to different rates of basis set convergence. HF/aug-cc-
pCV5Z $\mathcal{K}_{ii}$ were assumed to be at the CBS limit as it only has an
RMS deviation of $2.5\times 10^{-4}$ a.u vs aug-cc-pCVQZ results (which should
be sufficiently small, in light of the empirically observed exponential
convergence of HF energies vs cardinal number of the basis setJensen (2005);
Karton and Martin (2006)). Similarly, DFT/aug-pc-4Jensen (2001, 2002a, 2002b);
Jensen and Helgaker (2004); Jensen (2007, 2010) values were assumed to be at
the CBS limit for functionals from Rungs 1-4 of Jacob’s ladder. The virtual
orbital dependent correlation contribution $\bm{\mathcal{K}}^{\textrm{corr}}$
(in MP2/CC/double hybrids) was extrapolated to the CBS limit via the two point
extrapolation formulaHalkier  _et al._ (1999)
${\mathcal{K}^{\textrm{corr,}}_{ii}}^{\textrm{n}}={\mathcal{K}^{\textrm{corr,}}_{ii}}^{\textrm{CBS}}+\dfrac{A}{n^{3}}$
from aug-cc-pCVTZ ($n=3$) and aug-cc-pCVQZ ($n=4$) (which is adequate for
dipolesHait and Head-Gordon (2018b); Halkier  _et al._ (1999) and appears to
also be adequate for quadrupole momentsHalkier  _et al._ (1997); Zhang,
Calvin, and Valeev (2016); Bokhan, Trubnikov, and Bartlett (2016)).

Local exchange-correlation integrals for all DFT calculations were computed
over a radial grid with 99 points and an angular Lebedev grid with 590 points
for all atoms. Non-local correlation was evaluated on an SG-1 gridGill,
Johnson, and Pople (1993). Unrestricted (U) orbitals were employed for all CC
(except Be, where UHF breaks spatial symmetry) and non-double hybrid DFT
calculations. MP2 is known to yield non N-representable densities when spin-
symmetry breaksKurlancheek and Head-Gordon (2009), leading to poor dipoleHait
and Head-Gordon (2018b) and polarizabilityHait and Head-Gordon (2018d)
predictions. Consequently, MP2 with both restricted (R) and unrestricted
orbitals were carried out, and the restricted variant found to yield a
significantly smaller RMSE. Double hybrid calculations were subsequently
carried out with only R orbitals. Stability analysis was performed in the
absence of fields to ensure the orbitals correspond to a minimum, and the
resulting orbitals were employed as initial guesses for finite field
calculations (if any). The frozen-core approximation was not employed in any
calculation.

##  IV Results and discussions

Method | Class | RMSE | ME | MAX |  | Method | Class | RMSE | ME | MAX
---|---|---|---|---|---|---|---|---|---|---
|  | Full | NSP | SP |  |  |  |  |  | Full | NSP | SP |  |
CCSD | WFT | 0.003 | 0.002 | 0.003 | -0.001 | 0.022 |  |  SCAN0Hui and Chai (2016) | Rung 4 | 0.008 | 0.004 | 0.011 | 0.000 | 0.047
MP2 | WFT | 0.022 | 0.004 | 0.035 | 0.001 | 0.312 |  |  B97-2Wilson, Bradley, and Tozer (2001) | Rung 4 | 0.009 | 0.004 | 0.013 | 0.001 | 0.047
RMP2 | WFT | 0.010 | 0.004 | 0.015 | 0.002 | 0.074 |  |  TPSShStaroverov  _et al._ (2003) | Rung 4 | 0.009 | 0.004 | 0.013 | 0.002 | 0.053
HF | WFT | 0.027 | 0.015 | 0.039 | 0.007 | 0.301 |  |  PBE0Adamo and Barone (1999) | Rung 4 | 0.009 | 0.005 | 0.014 | 0.003 | 0.052
RHF | WFT | 0.020 | 0.015 | 0.026 | 0.006 | 0.091 |  |  B97Becke (1997) | Rung 4 | 0.010 | 0.006 | 0.013 | 0.005 | 0.049
|  |  |  |  |  |  |  |  HSE-HJSKrukau  _et al._ (2006) | Rung 4 | 0.010 | 0.005 | 0.014 | 0.003 | 0.051
SPW92Dirac (1931); Perdew and Wang (1992) | Rung 1 | 0.022 | 0.017 | 0.028 | 0.011 | 0.129 |  |  SOGGA11-XPeverati and Truhlar (2011a) | Rung 4 | 0.010 | 0.005 | 0.015 | 0.004 | 0.057
SlaterDirac (1931) | Rung 1 | 0.047 | 0.038 | 0.057 | 0.034 | 0.275 |  |  LRC-$\omega$PBEhRohrdanz, Martins, and Herbert (2009) | Rung 4 | 0.011 | 0.005 | 0.016 | 0.002 | 0.061
|  |  |  |  |  |  |  |  PW6B95Zhao and Truhlar (2005) | Rung 4 | 0.011 | 0.006 | 0.016 | 0.003 | 0.084
BPBEBecke (1988); Perdew, Burke, and Ernzerhof (1996) | Rung 2 | 0.013 | 0.008 | 0.017 | 0.004 | 0.062 |  |  $\omega$M05-DLin  _et al._ (2012) | Rung 4 | 0.011 | 0.007 | 0.015 | 0.005 | 0.064
mPW91Adamo and Barone (1998) | Rung 2 | 0.014 | 0.010 | 0.019 | 0.006 | 0.069 |  |  M06Zhao and Truhlar (2008a) | Rung 4 | 0.011 | 0.007 | 0.016 | 0.001 | 0.058
B97-D3Grimme (2006a) | Rung 2 | 0.016 | 0.010 | 0.022 | 0.007 | 0.103 |  |  $\omega$B97X-VMardirossian and Head-Gordon (2014) | Rung 4 | 0.012 | 0.008 | 0.016 | 0.006 | 0.066
PBEPerdew, Burke, and Ernzerhof (1996) | Rung 2 | 0.016 | 0.012 | 0.021 | 0.008 | 0.079 |  |  $\omega$B97X-DChai and Head-Gordon (2008) | Rung 4 | 0.012 | 0.006 | 0.017 | 0.003 | 0.064
N12Peverati and Truhlar (2012) | Rung 2 | 0.020 | 0.013 | 0.028 | 0.006 | 0.120 |  |  M06-2XZhao and Truhlar (2008a) | Rung 4 | 0.013 | 0.006 | 0.019 | 0.003 | 0.104
BLYPBecke (1988); Lee, Yang, and Parr (1988) | Rung 2 | 0.021 | 0.016 | 0.026 | 0.012 | 0.112 |  |  HFLYPLee, Yang, and Parr (1988) | Rung 4 | 0.014 | 0.013 | 0.015 | -0.004 | 0.051
SOGGA11Peverati, Zhao, and Truhlar (2011) | Rung 2 | 0.048 | 0.035 | 0.062 | 0.024 | 0.500 |  |  M11Peverati and Truhlar (2011b) | Rung 4 | 0.014 | 0.009 | 0.019 | 0.005 | 0.083
|  |  |  |  |  |  |  |  B3LYPBecke (1993) | Rung 4 | 0.014 | 0.009 | 0.019 | 0.007 | 0.068
SCANSun, Ruzsinszky, and Perdew (2015) | Rung 3 | 0.009 | 0.005 | 0.012 | 0.002 | 0.052 |  |  CAM-B3LYPYanai, Tew, and Handy (2004) | Rung 4 | 0.014 | 0.010 | 0.019 | 0.007 | 0.070
MS2Sun  _et al._ (2013) | Rung 3 | 0.010 | 0.004 | 0.014 | 0.000 | 0.072 |  |  M08-HXZhao and Truhlar (2008b) | Rung 4 | 0.015 | 0.009 | 0.020 | 0.009 | 0.079
TPSSTao  _et al._ (2003) | Rung 3 | 0.010 | 0.006 | 0.015 | 0.003 | 0.065 |  |  $\omega$B97M-VMardirossian and Head-Gordon (2016) | Rung 4 | 0.022 | 0.013 | 0.030 | 0.013 | 0.127
mBEEFWellendorff  _et al._ (2014) | Rung 3 | 0.011 | 0.005 | 0.017 | -0.001 | 0.074 |  |  MN15Haoyu  _et al._ (2016) | Rung 4 | 0.035 | 0.013 | 0.052 | 0.013 | 0.231
M06-LZhao and Truhlar (2006) | Rung 3 | 0.014 | 0.010 | 0.018 | -0.006 | 0.068 |  |  |  |  |  |  |  |
revM06-LWang  _et al._ (2017b) | Rung 3 | 0.015 | 0.013 | 0.018 | 0.002 | 0.091 |  |  DSD-PBEPBEKozuch and Martin (2013) | Rung 5 | 0.005 | 0.003 | 0.007 | 0.002 | 0.032
B97M-VMardirossian and Head-Gordon (2015) | Rung 3 | 0.026 | 0.013 | 0.037 | 0.011 | 0.157 |  |  XYGJ-OSZhang  _et al._ (2011) | Rung 5 | 0.006 | 0.003 | 0.008 | 0.002 | 0.028
M11-LPeverati and Truhlar (2011c) | Rung 3 | 0.038 | 0.023 | 0.052 | 0.005 | 0.332 |  |  PTPSSGoerigk and Grimme (2010) | Rung 5 | 0.006 | 0.003 | 0.009 | 0.002 | 0.032
MN15-LYu, He, and Truhlar (2016) | Rung 3 | 0.040 | 0.024 | 0.055 | 0.022 | 0.278 |  |  XYG3Zhang, Xu, and Goddard (2009) | Rung 5 | 0.006 | 0.003 | 0.009 | 0.002 | 0.032
|  |  |  |  |  |  |  |  $\omega$B97M(2)Mardirossian and Head-Gordon (2018) | Rung 5 | 0.008 | 0.005 | 0.010 | 0.003 | 0.037
|  |  |  |  |  |  |  |  B2GPPLYPKarton  _et al._ (2008) | Rung 5 | 0.008 | 0.005 | 0.011 | 0.004 | 0.042
|  |  |  |  |  |  |  |  $\omega$B97X-2(TQZ)Chai and Head-Gordon (2009) | Rung 5 | 0.008 | 0.006 | 0.010 | 0.005 | 0.044
|  |  |  |  |  |  |  |  B2PLYPGrimme (2006b) | Rung 5 | 0.010 | 0.007 | 0.013 | 0.005 | 0.047
Table 2: Errors in $\bm{\mathcal{K}}$ (in a.u) predicted by various DFAs for
the species in the dataset. Positive ME values indicate less compact (i.e.
more diffuse) densities, relative to the benchmark.

###  IV.1 Full dataset

The errors in DFA predictions for the full dataset are reported in Table 2,
along with errors for the wave function methods CCSD, MP2 and HF. Considering
the wave function theories (WFTs) first, we see that CCSD has the lowest RMSE
of 0.003 of all the methods considered, and gives fairly similar performance
across the NSP and SP subsets. MP2 gives good performance for the NSP dataset,
but N-representability failuresKurlancheek and Head-Gordon (2009) lead to
significantly worse results over the SP subset. Use of R orbitals ameliorates
this considerably, with RMP2 having a quite improved RMSE of 0.010 over the
full dataset. RMP2 is nonetheless still somewhat challenged by open-shell
systems (especially NF and NCl) where the artificial enforcement spin-symmetry
appears to be suboptimal. We further note that MP3 is also expected to have
the same N-representability failures as MP2 and is thus unlikely to lead to
significant improvements over MP2. This in fact has been observed for $\mu$,
although alternative orbital choices could lead to better MP
performance.Rettig  _et al._ (2020)

HF performs quite poorly due to lack of correlation, having an RMSE of 0.027
(which can be considered as a ceiling for reasonable DFA performance). Based
on the mean error (ME), HF has a tendency to make the variance too large (i.e.
HF densities are too diffuse). The SP species are much more challenging than
NSP, and RHF actually reduces RMSE substantially to 0.020. However, this is
almost solely on the account of two challenging alkali metal containing
species (Na2 and NaLi) and the RMSE for HF and RHF are quite similar upon
their exclusion. Nonetheless, this serves as a warning that spin-symmetry
breaking might compromise property predictions, despite being the natural
approach for improving the energy. In particular, spin-symmetry breaking in
diatomics like Na2 or NaLi results in densities similar to two independent
atoms vs a bonded molecule, leading to larger widths than the restricted
solution.

Coming to the DFAs, we observe that rung 1 local spin-density approximations
(LSDA) functionals fare considerably worse than CCSD or RMP2. Bare Slater
(LSDA) exchange has a rather large RMSE of 0.047, while inclusion of
correlation reduces RMSE substantially to 0.022 for SPW92. The SP subset is
significantly more challenging in both cases, with an RMSE that is almost
double the corresponding NSP value. The minimally parameterized nature of the
LSDA functionals mean that the SPW92 RMSE of 0.022 can also be considered as a
reasonable reference for judging DFAs against. In fact, the HF, RHF and SPW92
RMSEs collectively suggest that DFAs with RMSE larger than 0.02 are not really
accurate for predicting $\bm{\mathcal{K}}$.

Moving up Jacob’s ladder to rung 2 generalized gradient approximations (GGAs),
we observe that all but one of the functionals investigated improve upon
SPW92. SOGGA11 is the exception, being essentially as bad as bare Slater
exchange. On the other hand, BPBE and mPW91 have substantially lower RMSEs,
indicating considerable improvement over rung 1. However, the SP RMSE
continues to be nearly twice the NSP one for all GGAs.

Rung 3 meta-GGAs (mGGAs) further improve upon predictions, with the best mGGA
(SCAN) predicting a quite low RMSE of 0.009, while MS2 and TPSS are also
fairly reasonable. It is worth noting that all three were developed with a
heavy emphasis on non-empirical constraints and were specifically fit to the H
atom (one of the most challenging species in the dataset, as discussed later).
In contrast, modern, mGGAs developed principally by fitting to empirical
benchmark data appear to fare worse, with MN15-L, M11-L and B97M-V being
particularly disappointing as they fare worse than SPW92. Nonetheless, it is
worth noting that the modern, empirically fitted mBEEF functional yields a
respectable performance.

Hybrid functionals on rung 4 of Jacob’s ladder do not significantly improve
upon mGGAs. The best-performer is SCAN0, which only marginally improves upon
the parent SCAN functional. Similarly, TPSSh only slightly improves upon TPSS.
This general behavior stands in contrast to the case of dipole moments, where
hybrid functionals strongly improve upon lower rungsHait and Head-Gordon
(2018b). Some other decent performers are B97-2 (which was partially fitted to
densities) PBE0, B97, HSE-HJS and SOGGA11-X which also happen to be hybrid
GGAs. PW6B95 is the best-performing hybrid mGGA that was not already based on
an existing Rung 3 functional with good performance. Interestingly, HFLYP
improves significantly upon HF, mostly by virtue of dramatically reducing SP
errors (in large part because it does not break spin-symmetry for NaLi or
Na2). SP errors are roughly 2-3 times the NSP errors for nearly all other
hybrid functionals. It is also worth noting that that modern $\omega$B97M-V
and MN15 functionals fare particularly poorly, yielding performance worse than
SPW92. On the other hand, several empirically fitted functionals like M06 and
$\omega$B97X-V yield quite respectable performance (albeit significantly worse
than SCAN0).

Rung 5 double hybrid functionals reduce error significantly vs hybrids, with
DSD-PBEPBE approaching CCSD levels of RMSE for NSP species (although SP
species are much more challenging). Even B2PLYP (the poorest performer, and
coincidentally, the oldest) performs similarly to a good hybrid functional
like PBE0. Most rung 5 functionals also improve upon RMP2 (especially for the
SP subset), indicating perceptible benefit from the local exchange-correlation
contributions. It is also worth noting that the recently developed
$\omega$B97M(2) functional yields mediocre performance, despite being one of
the best-performers for energy predictionsMardirossian and Head-Gordon (2018)
(although it represents an enormous improvement over the parent $\omega$B97M-V
functional).

The positive ME values in Table 2 also indicates that most DFAs predict
slightly less compact densities (especially LSDA). This is perhaps not too
surprising in light of delocalization error present in the studied
functionalsHait and Head-Gordon (2018a). However, the connection between
density compactness and delocalization error is not always straightforward, as
shown later. Indeed, it can be seen that the local M06-L and mBEEF functionals
systematically predict too small $\bm{\mathcal{K}}$, which is contrary to
expectations based on delocalization error alone. HFLYP however has a negative
ME, consistent with overlocalization of density that is expected from 100% HF
exchange. On the other hand, HF systematically overestimates
$\bm{\mathcal{K}}$, which seems puzzling at a first glace as HF should
overlocalize if anythingLi and Yang (2017); Hait and Head-Gordon (2018a).
However, lack of correlation in HF leads to artificial symmetry breaking where
densities become more ‘atom-like’ than ‘bond-like’, leading to spurious
overestimation of $\bm{\mathcal{K}}$ for NaLi and Na2. In fact, lack of
correlation hinders electrons from occupying the same region of space and
could lead to wider ‘spread’ in the electron density.

Comparison to the results obtained earlier for dipole momentsHait and Head-
Gordon (2018b) indicates some similarities and differences. Ascending Jacob’s
ladder leads to improved predictions in both cases, and several modern,
empirically designed functionals are found to perform relatively poorly. In
fact, some of the functionals that were found to be amongst the top performers
for their rung in the dipole study (SCAN, DSD-PBEPBE, PBE0, SOGGA11-X etc.)
continue to perform well. The major differences are that hybrid functionals do
not significantly improve upon mGGAs and some functionals that are good for
$\bm{\mathcal{K}}$ predictions do not fare as well for dipoles. A clear
example of this is SCAN0 faring worse at predicting $\bm{\mu}$ than
$\omega$B97M-V, despite the former being the best hybrid for predicting
$\bm{\mathcal{K}}$ (while the latter is amongst the worst). Nonetheless, the
ability of several functionals to predict both $\bm{\mathcal{K}}$ and
$\bm{\mu}$ with low error is encouraging, as that indicates reasonable
performance for problems involving external electric fields.

It is also worthwhile to identify what species in the dataset are most
challenging for DFAs. Ordering the molecules in the dataset by the first
quartile of RMS error (over cartesian axes) across different DFAs reveals that
there is a break in the distribution after the first 10 species. These
difficult cases are: Be, H, Li, NaLi, Na2, BeH, LiH, H2, BeH2 and BF in
descending order of difficulty. Several of these cases are also known to be
challenging for dipole momentHait and Head-Gordon (2018b) and static
polarizabilityHait and Head-Gordon (2018d) predictions as well. NaLi in
particular was already known to be a very challenging case for many DFAs, and
it is unsurprising that the analogous Na2 is also challenging. Both of these
alkali metal dimers feature long bonds ($\sim$ 3 Å) and are SP at equilibrium
geometry, suggesting some multireference character (analogous to the isovalent
case of stretched H2). In addition, Na2 and NaLi are outliers with respect to
RMS deviation of CCSD $\mathcal{K}$ vs the CCSD(T) benchmark, having $\sim$
0.01 deviation (while the next largest deviation is about half of that, for
CH3Li). Near-exact calculations with the adaptive sampling configuration
interaction methodTubman  _et al._ (2016, 2020) (with the small cc-pVDZ
basisPrascher  _et al._ (2011)) however indicate that the species are not
particularly multireference (a single determinant has $\sim$ 90% weight in the
total wave function) and CCSD(T) is sufficiently accurate for $\mathcal{K}$
(as shown in the supporting information). Furthermore, the three atoms
represent the most challenging cases overall, with Li having a first quartile
RMSE $\sim$ 35% larger than NaLi.

###  IV.2 The case of challenging atoms

Figure 1: $\mathcal{K}_{ii}$ predicted for the H atom, by various DFAs
(ordered by their position in Table 2). The dashed line is the analytic value.

It is thus worthwhile to examine these atoms in greater detail for better
insight into functional behavior. H has only one electron and thus all DFA
errors are definitionally self-interaction driven. It is thus interesting to
observe that the error in $\bm{\mathcal{K}}$ does not appear to correlate well
with the delocalization error in many cases, as can be seen from Fig 1.
Perhaps the most clear example of this is a comparison between the local M06-L
functional and the range separated hybrid $\omega$B97M-V, with the former
surprisingly predicting a smaller H atom than exact quantum mechanics, while
the latter overestimates the size to essentially the same extent as SPW92
(contrary to the behavior seen for delocalization errorHait and Head-Gordon
(2018a)). Similarly, M06-2X predicts a larger H atom than M06 despite having
twice the HF exchange (54% vs 27%). It is thus apparent that the local
exchange-correlation components of several modern density functionals lead to
significant errors for the size of the H atom, which are counterintuitive from
the perspective of delocalization error (or fraction of HF exchange present).

On the other hand, if we only consider functionals that use the same local
exchange-correlation components hybridized with varying amounts of HF exchange
(such as PBE/PBE0), we find that error decreases with increase in the HF
exchange contribution. This is on account of the functional becoming strictly
more HF like. It is also worth noting that the double hybrid functionals
appear to have much lower error than hybrids, potentially due to relatively
smaller contributions from local exchange-correlation. While cancellation of
errors in larger systems is likely to make this less of an issue (as can be
seen from other species in the dataset), future density functional development
would likely benefit from including the size and polarizabilityHait and Head-
Gordon (2018d) of the H atom as soft constraints.

Figure 2: $\mathcal{K}_{ii}$ predicted for the Li atom, by various DFAs. The
dashed line is the benchmark (CCSD(T)/CBS) value and the HF value is 2.070.

Li superficially resembles H if the core electrons are neglected. However, Fig
2 shows that the error $\bm{\mathcal{K}}$ predictions are quite different.
Most functionals underestimate the size of the diffuse alkali metal atom,
predicting a more compact density than the benchmark. Delocalization errors
are normally expected to lead to less compact densities due to self-
interaction. This can indeed be seen from the overly repulsive mean-field
electron-electron repulsion potentials predicted by many DFAsGaiduk, Mizzi,
and Staroverov (2012); Gritsenko, Mentel, and Baerends (2016). The Li atom is
also not multireference, and thus the reason behind the size underestimation
by many functionals is not entirely clear.  HFLYP not being an outlier with
respect to overlocalization indicates that delocalization is unlikely to be a
major factor (and bare HF in fact overestimates the size). Nonetheless, most
of the modern functionals that fare poorly for H (B97M-V, MN15-L, MN15 and
$\omega$B97M-V) continue to fare poorly for Li and overestimate the size,
revealing considerable scope for improvement. The size of the Li atom is thus
another reasonable choice as a soft constraint for future functional fitting.

Figure 3: $\mathcal{K}_{ii}$ predicted for the Be atom, by various DFAs. The
dashed line is the benchmark (CCSD(T)/CBS) value and the HF value is 1.444.

In contrast to the preceding two cases, Fig 3 shows that most methods
overestimate the size of the Be atom by a similar amount, with some poor
performing outliers. In fact, the best mGGA and hybrid have $\sim 0.001$
variation in $\mathcal{K}_{ii}$ (while the best double hybrid is better by
$\sim 0.005$). This relatively uniform performance is likely a consequence of
the challenges faced by DFAs in accounting for the multireference character of
the Be atom (it has 0.37 effectively unpaired electronsHead-Gordon (2003)). It
is nonetheless worth noting that the worst performers are recently
parameterized modern functionals. It is not entirely clear whether the Be atom
should feature in datasets used for functional training because of the partial
multireference character of this problem, but it should definitely be employed
in test sets to gauge fit quality.

(a) Exact $\mathcal{K}_{zz}$ values, along with HF.

(b) Errors in $\mathcal{K}_{zz}$ predicted by various DFAs.

Figure 4: Behavior of $\mathcal{K}_{zz}$ for stretched H2.

##  V Stretched H2

The behavior of DFAs in predicting $\bm{\mathcal{K}}$ in non-equilibrium
configurations is potentially interesting. We only investigate the case of H2,
to avoid qualitative failures associated with delocalization error for polar
bondsDutoi and Head-Gordon (2006); Ruzsinszky  _et al._ (2006); Hait and Head-
Gordon (2018b, c) for general bond dissociations. We furthermore restrict
ourselves to a few well-behaved functional at Rung 4 and below, to avoid
incomplete spin-polarization based problemsHait, Rettig, and Head-Gordon
(2019b) . Double hybrid functionals are not considered due to
N-representability problems in UMP2,Kurlancheek and Head-Gordon (2009) and
divergence of RMP2 with bond stretching (indicating that they are quite
unsuitable for bond dissociation problems). Fig 4 shows that
$\mathcal{K}_{zz}$ grows quadratically with the internuclear distance $r$, and
Fig 4 thus only plots the error in $\mathcal{K}_{zz}$ for clarity. All of the
single-reference methods encounter a derivative discontinuity at the Coulson-
Fischer (CF) pointCoulson and Fischer (1949) (which is the onset of spin-
polarization) due to the orbital stability matrix being singular at that
point. The behavior right beyond this point is the most interesting, as the
two electrons are fairly strongly correlated in this regime. The DFAs other
than SPW92 yield remarkably similar behavior in this region, systematically
overestimating $\mathcal{K}_{zz}$ somewhat (and thus predicting more ‘atom’
like densities as opposed to a still partially bonded molecule). HF exhibits
the same qualitative behavior, but has much larger errors due to complete
absence of correlation (indeed, the systematic overestimation past the CF
point can be clearly observed in Fig 4 as well). Similar overestimation in
$\bm{\mathcal{K}}$ is seen for the Be atom, which also has partial
multireference character. Consistent systematic overestimation of
$\bm{\mathcal{K}}$ by several well-behaved functionals with varying levels of
delocalization error could in fact potentially be a signature of strong
correlation, as the electrons are forced to be less compact in the absence of
explicit (i.e. post mean-field) correlation that allows them to occupy the
same region of space.

##  VI Discussion and Conclusion

The objective of this work was to characterize the electron density of small
molecules using a scalar metric that goes beyond the first moment of the
electron density (i.e. dipole moments). Since the molecular quadrupole moment
is origin-dependent for molecules with non-zero charge or dipole, we elected
instead to characterize the $3\times 3$ matrix of second cumulants, or spatial
variances, of the electron density, $\bm{\mathcal{K}}$. The eigenvalues of
$\bm{\mathcal{K}}$ are thus a measure of the square of the characteristic
extent of the molecular density along each principal axis.

We produced a benchmark dataset, Var213, that contains 213 benchmark values of
$\bm{\mathcal{K}}$, from 100 small molecules, evaluated with CCSD(T)/CBS
correlation based on extrapolation towards the complete basis set limit from
aug-cc-pCVTZ and aug-cc-pCVQZ calculations, combined with HF/aug-cc-pCV5Z.
These reference values were then used to assess $\bm{\mathcal{K}}$ predictions
from 47 density functional approximations (DFAs). Broadly, the molecules
studied here are relatively straightforward in terms of their electronic
structure. In the language of ref. 5, this dataset should be considered “easy”
for DFAs, rather than “difficult” due to the absence of strong correlation or
delocalization effects.

The results show that it is possible to obtain quite reasonable
$\bm{\mathcal{K}}$ values in many cases, even with modern density functionals.
However, non-empirical functionals, especially at the mGGA level via SCAN, and
at the hybrid level via SCAN0, seem to fare better than more empirically
parameterized models. This constitutes useful independent validation of the
quality of electron densities from these functionals. Double hybrid rung 5
functionals yield the best overall performance: indeed all double hybrids
tested match or exceed the performance of SCAN, and only the early B2PLYP
double hybrid fails to outperform SCAN0.

Interestingly, the performance of the modern MN15-L, MN15, B97M-V and
$\omega$B97M-V functionals over the studied dataset is disappointing,
suggesting that the vastness of the meta-GGA spaceMardirossian and Head-Gordon
(2015) has offered these functionals the opportunity to predict reasonable
energies from fairly poor densities. There is an interesting lesson to be
drawn from this outcome, particularly for the minimally parametrized meta-GGA-
based functionals such as B97M-V and $\omega$B97M-V whose 12-16 parameter form
was inferred from the use of large nunbers of chemically relevant energy
differences. Specifically, such forms are evidently not fully constrained by
the data used to generate them. Independent data such as $\bm{\mathcal{K}}$
may indeed be very useful in attempting to develop improved data-driven
functional designs.

Overall, our $\bm{\mathcal{K}}$ benchmark tests show that the performance of
the best-performing functional strictly improves on ascending Jacob’s ladder,
indicating that extra complexity has the potential to improve behavior. Use of
this dataset and other density based information (and especially the H, Li and
Be atoms) for future functional development could thus provide sufficient soft
constraints to yield improved $\bm{\mathcal{K}}$ values, and by implication
electron densities as well. The results obtained here may thereby contribute
to offering a route to better computationally tractable approximations to the
exact functional.

## Acknowledgment

This research was initially supported by the Director, Office of Science,
Office of Basic Energy Sciences, of the U.S. Department of Energy under
Contract No. DE-AC02-05CH11231. Completion of the project was based upon work
performed by the Liquid Sunlight Alliance, a DOE Energy Innovation Hub,
supported by the U.S. Department of Energy, Office of Science, Office of Basic
Energy Sciences, under Award Number DE-SC0021266. Y.H.L. was funded via the UC
Berkeley College of Chemistry summer research program.

## Supporting Information

Geometries of species studied (zip). Computed values and analysis (xlxs).

## Data Availability

The data that supports the findings of this study are available within the
article and its supplementary material.

## References

  * Karton, Daon, and Martin (2011) A. Karton, S. Daon, and J. M. Martin, Chem. Phys. Lett. 510, 165 (2011).
  * Rezac and Hobza (2013) J. Rezac and P. Hobza, J. Chem. Theory Comput. 9, 2151 (2013).
  * Becke (2014) A. D. Becke, J. Chem. Phys. 140, 18A301 (2014).
  * Jones (2015) R. O. Jones, Rev. Mod. Phys. 87, 897 (2015).
  * Mardirossian and Head-Gordon (2017) N. Mardirossian and M. Head-Gordon, Mol. Phys. 115, 2315 (2017).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Levy (1979) M. Levy, Proc. Natl. Acad. Sci. U.S.A. 76, 6062 (1979).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Perdew and Schmidt (2001) J. P. Perdew and K. Schmidt, in  _AIP Conference Proceedings_ , Vol. 577 (AIP, 2001) pp. 1–20.
  * Goerigk  _et al._ (2017) L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi, and S. Grimme, Phys. Chem. Chem. Phys. 19, 32184 (2017).
  * Mardirossian and Head-Gordon (2018) N. Mardirossian and M. Head-Gordon, J. Chem. Phys. 148, 241736 (2018).
  * Najibi and Goerigk (2018) A. Najibi and L. Goerigk, J. Chem. Theory Comput. 14, 5725 (2018).
  * Perdew  _et al._ (1982) J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz Jr, Phys. Rev. Lett. 49, 1691 (1982).
  * Mori-Sánchez, Cohen, and Yang (2006) P. Mori-Sánchez, A. J. Cohen, and W. Yang, J. Chem. Phys. 125, 201102 (2006).
  * Hait and Head-Gordon (2018a) D. Hait and M. Head-Gordon, J. Phys. Chem. Lett. 9, 6280 (2018a).
  * Lonsdale and Goerigk (2020) D. R. Lonsdale and L. Goerigk, Phys. Chem. Chem. Phys. 22, 15805 (2020).
  * Cohen, Mori-Sánchez, and Yang (2008a) A. Cohen, P. Mori-Sánchez, and W. Yang, J. Chem. Phys. 129, 121104 (2008a).
  * Cohen, Mori-Sánchez, and Yang (2008b) A. J. Cohen, P. Mori-Sánchez, and W. Yang, Science 321, 792 (2008b).
  * Hait, Rettig, and Head-Gordon (2019a) D. Hait, A. Rettig, and M. Head-Gordon, J. Chem. Phys. 150, 094115 (2019a).
  * Dohm  _et al._ (2018) S. Dohm, A. Hansen, M. Steinmetz, Grimme, Stefan, and M. P. Checinski, J. Chem. Theor. Comput. 14, 2596 (2018).
  * Chan, Gill, and Kimura (2019) B. Chan, P. M. W. Gill, and M. Kimura, J. Chem. Theor. Comput. 15, 3610 (2019).
  * Hait and Head-Gordon (2018b) D. Hait and M. Head-Gordon, J. Chem. Theory Comput. 14, 1969 (2018b).
  * Rettig  _et al._ (2020) A. Rettig, D. Hait, L. W. Bertels, and M. Head-Gordon, J. Chem. Theory Comput. 16, 7473 (2020).
  * Medvedev  _et al._ (2017a) M. G. Medvedev, I. S. Bushmarinov, J. Sun, J. P. Perdew, and K. A. Lyssenko, Science 355, 49 (2017a).
  * Hammes-Schiffer (2017) S. Hammes-Schiffer, Science 355, 28 (2017).
  * Kepp (2017) K. P. Kepp, Science 356, 496 (2017).
  * Medvedev  _et al._ (2017b) M. G. Medvedev, I. S. Bushmarinov, J. Sun, J. P. Perdew, and K. A. Lyssenko, Science 356, 496 (2017b).
  * Graziano (2017) G. Graziano, Nat. Rev. Chem. 1, 0019 (2017).
  * Korth (2017) M. Korth, Angew. Chem. Int. Ed. 56, 5396 (2017).
  * Brorsen  _et al._ (2017) K. R. Brorsen, Y. Yang, M. V. Pak, and S. Hammes-Schiffer, J. Phys. Chem. Lett. 8, 2076 (2017).
  * Wang  _et al._ (2017a) Y. Wang, X. Wang, D. G. Truhlar, and X. He, J. Chem. Theory Comput. 13, 6068 (2017a).
  * Ranasinghe, Perera, and Bartlett (2017) D. S. Ranasinghe, A. Perera, and R. J. Bartlett, J. Chem. Phys. 147, 204103 (2017).
  * Mayer  _et al._ (2017) I. Mayer, I. Pápai, I. Bakó, and A. Nagy, J. Chem. Theory Comput. 13, 3961 (2017).
  * Gould (2017) T. Gould, J. Chem. Theory Comput. 13, 2373 (2017).
  * Leininger  _et al._ (2000) M. L. Leininger, W. D. Allen, H. F. Schaefer III, and C. D. Sherrill, J. Chem. Phys. 112, 9213 (2000).
  * Hait  _et al._ (2019) D. Hait, N. M. Tubman, D. S. Levine, K. B. Whaley, and M. Head-Gordon, J. Chem. Theory Comput. 15, 5370 (2019).
  * Hickey and Rowley (2014) A. L. Hickey and C. N. Rowley, J. Phys. Chem. A 118, 3678 (2014).
  * Verma and Truhlar (2017) P. Verma and D. G. Truhlar, Phys. Chem. Chem. Phys. 19, 12898 (2017).
  * Johnson  _et al._ (2019) A. I. Johnson, K. P. Withanage, K. Sharkas, Y. Yamamoto, T. Baruah, R. R. Zope, J. E. Peralta, and K. A. Jackson, J. Chem. Phys. 151, 174106 (2019).
  * Grotjahn  _et al._ (2020) R. Grotjahn, G. J. Lauter, M. Haasler, and M. Kaupp, J. Phys. Chem. A 124, 8346 (2020).
  * Hait and Head-Gordon (2018c) D. Hait and M. Head-Gordon, J. Chem. Phys. 148, 171102 (2018c).
  * Thakkar and Wu (2015) A. J. Thakkar and T. Wu, J. Chem. Phys. 143, 144302 (2015).
  * Hait and Head-Gordon (2018d) D. Hait and M. Head-Gordon, Phys. Chem. Chem. Phys. 20, 19800 (2018d).
  * Withanage  _et al._ (2019) K. P. Withanage, S. Akter, C. Shahi, R. P. Joshi, C. Diaz, Y. Yamamoto, R. Zope, T. Baruah, J. P. Perdew, J. E. Peralta, _et al._ , Phys. Rev. A 100, 012505 (2019).
  * Hait, Rettig, and Head-Gordon (2019b) D. Hait, A. Rettig, and M. Head-Gordon, Phys. Chem. Chem. Phys. 21, 21761 (2019b).
  * Cohen and Tantirungrotechai (1999) A. Cohen and Y. Tantirungrotechai, Chem. Phys. Lett. 299, 465 (1999).
  * De Proft, Tielens, and Geerlings (2000) F. De Proft, F. Tielens, and P. Geerlings, J. Mol. Struct. Theochem 506, 1 (2000).
  * Hofinger and Wendland (2002) S. Hofinger and M. Wendland, Int. J. Quantum Chem. 86, 199 (2002).
  * Matta (2010) C. F. Matta, J. Comput. Chem. 31, 1297 (2010).
  * Halkier  _et al._ (1997) A. Halkier, H. Koch, O. Christiansen, P. Jorgensen, and T. Helgaker, J. Chem. Phys. 107, 849 (1997).
  * Zhang, Calvin, and Valeev (2016) J. Zhang, J. A. Calvin, and E. F. Valeev, Mol. Phys. 114, 2894 (2016).
  * Bokhan, Trubnikov, and Bartlett (2016) D. Bokhan, D. N. Trubnikov, and R. J. Bartlett, J. Chem. Phys. 144, 234107 (2016).
  * Hollett, Kelly, and Poirier (2006) J. W. Hollett, A. Kelly, and R. A. Poirier, J. Phys. Chem. A 110, 13884 (2006).
  * Raghavachari  _et al._ (1989) K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, Chem. Phys. Lett. 157, 479 (1989).
  * Lee  _et al._ (2019) J. Lee, L. W. Bertels, D. W. Small, and M. Head-Gordon, Phys. Rev. Lett. 123, 113001 (2019).
  * Brakestad  _et al._ (2020) A. Brakestad, S. R. Jensen, P. Wind, M. D’Alessandro, L. Genovese, K. H. Hopmann, and L. Frediani, J. Chem. Theory Comput. 16, 4874 (2020).
  * Dutoi and Head-Gordon (2006) A. D. Dutoi and M. Head-Gordon, Chem. Phys. Lett. 422, 230 (2006).
  * Ruzsinszky  _et al._ (2006) A. Ruzsinszky, J. Perdew, G. Csonka, O. Vydrov, and G. Scuseria, J. Chem. Phys. 125, 194112 (2006).
  * Vydrov, Scuseria, and Perdew (2007) O. A. Vydrov, G. E. Scuseria, and J. P. Perdew, J. Chem. Phys. 126, 154109 (2007).
  * Shao  _et al._ (2015) Y. Shao, Z. Gan, E. Epifanovsky, A. T. B. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P. R. Horn, L. D. Jacobson, I. Kaliman, R. Z. Khaliullin, T. Kús, A. Landau, J. Liu, E. I. Proynov, Y. M. Rhee, R. M. Richard, M. A. Rohrdanz, R. P. Steele, E. J. Sundstrom, H. L. Woodcock III, P. M. Zimmerman, D. Zuev, B. Albrecht, E. Alguire, B. Austin, G. J. O. Beran, Y. A. Bernard, E. Berquist, K. Brandhorst, K. B. Bravaya, S. T. Brown, D. Casanova, C.-M. Chang, Y. Chen, S. H. Chien, K. D. Closser, D. L. Crittenden, M. Diedenhofen, R. A. DiStasio Jr., H. Dop, A. D. Dutoi, R. G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M. W. D. Hanson-Heine, P. H. P. Harbach, A. W. Hauser, E. G. Hohenstein, Z. C. Holden, T.-C. Jagau, H. Ji, B. Kaduk, K. Khistyaev, J. Kim, J. Kim, R. A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C. M. Krauter, K. U. Lao, A. Laurent, K. V. Lawler, S. V. Levchenko, C. Y. Lin, F. Liu, E. Livshits, R. C. Lochan, A. Luenser, P. Manohar, S. F. Manzer, S.-P. Mao, N. Mardirossian, A. V. Marenich, S. A. Maurer, N. J. Mayhall, C. M. Oana, R. Olivares-Amaya, D. P. O’Neill, J. A. Parkhill, T. M. Perrine, R. Peverati, P. A. Pieniazek, A. Prociuk, D. R. Rehn, E. Rosta, N. J. Russ, N. Sergueev, S. M. Sharada, S. Sharmaa, D. W. Small, A. Sodt, T. Stein, D. Stück, Y.-C. Su, A. J. W. Thom, T. Tsuchimochi, L. Vogt, O. Vydrov, T. Wang, M. A. Watson, J. Wenzel, A. White, C. F. Williams, V. Vanovschi, S. Yeganeh, S. R. Yost, Z.-Q. You, I. Y. Zhang, X. Zhang, Y. Zhou, B. R. Brooks, G. K. L. Chan, D. M. Chipman, C. J. Cramer, W. A. Goddard III, M. S. Gordon, W. J. Hehre, A. Klamt, H. F. Schaefer III, M. W. Schmidt, C. D. Sherrill, D. G. Truhlar, A. Warshel, X. Xua, A. Aspuru-Guzik, R. Baer, A. T. Bell, N. A. Besley, J.-D. Chai, A. Dreuw, B. D. Dunietz, T. R. Furlani, S. R. Gwaltney, C.-P. Hsu, Y. Jung, J. Kong, D. S. Lambrecht, W. Liang, C. Ochsenfeld, V. A. Rassolov, L. V. Slipchenko, J. E. Subotnik, T. Van Voorhis, J. M. Herbert, A. I. Krylov, P. M. W. Gill, and M. Head-Gordon, Mol. Phys. 113, 184 (2015).
  * Dunning Jr (1989) T. H. Dunning Jr, J. Chem. Phys. 90, 1007 (1989).
  * Woon and Dunning Jr (1995) D. E. Woon and T. H. Dunning Jr, J. Chem. Phys. 103, 4572 (1995).
  * Peterson and Dunning Jr (2002) K. A. Peterson and T. H. Dunning Jr, J. Chem. Phys. 117, 10548 (2002).
  * Prascher  _et al._ (2011) B. P. Prascher, D. E. Woon, K. A. Peterson, T. H. Dunning, and A. K. Wilson, Theor. Chem. Acc. 128, 69 (2011).
  * Jensen (2005) F. Jensen, Theor. Chem. Acc. 113, 267 (2005).
  * Karton and Martin (2006) A. Karton and J. M. Martin, Theor. Chem. Acc. 115, 330 (2006).
  * Jensen (2001) F. Jensen, J. Chem. Phys. 115, 9113 (2001).
  * Jensen (2002a) F. Jensen, J. Chem. Phys. 116, 7372 (2002a).
  * Jensen (2002b) F. Jensen, J. Chem. Phys. 117, 9234 (2002b).
  * Jensen and Helgaker (2004) F. Jensen and T. Helgaker, J. Chem. Phys. 121, 3463 (2004).
  * Jensen (2007) F. Jensen, J. Phys. Chem. A 111, 11198 (2007).
  * Jensen (2010) F. Jensen, J. Chem. Theory Comput. 6, 2726 (2010).
  * Halkier  _et al._ (1999) A. Halkier, W. Klopper, T. Helgaker, and P. Jørgensen, J. Chem. Phys. 111, 4424 (1999).
  * Gill, Johnson, and Pople (1993) P. M. Gill, B. G. Johnson, and J. A. Pople, Chem. Phys. Lett. 209, 506 (1993).
  * Kurlancheek and Head-Gordon (2009) W. Kurlancheek and M. Head-Gordon, Mol. Phys. 107, 1223 (2009).
  * Hui and Chai (2016) K. Hui and J.-D. Chai, J. Chem. Phys. 144, 044114 (2016).
  * Wilson, Bradley, and Tozer (2001) P. J. Wilson, T. J. Bradley, and D. J. Tozer, J. Chem. Phys. 115, 9233 (2001).
  * Staroverov  _et al._ (2003) V. N. Staroverov, G. E. Scuseria, J. Tao, and J. P. Perdew, J. Chem. Phys. 119, 12129 (2003).
  * Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
  * Becke (1997) A. D. Becke, J. Chem. Phys. 107, 8554 (1997).
  * Krukau  _et al._ (2006) A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, J. Chem. Phys. 125, 224106 (2006).
  * Dirac (1931) P. A. Dirac, in  _Proc. R. Soc. A_ , Vol. 133 (The Royal Society, 1931) pp. 60–72.
  * Perdew and Wang (1992) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
  * Peverati and Truhlar (2011a) R. Peverati and D. G. Truhlar, J. Chem. Phys. 135, 191102 (2011a).
  * Rohrdanz, Martins, and Herbert (2009) M. A. Rohrdanz, K. M. Martins, and J. M. Herbert, J. Chem. Phys. 130, 054112 (2009).
  * Zhao and Truhlar (2005) Y. Zhao and D. G. Truhlar, J. Phys. Chem. A 109, 5656 (2005).
  * Becke (1988) A. D. Becke, Phys. Rev. A 38, 3098 (1988).
  * Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Lin  _et al._ (2012) Y.-S. Lin, C.-W. Tsai, G.-D. Li, and J.-D. Chai, J. Chem. Phys. 136, 154109 (2012).
  * Adamo and Barone (1998) C. Adamo and V. Barone, J. Chem. Phys. 108, 664 (1998).
  * Zhao and Truhlar (2008a) Y. Zhao and D. G. Truhlar, Theor. Chem. Acc. 120, 215 (2008a).
  * Grimme (2006a) S. Grimme, J. Comput. Chem. 27, 1787 (2006a).
  * Mardirossian and Head-Gordon (2014) N. Mardirossian and M. Head-Gordon, Phys. Chem. Chem. Phys. 16, 9904 (2014).
  * Chai and Head-Gordon (2008) J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys. 10, 6615 (2008).
  * Peverati and Truhlar (2012) R. Peverati and D. G. Truhlar, J. Chem. Theory Comput. 8, 2310 (2012).
  * Lee, Yang, and Parr (1988) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
  * Peverati, Zhao, and Truhlar (2011) R. Peverati, Y. Zhao, and D. G. Truhlar, J. Phys. Chem. Lett. 2, 1991 (2011).
  * Peverati and Truhlar (2011b) R. Peverati and D. G. Truhlar, J. Phys. Chem. Lett. 2, 2810 (2011b).
  * Becke (1993) A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
  * Sun, Ruzsinszky, and Perdew (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
  * Yanai, Tew, and Handy (2004) T. Yanai, D. P. Tew, and N. C. Handy, Chem. Phys. Lett. 393, 51 (2004).
  * Sun  _et al._ (2013) J. Sun, R. Haunschild, B. Xiao, I. W. Bulik, G. E. Scuseria, and J. P. Perdew, J. Chem. Phys. 138, 044113 (2013).
  * Zhao and Truhlar (2008b) Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput. 4, 1849 (2008b).
  * Tao  _et al._ (2003) J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys. Rev. Lett. 91, 146401 (2003).
  * Mardirossian and Head-Gordon (2016) N. Mardirossian and M. Head-Gordon, J. Chem. Phys. 144, 214110 (2016).
  * Wellendorff  _et al._ (2014) J. Wellendorff, K. T. Lundgaard, K. W. Jacobsen, and T. Bligaard, J. Chem. Phys. 140, 144107 (2014).
  * Haoyu  _et al._ (2016) S. Y. Haoyu, X. He, S. L. Li, and D. G. Truhlar, Chem. Sci. 7, 5032 (2016).
  * Zhao and Truhlar (2006) Y. Zhao and D. G. Truhlar, J. Chem. Phys. 125, 194101 (2006).
  * Wang  _et al._ (2017b) Y. Wang, X. Jin, S. Y. Haoyu, D. G. Truhlar, and X. He, Proc. Natl. Acad. Sci. U.S.A. 114, 8487 (2017b).
  * Kozuch and Martin (2013) S. Kozuch and J. M. Martin, J. Comput. Chem. 34, 2327 (2013).
  * Mardirossian and Head-Gordon (2015) N. Mardirossian and M. Head-Gordon, J. Chem. Phys. 142, 074111 (2015).
  * Zhang  _et al._ (2011) I. Y. Zhang, X. Xu, Y. Jung, and W. A. Goddard, Proc. Natl. Acad. Sci. U.S.A. 108, 19896 (2011).
  * Peverati and Truhlar (2011c) R. Peverati and D. G. Truhlar, J. Phys. Chem. Lett. 3, 117 (2011c).
  * Goerigk and Grimme (2010) L. Goerigk and S. Grimme, J. Chem. Theory Comput. 7, 291 (2010).
  * Yu, He, and Truhlar (2016) H. S. Yu, X. He, and D. G. Truhlar, J. Chem. Theory Comput. 12, 1280 (2016).
  * Zhang, Xu, and Goddard (2009) Y. Zhang, X. Xu, and W. A. Goddard, Proc. Natl. Acad. Sci. U.S.A. 106, 4963 (2009).
  * Karton  _et al._ (2008) A. Karton, A. Tarnopolsky, J.-F. Lamére, G. C. Schatz, and J. M. Martin, J. Phys. Chem. A 112, 12868 (2008).
  * Chai and Head-Gordon (2009) J.-D. Chai and M. Head-Gordon, J. Chem. Phys. 131, 174105 (2009).
  * Grimme (2006b) S. Grimme, J. Chem. Phys. 124, 034108 (2006b).
  * Li and Yang (2017) C. Li and W. Yang, J. Chem. Phys. 146, 074107 (2017).
  * Tubman  _et al._ (2016) N. M. Tubman, J. Lee, T. Y. Takeshita, M. Head-Gordon, and K. B. Whaley, J. Chem. Phys. 145, 044112 (2016).
  * Tubman  _et al._ (2020) N. M. Tubman, C. D. Freeman, D. S. Levine, D. Hait, M. Head-Gordon, and K. B. Whaley, J. Chem. Theory. Comput. 16, 2139 (2020).
  * Gaiduk, Mizzi, and Staroverov (2012) A. P. Gaiduk, D. Mizzi, and V. N. Staroverov, Phys. Rev. A 86, 052518 (2012).
  * Gritsenko, Mentel, and Baerends (2016) O. Gritsenko, Ł. Mentel, and E. Baerends, J. Chem. Phys. 144, 204114 (2016).
  * Head-Gordon (2003) M. Head-Gordon, Chem. Phys. Lett. 372, 508 (2003).
  * Coulson and Fischer (1949) C. A. Coulson and I. Fischer, Philos. Mag. 40, 386 (1949).