Automatic Construction of the Initial Orbitals for Efficient Generalized

Subscriber access provided by University of Winnipeg Library

Quantum Electronic Structure

Automatic Construction of the Initial Orbitals for Efficient Generalized Valence Bond Calculations of Large Systems Qingchun Wang, Jingxiang Zou, Enhua Xu, Peter Pulay, and Shuhua Li J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00854 • Publication Date (Web): 27 Nov 2018 Downloaded from http://pubs.acs.org on November 29, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

1 2

Automatic Construction of the Initial Orbitals for Efficient Generalized Valence Bond Calculations of Large Systems

3

Qingchun Wang,† Jingxiang Zou,† Enhua Xu,‡ Peter Pulay,*,¶ Shuhua Li*,†

4 5 6 7 8 9 10

of Chemistry and Chemical Engineering, Key Laboratory of Mesoscopic Chemistry of Ministry of Education, Institute of Theoretical and Computational Chemistry, Nanjing University, 210023, P. R. China ‡Graduate School of Science, Technology, and Innovation, Kobe University, Nada-ku, Kobe 6578501, Japan ¶Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701, United States

11

KEYWORDS: generalized valence bond theory, initial orbitals, active space, large

12

systems

13

ABSTRACT: We propose an efficient general strategy for generating initial orbitals for

14

generalized valence bond (GVB) calculations which makes routine black-box GVB calculations on

15

large systems feasible. Two schemes are proposed, depending on whether the restricted Hartree-

16

Fock (RHF) wavefunction is stable (scheme I) or not (scheme II). In both schemes, the first step is

17

the construction of active occupied orbitals and active virtual orbitals. In scheme I, active occupied

18

orbitals are composed of the valence orbitals (the inner core orbitals are excluded) and the active

19

virtual orbitals are obtained from the original virtual space by requiring its maximum overlap with

20

the virtual orbital space of the same system at a minimal basis set. In scheme II, active occupied

21

orbitals and active virtual orbitals are obtained from the set of unrestricted natural orbitals (UNOs),

22

which are transformed from two sets of unrestricted HF spatial orbitals. In the next step, the active

23

occupied orbitals and active virtual ones are separately transformed to localized orbitals. Localized

24

occupied and virtual orbital pairs are formed using the Kuhn-Munkres (KM) algorithm and are used

25

as the initial guess for the GVB orbitals. The optimized GVB wavefunction is obtained using the

26

second-order self-consistent-field algorithm in the GAMESS program. With this procedure, GVB

27

energies have been obtained for the lowest singlet and triplet states of polyacenes (up to decacene

28

with 96 pairs) and the singlet ground state of two di-copper-oxygen-ammonia complexes. We have

29

also calculated the singlet-triplet gaps for some polyacenes and the relative energy between two di-

30

copper-oxygen-ammonia complexes with the block-correlated second-order perturbation theory

31

based on the GVB reference.

†School

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Page 2 of 28

I. INTRODUCTION

2

The development of theoretical methods capable of describing the electronic structures of

3

molecules at arbitrary structures (from equilibrium to dissociation) is still a very active field.1, 2 Near

4

equilibrium, the Hartree-Fock (HF) wavefunction is usually a good reference function, and HF-

5

based single reference (SR) electron correlation methods, such as second-order Møller-Plesset

6

perturbation theory (MP2),3 Coupled Cluster (CC) singles and doubles (CCSD),4, 5 and CC singles,

7

doubles, and perturbative triples (CCSD(T)),6 are very successful.7-9 However, the accuracy of

8

single reference electron correlation methods is not satisfactory for molecules far from the

9

equilibrium region, diradical or polyradical systems, and transition metal-containing compounds.10-

10

16

11

determinant. To give a qualitatively correct description of these systems, a multireference zeroth-

12

order wavefunction, which includes all important configurations, should be used.

In such strongly correlated systems other determinants may be nearly as important as the HF

13

The most common multireference (MR) wavefunction is the complete active space self-

14

consistent-field (CASSCF) wavefunction17 which contains all configurations within a limited orbital

15

space, the active space. The active space is composed of partially filled, usually near-degenerate

16

orbitals, e.g., bonding and anti-bonding orbitals of breaking bonds in chemical reactions.18, 19 Based

17

on the CASSCF reference function, various electron correlation methods have been proposed to

18

further incorporate the dynamic correlation energy in a similar way as post-HF methods. For

19

example, CASPT2 (CASSCF-based second-order perturbation theory)20, 21 is one of the most widely

20

used CASSCF-based methods. These CASSCF-based methods have been demonstrated to be quite

21

successful for a variety of chemical problems.22-25 However, these approaches are often limited to

22

systems with small or medium-sized active spaces (less than 20 electrons in 20 orbitals),26 since the

23

computational cost increases exponentially with the size of the active space. In recent years, several

24

new MR electron correlation methods have been developed and applied to a number of chemical

25

systems.27-31 For example, for linear or quasi-linear molecules with larger active spaces (such as 50

26

electrons in 50 active orbitals), the density matrix renormalization group (DMRG)29-31 can give quite

27

accurate energies for the ground-state and a few low-lying excited states. However, further

28

perturbative treatments on the top of the DMRG reference function are still needed to obtain

29

quantitative descriptions. For all these MR electron correlation methods, a key question is how to 2


Page 3 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

choose a suitable active space for a given system. In most cases, the selection of the active space is

2

done manually with the help of chemical intuition and experience. It should be mentioned that a

3

number of different techniques have been proposed for the automatic selection of active space.18, 19,

4

32-38

5

natural orbitals (UNOs)32 proposed by Pulay is the simplest one, which can automatically lead to

6

active orbitals. A number of MR electron correlation methods based on the UNOs have been

7

established.39-46 However, the active space selected using the UNOs criterion may change with the

8

geometry of the molecule, leading to discontinuities on the potential energy surface.32, 34

Among these techniques, the transformation from unrestricted HF orbitals to unrestricted

9

Other multi-reference wavefunctions, for instance the antisymmetric product of strongly

10

orthogonal geminals (APSG)47-51 and its special case, the perfect-pairing generalized valence bond

11

(GVB-PP, denoted as GVB in short) wavefunction,52 have also been considered as the reference

12

functions for strongly correlated systems. Although the APSG and GVB reference functions can

13

only recover a small portion of static correlation for general molecules, they are able to give correct

14

descriptions of single bond breaking processes. They are also inexpensive in comparison to

15

CASSCF and DMRG. A number of electronic structure methods based on the GVB reference

16

function have been developed.13,

17

valence bond (CCVB) approach,61, 62 which provides good descriptions for a number of strongly

18

correlated systems. We proposed a block-correlated second-order perturbation theory based on the

19

GVB

20

ground-state energies than the GVB method, but its accuracy is not high enough for quantitative

21

descriptions in some cases.13, 67 In comparison with CASSCF and DMRG, GVB has not been widely

22

used as a reference wavefunction. A possible reason may be ascribed to the fact that an efficient

23

optimization of the GVB wavefunction for medium-sized and large systems remains a challenging

24

problem. In fact, several methods have been developed for optimizing the GVB wavefunction of a

25

general molecule.49, 68-70 Among these methods, the second-order SCF (SOSCF) method generally

26

shows a good performance in comparison with other existing methods. However, this method also

27

suffers from the convergence problem in many cases, because it strongly depends on the initial guess

28

of the GVB orbitals (due to its second-order nature). Several approaches71, 72 of building an initial

29

guess of the GVB orbitals have been proposed. For example, Sano constructs the initial guess by

reference

function

53-66

The Head-Gordon group developed the coupled cluster

(GVB-BCPT2).13

This

method

produces

much

better

3



Page 4 of 28

1

using localized occupied molecular orbitals for the first (more strongly occupied or bonding) GVB-

2

NO orbital in a geminal pair, and maximizing the exchange integral between the 1st GVB-NOs and

3

the 2nd (usually antibonding) GVB-NOs in the virtual space.72 In spite of promising early results,

4

there has been little follow-up of this work, and the applications did not include strongly correlated

5

systems. Thus, it is still highly desirable to have a simple and effective strategy for constructing

6

automatically a good general initial GVB guess.

7

In this article, a robust and efficient procedure will be proposed to build the initial orbitals for

8

GVB calculations. Depending on whether the restricted HF (RHF) solution is stable or not, one or

9

two schemes based on the canonical HF orbitals will be adopted to construct initial orbitals of the

10

GVB orbitals. With these initial orbitals, the SOSCF algorithm implemented in the GAMESS73

11

program will be used to perform GVB calculations of a variety of systems including a series of

12

polyacenes and two transition-metal compounds. Our calculations demonstrated that converged

13

GVB orbitals can be easily obtained for large organic systems such as decacene (with 96 pairs). In

14

addition, we will show that with the GVB orbitals an automatic selection of the active space is

15

possible, and the selected active space may be used for subsequent multi-reference electron

16

correlation calculations.

17

The paper is arranged as follows: In Section. II, the GVB method, and our two schemes of

18

constructing initial guess GVB orbitals are introduced. In Section. III, we will perform GVB

19

calculations on a variety of molecules, and perform several multi-reference electron correlation

20

calculations with the GVB reference or the active space determined from the GVB orbitals to obtain

