See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/265649680

Simultaneous Visualization of Covalent and Non-Covalent Interactions Using

Regions of Density Overlap

Article in Journal of Chemical Theory and Computation · September 2014

DOI: 10.1021/ct500490b · Source: PubMed

CITATIONS READS

196 719

2 authors, including:

Piotr de Silva
Technical University of Denmark
80 PUBLICATIONS 1,303 CITATIONS

SEE PROFILE

All content following this page was uploaded by Piotr de Silva on 01 December 2017.

The user has requested enhancement of the downloaded file.

 Article

pubs.acs.org/JCTC
Terms of Use

Simultaneous Visualization of Covalent and Noncovalent

Interactions Using Regions of Density Overlap
Piotr de Silva* and Clémence Corminboeuf*
Laboratory for Computational Molecular Design, Institut des Sciences et Ingénierie Chimiques, Ecole Polytechnique Fédérale de
Lausanne, CH-1015 Lausanne, Switzerland
*
S Supporting Information

ABSTRACT: We introduce a density-dependent bonding

descriptor that enables simultaneous visualization of both
covalent and noncovalent interactions. The proposed quantity
is tailored to reveal the regions of space, where the total
electron density results from a strong overlap of shell, atomic,
or molecular densities. We show that this approach is
successful in describing a variety of bonding patterns as well
as nonbonding contacts. The Density Overlap Regions
Indicator (DORI) analysis is also exploited to visualize and
quantify the concept of electronic compactness in supramolecular chemistry. In particular, the scalar field is used to compare the
compactness in molecular crystals, with a special emphasis on quaterthiophene derivatives with enhanced charge mobilities.

1. INTRODUCTION well-defined energies. Also, depending on the method, LMOs

Visualization of bonding interactions between atoms and can be a combination of CMOs of different symmetry.
The Lewis picture of chemical bonding can be restored from
molecules is a long-standing quest in computational chemistry.
a one-particle density matrix in atomic orbitals basis by Natural
During its development many methods have been proposed to
Bond Orbitals (NBO) analysis.15,16 NBOs can be seen as an
serve this purpose. The main interest lies in creating a tool that
extension of LMOs, where the constraint on occupation
enables not only to see the interaction but also to interpret its
numbers has been released in favor of stricter localization.
character and properties. There is no general agreement on Therefore, NBOs yield maximum occupancy one- and two-
how to derive an optimal method and the existing ones are center orbitals. The bonding patterns beyond the two-center
based on many very different ideas based on orbital paradigm, can be studied with the Adaptive Natural Density
transformations, localization descriptors, topological density Partitioning (AdNDP)17,18 method, which is a generalization of
analysis, and others. This problem can be traced back to the the NBO concepts.
lack of a clear and unambiguous definition of a bond in The interpretation of chemical bonding in terms of molecular
quantum mechanics. Therefore, a chemical bond together with orbitals remains nonetheless ambiguous. For an N electron
other notions such as electron shells, lone pairs, aromaticity, closed-shell system, there are N/2 occupied orbitals, which
atomic charges, (hyper-) conjugation, strain, etc. constitute a extend over the entire space or a restricted region of space,
rich set of “fuzzy”, yet invaluably useful concepts.1−4 which means that orbital densities overlap with each other. As a
The fundamental model of chemical bonding is based on result, the do not offer an intuitive depiction and detailed
one-determinantal electronic structure methods such as information on the nature and location of electron pairs. A
Hartree−Fock or Kohn−Sham density functional theory different class of bonding analysis methods is based on scalar
(DFT). The fundamental ingredients of these methods, namely fields, which detect localized electrons in real space. These
canonical molecular orbitals (CMO), represent bonds, lone localization functions are usually computed from molecular
pairs, and core electrons. The analysis of bonding effects is orbitals or density matrices, such as Electron Localization
based on their symmetry as well as corresponding orbital Function (ELF),19−21 Localized Orbital Locator (LOL),22,23
energies. An undesired feature of CMOs is that they are Electron Localizability Indicator (ELI),24−26 or nonadditive
delocalized over large parts of a molecule, which is not Fisher information.27,28 Methods relying only on electron
compatible with a typically two-center character of an individual density have also been proposed recently, such as Localized
bond, that is, the Lewis picture. This character is restored by Electrons Detector (LED)29,30 or Single-Exponential Decay
localized molecular orbitals (LMO), which are obtained by Detector (SEDD).31,32 The starting point and rationale for
some unitary transformation of occupied CMOs. A number of different bonding detectors is based on different physical
localization procedures has been proposed, which differ by the
criterion used for the transformation.5−14 As LMOs are not Received: June 6, 2014
eigenfunction of a one-electron Hamiltonian, they do not have Published: June 30, 2014

© 2014 American Chemical Society 3745 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756
Journal of Chemical Theory and Computation Article

