Prachi Sharma‡
^{a},
Varinia Bernales‡^{a},
Stefan Knecht*^{b},
Donald G. Truhlar*^{a} and
Laura Gagliardi*^{a}
^{a}Department of Chemistry, Chemical Theory Center, Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455, USA. E-mail: truhlar@umn.edu; gagliard@umn.edu
^{b}Laboratory of Physical Chemistry, ETH Zürich, Vladimir-Prelog-Weg 2, CH-8093 Zürich, Switzerland. E-mail: stefan.knecht@phys.chem.ethz.ch
First published on 26th November 2018
The density matrix renormalization group (DMRG) is a powerful method to treat static correlation. Here we present an inexpensive way to calculate correlation energy starting from a DMRG wave function using pair-density functional theory (PDFT). We applied this new approach, called DMRG-PDFT, to study singlet–triplet gaps in polyacenes and polyacetylenes that require active spaces larger than the feasibility limit of the conventional complete active-space self-consistent field (CASSCF) method. The results match reasonably well with the most reliable literature values and have only a moderate dependence on the compression of the initial DMRG wave function. Furthermore, DMRG-PDFT is significantly less expensive than other commonly applied ways of adding additional correlation to DMRG, such as DMRG followed by multireference perturbation theory or multireference configuration interaction.
One example of strongly-correlated systems is the family of polyacenes (also called acenes), the ground states of which have polyradical character.^{9} The acenes are a series of compounds consisting of fused benzene rings; in the present article we only consider the case where they are arranged in a linear fashion as shown in Fig. 1. Polyacenes are especially interesting because their electronic properties make them useful for biodegradable and low-cost electronics.^{10–13} The smaller ones, such as naphthalene and anthracene, are used industrially to make various dyes, while larger acenes, like tetracene and pentacene along with their derivatives, are used as organic light-emitting diodes (OLEDs) and organic field-effect transistors (OFETs).^{14–25} Recently, tetracene and pentacene have been found to undergo singlet-fission where a high-energy singlet exciton converts to two low-energy triplet excitons, thereby increasing the efficiency of solar cells by up to 40%.^{26–33}
Polyenes are another example of strongly correlated systems that contains alternating single and double bonds; in the present article we will consider polyacetylenes as examples of this class of compounds. They are often used as model systems to understand the electronic properties of more complicated biological systems such as carotenoids, retinal, and various antibiotics such as amphotericin B, nystatin, etc.^{34,35} Polyacetylenes also show high electrical conductivity, which further increases on doping with p-type dopants such as Br_{2}, I_{2}, Cl_{2}, and AsF_{5}.^{36–40}
A theoretical method that plays a key role in many methods for treating multireference systems is CASSCF,^{41–43} which is a special case of the multiconfiguration self-consistent-field (MCSCF) method. In the CASSCF method, some orbitals are restricted to be doubly occupied in all configuration state functions; these are called inactive orbitals. An additional group of orbitals, called active orbitals, are allowed to have variable occupancy. The CASSCF wave function is then a configuration interaction wave function including all possible occupancies of the active orbitals; this is called full configuration interaction (FCI) in the active space. CASSCF is limited by the exponential increase in cost with size of the active space such that 20 electrons in 20 orbitals (20, 20) is the largest reported active space treated by the conventional CASSCF method.^{44} Although the CASSCF approach was originally designed as a way to recover the static correlation, inevitably it also includes some of the remaining electron correlation, which is called dynamic correlation, but usually only a small fraction of it. The dynamic correlation not included in a CASSCF calculation is called external correlation,^{45–47} and one of the main challenges in quantum chemistry is to design methods to include the external correlation energy efficiently.
The density matrix renormalization group (DMRG)^{48–51} algorithm is a replacement of the FCI solver in the CASSCF or CASCI methods.^{52,53} (CASCI is like CASSCF, but the orbitals are not self-consistently optimized for a given configuration list. In the present work we use self-consistently optimized orbitals.) The major advantage of DMRG comes from its polynomial scaling with respect to the size of the active space; this allows practical computation of numerically well-converged solutions for active spaces three to four times larger than standard CASSCF.^{54} Although multiconfigurational methods such as CASSCF or DMRG can be effectively applied to capture static correlation, and they in principle reduce the difficulty in choosing which orbitals to include in the active space (because they allow more orbitals to be included), they are not efficient at capturing all the correlation energy because of the slow convergence of external correlation with respect to the number of orbitals and configurations included.
A number of approaches have been advanced to make it more efficient to treat the dynamical correlation not included in an MCSCF-type calculation. For example, approaches that have been combined with DMRG include internally contracted multireference CI (MRCI),^{55,56} multireference perturbation theory (MRPT),^{57–68} and wave function theory-short range density functional theory (WFT-srDFT).^{69} In the present study we introduce the use of multiconfiguration pair-density functional theory (MC-PDFT)^{70,71} to include dynamic electron correlation beyond that captured within a DMRG wave function; this two-step procedure will be called DMRG-PDFT. The major advantage of DMRG-PDFT is that it treats both static and dynamic electron-correlation and gives accurate results at significantly lower memory and computational costs than multireference perturbation treatments and MRCI approaches. Therefore, it can attain high accuracy that would be too expensive to be practical with other methods.
We briefly review DMRG and MC-PDFT in Section 2, followed by a presentation of DMRG-PDFT. The computational methods are in Section 3. Results and discussion are in Section 4, and conclusions are in Section 5.
The non-relativistic electronic Hamiltonian may be written in second quantization as^{81}
(1) |
(2) |
(3) |
Each spatial orbital l (also referred to as a site in DMRG terminology) is associated with four possible occupations:
n_{l} = {|vac〉,|↑〉,|↓〉,|↑↓〉}. | (4) |
The eigenstate |Ψ〉 of the Hamiltonian (1) can be written in an occupation number basis as
(5) |
(6) |
The size of the matrices A^{nl} increases exponentially towards the middle of the product, then decreases again.^{53} Eqn (6) is called a matrix product state (MPS). Note, that the first and the last matrices in eqn (6) are in practice 1 × 4-dimensional row and 4 × 1-dimensional column vectors, respectively, such that the final contraction of the matrices A^{nl} yields a scalar CI coefficient C^{n1n2…nL}.
If no truncation is made in the matrix products, the MPS is equivalent to full configuration interaction (FCI). For practical work, one retains at most M terms in each of the matrix multiplication steps; this approximation is called compression, and M is called the bond dimension. Compression is the main approximation of DMRG as compared to FCI. The compressed MPS is then variationally optimized to give an upper bound to the ground-state energy. Thus, DMRG may be considered to be a way to calculate an approximate wave function that gives an upper bound to the FCI energy. We note though that we (and others) apply FCI only within the active space, so what we obtain is an upper bound to the CASSCF energy, which in turn is an upper bound to the FCI energy, which is an upper bound to the complete configuration interaction (CCI) energy, which is exact (CCI is FCI with a complete basis set).
As we have just explained, if one did not truncate the bond dimension in DMRG, the resulting wave function would correspond to FCI among the active orbitals, i.e., to CASSCF. In this context, DMRG may be considered to offer a possibility to calculate an approximate CASSCF wave function in a controlled manner such that it gives a variational upper bound to the CASSCF energy. However, its efficiency allows one to use many more active orbitals than in conventional CASSCF. The approximate solution to the large-active-space CASSCF problem will typically have a lower (and hence more accurate) energy than the uncompressed solution to the small-active-space CASSCF problem.
(7) |
The MC-PDFT energy can be written in terms of inactive and active orbitals as
(8) |
(9) |
(10) |
The active space is denoted as usual as (n, m) where n is the number of active electrons, and m is the number of active orbitals. The active spaces used here correspond to all the π electrons distributed in all the valence π orbitals; thus n = m. For all polyacene calculations, we constrained the wave functions to D_{2h} symmetry, where the ground state has symmetry ^{1}A_{g}, and the lowest triplet state has symmetry ^{3}B_{3u}. For the polyacetylene calculations, C_{2h} symmetry was imposed. The ground state has symmetry ^{1}A_{g}, and the lowest triplet state has symmetry ^{3}B_{u}.
The convergence with respect to the bond dimension M was tested from 100 to 2000 for naphthalene and anthracene, and M was set to 500 for higher acenes. The convergence for polyacetylenes was studied from M = 4 to M = 1000. The DMRG calculations were performed with QCMAQUIS.^{72,75,84,89}
M | DMRG | DMRG-PDFT | GAS-PDFT^{a} | Reference values | ||||
---|---|---|---|---|---|---|---|---|
Vert. | Ad. | Vert. | Ad. | Vert. | Ad. | Vert. | Ad.^{d} | |
a Generalized active-space pair-density functional theory with tPBE functional and 6-31G+(d,p) basis set from ref. 85.b DMRG-externally correlated multireference CI, DMRG-ec-MRCISD+Q at geometries optimized by UB3LYP/6-31G(d); see ref. 91 for details.c Restricted CCSD(T) with pV∞Z basis set at geometries optimized by B3LYP/cc-PVTZ; see ref. 90 for details.d Vibrationally corrected experimented values; see ref. 95 and 96 for experimental details. ΔZPE = −0.14 eV calculated by B3LYP/6-31G(d,p). | ||||||||
3.36 | 3.06 | 3.43,^{b} 3.30^{c} | 2.78, 2.79 | |||||
100 | 3.08 | 2.67 | 3.31 | 2.89 | ||||
200 | 3.05 | 2.66 | 3.35 | 2.91 | ||||
500 | 3.05 | 2.66 | 3.35 | 2.91 | ||||
1000 | 3.05 | 2.66 | 3.35 | 2.91 |
In Table 2, vertical and adiabatic S–T gaps for different bond dimensions are reported for anthracene. We see convergence at M = 1000, and we observe a significant improvement in both vertical and adiabatic singlet–triplet gaps when the MC-PDFT method is used to calculate correlation energy starting from a DMRG wave function. The S–T gap predicted by DMRG decreases when the bond dimension is increased, but the DMRG-PDFT gap increases along the same sequence of M values. A key result is that the DMRG-PDFT results show less dependence on M than do the plain DMRG results; in fact DMRG-PDFT gives results close to the best previously available values^{90–93,95} even for an M value as small as M = 100.
M | DMRG | DMRG-PDFT | GAS-PDFT^{a} | Other literature values | ||||
---|---|---|---|---|---|---|---|---|
Vert. | Ad. | Vert. | Ad. | Vert. | Ad. | Vert. | Ad.^{d} | |
a Generalized active-space pair-density functional theory with tPBE functional and 6-31G+(d,p) basis set from ref. 85.b DMRG-externally correlated multireference CI, DMRG-ec-MRCISD+Q at geometries optimized by UB3LYP/6-31G(d); see ref. 91.c Restricted CCSD(T) (FPA-5Z3) at geometries optimized by B3LYP/cc-PVTZ. FPA − 5Z3 = (E_{CCSD(T)/cc-pVTZ} − E_{MP4/cc-pVTZ}) + S_{MP2} − 5Z + (S_{MP4} − 4Z − S_{MP2} − 4Z), where S_{MP2} − 5Z are obtained as the sum of the HF energy and MP2 electron correlation energy, both extrapolated to the CBS limit using Schwartz extrapolations employing HF and MP2 energies obtained using the cc-pVTZ, cc-pVQZ, and cc-pV5Z basis sets; ref. 90.d Vibrationally corrected experimented values; see ref. 95 and 97 for experimental details. ΔZPE = −0.10 eV calculated by B3LYP/6-31G(d,p). | ||||||||
2.22 | 1.97 | 2.47,^{b} 2.46^{c} | 1.95–1.97 | |||||
100 | 2.46 | 2.05 | 2.28 | 2.00 | ||||
200 | 2.42 | 2.03 | 2.26 | 1.97 | ||||
500 | 2.34 | 1.97 | 2.33 | 2.00 | ||||
1000 | 2.31 | 1.96 | 2.36 | 2.02 | ||||
2000 | 2.30 | 1.95 | 2.38 | 2.04 |
The adiabatic S–T gaps for various systems, from naphthalene to heptacene, obtained with M = 500 are compared to literature values in Table 3. We observe that DMRG-PDFT gives very good agreement with the vibrationally corrected experimental S–T gaps for adiabatic singlet–triplet gaps with a mean unsigned deviation (MUD) of only 0.06 eV while DMRG without external correlation gives an MUD of 0.17 eV. While DMRG-CASPT2 and ACI-DSRG-MRPT2 underestimate the S–T gap with respect to the experiment, CCSD(T) tends to overestimate. The generalized active-space-pair density functional theory (GAS-PDFT) gives excellent agreement with the experimental results. In GAS^{101} based methods, the active space is further divided into smaller subspaces and only certain excitations are allowed within them, thereby making these methods more affordable. One difficulty with these approaches, though, is to choose the appropriate subspaces and appropriate restrictions on the kinds of included electron excitations because such choices are not systematic and often require subjective chemical intuition and specialized expertise.
DMRG^{a} | DMRG-PDFT^{a} | GAS-PDFT^{b} | pp-RPA^{c} | DMRG-CASPT2^{d} | DMRG-CASPT2^{e} | CCSD(T)^{f} | DMRG-ec-MRCISD+Q^{g} | ACI-DSRG-MRPT2^{h} | ΔZPE^{i} | Exp. | |
---|---|---|---|---|---|---|---|---|---|---|---|
a M = 500.b Generalized active-space pair-density functional theory with the tPBE functional and the 6-31+G(d,p) basis set; for details see ref. 85.c Particle–particle random phase approximation (pp-RPA) with the B3LYP functional at geometries optimized by UB3LYP/6-31G*. See ref. 13 for details.d DMRG-CASPT2 at geometries optimized by CASPT2-D/cc-pVTZ(-f); see ref. 98 for details.e DMRG-CASPT2 at geometries optimized by CAM-B3LYP/6-31G*; see ref. 98 for details.f Restricted CCSD(T) (FPA-5Z3) at geometries optimized by B3LYP/cc-PVTZ with added ZPE correction. FPA − 5Z3 = (E_{CCSD(T)/cc-pVTZ} − E_{MP4/cc-pVTZ}) + S_{MP2} − 5Z + (S_{MP4} − 4Z − S_{MP2} − 4Z), where S_{MP2} − 5Z are obtained as the sum of the HF energy and MP2 electron correlation energy, both extrapolated to the CBS limit using Schwartz extrapolations employing HF and MP2 energies obtained using the cc-pVTZ, cc-pVQZ, and cc-pV5Z basis sets.; see ref. 90 for details.g DMRG-ec-MRCISD+Q on geometry optimized by UB3LYP/6-31G(d); see ref. 91 for details.h Adaptive CI with density a density-fitted implementation of second-order perturbative multiconfiguration driven similarity renormalization group (ACI-DSRG-MRPT2) calculations; see ref. 99.i Zero-point energy correction; calculated by B3LYP/6-31G(d,p).j See ref. 95.k Ref. 96.l Ref. 97.m Ref. 100.n Ref. 94.o Mean unsigned deviation with respect to vibrationally corrected experimental S–T gap. | |||||||||||
Naphthalene (10, 10) | 2.66 | 2.91 | 3.06 | 2.87 | — | — | 2.85 | 2.71 | 2.70 | −0.14 | 2.64^{j}, 2.65^{k} |
Anthracene (14, 14) | 1.97 | 2.00 | 1.97 | 1.98 | 1.73 | 1.69 | 2.09 | 1.81 | 1.87 | −0.10 | 1.85^{j}, 1.87^{l} |
Tetracene (18, 18) | 1.54 | 1.37 | 1.46 | 1.39 | 1.29 | 1.18 | 1.45 | 1.23 | 1.23 | −0.08 | 1.28^{j} |
Pentacene (22, 22) | 1.24 | 0.98 | 1.10 | 0.98 | 0.86 | 0.82 | 1.10 | 0.92 | 0.78 | −0.06 | 0.86 ± 0.03^{m} |
Hexacene (26, 26) | 0.93 | 0.73 | 0.85 | 0.66 | 0.58 | 0.62 | 0.77 | 0.67 | 0.49 | −0.06 | 0.54 ± 0.05^{n} |
Heptacene (30, 30) | 0.67 | 0.62 | 0.72 | 0.39 | — | — | 0.58 | 0.48 | 0.33 | −0.05 | — |
MUD^{o} | 0.17 | 0.06 | 0.04 | 0.04 | 0.08 | 0.13 | 0.17 | 0.07 | 0.09 |
The vertical S–T gaps are reported in Table 4. Because it underestimates external correlation, DMRG overestimates the S–T gap for acenes larger than naphthalene. The adiabatic and vertical S–T gaps for DMRG and DMRG-PDFT are compared with experimental values for larger polyacenes in Fig. 2, showing reasonable agreement for both DMRG and DMRG-PDFT adiabatic S–T gaps.
DMRG | DMRG-PDFT^{a} | GAS-PDFT^{b} | CCSD(T)^{c} | DMRG-MRCISD+Q^{d} | |
---|---|---|---|---|---|
a M = 500.b GAS-PDFT with 6-31+G(d,p) basis set; for details see ref. 85.c Restricted CCSD(T) with pV∞Z basis set at geometries optimized by B3LYP/cc-PVTZ; see ref. 90 for details.d DMRG-ec-MRCISD+Q on geometry optimized by UB3LYP/6-31G(d); see ref. 91 for details. | |||||
Naphthalene (10, 10) | 3.05 | 3.35 | 3.36 | 3.30 | 3.43 |
Anthracene (14, 14) | 2.34 | 2.33 | 2.22 | 2.46 | 2.47 |
Tetracene (18, 18) | 1.88 | 1.58 | 1.69 | 1.75 | 1.81 |
Pentacene (22, 22) | 1.56 | 1.13 | 1.29 | 1.36 | 1.36 |
Hexacene (26, 26) | 1.19 | 0.79 | 0.99 | 0.99 | 0.98 |
Heptacene (30, 30) | 0.81 | 0.61 | 0.75 | 0.78 | 0.67 |
Fig. 2 Adiabatic and vertical singlet–triplet gaps for polyacenes. Experimental values are vibrationally corrected. |
Before we conclude this session, we would like to briefly comment on the active space dependence of DMRG and DMRG-PDFT. It is well known that all multi-reference methods, both standard ones and new ones, have some dependence on active space. Table S5 in ESI† shows the effect of active-space expansion in the naphthalene system. DMRG-PDFT has an equal or smaller dependence on active space than does DMRG alone. In going from a (10, 10) active space to a (10, 20) active space both the vertical and adiabatic DMRG-PDFT excitation energies change by less than 0.1 eV for M = 500. The active space dependence of PDFT has been discussed in several papers (see for example ref. 70 and 71). The effect is not aggravated by using DMRG as a reference wave function.
Number of monomers | Active space | DMRG | DMRG-PDFT | CASSCF | CASPT2 | Literature values | Exp.^{a} |
---|---|---|---|---|---|---|---|
a Experimental band maxima for ethylene,^{102–105} butadiene,^{106} and hexatriene.^{107}b CCSD(T)/cc-pVTZ result from ref. 108.c UCCSD result from ref. 109.d Multireference Møller–Plesset study corrected for basis-set and active-space effects, from ref. 110.e Mean unsigned deviation from experiment. | |||||||
1 | (2, 2) | 4.34 | 4.67 | 4.34 | 4.54 | 4.63^{b} | 4.3–4.6 |
2 | (4, 4) | 3.37 | 3.46 | 3.37 | 3.38 | 3.45^{b}, 3.20^{d} | 3.22 |
3 | (6, 6) | 2.80 | 2.79 | 2.80 | 2.73 | 2.80^{b}, 2.40^{d} | 2.61 |
4 | (8, 8) | 2.43 | 2.37 | 2.43 | 2.33 | 2.42^{c}, 2.10^{d} | 2.10 |
5 | (10, 10) | 2.29 | 1.99 | 2.19 | 2.07 | 2.20^{c}, 1.89^{d} | |
6 | (12, 12) | 2.20 | 1.79 | 2.01 | 1.88 | 2.00^{c} | |
7 | (14, 14) | 2.17 | 1.59 | 1.88 | 1.75 | 1.90^{c} | |
8 | (16, 16) | 2.20 | 1.52 | ||||
9 | (18, 18) | 1.03 | 0.07 | ||||
MUD^{e} | 0.20 | 0.18 | 0.20 | 0.16 |
We plotted the convergence of singlet–triplet gaps in DMRG and DMRG-PDFT with respect to conventional CASSCF and MC-PDFT in Fig. 3. DMRG-PDFT converges faster with respect to M and shows less dependence on M than does DMRG. In Fig. 4, we show the average time required for DMRG and DMRG-PDFT with M = 50 as compared to the time for CASPT2 calculations for polyacetylenes (the times are averaged over singlet and triplet). All the calculations were performed on a single processor with a maximum memory of 62 gigabytes. The DMRG-PDFT calculations take considerably less time than the corresponding CASPT2 calculations for higher polyacetylenes (five or more monomers). For example, for seven monomers, CASPT2 takes fifteen times longer than DMRG-PDFT.
Footnotes |
† Electronic supplementary information (ESI) available: Absolute energies in hartrees, timings, geometry comparison, and Cartesian coordinates. See DOI: 10.1039/c8sc03569e |
‡ These authors contributed equally. |
This journal is © The Royal Society of Chemistry 2019 |