21

accurate relative energies of two different states or two isomers. Finally, a brief summary will be

22

given in Section IV.

23

II. METHODOLOGY

24

A. A Summary of Generalized Valence Bond Theory

25 26 27

For a given system, the general form of the many-electron GVB wavefunction (in second quantized form) can be defined as

 GVB   core  pair  open vac

(1)

4


Page 5 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Nc

 core   ai ai ,

1

(2)

i

Np

 pair    u ,

2

(3)

u

No

 open   ak .

3

(4)

k

ak

4

In the equations above,

5

and  u stands for the creation operator of a singlet pair (or geminal) of electrons in the uth pair,

6

which is formulated as below,

represents the creation operator of a spin-up electron in the kth orbital,



 u  C1u au1 au1  C2u au2 au2 .

7

(5) u

u

8

Here u1 (or u2) is the 1st (or 2nd) natural orbital (NO) in the uth pair, respectively, and C1 (or C2 )

9

is the corresponding configuration interaction (CI) coefficient. For each pair, two NOs can be

10

conceptually considered as bonding and anti-bonding orbitals, and the pair function (derived from

11

the action of  u

12

dissociation of a single covalent bond. All orbitals (core orbitals, orbitals in all pairs, and singly

13

occupied orbitals) should be orthogonal, meeting the strong orthogonality condition. The ground-

14

state GVB energy can be expressed as



on the vacuum state) has the advantage of qualitatively describing the

occ

15

occ occ

E   GVB Hˆ  GVB =  wi i hˆ i    xij J ij  yij K ij  i

i

(6)

j

16

where the orbital occupation coefficient wi and the two-electron coupling coefficients xij, yij are

17

defined in the previous work,52 depending on the wavefunction.

18

one-electron integral, while Jij and Kij are the standard MO-based two-electron Coulomb (Jij) and

19

exchange (Kij) integrals.

i hˆ i is the standard MO-based

20

Both the CI coefficients and molecular orbitals (MOs) need to be optimized in a GVB calculation.

21

Usually the optimization of the CI coefficients and MOs can be done separately, since their

22

couplings are not strong. For the CI coefficients, with a set of initial MOs, one only needs to solve

23

a series of coupled local CI equations iteratively. For the optimizations of the MOs, the approximate 5



Page 6 of 28

1

second order SCF method69,

2

multiconfigurational calculations. In this approach, a new set of orbitals is obtained by transforming

3

the orbital coefficient matrix by a unitary matrix represented as the exponential of an antisymmetric

4

matrix

70

has shown better performance than other methods for

Cnew  ColdU  Cold exp( A),

5

exp( A)  1  A 

(18)

1 2 A L . The algorithm implemented in the GAMESS program 2

6

where

7

approximates the exact (and quite expensive) orbital Hessian matrix by an approximate diagonal

8

one in the first iteration, and updates the Hessian using a variable metric algorithm during

9

subsequent iteration steps, and it works well in practice. It has been employed in our calculations.

10

B. Automatic Construction of the Initial GVB Orbitals.

11

The first issue in the construction of the initial guess GVB orbitals is to determine the number of

12

GVB pairs for a given molecule. The inner core orbitals play a little role in chemical or physical

13

processes and are doubly occupied orbitals in the GVB wavefunction. Thus, the number of the GVB

14

pairs (n) is n = min(a, b), in which a is the number of unoccupied molecular orbitals in a minimum

15

basis set (say, STO-6G), and b is the number of valence occupied orbitals.

16

17 18 19 20

Figure 1. The flowchart of constructing GVB initial guess. To generate automatically the initial guess GVB orbitals, two schemes (scheme I and II) are 6


Page 7 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

proposed to deal with two different cases. The whole procedure of constructing the GVB initial

2

guess is illustrated in Figure 1. For a given molecule, the first step is to perform a restricted HF

3

(RHF) calculation for the singlet state at the given basis set. This is followed by a triplet

4

(RHFUHF) instability test, which corresponds to a CIS (configuration interaction singles)

5

procedure based on the RHF reference. Instability is indicated when one or more eigenvalues of the

6

triplet instability matrix are negative.74 There are several types of RHF instability75, 76 but we check

7

only the triplet (RHFUHF) instability. This indicates whether a more stable real UHF

8

wavefunction exists. Čížek and Paldus77 mention the mathematical possibility of a more stable UHF

9

solution even if the lowest RHF wavefunction is stable, but we are not aware of an actual example.

10

If the RHF wavefunction is stable, scheme I based on the RHF orbitals will be adopted. Otherwise,

11

an energetically more stable UHF wavefunction is obtained, and scheme II based on the UHF

12

orbitals will be adopted. The detailed procedure of each scheme will be discussed separately in the

13

following subsections.

14

15

7



1 2 3 4

Figure 2. (a) The schematic diagram of scheme I; (b) the schematic diagram of scheme II.

5

As shown in Figure 2(a), scheme I contains the following three steps:

Page 8 of 28

6

(1) The number of the active virtual orbitals should be determined first. Since a GVB pair is

7

generally believed to contain a “bonding” orbital and an “anti-bonding” orbital, the active

8

virtual orbitals will be used to generate the corresponding “anti-bonding” orbitals for all GVB

9

pairs. A simple strategy is to require the active virtual orbitals to have maximum overlap with

10

the RHF virtual orbitals of the same system in a minimum basis set calculation (STO-6G in this

11

work). We then construct the overlap matrix between the canonical virtual orbitals

12

a , a  1, 2,L

13

b , b  1, 2,L n obtained using the minimum basis (m  n). By performing a singular value

14

decomposition (SVD) of this mn overlap matrix, we obtain two unitary matrices, U and V,

15

which can transform the canonical virtual RHF orbitals and the reference (minimum basis)

16

virtual orbitals to sets of new orbitals

, m (obtained using the working basis set) and the reference virtual orbitals

% and % a

m

17

% a =  kU ka k

separately:

b

n

% b =  lVlb

(7)

l 1

% have non-zero overlap with the reference virtual ones.

18

The first n transformed virtual orbitals

19

These are defined as active virtual orbitals; the remaining m-n transformed virtual orbitals are

a

8


Page 9 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

inactive virtual orbitals. Only the active virtual orbitals will be used in constructing the initial GVB

2

orbitals but the inactive virtual orbitals are also required for the GVB orbital optimization. The

3

procedure described above has been applied previously to obtain a few active virtual orbitals, which

4

were used in CCSD(T)-h calculations.78, 79

5

(2) Canonical occupied orbitals and active virtual orbitals are transformed into localized occupied

6

and virtual orbitals, separately. We use the Boys80 or Pipek-Mezey (PM)81 localization method.

7

(3) Generate the initial GVB orbitals by grouping a set of localized occupied orbitals and a set of

8

localized virtual orbitals into all GVB pairs automatically. Here the task of partitioning

9

localized occupied and virtual orbitals into different pairs is similar to the previous personnel-

10

assignment problem introduced by Kuhn82 and later developed by Munkres.83 An optimal

11

assignment may be achieved by maximizing the following function

W=  i | rˆ | a

12

i , aP

2

=   i | x | a  i , aP

2

 i | y | a

2

 i | z | a

2

 

(9)

13

Here the transition-dipole-like integral between an occupied orbital and a virtual orbital is used

14

as a measure of spatial vicinity between two orbitals in the Pth GVB pair. This index is

15

consistent with the idea that two orbitals (an occupied and a virtual) in a given GVB pair should

16

be spatially close to each other in a certain region. The detailed procedure of the assignment

17

process between occupied orbitals and virtual orbitals is shown in the supporting information

18

(SI).

19

In additional to orbitals in all GVB pairs, the initial GVB orbitals also contain doubly

20

occupied orbitals, and open-shell orbitals for open-shell species. These open-shell orbitals are

21

the same as in ROHF calculations, and no unitary transformation (or subsequent localization)

22

on these orbitals are required. For most systems, the computational cost of obtaining the initial

23

GVB guess orbitals is insignificant, compared with the full optimization of GVB orbitals.

24

As mentioned above, scheme II should be applied for systems in which the UHF

25

wavefunction is more stable than the (lowest energy) RHF wavefunction. Scheme II is similar

26

to scheme I, except that the first step is different. Here the first step is to transform two sets of

27

UHF spatial orbitals (spin-up and spin-down) into one orthonormal set of unrestricted charge

28

natural orbitals. At first, one constructs the overlap matrix D between spin-up and spin-down 9



  

  , 

1

occupied spatial orbitals

2

Decomposition of D to generate the corresponding orbitals

3

can be mathematically formulated as below,

and

i

i

Page 10 of 28

and then performs a Singular Value

% 

i

and

% . This process 

i

Dij  i  j , D = UdV †

4

n

n

i 1

j 1

    % % i   k U ki， j   l Vlj

5

 %  %  dij  dii ij i j

6

(10)

(11)

(12)

7

Here the matrix d is a diagonal matrix with the singular values of D as the diagonal matrix elements,

8

and the unitary matrices U and V are left and right eigenvectors of D. The corresponding orbitals

9

% 

i

(or

% ) are unitary transformations of the spin-up (or spin-down) occupied orbitals. 

i