assumptions, nevertheless, usually they provide some local illustrations of covalent bonds of varying nature, π-electron
measure of Pauli repulsion, which is related directly to the delocalization and intermolecular interactions. We also use this
kinetic energy of electrons. approach to visualize and quantify electronic implications of
Another prominent method for bonding analysis is the compactness in organic molecular crystals. In particular, we
quantum theory of atoms in molecules (QTAIM).33,34 Within analyze the relation between so understood compactness and
this theory, the total electron density is partitioned into charge mobilities in a molecular crystal relevant to the field of
nonoverlapping atomic basins. The whole analysis is based on organic electronics.
topological properties of electron density, such as the character
of stationary points and existence of bonding paths defined in 2. THEORY
terms of the density gradient. The Laplacian of the density is Considering solid foundations and the remarkable success of
used as an indicator of local charge concentration,35 which density functional theory, today there is very little doubt that
yields also information about bonding in real space. Additional electron density contains all the useful information about
information about the character of the bond can be extracted electronic structure of molecular systems. However, there are
from various characteristics calculated at the bond critical point, not many well established methods to visualize bonding
such as the bond ellipticity36 or metallicity.37−39 patterns in molecules, that are based solely on electron density.
All the above-mentioned approaches have gained much The recently introduced single-exponential decay detec-
popularity and were applied to many computational studies of tor31,32,50 succeeds in revealing atomic shells, lone electron
molecular systems. Additionally, QTAIM is widely used in pairs and core electrons. Moreover, SEDD-based atomic shell
analysis of experimental charge densities,40−42 whereas populations50 stay in close quantitative agreement with the
application of orbital-based localization functions requires Aufbau principle, which fixes the integer number of electrons in
further approximations.43,44 QTAIM, in principle, allows to atomic shells. Initially, SEDD was introduced as an arbitrary
characterize both covalent and noncovalent interactions; dimensionless function of electron density, meant to discover
however, their representation in terms of stationary points single-exponential regions of the density,
and bond paths is not very intuitive. On the other hand,
localization functions provide a meaningful representation of ⎡ ⎛ 2 ⎞2 ⎤
⎢ ∇ ρ(r)
⎜ ∇ ρ(r) ( ) ⎟⎟ ⎥⎥
atomic shells, lone pairs, and covalent bonds, but usually, they ⎢
do not detect noncovalent interactions. This gap has been filled SEDD(r) = ln 1 + ⎜
⎢ ⎜ ρ(r) ⎟⎥
recently by the noncovalent interaction (NCI) index,45,46 which ⎢⎣ ⎝ ⎠ ⎦⎥
is based on the reduced density gradient (RDG) and uses (1)
concepts from QTAIM to distinguish the character of where a square of a vector refers to a scalar product. This
interactions. The dependence of NCI on electron density convention is used throughout this work and the explicit form
only enables a straightforward analysis of experimental charge of this equation and the following ones in terms of density
densities.47 NCI visualizes both intra- and intermolecular weak derivatives is given in the Supporting Information. The idea
interactions through RDG isosurfaces at low electron density behind SEDD was based on an Ansatz that, for localized
values, however, does not reveal covalent bonds. The interest in electrons, the density can be described locally by a single
studying simultaneously strong and weak interactions has led to orbital, which decays approximately exponentially from an
applications of combined ELF/NCI analysis.48,49 arbitrary point. Later, a physical interpretation was attributed to
In this work, we introduce a scalar field, which reveals both it in terms of the local wave vector k(r) = ((∇ρ(r))/(ρ(r)) and
covalent and noncovalent interactions in the same value range. the homogeneous electron gas (HEG) reference system.50 The
It should be seen as a modification of the Single-Exponential properties of SEDD are such that it has low values within
Decay Detector,31,32,50 which was proposed by one of us as a localized electron shells, bonds, and lone pairs and goes to
density-based bonding descriptor. Although, in principle, infinity far from the molecule. For noncovalent intermolecular
SEDD reveals both covalent and noncovalent interactions, interactions, a decrease in SEDD values was observed in the
the latter ones are often not well resolved due to the significant region between two molecules;31 so in principle, it discerns
numerical noise in the low density regions, stemming from the strong at weak interactions simultaneously. Unfortunately,
use of finite atomic basis sets. The approach presented here is SEDD suffers from relatively strong basis set dependence
free from this flaw and additionally allows for a convenient resulting in numerical noise in low density regions. Therefore,
transformation to the range of values restricted between 0 and analysis of noncovalent interactions is difficult in practice. We
1. Contrary to the majority of bonding descriptors, the propose a new descriptor, which is also based on the idea of
introduced quantity does not directly measure where do the detecting single-exponential parts of the density, however,
electrons locate but rather focuses on geometrical features of exhibits a more regular behavior and enables a simultaneous
the electron density. More specifically, one here reveals regions analysis of covalent and noncovalent interactions.
of space where the electron density between atoms, molecules For exponential densities, ρ(r) ∼ e−λ|r−r0|, the term ξ(r) =
or atomic shells clashes. As detailed in the next section, these ∇((∇ρ(r)/ρ(r))2/ρ(r))2 in eq 1 is proportional to η(r) =
overlapping regions can be identified by the deviation of their (∇k2(r)/(k3F))2, where k(r) is the local wave vector and kF =
densities from that of localized electrons, which decay (3π2ρ)1/3 is the Fermi wave vector in the HEG model.
approximately exponentially from an arbitrary point. A physical Therefore, SEDD has low values where the gradient of the
interpretation of the proposed descriptor has been provided in squared local wavenumber is small compared to the volume of
terms of the local wave vector51−53 (i.e., ∇ρ(r)/ρ(r)), which is the Fermi sphere. The square of the local wavenumber can be
constant for single-exponential densities and thus reflective of intuitively related to the local kinetic energy per particle;
the shape of the density. however, the latter is not a well-defined quantity.54−57 The
We apply the introduced scalar field to several illustrative reference of HEG is employed to obtain a dimensionless
studies of bonding effects in molecular systems. They include quantity, nevertheless, it is not the only possibility. Another
3746 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756
Journal of Chemical Theory and Computation Article

choice is to use self-reference and divide ▽k2(r) by the proper

power of k2(r) itself. This leads to the following dimensionless
quantity
⎛ 2 ⎞2

(∇k 2(r))2
⎜∇
⎝ ( ) ⎟⎠
∇ ρ(r)
ρ(r)
θ(r) = =
(k 2(r))3 6
( )
∇ ρ(r)
ρ(r) (2)

Both, ξ(r) and θ(r) are equal to 0 if ρ(r) is exactly single-

exponential; however, they have different properties for
approximately exponential densities, which can be found in
molecular systems. In vicinity of nuclei, the density behaves
similar to e−2Zα|r−rα|, where Z is the nuclear charge and rα is the
position of a nucleus.58 Therefore, in these regions θ ≈ 0, as the
nominator in eq 2 is approximately 0 and the denominator is a
constant. The same holds far from the molecule, where the
single-exponential decay is dictated by the ionization potential
1/2
I: ρ(r) ∼ e−2(2I) |r|.59
In bonding regions, the density is characterized by small
gradients, in particular ▽ρ(r) = 0 at the bond critical point
(BCP). In this case, both, the nominator and the denominator
go to zero. Since the denominator decays to zero faster while
approaching the BCP, θ(r) → + ∞. This contrasting behavior
in the atomic core regions and at bonds shows that θ(r) cannot
be treated as a localization function. Instead, it is a
renormalization of the density, allowing to distinguish its
character from geometrical properties and indicating bonding
regions wherever θ(r) → + ∞.
Considering that θ(r) is unbound from above, it is not very
convenient from the point of view of its visualization.
Therefore, we make use of the standard mapping to the [0,1]
range
θ(r)
γ(r) =
1 + θ(r) (3)

Now, γ(r) is close to 1 in bonding regions and close to 0 at

nuclei and far from the molecule. The asymptotic behavior of
γ(r) at nuclear positions, BCPs and infinity was deduced from
its analytic form. Nevertheless, the performance in detecting
bonding regions of realistic systems has to be verified
numerically.
The role of γ is to probe geometrical features of electron
density by measuring the deviations from single-exponentiality.
This idea is based on the fact that atomic densities are Figure 1. Density ρ(r) (blue) and γ(r) (purple) for the H2
approximately piecewise exponential,60−62 whereas in the promolecule at the internuclear separation of (a) 1.4 au, (b) 10 au,
interaction regions they are not, due to the overlap of two or (c) 10 au (close-up on low density values). Densities of two
more atomic densities. The same is true for interactions constituting hydrogen atoms are in green.
between whole molecules as the molecular density tails also
decay exponentially. This is a different approach compared to
the majority of bonding descriptors, which aim at revealing by the density of a single hydrogen. On the other hand, within
where electrons are localized. Therefore, it is instructive to the bond, the curvature of ρ is dictated by the overlap of two
show how γ performs in practice for a model system tails and γ > 0. In the bond midpoint, where the overlap is the
representing an overlap of two simple densities. To this end, strongest, γ reaches its upper bound.
we take a density which is a sum of two displaced hydrogen An analogous picture, but for a 10 au separation, is given in
densities, i.e. an H2 promolecule. Figure 1b. Here, the two densities are very well separated and
Figure 1a shows such density ρ (blue line) for a separation of their overlap is negligible on the scale of the density itself.
1.4 au between hydrogen atoms, which is very close to the Nevertheless, γ gives a similar picture as for the strong overlap
equilibrium bond length of the real H2 molecule. The case. Figure 1c shows the same density but only up to a 1000
decomposition of ρ into atomic contributions is marked with times smaller value. Now it is evident that γ reveals strong
green lines. It is evident that γ (purple line) is 0 outside the atomic density overlap regions regardless of the magnitude of
bonding region, where the total density is strongly dominated the density, which means that strong and weak interactions are
3747 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756
Journal of Chemical Theory and Computation Article

treated on an equal footing. It has to be noted, that the overlap

region between two densities is discovered just by analyzing the
shape of the total density, without any a priori knowledge of
how it was constructed.
Since the meaning of γ has been identified as a detector of
strong relative density overlap regions, we propose a name
Density Overlap Regions Indicator (DORI), which will replace
the γ symbol from now on.

3. COMPUTATIONAL DETAILS
Electronic structure of the investigated molecules was
calculated with the ADF63 package at the generalized Kohn−
Sham level of theory using B3LYP exchange-correlation
functional and the TZP basis set. Unless otherwise stated,
geometries were optimized at the same level. All presented
properties, including DORI, were calculated on a grid with a
locally modified version of DGrid program.64 ParaView65 was
used to visualize the results.
Patches to compute DORI with DGrid64 and NCIPLOT46
programs are available from the authors on request.

4. TYPICAL CHEMICAL INTERACTIONS

4.1. Intramolecular Interactions. The capability of DORI
to visualize typical bonding patterns is evaluated on a selection
of small molecular systems that are (a) H2, (b) O2, (c) N2, (d)
F2, (e) CO, (f) CO2, (g) LiH, (h) B2H6, (i) C2H6, (j) C2H4, (k)
C2H2, (l) C6H6. Figure 2 shows color-coded maps of DORI,
where blue color corresponds to DORI(r) ≈ 1 and the red
DORI(r) ≈ 0. The positions of nuclei and bonds are indicated
explicitly only when needed for the sake of clarity.
The hydrogen molecule (Figure 2a) is characterized by the
presence of a direct lenticular bond. This contrasts with the
picture given by bonding descriptors based on local Pauli Figure 2. DORI representation of typical bonds in a selection of
repulsion measures (e.g., ELF), which intrinsically cannot reveal molecular systems: (a) H2, (b) O2, (c) N2, (d) F2, (e) CO, (f) CO2,
any bonding for a closed-shell two-electron system. The double (g) LiH, (h) B2H6, (i) C2H6, (j) C2H4, (k) C2H2, (l) C6H6. Values in
the range from 0 (red) to 1 (blue).
bond in the triplet oxygen O2 (Figure 2b) is described by a
basin of DORI values close to 1, which is located around the
bond midpoint. Around nuclei, DORI is close to zero, side of the latter atom. Figure 2h shows the bonding pattern of
indicating the position of localized 1s electrons. Outside B2H6. From the valence perspective, diborane is an example of a
localized cores, the overlap region between core and valence three-center two-electron bond, where hydrogen acts as a
electrons results in high DORI values. Between these overlap bridging atom. From the perspective of DORI the bond results
regions and the bonding region DORI ≈ 0.5. In contrast with from a clash of atomic densities inside the BHBH four-member
ELF, the lone electron pairs on oxygen atoms, which do not ring. The interaction region is delocalized over the four atoms,
stem from the overlap of atomic densities, are not explicitly but still, one discerns a direct B−H bonding.
revealed by DORI. N2 (Figure 2c) and F2 (Figure 2d) shows The pictures of ethane, ethene, ethyne, and benzene are
essentially the same features as O2 but with larger (triple bond) given in Figures 2i−l. All these hydrocarbons display direct C−
and smaller (single bond) bonding regions, respectively. The C and C−H bonds. The shape of the single C−C bond is
bonding region of N2 merges with core/valence overlap lenticular, whereas multiple bonds become more cylindrical. π
domains, whereas in F2 DORI falls to 0 in between the atomic bonds are not explicitly discernible as atomic densities overlap
and bonding zones and exhibit a more discotic shape. This strongly along the bond axis. In the case of benzene (Figure 2l
distinct character of the DORI bonding domain in F2 is in line and 3a vide infra), the delocalized bonding pattern, manifested
with the fact that the F−F bond is one of the weakest of all in the lack of bond alternation, is nevertheless clearly visible
covalent bonds due to large electrostatic repulsion.66 along with a steric clash at the ring center. The latter feature,
Unlike their contrasting Lewis structures, the DORI pictures not captured by ELF, demonstrates DORI’s ability to reveal
of CO and CO2 (Figure 2e and f) are very similar. The bonds noncovalent interactions within molecule in addition to
merge with the carbon’s atomic region that is bigger than that covalent bonds. In this context, the π electronic structure of
of the oxygen atoms. This results from a larger spatial extent of C4H4 is certainly the most striking example. As emphasized by
core electrons density due to a smaller nuclear charge. In case Schleyer and Politzer67−69 the small 4π-electron annulene
of an ionic LiH (Figure 2g), the density on lithium is strongly should be regarded as a unique molecule rather than as the
polarized and forms a cavity inside the interaction region antiaromatic paradigm. Its uniqueness arises from its ring strain
between two atoms. The strong interaction region is present and from the presence of two strongly localized and dense π
not only between hydrogen and lithium, but also on the other bonds in close proximity. Consideration of isosurface of DORI
3748 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756
Journal of Chemical Theory and Computation Article

−
Figure 3. DORIπ = 0.7 isosurface for (a) benzene, (b) C4H4, (c) C4H2+ 2+ +
4 , (d) C4H4 in the geometry of C4H4, (e) C5H5 , (f) C5H5 , (g) 1,4-
benzoquinone.

calculated from the π electron density (Figure 3b) perfectly of their attractive or repulsive nature. However, the
illustrates the density clash arising from this proximity. Removal decomposition of ▽2ρ into the sum of Hessian eigenvalues
of two electrons releases the π localization constrain typical for ▽2ρ = λ1 + λ2 + λ3, (λ1 ≤ λ2 ≤ λ3) provides a much clearer
antiaromatic systems and restores the delocalized 4c−3e picture view. In particular, the second eigenvalue λ2 < 0 is known to
in D4h or even D2h C4H2+ 4 (Figure 3c and d, respectively). DORI identify bonding regions, while λ2 > 0 indicates nonbonding
provides a concrete evidence for the repulsive π−π interaction interactions. Along with its sign, the magnitude of the
(strong Pauli repulsion) between parallel double bonds. The interaction is estimated from the values of the density itself;
repulsion between the two localized double bonds is also therefore, we use sgn(λ2)ρ(r) as a complementary scalar field
apparent in C5H+5 (Figure 3e) or benzoquinone (Figure 3g) but (see ref 45 for more details). A valuable alternative for
to a lesser extent, due to broader angles or larger distance characterizing the interactions would be to use eigenvalues of
between them. Akin to C4H2+ 4 , the delocalization pattern is also the stress tensor70−74 instead of the electron density Hessian.
recovered in aromatic C5H−5 (Figure 3f). However, calculating the stress tensor is more computationally
4.2. Intermolecular Interactions. So far, DORI was demanding and requires knowledge of the one-particle reduced
shown to reveal myriad bonding patterns as well as steric density matrix. Since we aim at a solely density-based analysis,
clashes within molecules. However, the same tool can be used we do not exploit this possibility in the current work.
to untangle the nature of intermolecular interactions. Since The DORI capability of depicting intermolecular interactions
DORI identifies regions of overlapping density, it does not is demonstrated on a series of typical noncovalently bound
directly carry information on the strength of the interaction and dimers taken from the S22 set,75 which enable a direct
does not distinguish between attraction and repulsion. In line comparison with the NCI index (Supporting Information to ref
with the NCI index,45 this limitation is resolved by combining 45). Figure 4 shows DORI = 0.9 isosurfaces for the selected
the analysis of DORI with that of the electron density Laplacian dimers with color-coded values of sgn(λ2)ρ(r). As both NCI
(▽2ρ). The regions of noncovalent interactions are charac- and DORI probe the shape of the density, the character and
terized by positive values of the density Laplacian irrespective magnitude of the noncovalent interactions are captured in a
3749 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756
Journal of Chemical Theory and Computation Article

two-dimensional color-coded maps with superimposed electron

density isocontours for a water dimer at five intermolecular

distances taken from the S22 × 5 set.76

Figure 4. DORI = 0.9 isosurfaces for (a) water dimer, (b) ammonia
dimer, (c) methane dimer, (d) formic acid dimer, (e) π-stacked
benzene dimer, (f) T-shaped benzene dimer, (g) π-stacked adenine−
thymine, (h) hydrogen bonded adenine−thymine. Isosurfaces are
color-coded with sgn(λ2)ρ(r) in the range from −0.02 au (red) to 0.02
au (blue).

similar manner by their isosurfaces and sgn(λ2)ρ(r). It is

evident that DORI reveals all sorts of intermolecular
interactions going from hydrogen bonds, pure dispersion
interactions to steric clashes. However, what clearly distin-
guishes DORI is that covalent bonds and the core/valence
interface of the constituent atoms are also visible in the same
[0−1] range without imposing any arbitrary thresholds on the
electron density. All in all, DORI offers a coherent and
comprehensive description of all the chemically relevant
interactions present in complex molecular systems.
The interaction regions are detected based on the geo-
metrical deformation of the electron density, discovering where
densities of different entities clash in a molecular complex. As a
result, π-stacking interactions (e.g., Figure 4e) will typically
display larger intermolecular DORI domains than the more
localized H-bonds or edge-to-face interactions (e.g., Figure 4f Figure 5. DORI maps with superimposed electron density isocontours
and h). To better understand this size variation, it is instructive (white lines) for water dimer at (a) 90%, (b) 100%, (c) 120%, (d)
to determine how does DORI relate to the shape of the density 150%, (e) 200% of the equilibrium bond length. Isovalues for density
isocontours and how does the size of DORI domains depend contours are on logarithmic scale. Values in the range from 0 (red) to
on the distance between interacting species. Figure 5 shows 1 (blue).

3750 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756

Journal of Chemical Theory and Computation Article

The density isocontours can be divided into two categories: reveal information on electron density overlaps and is hence
those encompassing only one water molecule and those well suited to analyze the interacting region. In contrast to the
encompassing the whole dimer. The interacting region deformation densities, DORI is achieved without referring to
(shown in blue/green) coincides with the one where the promolecular densities.
density isosurfaces of individual molecules collide and merge Figure 6 shows color-coded maps for [1.1.1]-, [2.1.1]-,
into a single isosurface. This contact region is characterized by a [2.2.1]-, and [2.2.2]-propellanes. Each cut-plane contains one
stronger curvature of isocontours, which is a result of the
overlap of densities attributed to distinct molecules. At 90% of
the water dimer equilibrium hydrogen bond length (Figure 5a),
DORI identifies the region where the density isocontours
distort due to the interactions.
At larger bond distances (Figures 5b−e), the merging
between the water molecules occurs with contours of lower
isovalues, which results in larger contact area. This means that
for a given system, the volume of a DORI domain increases
with decreasing the interaction strength. From another
perspective, DORI identifies where the molecular densities
merge into a strongly overlapping supermolecular density
irrespective of how large this density is. One should remind,
however, that the magnitude and sign of the interactions can be
brought by the analysis of sgn(λ2)ρ(r). For the largest distance
(Figure 5e), the interaction region splits into two domains due
to a numerical artifact, stemming from the finite atomic orbital
basis set expansion, which is know to have a nonphysical
behavior at density tails.77 Figure 6. DORI maps for (a) [1.1.1]-propellane, (b) [2.1.1]-
propellane (three-carbon ring plane), (c) [2.1.1]-propellane (four-
5. ILLUSTRATIVE APPLICATIONS carbon ring plane), (d) [2.2.1]-propellane (three-carbon ring plane),
(e) [2.2.1]-propellane (four-carbon ring plane), (f) [2.2.2]-propellane.
5.1. Propellanes. DORI is evidently capable of probing Values in the range from 0 (red) to 1 (blue).
covalent bonds and noncovalent interactions, but how does it
perform for cases where the existence of a bond is ambiguous?
Propellanes are good examples of such situations. The character of the three rings of a given propellane. The planes for both the
of the so-called inverted C−C bond between bridgehead three- and four-membered rings are represented for [2.1.1]-
carbon atoms has been investigated with various method- and [2.2.1]-propellanes. In [1.1.1]-propellane, the interaction
ologies,27,78−86 including analysis of experimental densities.87,88 between bridgehead carbons is clearly different from the other
The general agreement is that this bridgehead bond exhibits covalent C−C bonds (Figure 6a). The bonding region does not
some degree of covalency in small propellanes although the connect the two carbons directly but merges with the two other
pure covalent bond is only achieved in [2.2.2]-propellane. The bonds through the ring center. The feature occurs for all the
conceptual explanations for the peculiar nature of the bond in three-membered rings as well as for the four-membered rings of
the smallest [1.1.1]-propellane are specific to the employed [2.1.1]-propellane (Figure 6b,c). In contrast, the bridgehead
method. For instance, the valence bond (VB) analysis85,86 uses C−C bonds of the four-membered rings in the larger [2.2.1]
the “charge-shift” bond terminology due to a dominant energy and [2,2,2] polycyclic analogues are separated from the DORI
contribution coming from the covalent-ionic resonance domain at the ring center akin to rest of the bonds.
structures. On the other hand, the analysis of the deformation The DORI isosurfaces given in Figure 7 help connecting the
density and entropy displacement83 leads to an interpretation local character of the central bond with the overall bonding
in terms of “through-bridge” interaction, whereas “through- pattern of the same propellane series. They reveal a bonding
space” interactions are considered as more important in larger domain delocalized among the three rings in the smallest
propellanes. This incoherent description and terminology system, which contrasts with [2.2.2]-propellane that shows 2c-
results in an apparent controversy in the interpretation of the 2e covalent bonds only. [2.1.1]- and [2.2.1]-propellanes
nature of the central bond. represent intermediate situations in which the 2c−2e bonding
The unusual character of the inverted bond in [1.1.1]- pattern is gradually strengthened through the substitution of a
propellane can be traced back to the unique properties of three- by a four-membered ring. The consideration of the
cyclopropane. The smallest cycloalkane has a Baeyer strain far combined DORI-Laplacian analysis (Figure 8) further amplifies
lower than expected based simply on angle deformation.89 This and clarifies this contrast between the large and small polycyclic
unexpected stability was originally attributed to σ-aromaticity,90 systems: the interaction between bridgehead atoms in [1.1.1]-
although most recent interpretations re-emphasize alternative and [2.1.1]-propellanes (Figures 8a, b) is noncovalent (blue,
electronic effects such as rehybridization and strong geminal ▽2ρ(r) > 0), whereas, the red ▽2ρ(r) < 0 zone in [2.2.2]-
hyperconjugation (see refs 16, 91, and 92). While the strain propellane indicates covalence (Figure 8d). In line with other
affects mostly the energy, electronic effects manifest themselves investigations, the DORI analysis distinguishes between the
in qualitative properties of the wave function and the density. It classical 2c-2e C−C bond between bridgehead carbon atoms of
is therefore not surprising that three fused cyclopropane rings large propellane and the atypical noncovalent interactions
with a relatively short inverted bridgehead−bridgehead bond characteristics of [1.1.1]-propellane. It is important to stress
(ca. 1.60 Å93) such as in [1.1.1]-propellane may results in that the seemingly delocalized pattern of the latter is not
exceptional type of interactions. DORI is specifically devised to reflective of a three-dimensional σ- delocalization but rather of a
3751 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756
Journal of Chemical Theory and Computation Article

noncovalent interaction resembling a steric clash imposed by

the small triangular framework.
5.2. Compactness in Supramolecular Chemistry:
Quaterthiophene Derivatives Case. Noncovalent interac-
tions govern various structural and energetic phenomena in
biology, chemistry, and materials science. A highly relevant
example is the performance of organic electronic devices, which
depends strongly on the supramolecular organization of the π-
conjugated units. On this basis, we now push applications of
DORI one step forward and draw a quantitative relationship
between the mutual arrangement of oligothiophene derivatives
in the condensed phase and charge mobility.
Charge transport in organic semiconductors relies upon
numerous factors, including nuclear dynamics and reorganiza-
tion energy with the main prerequisite being a large π-electron
overlap.94 A practical means to enhance charge carrier mobility
is to design molecular crystals in which constituent units are
more densely packed.95 This compactness paradigm is driven
by the relationship between increased molecular compactness
and electronic coupling.
While compactness is easily defined in terms of unit cell size
or more generally, in terms of atom pairwise distances, its
electronic implication is not as straightforward. Surely, more
compact materials exhibit larger overlap of electron densities
but there is no unique way to quantify, visualize and validate
this assumption. Since the electron overlap is the major factor
influencing semiconducting properties, it would be useful to
exploit DORI as a direct measure of “electronic compactness”.
Figure 7. DORI = 0.9 isosurfaces for (a) [1.1.1]-propellane, (b) A well-suited application is to compare the compactness of
[2.1.1]-propellane, (c) [2.2.1]-propellane, (d) [2.2.2]-propellane. bare quaterthiophene molecules in a unit cell with that of
quaterthiophene substituted with terminal hydrogen-bonded
side groups. As recently demonstrated,96 the realization of
unprecedented motifs that include terminal acetamide
functions offers the possibility of intermolecular NH···OC
hydrogen-bonding. The hydrogen-bonded acetamides provide
an ideal way to guide crystallization and improve the
performance of π-type organic semiconductors through the
reinforcement of π- stacking and edge-to-face interactions. The
crystalline quaterthiophene diacetamide is believed to exhibit
denser packing compared to bare quaterthiophene crystal
structure. In fact, the reported volume of a quaterthiophene
diacetamide unit cell is 291 Å3, which is noticeably smaller than
for α-quaterthiophene (307 Å3), α,ω-dimethylquaterthiophene
(297 Å3) and α,ω-dihexylquaterthiophene (307 Å3). It was
proven experimentally that the apparent denser atomic packing
results in enhanced field-effect mobility.
Thus, DORI will serve to visualize and compare the strength
of the electron density overlap between neighboring molecules
in individual crystals. DORI was computed for all nearest-
neighbor dimers in a unit cell of quaterthiophene diacetamide,
α-quaterthiophene, α,ω-dimethylquaterthiophene, and α,ω-
dihexylquaterthiophene. Since, the charge transport properties
depend upon the coupling of π-electrons between the
thiophene moieties, the side chains were replaced by hydrogen
atoms in order to ensure that only the overlap between
quaterhiophene cores’ densities was probed. The X-ray
structures were taken from ref 96, and the positions of the
hydrogen atoms were optimized computationally for a single
unit cell. Figure 9a displays a unit cell of the quaterthiophene
Figure 8. DORI = 0.995 isosurfaces for (a) [1.1.1]-propellane, (b)
[2.1.1]-propellane, (c) [2.2.1]-propellane, (d) [2.2.2]-propellane, with diacetamide without the terminal chains. The molecular
color-coded ▽2ρ(r) in the range from −0.1 au (red) to 0.1 au (blue). labeling is used consistently for all the quaterthiophene
derivatives. Note that the only important structural difference
between the investigated crystals is the symmetry breaking
3752 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756
Journal of Chemical Theory and Computation Article

Figure 9. (a) Structure of a unit cell of the quatertiophene diacetamide; DORI = 0.9 isosurfaces for (b) 1−9, (c) 6−9, α-quaterthiophene dimer, and
(d) 1−9, e) 6−9, quaterthiophene diacetamide dimer. Isosurfaces are color-coded with sgn(λ2)ρ(r) in the range from −0.01 au (red) to 0.01 au
(blue).

Table 1. Electron Density Integrals Over VDORI=0.8

int Intermolecular Interaction Domains
system dimer 1−9 dimer 8−9 dimer 6−9 dimer 1−7
α-quaterthiophene 1.07 1.07 0.78 0.78
α,ω-dimethylquaterthiophene 1.12 1.12 0.80 0.80
α,ω-dihexylquaterthiophene 1.03 1.03 0.77 0.77
quaterthiophene diacetamide 1.31 1.31 0.85 0.84

occurring in the presence of hydrogen bonds, which results in compactness of the diacetamide-containing crystal is also
symmetry distinct (1−7, 6−9) and (1−9, 8−9) dimer pairs in noticeable by the more intense coloring of the respective
quatertiophene diacetamide. The latter crystallizes in the DORI isosurfaces.
triclinic P-1 instead of monoclinic or orthorhombic space These visual indicators are a direct consequence of a stronger
group. density overlap. However, the electronic overlap between
Figures 9b−e display DORI = 0.9 isosurfaces for 1−9 and 6− quaterthiophene units can be placed on a more quantitative
9 dimers of α-quatherthiophene and quaterthiophene diac- ground by exploiting the intermolecular regions identified by
etamide. Color-coded sgn(λ2)ρ(r) is visible on the surfaces. For DORI. The integral of the density over this interaction region
each of the dimers, DORI visualizes four types of interactions: gives the number of overlapping electrons. Obviously, the
covalent bonds, steric clashes at the thiophene ring centers, choice of the isovalue for the surface enclosing the integration
intramolecular noncovalent interactions between the sulfur and volume is arbitrary, but for a fixed and carefully chosen value
hydrogen atoms and the intermolecular interactions between the overlap between analogous dimers in the crystal lattice can
the quaterthiophene units, which are most relevant to the be directly compared. Since the intermolecular interaction
present purpose. A clear manifestation of electronic compact- region may merge with other domains representing intra-
ness is that the DORI domains associated with intermolecular molecular interactions, the isovalue should be either large
interactions are affected by the mutual arrangement of the enough to fully disconnect the distinct domains or relatively
aromatic cores. Even though the changes are subtle, a close small to account for the same interactions in all the systems.
inspection reveals that the domains become more elongated in We have investigated both options, setting the isovalues to
quaterthiophene diacetamide. Furthermore, due to the shorter DORI = 0.8 and DORI = 0.95, which fulfilled the above-
intermolecular distances, a larger number of sulfur−hydrogen mentioned conditions.
domains merge with those arising from intermolecular We define a DORI-based compactness index for the systems
interactions at the chosen isovalue. Finally, the enhanced under investigation as an integral of the electron density over
3753 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756
Journal of Chemical Theory and Computation Article