10

 (1) % Then, for each pair of the corresponding orbitals ( % i and i ), we can obtain two UNOs ( i and

11

i(2) ), which are shown below,

12

i(1) 

1  % (% i  i ) 2(1  dii )

(13)

13

i(2) 

1  % (% i  i ) 2(1  dii )

(14)

i(1) and i(2) are 1  dii and 1  dii , respectively.84,

14

The corresponding occupation numbers for

15

85

16

denominator of equation (14) will approach to zero, meaning that the construction of

17

 % needed. This result indicates that two corresponding orbitals ( % i and i ) have almost identical

18

 spatial components, and only one UNO (for example, the corresponding % orbital) should be i

19

reserved. In fact, the corresponding UNO can be considered as an RHF-like doubly occupied orbital.

20

In addition, for a high-spin UHF wavefunction some singular values will be exactly zero. For

21

example, in the triplet UHF solution, two singular values will be zero. In such cases, two

However, if a singular value is very close to 1.0 (larger than 0.99999 in this work), the

i(2) is not

10


Page 11 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

 corresponding orbitals % i are considered to be two UNOs and correspond to two singly occupied

2

orbitals in the subsequent GVB calculation. According to the procedure described above, one may

3

transform two sets of UHF orbitals (spin-up and spin-down) into only one orthonormal set of UNOs.

4

The UNOs with large (larger than 0.99999) singular values will be doubly occupied in the initial

5

GVB wavefunction, and UNOs with zero singular values will be singly occupied high-spin GVB

6

orbitals. The remaining UNOs are divided into two subsets,

7



8

orbitals, and the latter subset corresponds to active virtual orbitals, as discussed in the first step of

9

scheme I. Two other steps in scheme II are the same as those in scheme I. A schematic diagram of

(2) i



(1) i

, i  1, 2,K 

and

, i  1, 2,K  . The former (more strongly occupied) subset corresponds to active occupied

10

scheme II is shown in Figure 2(b).

11

C. Computational Details. In the current implementation, the RHF and UHF calculations along

12

with the RHF stability check are performed with the Gaussian1686 package. Our own program is

13

used to generate the GVB guess orbitals for both scheme I (for the RHF wavefunction) and scheme

14

II (for the UHF wavefunction). The Boys or PM localization method is used in our program to get

15

localized MOs. With the GVB guess orbitals from our program, a slightly modified GAMESS

16

program (with the SOSCF method) is employed to carry out all GVB calculations in this study. For

17

some systems, a GVB-based block-correlated second-order perturbation theory (GVB-BCPT2)

18

developed by our group13 is used to get more accurate energies. For some systems, we also perform

19

CASSCF calculations with the Gaussian16 package. Besides, the PySCF87 program is applied to

20

perform electron correlation calculations with methods like CASSCF-NEVPT2.

21

III. RESULTS AND DISCUSSION

22

In this section, the GVB method will be applied to investigate a series of polyacenes (up to

23

decacene) and two isomers of the di-copper-oxygen-ammonia complex. For polyacenes, their

24

singlet and triplet states will be studied, and the singlet-triplet gaps calculated at the GVB and GVB-

25

BCPT2 levels will be reported. For two di-copper-oxygen-ammonia complexes, their relative

26

energies at the GVB and GVB-BCPT2 levels will be discussed. Moreover, we will show that the

27

GVB orbitals may be used to construct a reasonable minimum active space for CASSCF and 11



Page 12 of 28

1

CASSCF-NEVPT2 calculations. These more accurate energies may be used to benchmark the GVB

2

and GVB-BCPT2 energies.

3

A. Polyacenes

4

Polyacenes are a series of linear fused benzene oligomers. For longer polyacenes, a number of

5

experimental88-95 and theoretical96-109 studies have shown that their ground states are of increasingly

6

polyradical character. This feature makes polyacenes as good test systems for MR electronic

7

structure methods. Here the lowest singlet (ground-state) and the lowest triplet state of polyacenes

8

(from benzene to decacene) will be studied with GVB and GVB-BCPT2 methods, and the calculated

9

vertical and adiabatic singlet-triplet (S-T) gaps will be compared with those from other theoretical

10

methods or experimental data.88-90, 104, 105, 108 To the best of our knowledge, no GVB calculations for

11

systems with more than 50 GVB pairs have been reported before, except full-valence calculations

12

with a variant GVB method, CCVB-SD,107 for the polyacene with 17 benzene rings. In our

13

calculations, the largest polyacene is decacene, in which 192 electrons in 96 GVB pairs are

14

correlated. With our strategy, GVB calculations for all studied systems can be routinely done in a

15

black-box way.

16

The ground-state (singlet) geometries of various polyacenes were optimized at the UB3LYP/6-

17

31G(d) level. All GVB calculations and GVB-based calculations were carried out at the Cartesian

18

cc-pVDZ basis set. At the optimized geometry of each species, we find that the UHF wavefunction

19

is always more stable than the corresponding RHF wavefunction using the cc-pVDZ basis set.

20

Detailed electronic energies at the HF level can be found in the SI. For polyacenes, scheme II is

21

preferred because the UHF solution is more stable. In constructing the initial GVB orbitals, the Boys

22

localization method will be used, unless stated elsewhere.

23 24 25 26

Table 1. Singlet and triplet GVB energies (Eh) and vertical S-T gaps (kcal/mol) of polyacenes calculated with various theoretical methods. Size of polyacenes

GVB energies (S=0)a

GVB energies (S=1)a

GVB gapa

GVB gapb

GVBBCPT2 gapb

SFCCSD gapc

GASPDFT gapd

1

-230.919494

-230.740319

112.4

112.4

103.7

106.7

2

-383.706239

-383.556277

94.10

93.84

95.81

78.14

74.9

3

-536.480241

-536.352282

80.30

80.09

67.18

56.76

50.4

4

-689.252703

-689.133491

74.81

74.64

50.06

39.85

34.7 12


Page 13 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 2 3 4 5

5

-842.018728

-841.931731

54.59

54.70

17.62

28.62

25.0

6

-994.781710

-994.702437

49.74

49.89

3.31

7

-1147.537575

-1147.483855

33.71

10.6

8

-1300.294400

-1300.250415

27.60

5.4

9

-1453.051865

-1453.013985

23.77

4.4

10

-1605.812021

-1605.782004

18.84

3.7

17.3

aFull-valence

GVB calculations for the singlet and triplet states, respectively. GVB calculations for the singlet and triplet states, as discussed in the text. cData taken from Ref 105, using the 6-31+G(d) basis set. dData taken from Ref 108, using the 6-31+G(d,p) basis set. bSimplified

6

For polyacenes up to decacene, the lowest singlet and triplet energies calculated with GVB are

7

listed in Table 1. The GVB energies calculated with scheme II are shown here. Scheme I gives

8

converged GVB ground-state energies for these systems which differ from Scheme II by at most

9

1.68 mEh (see Table S2). These small errors can be traced back to the different components of GVB

10

orbitals obtained with two schemes because the Boys localization method (leading to - mixing

11

orbitals) was used. The GVB energies calculated with scheme II are usually somewhat lower than

12

those with scheme I. For the triplet state, the UNOs obtained from the corresponding UHF orbitals

13

are used to generate the GVB guess orbitals.

14

With all singlet and triplet GVB energies, we can predict the vertical S-T gaps, as also collected

15

in Table 1. As expected, the GVB S-T gaps decrease with increasing the system size. This tendency

16

is in qualitative agreement with other theoretical results from spin-flip CCSD (SF-CCSD)105 and

17

general active space pair density functional theory (GAS-PDFT).108 However, for the same system

18

the magnitude of the GVB S-T gap is much larger than the results from these other two methods.

19

For several small members of polyacenes, we have also carried out block-correlated GVB-BCPT2

20

calculations on both the singlet and triplet states. Since a GVB-BCPT2 calculation is relatively

21

expensive, we excluded the GVB pairs corresponding to C-H bonds, and reoptimized the

22

wavefunction to simplify the GVB calculation. This reduces the number of GVB pairs in benzene

23

from 15 in the full valence GVB calculations to 9 in the simplified calculation. Using the full valence

24

GVB orbitals as the initial guess, we have easily obtained the singlet and triplet energies of the

25

simplified GVB calculations. Table 1 shows that the vertical S-T gaps from full valence GVB and

26

simplified GVB calculations are very close. The deviations for all six compounds are all less than 13



Page 14 of 28

1

0.3 mEh. Therefore the GVB-BCPT2 S-T gaps based on simplified GVB wavefunctions are

2

expected to be good approximations to those based on full valence GVB wavefunctions. On the

3

basis of simplified GVB calculations, the calculated GVB-BCPT2 S-T gaps for systems up to

4

hexacene are also shown in Table 1. Obviously, the GVB-BCPT2 S-T gaps are much closer to the

5

results from two other methods than the GVB results for longer polyacenes (with three or more

6

benzene rings). This indicates that the inclusion of dynamic correlation from GVB-BCPT2

7

calculations significantly reduces the S-T gaps for longer polyacenes.

8

We also calculated the adiabatic GVB and GVB-BCPT2 S-T gaps for some polyacenes (from

9

naphthalene to hexacene). The corresponding triplet geometries are also optimized at the

10

UB3LYP/6-31G(d) level. The calculated adiabatic S-T gaps, together with those from other

11

theoretical methods104, 108 and the experimental data,88-90 are collected in Table 2. Again, one can

12

see that the S-T gaps from GVB-BCPT2 are significant improvements over the GVB results. The

13

dependence of the adiabatic S-T gaps upon the system size from GVB-BCPT2 is in qualitative

14

agreement with those from MR-AQCC104 and GAS-PDFT108 methods.

15 16 17 18

Table 2. Triplet GVB energies (Eh) and adiabatic S-T gaps (kcal/mol) of polyacenes from various theoretical methods and experiments Size of polyacenes

GVB energies (S=1)a

GVB gapa

GVB gapb

GVBBCPT2 gapb

MRCISD+Q gapc

MRAQCC gapc

GASPDFT gapd

Exp.e

2

-383.581840

78.1

77.5

85.3

65.5

67.6

64.7

60.9

3

-536.380968

62.3

61.9

60.2

48.4

50.0

43.1

42.6

4

-689.148205

65.6

65.4

44.1

38.5

36.9

28.8

29.4

5

-841.943790

47.0

47.0

23.5

27.7

24.7

20.5

19.8

6

-994.717194

40.5

40.6

2.7

24.2

17.3

15.0

12.4

19 20 21 22 23 24 25 26 27

aFull-valence

28

orbital C2  0.1 ) are shown in bold text. One calculation was started with Boys localized UHF

29 30

natural orbitals as initial GVB orbitals, and another with PM localized orbitals.

GVB calculations for the singlet and triplet states, respectively. GVB calculations for the singlet and triplet states, as discussed in the text. cData taken from Ref 104, using a -RAS/CAS(4,4)/AUX/6-31G reference space. dData taken from Ref 108, using the 6-31+G(d,p) basis set. eRef 88 for size of polyacenes n = 2-4, Ref 89 for n = 5, and Ref 90 for n = 6, respectively. bSimplified