Table 2. Electron Density Integrals Over VDORI=0.95

int Intermolecular Interaction Domains
system dimer 1−9 dimer 8−9 dimer 6−9 dimer 1−7
α-quaterthiophene 0.20 0.20 0.12 0.12
α,ω-dimethylquaterthiophene 0.23 0.23 0.13 0.13
α,ω-dihexylquaterthiophene 0.23 0.23 0.12 0.12
quaterthiophene diacetamide 0.26 0.24 0.18 0.21

the intermolecular interaction region determined by a DORI Our approach can probe compactness not in terms of nuclear
isosurface arrangement but rather in terms of overlap of electron densities.
The capability of the index was demonstrated on a
IDORI = ∫V DORI
ρ(r)dr
(4)
quaterthiophene derivative designed to crystallize in denser
int packing environment and exhibit better charge transport
This concept is general and can serve to quantify overlap properties.
effects between any two fragments as long as the isovalue is Several attractive features distinguish DORI from other
chosen in a way that produces well-separated domains. It is bonding detectors. First the introduced scalar field depends
important to stress that the DORI-based compactness carries only on the electron density; thus, it is well-defined at any level
different information than the volume of a unit cell, which, of theory. In particular, the DORI analysis is easily applicable to
unlike the present index, also reflects the volume of the densities obtained from post-Hartree−Fock methods, orbital-
terminal chains. free approaches as well as to experimental densities. Second, the
The computed integrals are given in Tables 1 and 2. For the values of the descriptor are system-independent. Due to the
quaterthiophene diacetamide, they are consistently larger than effective [0−1] mapping, bonding patterns and noncovalent
those of other quaterthiophenes, irrespective of the chosen interactions can simultaneously be visualized on equal footing
isovalue and of the dimer considered. This strongly indicates for every system. Finally, the ability of DORI to reveal the local
that insertion of hydrogen-bonded substituents leads to character of electron density is also a promising prerequisite for
its use in the development of approximate density functionals.

■
enhanced density overlap that is at the origin of the measured
increased field-effect mobility. More subtle effects such as the
symmetry breaking expected for the dimer pairs (1−9, 8−9) ASSOCIATED CONTENT
and (6−9, 1−7) or the relative ordering of the three other * Supporting Information
S
quatertiophene crystals are difficult to capture and, con- Expansion of equations in the form explicitly dependent on
sequently, rely more upon the chosen DORI isovalue. The electron density derivatives. This material is available free of
consideration of two distinct isovalues is thus generally charge via the Internet at http://pubs.acs.org/.

■
recommended but one expects that meaningful differences in
density overlap should lead to robust trends independent of the AUTHOR INFORMATION
chosen isovalue.
Corresponding Authors
6. CONCLUSIONS *Email: piotr.desilva@epfl.ch.
*Email: clemence.corminboeuf@epfl.ch.
In this work, we introduced a density-dependent scalar field
Notes
designed to simultaneously identify covalent and noncovalent
The authors declare no competing financial interest.

■
interactions in molecular systems. The proposed quantity,
DORI, is a modification of a previously introduced SEDD
detector,31,32,50 which uses a different reference to obtain the ACKNOWLEDGMENTS
dimensionless quantity. This modification results in appealing Funding from EPFL, from the European Research Council
properties enabling the use of DORI as a universal indicator of (ERC Grants 306528, COMPOREL) and from the Swiss NSF
intra- and intermolecular interactions. Grant no. 200021_137529 is gratefully acknowledged. We
DORI carries information about regions in space where the thank Prof. Holger Frauenrath and his group members for their
total density arises from a strong overlap of individual atomic or inspiring experimental work on quaterthiophene diacetamide.
molecular densities. As such, it should be seen as a scalar field
discovering particular geometrical features of the density, which
can be rationalized in terms of local wavenumber. The
■ REFERENCES
(1) Schleyer, P. v. R. Chem. Rev. 2005, 105, 3433−3435.
analytical properties ensure that DORI is a versatile tool to (2) Popelier, P. L. A. Faraday Discuss. 2007, 135, 3−5.
study bonding patterns and to visualize myriad intra- and (3) Alabugin, I. V.; Gilmore, K. M.; Peterson, P. W. Wiley Interdiscip.
intermolecular interactions, understood as regions of large Rev. Comput. Mol. Sci. 2011, 1, 109−141.
density deformations. The combination of DORI with the (4) Gonthier, J. F.; Steinmann, S. N.; Wodrich, M. D.; Corminboeuf,
analysis of the quantity sgn(λ2)ρ(r) allows for further C. Chem. Soc. Rev. 2012, 41, 4671−4687.
differentiations between bonding and nonbonding interactions (5) Boys, S. F. Rev. Mod. Phys. 1960, 32, 296.
and for an estimate of their magnitude. (6) Edmiston, C.; Ruedenberg, K. Rev. Mod. Phys. 1963, 35, 457.
(7) Pipek, J.; Mezey, P. G. J. Chem. Phys. 1989, 90, 4916.
The utility of DORI was illustrated on various intramolecular (8) Subotnik, J. E.; Dutoi, A. D.; Head-Gordon, M. J. Chem. Phys.
phenomena involving visualization of covalent bonding 2005, 123, 114108.
patterns, steric clashes as well as of typical noncovalent (9) Auqilante, F.; Pedersen, T. B.; de Merás, A. S.; Koch, H. J. Chem.
interactions occurring between and within molecules. The Phys. 2006, 125, 174101.
analysis has also been exploited to visualize and quantify the (10) Jansík, B.; Høst, S.; Kristensen, K.; Jørgensen, P. J. Chem. Phys.
concept of electronic compactness in supramolecular chemistry. 2011, 134, 194104.