Table 3. CI coefficients for the first 12 pairs in two singlet GVB calculations on hexacene. Both calculations include 60 GVB pairs, corresponding to a UNO occupation number threshold of 0.00001. The coefficients of strongly correlated pairs (pairs with the coefficient of the second natural u

No.

Boys Localization

PM Localization 14


Page 15 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


pair 1 2 3 4 5 6 7 8 9 10 11 12

1st NO 0.9836 0.9836 0.9836 0.9836 0.9958 0.9958 0.9958 0.9958 0.9958 0.9958 0.9958 0.9958

2nd NO -0.1801 -0.1801 -0.1801 -0.1801 -0.0920 -0.0920 -0.0919 -0.0919 -0.0918 -0.0917 -0.0916 -0.0916

1st NO 0.9843 0.9843 0.9843 0.9843 0.9883 0.9883 0.9883 0.9884 0.9915 0.9915 0.9915 0.9915

2nd NO -0.1763 -0.1763 -0.1763 -0.1763 -0.1526 -0.1526 -0.1526 -0.1526 -0.1301 -0.1301 -0.1301 -0.1301

1 2

3

4 5

Figure 3. (a) GVB orbitals in two GVB pairs for singlet hexacene. The indices to the left of the 15



Page 16 of 28

1 2 3 4 5 6 7 8

orbitals are the same as those in Table 3, i.e. pairs 2, 3, 4 are symmetry images of the pair 1, and pairs 5, 6, 7, 8, 9, 11, 12 are symmetry images of the pair 10. The initial GVB orbitals were obtained with scheme II and the Boys localization method; (b) GVB orbitals of three GVB pairs for singlet hexacene. Pairs 2, 3, 4 are symmetry images of the pair 1; pairs 6, 7, 8 are symmetry images of pair 5; and pairs 10, 11, 12 are symmetry images of pair 9. The initial GVB orbitals were obtained with scheme II and the PM localization method. CI coefficients for the first and second natural orbitals of the GVB pairs are included in parentheses.

9

Another important issue is to examine the effect of the localization method employed in

10

constructing the initial GVB orbitals on the obtained GVB energies and GVB orbitals. We will

11

take hexacene as an example to illustrate this point. Scheme II with two different localization

12

methods (Boys and PM) was used to build the initial GVB orbitals. The fully optimized GVB

13

energies for the ground state are -994.781710 and -994.775987 Eh, respectively. Thus, the GVB

14

energy obtained with the PM localization method is about 5.7 mEh higher than that obtained with

15

the Boy localization method. When the PM localized orbitals are used in constructing the initial

16

GVB orbitals, both scheme I and II lead to exactly the same GVB energies (see Table S3 for

17

hexacene). It is interesting to analyze the components of orbitals in all GVB pairs and the CI

18

coefficients for all GVB pairs. A GVB pair may be considered to be a strongly correlated pair (with

19

some diradical character) if the magnitude of the CI coefficient of the second NO is larger than a

20

given value (say 0.1). The CI coefficients for the first 12 of 60 pairs in two GVB calculations are

21

listed in Table 3, with coefficients of strongly correlated pairs in bold text. One can see that the

22

number of strongly correlated pairs ( C2  0.1 ) is dependent on the type of localized orbitals

23

used in the construction of the initial GVB guess. Using this definition, there are 4 strongly

24

correlated GVB pairs if the Boys localized orbitals are used, and 12 strongly correlated pairs if the

25

PM localized orbitals are employed. It appears that the GVB wavefunction obtained with the PM

26

localization method has stronger multi-reference character than that obtained with the Boys method.

27

Orbitals of the most strongly correlated GVB pairs in both cases are shown schematically in Figure

28

3. In both GVB solutions, two GVB orbitals in each strongly correlated pair are bonding and

29

antibonding  orbitals, respectively. It is well known that the localized orbitals from the PM

30

localization method are either  or  orbitals, but the localized orbitals from the Boys method can

31

be of mixed - character.81 Our calculations show that with PM localized orbitals, the GVB initial

32

orbitals are either  or  orbitals, and the resulting GVB orbitals are also of pure  or  character,

u

16


Page 17 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

indicating no mixing of  and  orbitals in the orbital optimization. However, with Boys localized

2

orbitals, the initial GVB orbitals have always mixed - character, and the resulting GVB orbitals

3

can have either pure  or , or mixed - character. This result indicates that the correlation

4

between mixed - initial GVB orbitals leads to three types of GVB orbitals in the orbital

5

optimization. In general, the Boys localization method is recommended since slightly lower GVB

6

energies will be obtained with the corresponding initial GVB orbitals.

7

It can also be seen from Table 3 that, for hexacene, the GVB wavefunction keeps the symmetry

8

of the molecule (D2h). For example, the orbitals in pairs 2, 3, 4 are symmetrical images of those in

9

GVB pair 1 (shown in Figure 3(a)). For longer polyacenes, the GVB wavefunction exhibits slight

10

symmetry breaking (see Figure S1(a) for heptacene). Symmetry breaking of GVB pairs was also

11

observed in Co(CH3)2+ by Perry et al.110

12

B. Two Isomers of the Di-copper-oxygen-ammonia Complex

13

Transition metal complexes containing [Cu2O2]2+ have high catalytic activity in biochemistry.

14

However, several low-lying isomers may co-exist, depending on the ligands coordinated to two

15

copper atoms. Previous studies have demonstrated that the electronic structures of these systems

16

have strong multi-reference characters, and the evaluation of the relative energies among different

17

isomeric species is quite challenging for many theoretical methods.15, 28, 111-114 Here we are interested

18

in the relative energy between two isomers of a di-copper-oxygen-ammonia complex, labeled as bis

19

and per, as shown in Figure 4. The relative energies of these two isomers have been investigated by

20

many theoretical methods.15, 28 Here we will first investigate whether the use of the initial GVB

21

orbitals constructed with our procedure can lead to the converged GVB energies for two isomers.

22

Then, the relative energies of two isomers calculated with GVB and GVB-BCPT2 will be compared

23

with other theoretical methods. In addition, we will show that the GVB orbitals may be used to

24

automatically select the minimum active space for other multi-reference methods. The geometries

25

of two species were taken from Ref 111. The Cartesian cc-pVTZ basis set is employed for Cu and

26

O atoms, and cc-pVDZ for N and H atoms. Relativistic effects are not considered, although there is

27

indication that they may play an important role in these two species.113, 114

28 17



1 2 3 4

Page 18 of 28

Figure 4. Molecular structures (bis and per) for two isomers of {[Cu(NH3)]2O2}2+

5 6 7 8 9 10

Figure 5. (a) GVB orbitals of 2 strongly correlated pairs in the bis structure; (b) GVB orbitals of 2 strongly correlated pairs in the per structure.

11 12 13 14 15

Figure 6. CASSCF(4,4) natural orbitals and corresponding natural orbital occupation number (NOON) of 2 strongly correlated pairs in the bis structure (a) and the per structure (b)

16

For the studied compounds, the UHF solution is more stable than the RHF solution, so scheme

17

II was used. Two UHF solutions exist for these species; the lowest UHF solution was employed for

18

constructing the initial GVB orbitals. In the per structure, two copper atoms are found to be

19

antiferromagnetically coupled in the lowest UHF solution. Detailed electronic energies at the HF

20

level can be found in the SI. The Boys localization method is used in constructing the initial GVB

21

orbitals. With the initial GVB orbitals, the SOSCF algorithm in the GAMESS program leads to the 18


Page 19 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

converged GVB(16) energies -3539.779994 Eh (for bis) and -3539.800928 Eh (for per), respectively.

2

We also conducted GVB-BCPT2 calculations on these two isomers to obtain the corresponding

3

relative energy. The calculated energy differences, E = E(bis)–E(per), obtained with GVB, GVB-

4

BCPT2, and several other methods are listed in Table 4. The most accurate method applied

5

previously is probably the (24,24) (24 electrons in 24 orbitals) DMRG-CASPT2 study of Pierlot

6

and coworkers.15 Our GVB-BCPT2 relative energy (21.6 kcal/mol) is quite close to the DMRG-

7

CASPT2 value (22.6 kcal/mol), being much better than the GVB value. This shows that the

8

inclusion of the post-GVB dynamic correlation is essential for estimating the relative energies of

9

these two isomers.

10

By analyzing the components of the GVB orbitals, the electronic structure of these two isomers

11

can be understood. For both isomers, GVB calculations show that there are two strongly correlated

12

pairs. Four GVB orbitals in these two pairs are illustrated in Figure 5. For the bis structure, two

13

GVB pairs correspond to two Cu-O bonds mainly formed by one 3d (on Cu) and one 2p (on O)

14

orbitals. In the per structure, one GVB pair corresponds to the O-O σ bond and the other pair mainly

15

consists of two weakly coupled 3d orbitals on two Cu atoms. By analyzing the CI coefficients for

16

two GVB pairs of each species, we may conclude that the bis isomer has a strong tetraradical

17

character, and the per isomer has a strong diradical character. On the basis of two strongly correlated

18

GVB pairs for each structure, one may consider a minimum active space to consist of four electrons

19

in four GVB orbitals in these two GVB pairs. Using these active GVB orbitals as the initial orbitals,

20

we have carried out CASSCF(4,4) and CASSCF(4,4)-NEVPT2 calculations. The relative energies

21

obtained at CASSCF and CASSCF-NEVPT2 levels are also listed in Table 4. The CASSCF(4,4)-

22

NEVPT2 energy difference is 24.4 kcal/mol, being close to the RASPT2(24, 28)28 result (25.2

23

kcal/mol) and the DMRG(24,24)-CASPT215 value (22.6 kcal/mol) reported previously. This result

24

suggests that the minimum active space, which can be automatically derived from GVB calculations,

25

may be suitable for CASSCF and CASSCF-based multi-reference electron correlation calculations.

26

To check whether the GVB orbitals in two strongly correlated pairs are good approximations to

27

active orbitals in CASSCF(4,4) calculations, we have shown four natural orbitals from CASSCF(4,4)

28

calculations in Figure 6. For the per structure, its four natural orbitals are almost identical to the 19



Page 20 of 28

1

corresponding GVB orbitals. For the bis isomer, two natural orbitals with larger (smaller)

2

occupation numbers may be considered as a linear combination of the first (second) natural orbitals

3

of two GVB pairs. Thus, the GVB orbitals are very good initial guess orbitals for active orbitals in

4

CASSCF calculations.

5 6 7 8

Table 4. Relative energies (in kcal/mol) between two isomers of {[Cu(NH3)]2O2}2+calculated with different methods.

9 10 11 12 13 14

Method E = E(bis) – E(per) GVB(16) 13.1 GVB(16)-BCPT2 21.6 CASSCF(4,4) 8.5 CASSCF(4,4)-NEVPT2 24.4 a RASPT2(24, 28) 25.2 DMRG(24,24)b 19.5 b DMRG(24,24)-CASPT2 22.6 a: Results taken from Ref 28. b: Results taken from Ref 15. The number of renormalized states M=1500. The Douglas-Kroll-Hess (DKH) Hamiltonian and the extensive relativistic atomic natural orbital (ANO-RCC) type basis sets are employed.

IV. CONCLUSIONS

15

In this paper, we propose an efficient strategy of automatically building the initial orbitals to allow

16

calculations of generalized valence bond (GVB) wavefunctions feasible for large systems. Two

17

different schemes are proposed to deal with different systems. The first scheme (scheme I) is

18

suitable for systems in which the RHF wavefunction is stable, while the second scheme (scheme II)

19

is for systems where the UHF wavefunction is more stable than the RHF one. In scheme I, the active

20

occupied orbitals consist of the valence orbitals of molecules, while the active virtual orbitals are

21

constructed from the original virtual space by requiring maximum overlap with the virtual orbitals

22

of the same system in an auxiliary minimum basis set. Scheme II differs from scheme I, in that the

23

active occupied and active virtual orbitals are obtained from the transformation of two sets of UHF

24

spatial orbitals into one set of Unrestricted Natural Orbitals. Next, these active occupied orbitals and

25

active virtual orbitals are separately transformed into localized orbitals using the Boys or PM

26

localization method. The GVB initial guess is then obtained by grouping active occupied and virtual

27

localized orbitals into different pairs using the Kuhn-Munkres algorithm. With the initial GVB guess

28

constructed above and the existing SOSCF algorithm for optimizing GVB orbitals, efficient black20


Page 21 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

box GVB calculations for large systems have been achieved. Illustrative calculations have shown

2

that GVB energies can be easily obtained for the lowest singlet and triplet states of polyacenes (up

3

to decacene) and the singlet ground state of two di-copper-oxygen-ammonia complexes. We have

4

also calculated the singlet-triplet gaps for some polyacenes and the relative energy between two di-

5

copper-oxygen-ammonia complexes at the GVB-BCPT2 levels. In addition, we have shown that

6

analysis of the GVB orbitals in strongly correlated pairs may be used to achieve an automatic

7

definition of the minimum active space for multi-reference electronic structure calculations like

8

CASSCF and DMRG.

9

The present strategy of building the initial GVB guess has been tested on a number of molecules

10

in the singlet and triplet states. More chemical systems with different spin states will be tested in

11

the near future. We are implementing the block-correlated third-order perturbation theory based on

12

the GVB reference (GVB-BCPT3) to get more accurate energies. Work along these two directions

13

is in progress.

14 15

ASSOCIATED CONTENT

16

Supporting Information

17

Details about Kuhn-Munkres (KM) algorithm. Cartesian coordinates and UB3LYP, UHF, GVB

18

energies, part CI coefficients and orbitals of singlets and triplet polyacenes. Cartesian coordinates

19

and Singlet energies of two forms of {[Cu(NH3)]2O2}2+.

20 21

AUTHOR INFORMATION

22 23 24 25 26 27 28 29 30 31 32

Corresponding Author *E-mail: [email protected]. *E-mail: [email protected]. ORCID Qingchun Wang: 0000-0003-0966-9782 Peter Pulay: 0000-0003-4410-4373 Shuhua Li: 0000-0001-6756-057X

33

ACKNOWLEDGMENTS

Notes: The authors declare no competing financial interest.

21



Page 22 of 28

1 2 3 4 5 6 7 8 9 10

This work was supported by the National Natural Science Foundation of China (Grants No. 21333004, NO. 21673110, NO. 21833002), the Natural Science Foundation of JiangSu Province (Grant No. SBK2015041570), and the “Fundamental Research Funds for the Central Universities”. The numerical calculations in this paper have been done on the computing facilities in the High Performance Computing Center (HPCC) of Nanjing University. It was also supported by the U. S. National Science Foundation (NSF), Grant No. DMR-1609650 and the U. S. National Institutes of Health (NIH), Grant No. NIGMS 1R01GM120578, and the Arkansas Biosciences Institute, the major research component of the Arkansas Tobacco Settlement Proceeds Act of 2000.

11

REFERENCES

12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54

(1) Andrade, X.; Strubbe, D.; De Giovannini, U.; Larsen, A. H.; Oliveira, M. J. T.; Alberdi-Rodriguez, J.; Varas, A.; Theophilou, I.; Helbig, N.; Verstraete, M. J.; Stella, L.; Nogueira, F.; Aspuru-Guzik, A.; Castro, A.; Marques, M. A. L.; Rubio, A., Real-space grids and the Octopus code as tools for the development of new simulation approaches for electronic systems. Phys. Chem. Chem. Phys. 2015, 17, 31371-31396. (2) Grimme, S.; Schreiner, P. R., Computational Chemistry: The Fate of Current Methods and Future Challenges. Angew. Chem. Int. Ed. 2018, 57, 4170-4176. (3) Møller, C.; Plesset, M. S., Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618-622. (4) Purvis, G. D.; Bartlett, R. J., A full coupled ‐ cluster singles and doubles model: The inclusion of disconnected triples. J. Chem. Phys. 1982, 76, 1910-1918. (5) Bartlett, R. J., Coupled-cluster approach to molecular structure and spectra: a step toward predictive quantum chemistry. J. Phys. Chem. 1989, 93, 1697-1708. (6) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M., A fifth-order perturbation comparison of electron correlation theories. Chem. Phys. Lett. 1989, 157, 479-483. (7) Lee, T. J.; Scuseria, G. E., Achieving Chemical Accuracy with Coupled-Cluster Theory. In Quantum Mechanical Electronic Structure Calculations with Chemical Accuracy, Langhoff, S. R., Ed. Springer Netherlands: Dordrecht, 1995; pp 47-108. (8) Cremer, D., Møller–Plesset perturbation theory: from small molecule methods to methods for thousands of atoms. WIREs Comp. Mol. Sci. 2011, 1, 509-530. (9) Řezáč, J.; Šimová, L.; Hobza, P., CCSD[T] Describes Noncovalent Interactions Better than the CCSD(T), CCSD(TQ), and CCSDT Methods. J. Chem. Theory Comput. 2013, 9, 364-369. (10) Bartlett, R. J.; Musiał, M., Coupled-cluster theory in quantum chemistry. Rev. Mod. Phys. 2007, 79, 291-352. (11) Tarumi, M.; Kobayashi, M.; Nakai, H., Generalized Møller–Plesset Multiconfiguration Perturbation Theory Applied to an Open-Shell Antisymmetric Product of Strongly Orthogonal Geminals Reference Wave Function. J. Chem. Theory Comput. 2012, 8, 4330-4335. (12) Li Manni, G.; Ma, D.; Aquilante, F.; Olsen, J.; Gagliardi, L., SplitGAS Method for Strong Correlation and the Challenging Case of Cr2. J. Chem. Theory Comput. 2013, 9, 3375-3384. (13) Xu, E.; Li, S., Block correlated second order perturbation theory with a generalized valence bond reference function. J. Chem. Phys. 2013, 139, 174111. (14) Yang, K. R.; Jalan, A.; Green, W. H.; Truhlar, D. G., Which Ab Initio Wave Function Methods Are Adequate for Quantitative Calculations of the Energies of Biradicals? The Performance of CoupledCluster and Multi-Reference Methods Along a Single-Bond Dissociation Coordinate. J. Chem. Theory Comput. 2013, 9, 418-431. (15) Phung, Q. M.; Wouters, S.; Pierloot, K., Cumulant Approximated Second-Order Perturbation Theory Based on the Density Matrix Renormalization Group for Transition Metal Complexes: A Benchmark Study. J. Chem. Theory Comput. 2016, 12, 4352-4361. (16) Tsuchimochi, T.; Ten-no, S., Black-Box Description of Electron Correlation with the Spin-Extended Configuration Interaction Model: Implementation and Assessment. J. Chem. Theory Comput. 2016, 12, 1741-1759. (17) Roos, B. O., The Complete Active Space Self ‐ Consistent Field Method and its Applications in Electronic Structure Calculations. In Adv. Chem. Phys., 1987; pp 399-445. (18) Rogers, D. M.; Wells, C.; Joseph, M.; Boddington, V. J.; McDouall, J. J. W., On the choice of active 22


Page 23 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56


space orbitals in MCSCF calculations. J. Mol. Struct.: THEOCHEM 1998, 434, 239-245. (19) Sayfutyarova, E. R.; Sun, Q.; Chan, G. K.-L.; Knizia, G., Automated Construction of Molecular Active Spaces from Atomic Valence Orbitals. J. Chem. Theory Comput. 2017, 13, 4063-4078. (20) Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K., Second-order perturbation theory with a CASSCF reference function. J. Phys. Chem. 1990, 94, 5483-5488. (21) Andersson, K.; Malmqvist, P. Å.; Roos, B. O., Second‐order perturbation theory with a complete active space self‐consistent field reference function. J. Chem. Phys. 1992, 96, 1218-1226. (22) Werner, H.-J., Third-order multireference perturbation theory The CASPT3 method. Mol. Phys. 1996, 89, 645-661. (23) Pierloot, K., The CASPT2 method in inorganic electronic spectroscopy: from ionic transition metal to covalent actinide complexes∗. Mol. Phys. 2003, 101, 2083-2094. (24) Fang, T.; Li, S., Block correlated coupled cluster theory with a complete active-space self-consistentfield reference function: The formulation and test applications for single bond breaking. J. Chem. Phys. 2007, 127, 204108. (25) Vancoillie, S.; Malmqvist, P. Å.; Veryazov, V., Potential Energy Surface of the Chromium Dimer Rere-revisited with Multiconfigurational Perturbation Theory. J. Chem. Theory Comput. 2016, 12, 16471655. (26) Aquilante, F.; Autschbach, J.; Carlson, R. K.; Chibotaru, L. F.; Delcey, M. e. G.; Vico, L. D.; Galván, I. F.; Ferré, N.; Frutos, L. M.; Gagliardi, L.; Garavelli, M.; Giussani, A.; Hoyer, C. E.; Manni, G. L.; Lischka, H.; Ma, D.; Malmqvist, P. Å.; M€uller, T.; Nenov, A.; Olivucci, M.; Pedersen, T. B.; Peng, D.; Plasser, F.; Pritchard, B.; Reiher, M.; Rivalta, I.; Schapiro, I.; Javier, S.-M.; Stenrup, M.; Truhlar, D. G.; Ungur, L.; Valentini, A.; Vancoillie, S.; Veryazov, V.; Vysotskiy, V. P.; Weingart, O.; Zapata, F.; Lindh, R., Molcas 8: New capabilities for multiconfigurational quantum chemical calculations across the periodic table. J. Comput. Chem. 2016, 37, 506-541. (27) Olsen, J.; Roos, B. O.; Jo/rgensen, P.; Jensen, H. J. r. A., Determinant based configuration interaction algorithms for complete and restricted configuration interaction spaces. J. Chem. Phys. 1988, 89, 21852192. (28) Malmqvist, P. Å.; Pierloot, K.; Shahi, A. R. M.; Cramer, C. J.; Gagliardi, L., The restricted active space followed by second-order perturbation theory method: Theory and application to the study of CuO2 and Cu2O2 systems. J. Chem. Phys. 2008, 128, 204109. (29) Chan, G. K.-L.; Sharma, S., The Density Matrix Renormalization Group in Quantum Chemistry. Annu. Rev. Phys. Chem. 2011, 62, 465-481. (30) Takeshi, Y.; Yuki, K.; Wataru, M.; Jakub, C.; Nguyen, L. T.; Masaaki, S., Density matrix renormalization group for ab initio Calculations and associated dynamic correlation methods: A review of theory and applications. Int. J. Quantum Chem. 2015, 115, 283-299. (31) Szilárd, S.; Max, P.; Valentin, M.; Gergely, B.; Frank, V.; Reinhold, S.; Örs, L., Tensor product methods and entanglement optimization for ab initio quantum chemistry. Int. J. Quantum Chem. 2015, 115, 1342-1391. (32) Pulay, P.; Hamilton, T. P., UHF natural orbitals for defining and starting MC‐SCF calculations. J. Chem. Phys. 1988, 88, 4926-4933. (33) Anglada, J. M.; Bofill, J. M., Practical remarks on the selection of the active space in the CAS-SCF wavefunction. Chem. Phys. Lett. 1995, 243, 151-157. (34) Keller, S.; Boguslawski, K.; Janowski, T.; Reiher, M.; Pulay, P., Selection of active spaces for multiconfigurational wavefunctions. J. Chem. Phys. 2015, 142, 244104. (35) Stein, C. J.; Reiher, M., Automated Selection of Active Orbital Spaces. J. Chem. Theory Comput. 2016, 12, 1760-1771. (36) Dutta, A. K.; Nooijen, M.; Neese, F.; Izsák, R., Automatic active space selection for the similarity transformed equations of motion coupled cluster method. J. Chem. Phys. 2017, 146, 074103. (37) Alexander, B. C.; Andreas, H.; Stefan, G., The Fractional Occupation Number Weighted Density as a Versatile Analysis Tool for Molecules with a Complicated Electronic Structure. Chem. Eur. J. 2017, 23, 6150-6164. (38) Bao, J. J.; Dong, S. S.; Gagliardi, L.; Truhlar, D. G., Automatic Selection of an Active Space for Calculating Electronic Excitation Spectra by MS-CASPT2 or MC-PDFT. J. Chem. Theory Comput. 2018, 14, 2017-2025. (39) Yamaguchi, K.; Ohta, K.; Yabushita, S.; Fueno, T., DODS natural orbital (NO) CI investigations of 1,3 ‐diradicals: CH2NHO, CH2OO, and CH2CH2O. J. Chem. Phys. 1978, 68, 4323-4325. 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56

Page 24 of 28

(40) Kizashi, Y., Ab initio unrestricted Hartree–Fock (UHF) and UHF–natural orbital CI studies of ozone. Int. J. Quantum Chem. 1980, 18, 101-106. (41) Yamaguchi, K.; Yabushita, S.; Fueno, T.; Kato, S.; Morokuma, K.; Iwata, S., Ab initio UHF and UHF NO CI approaches for quasi-degenerate systems: methylene peroxide (CH2OO). Chem. Phys. Lett. 1980, 71, 563-568. (42) Bofill, J. M.; Pulay, P., The unrestricted natural orbital–complete active space (UNO–CAS) method: An inexpensive alternative to the complete active space–self‐consistent‐field (CAS–SCF) method. J. Chem. Phys. 1989, 90, 3637-3646. (43) Nishihara, S.; Saito, T.; Yamanaka, S.; Kitagawa, Y.; Kawakami, T.; Okumura, M.; Yamaguchi, K., MkMRCC, APUCC and APUBD approaches to 1,n-didehydropolyene diradicals: the nature of throughbond exchange interactions. Mol. Phys. 2010, 108, 2559-2578. (44) Saito, T.; Nishihara, S.; Yamanaka, S.; Kitagawa, Y.; Kawakami, T.; Yamada, S.; Isobe, H.; Okumura, M.; Yamaguchi, K., Singlet–triplet energy gap for trimethylenemethane, oxyallyl diradical, and related species: single- and multireference computational results. Theor. Chem. Acc. 2011, 130, 739-748. (45) Kawakami, T.; Saito, T.; Sharma, S.; Yamanaka, S.; Yamada, S.; Nakajima, T.; Okumura, M.; Yamaguchi, K., Full-valence density matrix renormalisation group calculations on meta-benzyne based on unrestricted natural orbitals. Revisit of seamless continuation from broken-symmetry to symmetryadapted models for diradicals. Mol. Phys. 2017, 115, 2267-2284. (46) Kawakami, T.; Sano, S.; Saito, T.; Sharma, S.; Shoji, M.; Yamada, S.; Takano, Y.; Yamanaka, S.; Okumura, M.; Nakajima, T.; Yamaguchi, K., UNO DMRG CASCI calculations of effective exchange integrals for m-phenylene-bis-methylene spin clusters. Mol. Phys. 2017, 115, 2154-2167. (47) Parks, J. M.; Parr, R. G., Theory of Separated Electron Pairs. J. Chem. Phys. 1958, 28, 335-345. (48) Kutzelnigg, W., Direct Determination of Natural Orbitals and Natural Expansion Coefficients of Many‐Electron Wavefunctions. I. Natural Orbitals in the Geminal Product Approximation. J. Chem. Phys. 1964, 40, 3640-3647. (49) Miller, K. J.; Ruedenberg, K., Electron Correlation and Separated ‐ Pair Approximation. An Application to Berylliumlike Atomic Systems. J. Chem. Phys. 1968, 48, 3414-3443. (50) Surján, P. R., An Introduction to the Theory of Geminals. In Correlation and Localization, Surján, P. R.; Bartlett, R. J.; Bogár, F.; Cooper, D. L.; Kirtman, B.; Klopper, W.; Kutzelnigg, W.; March, N. H.; Mezey, P. G.; Müller, H.; Noga, J.; Paldus, J.; Pipek, J.; Raimondi, M.; Røeggen, I.; Sun, J. Q.; Surján, P. R.; Valdemoro, C.; Vogtner, S., Eds. Springer Berlin Heidelberg: Berlin, Heidelberg, 1999; pp 63-88. (51) Surján, P. R.; Szabados, Á.; Jeszenszki, P.; Zoboki, T., Strongly orthogonal geminals: size-extensive and variational reference states. J. Math. Chem. 2012, 50, 534-551. (52) Bobrowicz, F. W.; Goddard, W. A., The Self-Consistent Field Equations for Generalized Valence Bond and Open-Shell Hartree—Fock Wave Functions. In Methods of Electronic Structure Theory, Schaefer, H. F., Ed. Springer US: Boston, MA, 1977; pp 79-127. (53) de Albuquerque Lins, J. O. M.; Nascimento, M. A. C., Theoretical investigation of the methane activation reaction on protonated zeolite from generalized valence-bond plus configuration interaction calculations. J. Mol. Struct.: THEOCHEM 1996, 371, 237-243. (54) Cullen, J., Generalized valence bond solutions from a constrained coupled cluster method. Chem. Phys. 1996, 202, 217-229. (55) Murphy, R. B.; Pollard, W. T.; Friesner, R. A., Pseudospectral localized generalized Mo/ller–Plesset methods with a generalized valence bond reference wave function: Theory and calculation of conformational energies. J. Chem. Phys. 1997, 106, 5073-5084. (56) Francesco, F.; A., G. W., GVB–RP: A reliable MCSCF wave function for large systems. Int. J. Quantum Chem. 1999, 73, 1-22. (57) Friesner, R. A.; Murphy, R. B.; Beachy, M. D.; Ringnalda, M. N.; Pollard, W. T.; Dunietz, B. D.; Cao, Y., Correlated ab Initio Electronic Structure Calculations for Large Molecules. J. Phys. Chem. A 1999, 103, 1913-1928. (58) Sejpal, M.; Messmer, R. P., Calculations using generalized valence bond based Møller–Plesset perturbation theory. J. Chem. Phys. 2001, 114, 4796-4804. (59) Rosta, E.; Surján, P. R., Two-body zeroth order Hamiltonians in multireference perturbation theory: The APSG reference state. J. Chem. Phys. 2002, 116, 878-890. (60) Anderson, A. G.; III, W. A. G., Generalized valence bond wave functions in quantum Monte Carlo. J. Chem. Phys. 2010, 132, 164110. (61) Small, D. W.; Head-Gordon, M., Post-modern valence bond theory for strongly correlated electron 24


Page 25 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56


spins. Phys. Chem. Chem. Phys. 2011, 13, 19285-19297. (62) Small, D. W.; Head-Gordon, M., A fusion of the closed-shell coupled cluster singles and doubles method and valence-bond theory for bond breaking. J. Chem. Phys. 2012, 137, 114103. (63) Fracchia, F.; Filippi, C.; Amovilli, C., Size-Extensive Wave Functions for Quantum Monte Carlo: A Linear Scaling Generalized Valence Bond Approach. J. Chem. Theory Comput. 2012, 8, 1943-1951. (64) Chatterjee, K.; Pernal, K., Excitation energies from time-dependent generalized valence bond method. Theor. Chem. Acc. 2015, 134, 118. (65) Chatterjee, K.; Pastorczak, E.; Jawulski, K.; Pernal, K., A minimalistic approach to static and dynamic electron correlations: Amending generalized valence bond method with extended random phase approximation correlation correction. J. Chem. Phys. 2016, 144, 244111. (66) Filatov, M.; Martinez, T. J.; Kim, K. S., Using the GVB Ansatz to develop ensemble DFT method for describing multiple strongly correlated electron pairs. Phys. Chem. Chem. Phys. 2016, 18, 21040-21050. (67) Xu, E.; Zhao, D.; Li, S., Multireference Second Order Perturbation Theory with a Simplified Treatment of Dynamical Correlation. J. Chem. Theory Comput. 2015, 11, 4634-4643. (68) Hunt, W. J.; Dunning, T. H.; Goddard, W. A., The orthogonality constrained basis set expansion method for treating off-diagonal lagrange multipliers in calculations of electronic wave functions. Chem. Phys. Lett. 1969, 3, 606-610. (69) Fischer, T. H.; Almlof, J., General methods for geometry and wave function optimization. J. Phys. Chem. 1992, 96, 9768-9774. (70) Chaban, G.; Schmidt, M. W.; Gordon, M. S., Approximate second order method for orbital optimization of SCF and MCSCF wavefunctions. Theor. Chem. Acc. 1997, 97, 88-95. (71) Langlois, J.-M.; Yamasaki, T.; Muller, R. P.; Goddard, W. A., III, Rule-Based Trial Wave Functions for Generalized Valence Bond Theory. J. Phys. Chem. 1994, 98, 13498-13505. (72) Sano, T., Elementary Jacobi rotation method for generalized valence bond perfect-pairing calculations combined with simple procedure for generating reliable initial orbitals. J. Mol. Struct.: THEOCHEM 2000, 528, 177-191. (73) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Jr, J. A. M., General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347-1363. (74) Seeger, R.; Pople, J. A., Self‐consistent molecular orbital methods. XVIII. Constraints and stability in Hartree–Fock theory. J. Chem. Phys. 1977, 66, 3045-3050. (75) Paldus, J.; Čĺžzek, J., Stability Conditions for the Solutions of the Hartree–Fock Equations for Atomic and Molecular Systems. II. Simple Open‐Shell Case. J. Chem. Phys. 1970, 52, 2919-2936. (76) Schlegel, H. B.; McDouall, J. J. W., Do You Have SCF Stability and Convergence Problems? In Computational Advances in Organic Chemistry: Molecular Structure and Reactivity, Ögretir, C.; Csizmadia, I. G., Eds. Springer Netherlands: Dordrecht, 1991; pp 167-185. (77) Čížek, J.; Paldus, J., Stability Conditions for the Solutions of the Hartree—Fock Equations for Atomic and Molecular Systems. Application to the Pi‐Electron Model of Cyclic Polyenes. J. Chem. Phys. 1967, 47, 3976-3985. (78) Shen, J.; Kou, Z.; Xu, E.; Li, S., The coupled cluster singles, doubles, and a hybrid treatment of connected triples based on the split virtual orbitals. J. Chem. Phys. 2012, 136, 044101. (79) Kou, Z.; Shen, J.; Xu, E.; Li, S., Hybrid Coupled Cluster Methods Based on the Split Virtual Orbitals: Barrier Heights of Reactions and Spectroscopic Constants of Open-Shell Diatomic Molecules. J. Phys. Chem. A 2013, 117, 626-632. (80) Foster, J. M.; Boys, S. F., Canonical Configurational Interaction Procedure. Rev. Mod. Phys. 1960, 32, 300-302. (81) Pipek, J.; Mezey, P. G., A fast intrinsic localization procedure applicable for ab initio and semiempirical linear combination of atomic orbital wave functions. J. Chem. Phys. 1989, 90, 4916-4926. (82) Kuhn, H. W., The Hungarian method for the assignment problem. Naval Research Logistics Quarterly 1955, 2, 83-97. (83) Munkres, J., Algorithms for the Assignment and Transportation Problems. J. Soc. Ind. Appl. Math. 1957, 5, 32-38. (84) King, H. F.; Stanton, R. E.; Kim, H.; Wyatt, R. E.; Parr, R. G., Corresponding Orbitals and the Nonorthogonality Problem in Molecular Quantum Mechanics. J. Chem. Phys. 1967, 47, 1936-1941. (85) Xu, E.; Shen, J.; Kou, Z.; Li, S., Coupled cluster with singles, doubles, and partial higher-order excitations based on the corresponding orbitals: The formulation and test applications for bond 25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56

Page 26 of 28

breaking processes. J. Chem. Phys. 2010, 132, 134110. (86) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A. V.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J. J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Keith, T. A.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian 16 Rev. A.03, Wallingford, CT, 2016. (87) Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S.; Wouters, S.; Chan, G. K.-L., PySCF: the Python-based simulations of chemistry framework. WIREs Comp. Mol. Sci. 2018, 8, e1340. (88) Siebrand, W., Radiationless Transitions in Polyatomic Molecules. II. Triplet ‐ Ground ‐ State Transitions in Aromatic Hydrocarbons. J. Chem. Phys. 1967, 47, 2411-2422. (89) Burgos, J.; Pope, M.; Swenberg, C. E.; Alfano, R. R., Heterofission in pentacene-doped tetracene single crystals. physica status solidi (b) 1977, 83, 249-256. (90) Angliker, H.; Rommel, E.; Wirz, J., Electronic spectra of hexacene in solution (ground state. Triplet state. Dication and dianion). Chem. Phys. Lett. 1982, 87, 208-212. (91) Kamada, K.; Ohta, K.; Shimizu, A.; Kubo, T.; Kishi, R.; Takahashi, H.; Botek, E.; Champagne, B.; Nakano, M., Singlet Diradical Character from Experiment. J. Phys. Chem. Lett. 2010, 1, 937-940. (92) Balaji, P.; Matthew, B.; R., P. S.; Anne ‐ Frances, M.; E., A. J., Synthesis and Structural Characterization of Crystalline Nonacenes. Angew. Chem. Int. Ed. 2011, 50, 7013-7017. (93) Watanabe, M.; Chang, Y. J.; Liu, S.-W.; Chao, T.-H.; Goto, K.; Islam, M. M.; Yuan, C.-H.; Tao, Y.-T.; Shinmyozu, T.; Chow, T. J., The synthesis, crystal structure and charge-transport properties of hexacene. Nat. Chem. 2012, 4, 574. (94) Huang, R.; Phan, H.; Herng, T. S.; Hu, P.; Zeng, W.; Dong, S.-q.; Das, S.; Shen, Y.; Ding, J.; Casanova, D.; Wu, J., Higher Order π-Conjugated Polycyclic Hydrocarbons with Open-Shell Singlet Ground State: Nonazethrene versus Nonacene. J. Am. Chem. Soc. 2016, 138, 10323-10330. (95) Zuzak, R.; Dorel, R.; Krawiec, M.; Such, B.; Kolmer, M.; Szymonski, M.; Echavarren, A. M.; Godlewski, S., Nonacene Generated by On-Surface Dehydrogenation. ACS Nano 2017, 11, 9321-9329. (96) Bendikov, M.; Duong, H. M.; Starkey, K.; Houk, K. N.; Carter, E. A.; Wudl, F., Oligoacenes: Theoretical Prediction of Open-Shell Singlet Diradical Ground States. J. Am. Chem. Soc. 2004, 126, 74167417. (97) Hachmann, J.; Dorando, J. J.; Avilés, M.; Chan, G. K.-L., The radical character of the acenes: A density matrix renormalization group study. J. Chem. Phys. 2007, 127, 134309. (98) Jiang, D.-e.; Dai, S., Electronic Ground State of Higher Acenes. J. Phys. Chem. A 2008, 112, 332-335. (99) Hajgató, B.; Szieberth, D.; Geerlings, P.; Proft, F. D.; Deleuze, M. S., A benchmark theoretical study of the electronic ground state and of the singlet-triplet split of benzene and linear acenes. J. Chem. Phys. 2009, 131, 224321. (100) Qu, Z.; Zhang, D.; Liu, C.; Jiang, Y., Open-Shell Ground State of Polyacenes: A Valence Bond Study. J. Phys. Chem. A 2009, 113, 7909-7914. (101) Hajgató, B.; Huzak, M.; Deleuze, M. S., Focal Point Analysis of the Singlet–Triplet Energy Gap of Octacene and Larger Acenes. J. Phys. Chem. A 2011, 115, 9282-9293. (102) Felix, P.; Hasan, P.; H., G. M.; Florian, L.; Rafael, R.; Joachim, B.; Thomas, M.; Ron, S.; Hans, L., The Multiradical Character of One‐ and Two‐Dimensional Graphene Nanoribbons. Angew. Chem. Int. Ed. 2013, 52, 2581-2584. (103) Small, D. W.; Lawler, K. V.; Head-Gordon, M., Coupled Cluster Valence Bond Method: Efficient Computer Implementation and Application to Multiple Bond Dissociations and Strong Correlations in the Acenes. J. Chem. Theory Comput. 2014, 10, 2027-2040. (104) Horn, S.; Plasser, F.; Müller, T.; Libisch, F.; Burgdörfer, J.; Lischka, H., A comparison of singlet and triplet states for one- and two-dimensional graphene nanoribbons using multireference theory. Theor. Chem. Acc. 2014, 133, 1511. (105) Ibeji, C. U.; Ghosh, D., Singlet-triplet gaps in polyacenes: a delicate balance between dynamic and 26


Page 27 of 28 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

static correlations investigated by spin-flip methods. Phys. Chem. Chem. Phys. 2015, 17, 9849-9856. (106) Bettinger, H. F.; Tönshoff, C.; Doerr, M.; Sanchez-Garcia, E., Electronically Excited States of Higher Acenes up to Nonacene: A Density Functional Theory/Multireference Configuration Interaction Study. J. Chem. Theory Comput. 2016, 12, 305-312. (107) Lee, J.; Small, D. W.; Epifanovsky, E.; Head-Gordon, M., Coupled-Cluster Valence-Bond Singles and Doubles for Strongly Correlated Systems: Block-Tensor Based Implementation and Application to Oligoacenes. J. Chem. Theory Comput. 2017, 13, 602-615. (108) Ghosh, S.; Cramer, C. J.; Truhlar, D. G.; Gagliardi, L., Generalized-active-space pair-density functional theory: an efficient method to study large, strongly correlated, conjugated systems. Chem. Sci. 2017, 8, 2741-2750. (109) Lehtola, S.; Parkhill, J.; Head-Gordon, M., Orbital optimisation in the perfect pairing hierarchy: applications to full-valence calculations on linear polyacenes. Mol. Phys. 2018, 116, 547-560. (110) Perry, J. K.; III, W. A. G.; Ohanessian, G., Inequivalence of equivalent bonds: Symmetry breaking in Co(CH3)2+. J. Chem. Phys. 1992, 97, 7560-7572. (111) Cramer, C. J.; Włoch, M.; Piecuch, P.; Puzzarini, C.; Gagliardi, L., Theoretical Models on the Cu2O2 Torture Track: Mechanistic Implications for Oxytyrosinase and Small-Molecule Analogues. J. Phys. Chem. A 2006, 110, 1991-2004. (112) Gherman, B. F.; Cramer, C. J., Quantum chemical studies of molecules incorporating a Cu2O22+ core. Coord. Chem. Rev. 2009, 253, 723-753. (113) Liakos, D. G.; Neese, F., Interplay of Correlation and Relativistic Effects in Correlated Calculations on Transition-Metal Complexes: The (Cu2O2)2+ Core Revisited. J. Chem. Theory Comput. 2011, 7, 15111523. (114) Witte, M.; Herres-Pawlis, S., Relativistic effects at the Cu2O2 core - a density functional theory study. Phys. Chem. Chem. Phys. 2017, 19, 26880-26889.

26 27 28 29

For Table of Contents Only

30 31

Automatic Construction of the Initial Orbitals for Efficient Generalized Valence Bond Calculations of Large Systems

32 33 34 35 36 37 38 39

Qingchun Wang,† Jingxiang Zou,† Enhua Xu,‡ Peter Pulay,*,¶ Shuhua Li*,† †School

of Chemistry and Chemical Engineering, Key Laboratory of Mesoscopic Chemistry of Ministry of Education, Institute of Theoretical and Computational Chemistry, Nanjing University, 210023, P. R. China ‡Graduate School of Science, Technology, and Innovation, Kobe University, Nada-ku, Kobe 6578501, Japan ¶Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701, United States

40

41 42 27



Page 28 of 28

1 2

28


Automatic Construction of the Initial Orbitals for Efficient Generalized

Recommend Documents