3754 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756

Journal of Chemical Theory and Computation Article

(11) de Silva, P.; Giebułtowski, M.; Korchowiec, J. Phys. Chem. Chem. (47) Saleh, G.; Gatti, C.; Lo Presti, L.; Contreras-García, J. Chem.
Phys. 2012, 14, 546−552. Eur. J. 2012, 18, 15523−15536.
(12) de Silva, P.; Makowski, M.; Korchowiec, J. Chimia 2012, 66, (48) Gillet, N.; Chaudret, R.; Contreras-García, I.̀ ; Yang, W.; Silvi, B.;
178−181. Piquemal, J.-P. J. Chem. Theory Comput. 2012, 8, 3993−3997.
(13) Høyvik, I.-M.; Jansik, B.; Jørgensen, P. J. Chem. Phys. 2012, 137, (49) Fang, D.; Chaudret, R.; Piquemal, J.-P.; Cisneros, G. A. J. Chem.
224114. Theory Comput. 2013, 9, 2156−2160.
(14) Høyvik, I.; Jansik, B.; Jørgensen, P. J. Comput. Chem. 2013, 34, (50) de Silva, P.; Korchowiec, J.; Wesolowski, T. A. J. Chem. Phys.
1456−1462. 2014, 140, 164301.
(15) Foster, J.; Weinhold, F. J. Am. Chem. Soc. 1980, 102, 7211− (51) Nagy, A.; March, N. Mol. Phys. 1997, 90, 271−276.
7218. (52) Bohórquez, H. J.; Boyd, R. J. J. Chem. Phys. 2008, 129, 024110.
(16) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. Rev. 1988, 88, (53) Nagy, A.; Liu, S. Phys. Lett. A 2008, 372, 1654−1656.
899−926. (54) Cohen, L. J. Chem. Phys. 1979, 70, 788−789.
(17) Zubarev, D. Y.; Boldyrev, A. I. Phys. Chem. Chem. Phys. 2008, 10, (55) Cohen, L. J. Chem. Phys. 1984, 80, 4277−4279.
5207−5217. (56) Ayers, P. W.; Parr, R. G.; Nagy, A. Int. J. Quantum Chem. 2002,
(18) Zubarev, D. Y.; Boldyrev, A. I. J. Org. Chem. 2008, 73, 9251− 90, 309−326.
9258. (57) Anderson, J. S.; Ayers, P. W.; Hernandez, J. I. R. J. Phys. Chem. A
(19) Becke, A. D.; Edgecombe, K. E. J. Chem. Phys. 1990, 92, 5397. 2010, 114, 8884−8895.
(20) Savin, A.; Nesper, R.; Wengert, S.; Fässler, T. F. Angew. Chem., (58) Kato, T. Commun. Pure Appl. Math. 1957, 10, 151−177.
Int. Ed. Engl. 1997, 36, 1808−1832. (59) Morrell, M. M.; Parr, R. G.; Levy, M. J. Chem. Phys. 1975, 62,
(21) Savin, A.; Becke, A.; Flad, J.; Nesper, R.; Preuss, H.; Von 549.
Schnering, H. Angew. Chem., Int. Ed. Engl. 1991, 30, 409−412. (60) Wang, W.-P.; Parr, R. G. Phys. Rev. A 1977, 16, 891.
(22) Schmider, H.; Becke, A. J. Mol. Struct.: THEOCHEM 2000, 527, (61) Kohout, M.; Savin, A.; Preuss, H. J. Chem. Phys. 1991, 95,
51−61. 1928−1942.
(23) Schmider, H. L.; Becke, A. D. J. Chem. Phys. 2002, 116, 3184− (62) Sperber, G. Int. J. Quantum Chem. 1971, 5, 189−214.
3193. (63) Te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca
(24) Kohout, M. Int. J. Quantum Chem. 2004, 97, 651−658. Guerra, C.; van Gisbergen, S. J.; Snijders, J. G.; Ziegler, T. J. Comput.
(25) Kohout, M.; Pernal, K.; Wagner, F. R.; Grin, Y. Theor. Chem. Chem. 2001, 22, 931−967.
Acc. 2004, 112, 453−459. (64) Kohout, M. DGrid, version 4.6, Radebeul. 2011.
(26) Kohout, M.; Pernal, K.; Wagner, F.; Grin, Y. Theor. Chem. Acc. (65) Henderson, A.; Ahrens, J.; Law, C. The ParaView Guide; Kitware
2005, 113, 287−293. Clifton Park, NY, 2004.
(27) Nalewajski, R. F.; de Silva, P.; Mrozek, J. J. Mol. Struct.: (66) Forslund, L. E.; Kaltsoyannis, N. New J. Chem. 2003, 27, 1108−
THEOCHEM 2010, 954, 57−74. 1114.
(28) Nalewajski, R. F.; de Silva, P.; Mrozek, J. In Theoretical and (67) Politzer, P.; Grice, M. E.; Murray, J. S.; Seminario, J. M. Can. J.
Computational Developments in Modern Density Functional Theory; Roy, Chem. 1993, 71, 1123−1127.
A., Ed.; Nova Science Publishers, Inc., 2012; pp 561−587. (68) Mo, Y.; Schleyer, P. v. R. Chem.Eur. J. 2006, 12, 2009−2020.
(29) Bohórquez, H. J.; Boyd, R. J. Theor. Chem. Acc. 2010, 127, 393− (69) Wu, J. I.; Mo, Y.; Evangelista, F. A.; von Ragué Schleyer, P.
400. Chem. Commun. 2012, 48, 8437−8439.
(30) Bohórquez, H. J.; Matta, C. F.; Boyd, R. J. Int. J. Quantum Chem. (70) Bader, R. J. Chem. Phys. 1980, 73, 2871−2883.
2010, 110, 2418−2425. (71) Tao, J.; Vignale, G.; Tokatly, I. Phys. Rev. Lett. 2008, 100,
(31) de Silva, P.; Korchowiec, J.; Wesolowski, T. A. ChemPhysChem 206405.
2012, 13, 3462−3465. (72) Tachibana, A. J. Mol. Struct.: THEOCHEM 2010, 943, 138−151.
(32) de Silva, P.; Korchowiec, J.; Ram J. S, N.; Wesolowski, T. A. (73) Guevara-García, A.; Echegaray, E.; Toro-Labbe, A.; Jenkins, S.;
Chimia 2013, 67, 253−256. Kirk, S. R.; Ayers, P. W. J. Chem. Phys. 2011, 134, 234106.
(33) Bader, R. F. Atoms in Molecules: A Quantum Theory; Oxford (74) Jenkins, S.; Kirk, S. R.; Guevara-García, A.; Ayers, P. W.;
University Press, 1994. Echegaray, E.; Toro-Labbe, A. Chem. Phys. Lett. 2011, 510, 18−20.
(34) Bader, R. F. Chem. Rev. 1991, 91, 893−928. (75) Jurec̆ka, P.; Šponer, J.; C̆ ernỳ, J.; Hobza, P. Phys. Chem. Chem.
(35) Bader, R.; MacDougall, P.; Lau, C. J. Am. Chem. Soc. 1984, 106, Phys. 2006, 8, 1985−1993.
1594−1605. (76) Gráfová, L.; Pitonák, M.; Rezac, J.; Hobza, P. J. Chem. Theory
(36) Bader, R.; Slee, T.; Cremer, D.; Kraka, E. J. Am. Chem. Soc. Comput. 2010, 6, 2365−2376.
1983, 105, 5061−5068. (77) de Silva, P.; Wesolowski, T. A. Phys. Rev. A 2012, 85, 032518.
(37) Jenkins, S. J. Phys.: Condens. Matter 2002, 14, 10251. (78) Jackson, J. E.; Allen, L. C. J. Am. Chem. Soc. 1984, 106, 591−
(38) Jenkins, S.; Ayers, P. W.; Kirk, S. R.; Mori-Sánchez, P.; Martn 599.
Pendás, A. Chem. Phys. Lett. 2009, 471, 174−177. (79) Wiberg, K. B.; Bader, R. F.; Lau, C. D. J. Am. Chem. Soc. 1987,
(39) Seriani, N. J. Phys.: Condens. Matter 2010, 22, 255502. 109, 985−1001.
(40) Koritsanszky, T. S.; Coppens, P. Chem. Rev. 2001, 101, 1583− (80) Kar, T.; Jug, K. Chem. Phys. Lett. 1996, 256, 201−206.
1628. (81) Adcock, W.; Brunger, M.; Clark, C.; McCarthy, I.; Michalewicz,
(41) Espinosa, E.; Molins, E.; Lecomte, C. Chem. Phys. Lett. 1998, M.; Von Niessen, W.; Weigold, E.; Winkler, D. J. Am. Chem. Soc. 1997,
285, 170−173. 119, 2896−2904.
(42) Volkov, A.; Gatti, C.; Abramov, Y.; Coppens, P. Acta Crystallogr., (82) Ebrahimi, A.; Deyhimi, F.; Roohi, H. J. Mol. Struct.:
Sect. A 2000, 56, 252−258. THEOCHEM 2003, 626, 223−229.
(43) Tsirelson, V.; Stash, A. Chem. Phys. Lett. 2002, 351, 142−148. (83) Nalewajski, R. F.; Broniatowska, E. J. Phys. Chem. A 2003, 107,
(44) Tsirelson, V.; Stash, A. Acta Crystallogr. Sect. B 2002, 58, 780− 6270−6280.
785. (84) Polo, V.; Andres, J.; Silvi, B. J. Comput. Chem. 2007, 28, 857−
(45) Johnson, E. R.; Keinan, S.; Mori-Sanchez, P.; Contreras-García, 864.
J.; Cohen, A. J.; Yang, W. J. Am. Chem. Soc. 2010, 132, 6498−6506. (85) Wu, W.; Gu, J.; Song, J.; Shaik, S.; Hiberty, P. C. Angew. Chem.,
(46) Contreras-García, J.; Johnson, E. R.; Keinan, S.; Chaudret, R.; Int. Ed. Engl. 2009, 121, 1435−1438.
Piquemal, J.-P.; Beratan, D. N.; Yang, W. J. Chem. Theory Comput. (86) Shaik, S.; Chen, Z.; Wu, W.; Stanger, A.; Danovich, D.; Hiberty,
2011, 7, 625−632. P. C. ChemPhysChem 2009, 10, 2658−2669.

3755 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756

 Journal of Chemical Theory and Computation Article

(87) Seiler, P.; Belzner, J.; Bunz, U.; Szeimies, G. Helv. Chim. Acta
1988, 71, 2100−2110.
(88) Messerschmidt, M.; Scheins, S.; Grubert, L.; Pätzel, M.;
Szeimies, G.; Paulmann, C.; Luger, P. Angew. Chem., Int. Ed. Engl.
2005, 44, 3925−3928.
(89) Schleyer, P. v. R.; McKee, W. C. J. Phys. Chem. A 2010, 114,
3737−3740.
(90) Dewar, M. J. Bull. Soc. Chim. Belg. 1979, 88, 957−967.
(91) Wu, W.; Ma, B.; Wu, J. I.; von Ragué Schleyer, P.; Mo, Y.
Chem.Eur. J. 2009, 15, 9730−9736.
(92) Wu, J. I.; von Ragué Schleyer, P. Pure Appl. Chem. 2013, 85.
(93) Wiberg, K. B.; Dailey, W. P.; Walker, F. H.; Waddell, S. T.;
Crocker, L. S.; Newton, M. J. Am. Chem. Soc. 1985, 107, 7247−7257.
(94) Coropceanu, V.; Cornil, J.; da Silva Filho, D. A.; Olivier, Y.;
Silbey, R.; Brédas, J.-L. Chem. Rev. 2007, 107, 926−952.
(95) Brédas, J.-L.; Calbert, J. P.; da Silva Filho, D.; Cornil, J. Proc.
Natl. Acad. Sci. U.S.A. 2002, 99, 5804−5809.
(96) Gebers, J.; Hartmann, L.; Schaer, M.; Suàrez, S.; Bugnon, P.;
Steinrü c k, H.-G.; Magerl, A.; Brinkmann, M.; Petraglia, R.;
Corminboeuf, C.; Frauenrath, H. (submitted) 2014,.

3756 dx.doi.org/10.1021/ct500490b | J. Chem. Theory Comput. 2014, 10, 3745−3756

View publication stats