Difficulties and Virtues in Assessing the Potential Energy Surfaces of

Jan 11, 2019 - The main difficulties and virtues in assessing Cκ's via DMBE theory are .... on the possibility of accounting by built-in construction...
2 downloads 0 Views 6MB Size
Subscriber access provided by EKU Libraries

A: Spectroscopy, Molecular Structure, and Quantum Chemistry

Difficulties and Virtues in Assessing the Potential Energy Surfaces of Carbon Clusters via DMBE Theory: Stationary Points of C# (#=2-10) at the Focal Point Carlos Murilo R. Rocha, Jing Li, and António J.C. Varandas J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b11632 • Publication Date (Web): 11 Jan 2019 Downloaded from http://pubs.acs.org on January 13, 2019

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 27 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

The Journal of Physical Chemistry

Difficulties and Virtues in Assessing the Potential Energy Surfaces of Carbon Clusters via DMBE Theory: Stationary Points of Cκ (κ=2-10) at the Focal Point C. M. R. Rocha,† Jing Li,‡ and A. J. C. Varandas∗,‡,† †Department of Chemistry and Coimbra Chemistry Centre, University of Coimbra 3004-535 Coimbra, Portugal. ‡School of Physics and Physical Engineering, Qufu Normal University, 273165, P.R. China E-mail: [email protected]

Abstract The previously reported potential energy surfaces (PESs) of ground-state singlet C3 and triplet C4 are here utilized as input for the construction of approximate cluster expansions for larger Cκ (κ=5-10) species. Relying primarily on the double many-body expansion (DMBE) approach, global potentials are obtained by summing the total interaction energies of all atomic sub-clusters up to four-body terms. The notable capability of the final forms in predicting good estimates of the linear global minima and their thermochemical/structural properties may provide important insights into the structure-determining nature of the (2+3+4) terms. The main difficulties and virtues in assessing Cκ s via DMBE theory are analyzed and new prospects given for the construction of global reliable PESs for the target species.

1

Introduction

New avenues on the field of carbon chemistry have been traced after the isolation and structural characterization of the soccer-ballshaped C60 molecule. 2 Such carbon allotrope consists of 60 sp2 -hybridized C atoms which are arranged in 12 pentagons and 20 hexagons to form an Ih -symmetrical structure. 2 Indeed, several other fullerene structures such as C20 , C36 , C70 , C76 , C78 , C84 and C90 have soon been recognized. 1 The astonishing properties and diverse potential applications of fullerenes motivated a huge research to identify new allotropes. These include single- and multi-walled carbon nanotubes, 3 single-layer graphene sheets 4 and the introduction of various elusive synthetic ones, 1 all with remarkable electrical, mechanical, chemical and optical properties. Undoubtedly, a detailed knowledge of the

Elemental carbon exists in nature as two common allotropes, diamond and graphite. They differ in the way C atoms are attached together to form extended crystalline networks. Diamond is made of sp3 -hybridized (tetrahedral) carbon atoms, while graphite comprises stacked sp2 carbon sheets (or graphene layers) arranged in a hexagonal manner. As expected from such distinct bonding patterns, such allotropic forms exhibit quite remarkable differences in their observable properties. Diamond is a transparent insulator material which is widely appreciated as the hardest known natural substance. On the other hand, graphite is an opaque, fragile dark material with exceptional conductivity properties due to its delocalized π electrons. 1

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 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

fascinating (often unexpected) complexities of large carbon aggregates is only attainable once the chemical properties of their involved precursors (i.e., the small pure carbon molecules) have been clarified. In fact, the small Cκ clusters have been the subject of an intense research effort and continue to be target prototypes for understanding a variety of complex chemical environments. 5,6 Molecules of this kind are known to appear in a myriad of astrophysical objects, particularly in circumstellar envelopes of carbon stars 7–9 and in the sun 10 as well as in comets 11 and interstellar molecular clouds. 12,13 Carbon is formed by nuclear fusions in the cores of stars (a phenomenon often referred to as 3α process 14 ) which, in the late stages of stellar evolution, is ejected into their outer layers and subsequently to the interstellar medium. In these environments, chemical reactions are prone to occur that transform atomic C gas into complex organic molecules and dust grains. 15 Small carbon aggregates are then thought to be fundamental building blocks for formation of interstellar C-rich compounds like HC2κ+1 N, Cκ H, Cκ H2 , Cκ N, Cκ O, Cκ S, polycyclic aromatic hydrocarbons and fullerenes (C60 and C70 ). 15,16 Apart from their astrochemical ubiquity, such clusters are prominent species in equilibrium hot carbon vapors and plasmas generated through energetic processing of Ccontaining materials, hydrocarbon flames and in soot-forming terrestrial environments. 5,6 Clearly, the elucidation of the structural, spectroscopic, and energetic properties of such species could provide valuable tools for unraveling their underlying chemistry and also for clarifying the unknown growth processes leading to the formation of larger aggregates. Yet, small Cκ clusters are highly reactive (shortlived) molecules under most laboratory conditions which makes their experimental characterization extremely involved. 5,6 Moreover, the existence of numerous close-in-energy isomeric forms (ranging from linear chains to monocyclic/polycyclic rings), high-density of lowlying singlet/triplet/quintet electronic states, pronounced biradical character and (pseudo-) Jahn-Teller effects all complicate their theoretical description. 17–26 Definitely, a detailed study

Page 2 of 27

of these molecules requires a fruitful interplay of state-of-the-art experimental techniques and high-level ab initio approaches. What bridges the gap between the latter two is unavoidably the potential energy surface(s), PES(s). Within the adiabatic approximation, 27 these describe the potential energy of a system in terms of the relative positions of the atoms that make up the system. Inevitably, such potentials carry precise information about the underlying species, and are therefore conceptually important tools for the analysis of all structural isomers, spectroscopy and chemical reaction dynamics. 28 From the theoretical perspective, PESs are obtained by solving the electronic Schrödinger equation at sufficiently many (fixed) nuclear configurations whose energies may subsequently be modelled by an analytic form. Amongst the most reliable global modeling schemes, the many-body expansion (MBE) 28 and, its improved variant, the double many-body expansion (DMBE), 29,30 play a prominent role, having acquired a high popularity. In both methodologies, the total interaction potential is defined by making an expansion of the energy in terms of all involved subclusters of atoms. 28 Besides accurately describing valence interactions, one of the main advantages of such approaches relies on the possibility of accounting by built-in construction for the correct asymptotic behavior of each nbody term so that all dissociation limits (as well as long-range interactions in DMBE) are in principle warranted. Of course, if all potentials of all n-body fragments are added, the manybody expansion enables an exact description of the PES of the target polyatomic; 31–33 and references therein. Generally, even if the series shows a fast convergence, high accuracy is attainable only when including the highest-order contributions to the potential. 28,30 In the present work we assess the reliability of the DMBE method in providing first estimates of global PESs for carbon clusters, Cκ (κ=5-10). The potentials are constructed using only the previously reported PESs of ground-state singlet C3 and triplet C4 , 33,34 thence by summing the total interaction energies of all atomic subclusters up to four-body terms. The quality of

ACS Paragon Plus Environment

2

Page 3 of 27 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

The Journal of Physical Chemistry

system assumes the form 28

the DMBE forms and their predictive nature are then judged on the basis of accurate theoretical and experimental results. The main difficulties and virtues in assessing Cκ s via DMBE method are discussed and new prospects given for constructing the global accurate PESs of the title clusters. Several empirical potentials have been proposed over the past few years with the aim at studying the structure and energetics of carbon clusters. 35–43 Of them, we highlight the Reactive Empirical Bond Order (REBO) potential of Brenner et al., 35 which is a generalization of the previous Brenner 36 and Abell-Tersoff 37 forms, and the most recent Reactive Force Field (ReaxFF) of van Duin et al. 38 Other simple and reliable potentials include the TLHT form of Takai et al. 39 and the Murrell-Mottram (MM) functions. 40–43 All these (generally written as a sum of two- and three-body terms only) share the common feature that their parameters have been calibrated from bulk graphite and diamond solid properties 39–42 and/or via hydrocarbon bond parameters. 35,36,38 Although reproducing the bonding patterns of fullerenes, graphene layers, nanotubes and other mediumto-large sized clusters, 43,44 they cannot be expected to be reliable for small Cκ s. 42 Our approach follows the opposite direction (“bottomup” instead of “top-down”) by not using any parameter but those inherent to the PESs of the involved elemental clusters. With this in mind, we expect to obtain predictions exempt of any input information but the number of C atoms.

Vabc...κ (R) = +

Va(1) +

X

(2)

Vab (Rab )

a ab (3) (n) Vabc (Rab , Rac , Rbc )+. . .+Vabc...κ (R)

(1)

abc (1)

where the one-body term Va is the energy of atom a in the state produced by adiabatically removing it from the cluster. If the reference is (1) chosen as all atoms in their ground states, Va is only non-zero when, on dissociation, atom a (2) is left in an excited state. In turn, Vab (Rab ) is a two-body interaction potential for the sub(3) system ab, Vabc (Rab , Rac , Rbc ) is the abc threebody energy, and so on. Note that the summations in Eq. (1) run over all n-body terms   κ! κ , (2) = n n!(κ − n)! each depending on Rn ≡ {R1 , R2 , . . . Rn(n−1)/2 }, the set of interparticle coordinates which is a subset of the total number of interatomic separations Rκ ≡ {R1 , R2 , . . . Rκ(κ−1)/2 }. Thus, we can rewrite Eq. (1) in a more compact form as κ

V (R ) =

κ X X

V (n) (Rn ).

(3)

n=1 Rn ⊂Rκ

As already noted, in the MBE, each n-body term vanishes as any of its constituent atoms is removed to infinity. Such a requirement is generally satisfied by writing V (n) (Rn ) = P (Rn )T (Rn ),

The paper is organized as follows. Section 2 surveys the MBE and DMBE methods. The details of the molecular potentials and global optimization algorithm are in sections 3 and 4, while the results and discussion are in section 5. Section 6 gathers the conclusions.

2

X

X

(4)

where P (Rn ) is an n-body polynomial in the interparticle coordinates and T (Rn ) is a range determining factor that vanishes when any of the coordinates in Rn becomes infinite. 28 A commonly accepted and conceptually convenient approach in obtaining realistic functional forms for molecular PESs is to make a further partition of the potential by splitting each n-body term of Eq. (3) into extended Hartree-fock (EHF) and dynamical correlation (dc) energy-type components. 30 Thence, in

MBE and DMBE methods: a synopsis

In MBE theory, the molecular PES of a κ-atom

ACS Paragon Plus Environment

3

The Journal of Physical Chemistry

DMBE theory, 29 the PES assumes the form

should reflect the intrinsic complexities of the molecule at hand, and hence must be modeled accordingly (i.e., be system-specific). 31,47–49 In fact, for small polyatomic systems, it is a longstanding belief that the main features of their PESs are contained in the two- and three-body terms, with the four-body (and higher-order) terms being regarded as a fine-tuning to give chemical accuracy. 28 If this is the case, only a simple and approximate representation of the latter would generally be required.

κ i X X h (n) (n) V (R ) = VEHF (Rn )+Vdc (Rn ) . κ

n=1 Rn ⊂Rκ

(5) Clearly, one of the major advantages of DMBE over the conceptually simpler MBE approach relies in the allowance for the different convergence rates of the n-body terms at short separations, where the EHF energy is the dominant contribution, and at large distances, where the dispersion (dynamical correlation contributions, dc, or even other long-range energy components such as due to simple electrostatic interactions) dominate. It is also fair to say that by performing such a double decomposition of the interaction energy, one can add physical insight into the origin of the potential. 30

3

# of n-body terms

40

180

(2+3+4)

VCκ

120

2-body 4-body

20

10

0 2

4

5

6

7

cluster size (k)

60

0

3

2

3

4

5

6

7

8

9

(2)

(3)

(4)

(R) = VC2 (R) + VC3 (R) + VC4 (R),

(6) where R is a collective variable defining the (2) κ(κ − 1)/2 interparticle distances; VC2 (R), (3) (4) VC3 (R) and VC4 (R) represent the sum of all two-, three- and four-body potentials [second summation of Eq. (5)]. As in this equation, each n-body term is split into EHF and dc energy contributions; for brevity, their functional forms will not be given here, with the reader being addressed to the original papers 33,34,49,50 for details. Note that we assume that all the fragments dissociate into ground-state C atoms, and hence there is no need for one-body terms in Eq. (6); the zero of energy is set as κ C(3 P ). To avoid numerical overflows in the evaluation of the factorials in Eq. (2), the number of twobody terms have been obtained by calculating the corresponding binomial coefficient as 51   κ = nint {exp [ln (κ) − ln (2) − ln (κ − 2)]} , 2

3-body 30

DMBE PESs

Following previous work on C4 , 33 the approximate cluster expansions for carbon clusters Cκ (κ=5-10) assume the form

240

# of n-body terms

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

Page 4 of 27

10

cluster size (k)

Figure 1: Number of n-body terms with n=2, 3, 4 as a function of the carbon cluster size κ; see Eq. (2).

Varandas and coworkers proposed simple, yet reliable, physically motivated forms for the various n-body terms of Eq. (5). For the simplest case of a diatomic potential energy curve, the DMBE method finds its usefulness in the well-known extended Hartree-Fock approximate correlation energy for two-body interactions (EHFACE2) and the EHFACE2U (EHFACE2 including united atom limit) models. 45 For the three-body energy terms, several approaches have also been advocated in obtaining realistic and flexible functional representations. 30,46 Note that a suitable analytical form

(7) with the number of three and four bodies being evaluated by the recurrence relation 51     κ−n κ κ . (8) = n+1 n+1 n

ACS Paragon Plus Environment

4

Page 5 of 27 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

The Journal of Physical Chemistry

Metropolis criterion: 54 if the resulting energy change ∆V = V (c0 )−V (c) < 0, the point is accepted (downhills), and the configuration with the displaced atom retained, i.e, c is set to c0 . Instead, if ∆V ≥ 0, a probabilistic criterion is utilized: the acceptance probability of the generated configuration (uphills) is p(∆V ) = exp (−∆V /T ). Random numbers are then uniformly distributed in the interval [0,1] and a selected number is compared with p(∆V ). If it is less than p(∆V ), the new point is accepted; otherwise, c0 is rejected. All the above steps

These are plotted in Figure 1 as a function of the cluster size κ; see Table 1 for the total number of n-body terms up to n = 10. Table 1: Number of n-body terms up to n=10 as a function of Cκ ; see Eq. (2).

Cluster C2 C3 C4 C5 C6 C7 C8 C9 C10

# of n-body terms κ=2 3 4 5 6 7 8 1 3 1 6 4 1 10 10 5 1 15 20 15 6 1 21 35 35 21 7 1 28 56 70 56 28 8 1 36 84 126 126 84 36 9 45 120 210 252 210 120 45

9

10

Cκ Table 2: Zero point (EZPE ) and total interac-

1 10

tion energies (ECκ ) of the global minima predicted from the V (2+3+4) PESs of Cκ (κ=2-10).

1

Cluster

For a fixed molecular geometry, the determination of such n-body terms (with n=2,3,4) has been accomplished by attributing indexes to each C atom (up to κ) in the Cκ cluster and, following the combinatorial algorithm of Nijenhuis and Wilf, 52 generating the total number of n-subsets that can be formed from an κ-set {1,2,...,κ}. Finally, to each of such subsets, n(n−1)/2 interparticle coordinates are assigned and the corresponding two-, three- and four-body potentials evaluated. Note that the coordinate system shown in Figure 15(b) of Ref. 33 is utilized for the latter.

4

Cκ a EZPE /kJ mol−1

ECκ b /kJ mol−1

9.8 (9.8c ) 20.3 (21.2c ) 31.3 (33.9c ) 38.9 (49.7c ) 49.2 (60.7c ) 58.9 (76.1c ) 67.2 (86.6c ) 75.9 (101.9c ) 85.7 (112.2c )

−610.0 (−605.8c ; −610.0±2.1d ) −1372.5 (−1342.6c ; −1324.1±13.0d ) −1860.5 (−1829.5c ; −1826.8±17.0d ) −2228.6 (−2546.1c ; −2524.6±17.0d ) −3155.9 (−3061.4c ; −3015.8±20.0d ) −3732.7 (−3747.3c ; −3729.1±20.0d ) −4310.3 (−4278.4c ; −4203.7±43.5c ) −4884.5 (−4952.0c ; −4932.7±43.5c ) −5459.4 (−5502.1c ; −5488.0±43.5e )

C2 `-C3 `-C4 `-C5 `-C6 `-C7 `-C8 `-C9 `-C10

Optimization method

a

Since we are primarily interested in the determination of the global minima of the PESs that outcome from the cluster expansions of Eq (6), a simulated annealing (SA) optimization algorithm as described and implemented by Corana et al. 53 has been utilized. Starting at some “high” computational temperature T (≈ 1000 K) and a user supplied Cartesian coordinate vector c ≡ {c1 , c2 , . . . c3κ }, a function evaluation is made and its initial value V (c) recorded. New candidates are next created around the current point c by applying random moves along each coordinate direction. 53 The new geometry, say c0 (displaced along ci ), is then accepted or rejected according to the

This work. Zero point energies calculated as

P

wi are 3κ−5 harmonic vibrational frequencies.

b

1 i 2 hwi ,

where

This work. The

zero of energy is set to the infinitely separated κ C(3 P ) atoms. c Ref. 55. Energies predicted from the W4 theory via Eq. (14). Zero point energies obtained at B3LYP/pc-2 level of d theory and scaled by 0.985. Experimental energies retrieved e from Ref. 56 via Eq. (14). This work. Total interaction energy predicted from a fit to the experimental data using a straight line; see text and Figure 6.

(through all elements of c) are repeated for a T reduction schedule (T 0 = r T , where r is between 0 and 1). The optimization process is stopped at a temperature low enough (T ≈ 0) that no more useful improvements can be expected, according to some termination criterion. This is

ACS Paragon Plus Environment

5

The Journal of Physical Chemistry 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

V(2+3)

V(2) R1 2.483

C2

Page 6 of 27

V(2+3+4)

R1 2.445

C3

R1 R2 2.524 2.435

C4

R1 R2 2.575 2.539

C5 R1 R2 R3 2.483 2.485 2.536

C6 R1 R2 R3 2.478 2.459 2.584

C7

R1 R2 R3 R4 2.488 2.456 2.534 2.668

C8

R1 R2 R3 R4 2.488 2.467 2.530 2.607

C9

R1 R2 R3 R4 R5 2.484 2.465 2.549 2.596 2.554

C10

Figure 2: Cκ (κ=2-10) global minima predicted from truncated DMBE expansions of Eq. (6). Also shown for comparison are the corresponding structures obtained from PESs in which only two- [V (2) ] and two- plus three-body potentials [V (2+3) ] are considered. Bond distances are in atomic units (a0 ); see also Tables 2 and S1 of Supplementary Material (SM) for zero point energies, harmonic vibrational frequencies and total interaction energies. set if the final potential energy obtained from the last temperature iteration differs from the corresponding value at the current temperature by < 1.0 × 10−6 Eh . To test (and/or accelerate) convergence toward true stationary points, the optimal SA structures were subsequently subjected to an iterative local optimization conjugate gradient (CG) algorithm as implemented by Gilbert and Nocedal. 57 The CG starts with the definition of a descent direction (dk ) as ( −gk for k = 1 dk = (9) −gk + βk dk−1 for k ≥ 2,

with the dot indicating a scalar product. For every interaction k, a step length αk is defined by means of a one-dimensional line search so that the new position vector ck+1 = ck + αk dk satisfies the “sufficient descent” (V (ck+1 ) ≤ V (ck ) + ξαk gkT dk ) and the “directional derivaT tive” (kgk+1 dk k ≤ µkgkT dk k) conditions; 0 < ξ < µ < 1 (typically, ξ = 10−4 and µ = 0.9). 57 The above steps are repeated until the convergence tolerance is reached. This is maintained until the largest component of the gradient vector is < 1.0 × 10−4 Eh a−1 0 (and the maximum energy change is < 1.0 × 10−6 Eh ). Given the highquality of the starting geometries, all global Cκ minima are typically obtained after 1 to 4 CG optimization steps. All geometry optimizations were followed by a harmonic vibrational analysis to confirm the nature of the stationary points (minima: only

where gk = g(ck ) is the 3κ gradient vector calculated at the molecular geometry ck and βk is a scalar chosen according to the formula 57 βk = [~gk · (~gk − ~gk−1 )] /kgk−1 k2 ,

(10)

ACS Paragon Plus Environment

6

Page 7 of 27

real frequencies). Normal modes and frequencies (wi ) are then obtained by diagonalizing the 3κ × 3κ mass-weighted Cartesian force constant matrix F , 58 LT F L = Λ, (11)

obtained from expansions in which only two[V (2) ] and two- plus three-body terms [V (2+3) ] are considered. The structural parameters, har2

Q = LT q,

Triplet

1

where Λ is a diagonal matrix consisting of 3κ eigenvalues λi = 4π 2 c2 wi2 (c is the speed of light) and L is an orthogonal matrix, with T denoting its transpose. Note that, for bent (linear) clusters, there are 3κ−6 (3κ−5) λs which are not zero and hence 3κ−6 (3κ−5) eigensolutions (modes of vibration and frequencies). The remaining six (five) roots are null and correspond to translations and rotations. The transformation between mass-weighted Cartesian displacement coordinates q = M1/2 ξ to normal coordinates Q is thus given by 58

ΔEST/102 kJ mol-1

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

The Journal of Physical Chemistry

k

ref.

1 3 5 7 9

(33) 2 (55) (33) 4 (20) (21) 6 (22) (23) 8 (6) (6) 10 (6)

k

ref.

0

-1

-2

Singlet -3

1

2

3

4

5

6

7

8

9

10

cluster size (k)

Figure 3: Literature survey of the energy differences between the lowest singlet and triplet states of linear carbon chains `-Cκ (κ=1-10). Blue and red shaded areas illustrate the regions where triplet and singlet electronic configurations tend to prevail for a given κ.

(12)

or, in terms of Cartesian displacements ξ, T

Q=L M

1/2

¯T

ξ = L ξ,

monic frequencies, zero point and total interaction energies [with respect to the infinitely separated κ C(3 P ) atoms] for the linear `-Cκ clusters so obtained are given in Tables 2 and S1 of the Supplementary Material (SM), together with available experimental results 6,56,59 as well as theoretical W4 ones reported by Karton et al. 55 Suffice to say that such composite quantum chemistry protocol relies on explicit extrapolations to the complete basis set limit for both Hartree-Fock (HF) and valence Coupled-Cluster (CC) correlation energies (up to iterative pentuples CCSDTQ5), in addition to correct additively for core-core/corevalence correlation, scalar relativistic and diagonal Born-Oppenheimer effects. 60 All reference geometries are obtained at the CCSD(T)/ccpVQZ (for brevity, CC/VQZ) level of theory. 60 For all atomization energies at 0 K P comparison, Cκ ( D0 ) reported in Refs. 56 and 55 were converted to total interaction energies (ECκ ) following the relation X  X Cκ ECκ = − D0Cκ + EZPE =− DeCκ ,

(13)

where M is a diagonal matrix whose elements ¯ are are nuclear masses and the columns of L representations of the normal modes in terms of Cartesian displacements. Note that all first and second-order derivatives of the potentials with respect to the latter were determined numerically using standard finite-difference schemes. 51 One should comment on passing that PESs for clusters as large as C20 were also obtained via Eq. (6). However, for κ > 10, the predicted optimal structures (typically of bent-type) either attain gradient components larger than the usual threshold (1.0 × 10−4 Eh a−1 0 ) or show imaginary frequencies along the bending modes. Thence, they are predicted to be not minima, which in this case too appear to agree to the best of our knowledge with observation.

5

Results and discussion

Figure 2 shows all global minima predicted from the truncated V (2+3+4) PESs of Eq. (6) using the SA/CG methodology. Also shown for comparison are the corresponding geometries

(14) Cκ where EZPE are scaled zero-point P Cκ vibrational 55 energies (Table 2) and De the corre-

ACS Paragon Plus Environment

7

The Journal of Physical Chemistry

sponding dissociation energies (at the Cκ s equilibrium geometries) with respect to the infinitely separated κ C(3 P ) atoms.

and   C2 (a 3 Πu ) + C2 (a 3 Πu )     3 e 1 Σ+ `-C3 (X g ) + C( P ) 3 − `-C4 ( Σg ) −→  C2 (a 3 Πu ) + 2C(3 P )     4C(3 P ),

0.2

2-body energy 0.1

Energy/Eh

0.3

-0.1

C2 - ref. (42)

(16b) (16c) (16d)

C2(a3Πu) - this work

-0.2

0.1

C2(X1Σ+g) - ref. (19) 1

2

3

4

5

3-body energy

0.2

C2 - ref. (43)

-0.3

(16a)

which further restricts any dissociation of the

0.0

6

7

Energy/Eh

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

Page 8 of 27

8

R/a0

Figure 4: Potential energy curves of C2 as a func-

0.0 -0.1

C3 - ref. (42)

TS

-0.2

C3 - ref. (43)

-0.3

Min

C3(1A’) - this work

-0.4

tion of the internuclear distance. The curve employed in the present work is compared with the corresponding two-body energy terms due to Eggen et al. 42 and Tailor et al. 43 implemented in their Murrell-Mottram empirical functions. For comparison, the most accu19 rate ground-state potential curve of C2 (X1 Σ+ is also g) shown. Total interaction energies are given with respect to the infinitely separated C atoms.

-0.5 -0.6

20

40

60

80

100

120

140

160

180

∠CCC/deg

Figure 5: Optimized bending potential for groundstate C3 obtained from the PES utilized in this work 34 and the Murrell-Mottram empirical functions due to Eggen et al. 42 and Tailor et al. 43 Also shown, for comparison (in the grayish part) are the corresponding three-body energy terms. Global minima and isomerization transition states are indicated by arrows. Total interaction energies are given with respect to the infinitely separated C atoms.

Following seminal work due to Pitzer and Clementi 61 who first recognized `-Cκ forms as low energy structures on their corresponding PESs, odd- and even-numbered chains (for κ > 3 − 2) correlate with 1 Σ+ g and Σg ground electronic states, respectively. Thence, the comparisons with the theoretical results will be done accordingly. Although recent ab initio calculations support the fact that, for even-numbered structures, monocyclic (singlet) isomers are almost isoenergetic or even more stable than linear (triplet) arrangements, 6,26 no cyclic stationary points have been predicted for Cκ (κ=4,6,8,10). This can be mostly attributed to the highly repulsive nature of the four-body terms at ring forms and/or the lack of higher-order (attractive) non-pairwise contributions. 33 Following Ref. 33, the global PESs of C3 and C4 dissociate adiabatically into ( C2 (a 3 Πu ) + C(3 P ) (15a) ˜ 1 Σ+ `-C3 (X ) −→ g 3C(3 P ), (15b)

Cκ (κ=5-10) PESs to occur according to the above fragments. Note that the above correlations are the only ones used for the analysis below, thus assuming that the difference in energy between the singlet and triplet states of C2 are small enough (even if magnified by its occurrence; see Table 1) to leave the possibility of using either its singlet or triplet state to satisfy the required spin correlations. It is then postulated that the odd linear Cκ are predicted in the singlet state, while the corresponding even ones are in the triplet: thence, all linear chains are in their most stable forms shown in Figure 3. It is appropriate to comment at this point on the reliability of the n-body terms here utilized when compared with other commonly used (semi-)empirical forms. 35–43 As an illustration, Figures 4 and 5 show two- and threebody energy terms for the C2 and C3 clusters,

ACS Paragon Plus Environment

8

Page 9 of 27

The Journal of Physical Chemistry 0

4 W4-expt.

kC(3P)

-1

DMBE-expt.

3 2 -1 ΔEC /10 kJ mol

-2 -3 -4

(2)

2 1

k

V

k

3 -1 EC /10 kJ mol

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

-5

V(2+3)

-6

V W4

(2+3+4)

0 -1

expt. -7

2

3

4

5

6

7

8

9

-2

10

cluster size (k)

2

3

4

5

6

7

8

9

10

cluster size (k)

Figure 6: Total interaction energies [with respect to

Figure 7: Deviations between the total interaction en-

the infinitely separated κ C(3 P ) atoms] of the linear global minima predicted from the V (2+3+4) PESs. W4 results are taken from Ref. 55; see Eq. (14). Experimental energies retrieved from Refs. 56 and 55 via Eq. (14). Black dashed curve represents a fit to the experimental data using a straight line; see text.

ergies predicted from the V (2+3+4) PESs of `-Cκ and the ones from W4 theory 55 as well as the experimental estimates of Gingerich et al. 56 Experimental uncertainties are shown by red/blue solid lines.

Figures 6 and 7 assess the accuracy of the predicted thermochemistry of the `-Cκ against available literature results. 55,56 Also shown, in Figure 8, are the deviations between averaged ¯ and those reported at C – C bond lengths (R) CC/VQZ level 55 as well as experimental estimates. 6,59 As shown, the total interaction energies so obtained agree reasonably well with the available measured data of Gingerich and coworkers 56 and the W4 55 results. This is interesting given the lack of higher-order n-body terms, and hence the approximate nature of the PESs. One exception to this trend is perceived in `-C5 . For such species, the V (2+3+4) form underestimates the interaction energy by about 300 kJ mol−1 when compared with experimental and ab initio data, which is of nearly the magnitude of the lacking (attractive) five-body term. Note, however, that for larger aggregates, discrepancies between the predicted and available interaction energies tend to decrease with increasing cluster size in a systematic way; see Figure 7. Such a behavior is likely to stem from the fact that higher-order terms tend to cancel each other 28,33 so that their final contribution to the total energy is minor. This is equivalent to assume that the expansion in Eq. (5) is rapidly convergent after n ' 4. As seen in Figure 8, the current DMBE PESs not only predict reasonably well the thermochemistry of the relevant global min-

from which the MM empirical functions 42,43 are built. Recall that such potentials in the former reference 42 were obtained from a fit to phonon frequencies and elastic constants of diamond as well as to the cohesive energy and intralayer spacing of graphite. 42 In the other, 43 they were reparametrized to reproduce the vibrational spectroscopy of carbon nanotubes and graphene-based materials. 43 As Figure 4 shows, the C2 curves of both PESs are strongly repulsive near their equilibrium distances (besides having the minimum too far in the case of Eggen et al. 42 ) when compared with our own ab initio potential. As a result, the three-body energies become rather more attractive for bent C3 arrangements and rather more repulsive for linear ones. Although the MM potentials are conceptually very simple and reliable, the “topdown” approach on which they (and others) are based favor, for small clusters, primarily cyclic (in some cases, unrealistic) structures as global minima, 42 starting at the trimer in Figure 5. They read from top to bottom in Table 2 (in kJ mol−1 , with the symmetry in brackets 42 ): −607.6 (D∞h ), −1156.1 (C2v ), −1776.1 (D4h ), −2683.0 (D5h ), −3246.2 (D6h ), −3779.2 (Cs ), −4596.6 (C2v ), −5172.2 (C2v ), −5994.3 (Td ). Clearly, we slightly enhance the accuracy while attributing the correct symmetries.

ACS Paragon Plus Environment

9

The Journal of Physical Chemistry 14

found to differ by about 300 cm−1 (in the worst cases) from the ones reported at CASPT2 21,22 and CCSD(T)-F12 23 levels of theory; of course, such differences also reflect the finite nature of the one-electron and N -electron basis sets. Interestingly, Figure 6 also suggests that the total interaction energies (either theoretical or experimental) behave as a linear function with respect to the size κ of the `-Cκ clus3 − Σg states). ters (for alternating 1 Σ+ g and As shown by the black dashed curve, the third-law measurements of Gingerich et al. 56 can be accurately represented by the relation ECκ (κ) = a + bκ, where a = 580.6 kJ mol−1 and b = −606.8 kJ mol−1 . From this and considC10 ering 55 EZPE = 112.2 kJ mol−1 in Eq. (14), we are able to estimate the atomization energy of `-C10 (3 Σ− g ) at 0 K (not reported so far) as being 5375.7±43.5 kJ mol−1 . This is to be compared with the ab initio value of 5389.9 kJ mol−1 reported by Karton et al. 55 and 5373.7 kJ mol−1 predicted from the V (2+3+4) form. So, despite the fact that monocyclic singlet structures were conjectured the most stable forms for the evennumbered Cκ [κ=4(D2h ), 6(D3h ), 8(C4h ), and 10(D5h )] series, 6,17,18,25,26,55 the good agreement between the estimated values of W4 55 for `-Cκ (3 Σ− g ) and the experimental atomization energies 56 is notable. Figure 9 shows sample contour plots for the collinear reactions `-C(κ−2) +C2 → `-C(κ−1) +C (with κ ≥ 5) obtained from the truncated DMBE expansions of Eq. (6). As expected, the global V (2+3+4) forms describe in a physically reasonable manner both valence and longrange features of the potential, including the correct asymptotic behavior at dissociation, a long-standing asset of DMBE theory. 29,30 In fact, based on the results here shown, we can definitely judge the current methodology and the DMBE PESs obtained therefrom as quite satisfactory in providing first estimates of the Cκ s true potentials. Actually, the excellent predictive capability of these V (2+3+4) forms gave us important insights into the structuredetermining nature of the (2+3+4) terms and also opened new avenues toward the construction of global accurate PESs of carbon clusters (and possibly other atomic clusters). Such a

12 10

Δ− R/10-2 a0

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

8 6

CC/VQZ-expt. DMBE-expt.

4

DMBE-CC/VQZ 2 0 -2

2

3

4

5

6

7

8

9

Page 10 of 27

10

cluster size (k)

Figure 8: Deviations between averaged C – C bond lengths predicted from the V (2+3+4) PESs of `-Cκ and the ones reported at CCSD(T)/VQZ level 55 as well as the corresponding experimental estimates. 6,59

ima, but also their corresponding structural parameters. Accordingly, the predicted average bond lengths are overestimated by only ≈ 0.1 a0 when compared with experimental and CC/VQZ optimized ones. Figure 2 also unravels the structure-determining nature of the (2+3+4)-body terms, as one goes from highly packed to completely unpacked (linear) structures by the gradual addition of two-, threeand four-body potentials in Eq. (6). It further shows that the inclusion of two- and three-body terms alone does not seem sufficient to provide reasonable geometries for elemental carbon clusters (this is actually the case for the MM potentials discussed above). For example, although the square planar C4 and cubic C8 structures are predicted as global minima on the corresponding V (2+3) expansions, they are confirmed by ab initio calculations to be highlying stationary structures. 20,24 A clearer picture on the convergent nature of the n-body terms can also be perceived from Figure 6. As already noted, the V (4) (V (3) ) terms act in such a way as to counterbalance the increased attraction of the V (3) (V (2) ) contributions. Although the V (n>4) terms do contribute to the final interaction potentials, their role should be of fine-tuning structural parameters, frequencies, and bond energies rather than a structure-determining one. For example, the harmonic frequencies for `-C5 , `-C6 and `-C7 are

ACS Paragon Plus Environment

10

Page 11 of 27 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

The Journal of Physical Chemistry

1

3

R2

7

5

12

8

2

(a)

R1

R2

4

1

9

6

15

10

16

11

R2

5

8

6

12 1

3 4

11 6

7

9

13

1

10

7

5

11

6

14

15

(f) R1

R2

7 6

11

8 9

12

2 3 4 5

12

9

(e)

R1

10

10 13

4

7 8

8

2 3

15

11

R2

5

13 14

(d)

R1

2 1

9

4 3

12 7

2

(c)

R1

R2

8

6

14

13

(b)

R1

1 2

13

14

10

11 3 4

5

12 13 10

9

Figure 9: Partially relaxed contour plots for the collinear reactions: (a) `-C3 + C2 → `-C4 + C, with contours equally spaced by 0.02 Eh , starting at −1.1 Eh ; (b) `-C4 +C2 → `-C5 +C, with spacing 0.03 Eh , starting at −1.5 Eh ; (c) `-C5 +C2 → `-C6 +C, with spacing 0.03 Eh , starting at −1.7 Eh ; (d) `-C6 +C2 → `-C7 +C, with spacing 0.03 Eh , starting at −1.9 Eh ; (e) `-C7 +C2 → `-C8 +C, with spacing 0.03 Eh , starting at −2.1 Eh . (f) `-C8 +C2 → `-C9 +C, with spacing 0.03 Eh , starting at −2.3 Eh . All hidden coordinates are partially relaxed in the interval 2.4 6 Ri /a0 6 2.8. task can then be accomplished by employing relatively simple functional forms for higherorder terms, say n = 5 or 6 (i.e., assuming that the MBE/DMBE series is rapidly converging). Unfortunately, this situation is difficult to judge since for, e.g, C10 , there are 210 and 120 6and 7-body terms that we might be tempted to ignore (Table 1), even knowing that their individual contributions can be minor (see, e.g., Figure 7 and Table 2). For larger clusters, however, it is fair to think that many of these lacking terms approximately vanish since the atoms in these sub-clusters may be far away from one another. Moreover, we cannot exclude the possibility of mutual cancelations as noted above. Indeed, due to the typical smallness of such terms and factorial nature of their number, 62 it may be even recommendable (despite the suggestion above) to ignore them rather than include crude estimates to avoid error proliferation. On the other hand, experience shows

that, for small polyatomics (those with three, four or even five atoms), it is necessary to retain all terms if chemical accuracy is desired. This was clearly seen for `-C5 and (to a lesser extent) `-C6 where the predictive capability of their PESs are largely affected by the lacking higher-order terms. In the worst case, if neither the MBE/DMBE is rapidly convergent nor higher-order terms are functionally simple, there are several benefits in performing truncated expansions of the molecular potentials. One of the most advantageous relates to their excellent predictive capability as first demonstrated here. Given the relatively easiness in construction and the high computational costs in evaluating ab initio gradients and hessians for such large species (such as the ones here envisaged), not only global minima but other extrema predicted from truncated PESs provide good guesses to guide for stationary point searches in actual electronic structure cal-

ACS Paragon Plus Environment

11

The Journal of Physical Chemistry 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

culations. Such an approach proved to be very useful in our previous study on C4 . 33 Of course, the quality and predictive nature of the truncated molecular potentials hinge intimately on: (i) accuracy, and (ii) type of n-body terms from which they have been built. Assuming that (i) has been assured in our previous PESs, we are left with the correct definition of the spin and spatial symmetries of the fragments. The present work also unraveled that inclusion of the (1 + 2 + 3 + 4) terms from the (lowest) triplet and singlet states of C3 and C4 , respectively, may be as important as the ones here considered. Such PESs are amongst our goals and are actually under completion.

6

Page 12 of 27

veis” (CENTRO-01-0145-FEDER-000014) cofinanced by the European Regional Development Fund (FEDER), through “Programa Operacional Regional do Centro” (CENTRO2020). C.M.R.R. thanks also CAPES Foundation (Ministry of Education of Brazil) for a scholarship (BEX 0417/13-0). J.L. thanks the support from National Natural Science Foundation of China (Grant No.11604179) and Shandong Natural Science Foundation (Grant No. ZR2016AQ18). A.J.C.V. gratefully acknowledges China’s Shandong Province “DoubleHundred Talent Plan” (2018).

References

Conclusions

(1) Georgakilas, V.; Perman, J. A.; Tucek, J.; Zboril, R. Broad Family of Carbon Nanoallotropes: Classification, Chemistry, and Applications of Fullerenes, Carbon Dots, Nanotubes, Graphene, Nanodiamonds, and Combined Superstructures. Chem. Rev. 2015, 115, 4744–4822.

The previously reported PESs of ground-state singlet C3 and triplet C4 have here been utilized as input for the construction of approximate cluster expansions for Cκ (κ=5-10) clusters. Relying primarily on the DMBE approach, global molecular potentials are obtained by summing the total interaction energies of all sub-clusters of atoms up to four-body terms. The predictive capabilities of the V (2+3+4) forms in providing reasonable estimates of Cκ s global minima and their associated thermochemical properties gave us important insights into the structuredetermining nature of the (2+3+4) terms. The main difficulties and virtues in assessing Cκ s via the DMBE method were then fully analyzed and new prospects given toward the construction of global accurate PESs for larger carbon clusters.

(2) Kroto, H. W.; Heath, J. R.; O’Brien, S. C.; Curl, R. F.; Smalley, R. E. C60 : Buckminsterfullerene. Nature 1985, 318, 162–163. (3) Iijima, S. Helical Microtubules of Graphitic Carbon. Nature 1991, 354, 56–58. (4) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. Science 2004, 306, 666–669.

Supporting Information

(5) Weltner, W.; van Zee, R. J. Carbon Molecules, Ions, and Clusters. Chem. Rev. 1989, 89, 1713–1747.

List of all structural parameters, harmonic frequencies and total interaction energies at the global minima of the Cκ DMBE forms.

(6) van Orden, A.; Saykally, R. J. Small Carbon Clusters: Spectroscopy, Structure, and Energetics. Chem. Rev. 1998, 98, 2313–2358.

Acknowledgement This work is supported by Fundação para a Ciência e a Tecnologia and Coimbra Chemistry Centre, Portugal, in the framework of the project “MATISMateriais e Tecnologias Industriais Sustentá-

(7) Hinkle, K. W.; Keady, J. J.; Bernath, P. F. Detection of C3 in the Circumstellar Shell

ACS Paragon Plus Environment

12

Page 13 of 27 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

The Journal of Physical Chemistry

of IRC+10216. Science 1988, 241, 1319– 1322.

Tunable VUV Experiments and Theory. J. Am. Chem. Soc. 2007, 129, 10229–10243.

(8) Bernath, P. F.; Hinkle, K. W.; Keady, J. J. Detection of C5 in the Circumstellar Shell of IRC+10216. Science 1989, 244, 562– 564.

(18) Yousaf, K. E.; Taylor, P. R. On the Electronic Structure of Small Cyclic Carbon Clusters. Chem. Phys. 2008, 349, 58–68. (19) Varandas, A. J. C. A Simple, Yet Reliable, Direct Diabatization Scheme. The 1 + Σg States of C2 . Chem. Phys. Lett. 2009, 471, 315–321.

(9) Hargreaves, R. J.; Hinkle, K.; Bernath, P. F. Small Carbon Chains in Circumstellar Envelopes. Mon. Not. R. Astron. Soc. 2014, 444, 3721–3728.

(20) Massó, H.; Senent, M. L.; Rosmus, P.; Hochlaf, M. Electronic structure calculations on the C4 cluster. J. Chem. Phys. 2006, 124, 234304–234311.

(10) Brault, J. W.; Testerman, L.; Grevesse, N.; Sauval, A. J.; Delbouille, L.; Roland, G. Infrared bands of C2 in the Solar Photospheric Spectrum. Astron. Astrophys. 1982, 108, 201–205.

(21) Massó, H.; Veryazov, V.; Malmqvist, P. A.; Roos, B. O.; Senent, M. L. Ab Initio characterization of C5 . J. Chem. Phys. 2007, 127, 154318–154326.

(11) Huggins, W. Preliminary Note on the Photographic Spectrum of Comet b 1881. Proc. R. Soc. Lond. 1881, 33, 1–3. (12) Lambert, D. L.; Sheffer, Y.; Federman, S. R. Hubble Space Telescope Observations of C2 Molecules in Diffuse Interstellar Clouds. Astrophys. J. 1995, 438, 740–749.

(22) Massó, H.; Senent, M. L. Ab Initio characterization of C6 . J. Phys. Chem. A 2009, 113, 12404–12410. (23) Al-Mogren, M. M.; Senent, M. L.; Hochlaf, M. Theoretical Characterization + of C7 , C− 7 , and C7 . J. Chem. Phys. 2013, 139, 064301–064309.

(13) Galazutdinov, G. A.; Musaev, F. A.; Krelowski, J. On the Detection of the Linear C5 Molecule in the Interstellar Medium. Mon. Not. R. Astron. Soc. 2001, 325, 1332.

(24) Sharapa, D.; Hirsch, A.; Meyer, B.; Clark, T. Cubic C8 : An Observable Allotrope of Carbon? ChemPhysChem 2015, 16, 2165–2171.

(14) Burbidge, E. M.; Burbidge, G. R.; Fowler, W. A.; Hoyle, F. Synthesis of the Elements in Stars. Rev. Mod. Phys. 1957, 29, 547–650.

(25) Wang, Z.-Q.; Hu, C.-E.; Chen, X.-R.; Cheng, Y.; Chen, Q.-F. Ab Initio Investigation of Structure, Spectrum, Aromaticity and Electronic Properties of C10 Carbon Cluster. Comput. Theor. Chem. 2017, 1118, 94–106.

(15) Henning, T.; Salama, F. Carbon in the Universe. Science 1998, 282, 2204–2210. (16) Cami, J.; Bernard-Salas, J.; Peeters, E.; Malek, S. E. Detection of C60 and C70 in a Young Planetary Nebula. Science 2010, 329, 1180–1182.

(26) Varandas, A. J. C. Even Numbered Carbon Clusters: Cost-Effective Wavefunction-Based Method for Calculation and Automated Location of Most Structural Isomers. Eur. Phys. J. D. 2018, 72, 134–140.

(17) Belau, L.; Wheeler, S. E.; Ticknor, B. W.; Ahmed, M.; Leone, S. R.; Allen, W. D.; Schaefer, H. F.; Duncan, M. A. Ionization Thresholds of Small Carbon Clusters:

ACS Paragon Plus Environment

13

The Journal of Physical Chemistry 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

(27) Born, M.; Huang, K. Dynamical Theory of Crystal Lattices; Clarendon Press: Oxford, 1954.

Page 14 of 27

Chemical Vapor Deposition of Diamond Films. Phys. Rev. B 1990, 42, 9458–9471. (37) Tersoff, J. Empirical Interatomic Potential for Carbon, with Applications to Amorphous Carbon. Phys. Rev. Lett. 1988, 61, 2879–2882.

(28) Murrell, J. N.; Carter, S.; Farantos, S. C.; Huxley, P.; Varandas, A. J. C. Molecular Potential Energy Functions.; John Wiley & Sons: Chichester, 1984.

(38) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396–9409.

(29) Varandas, A. J. C. A Double Many-Body Expansion of Molecular Potential Energy Functions. Mol. Phys. 1984, 53, 1303– 1325.

(39) Takai, T.; Lee, C.; Halicioglu, T.; Tiller, W. A. A Model Potential Function for Carbon Systems: Clusters. J. Phys. Chem. 1990, 94, 4480–4482.

(30) Varandas, A. J. C. Intermolecular and Intramolecular Potentials: Topographical Aspects, Calculation, and Functional Representation via A Double Many-Body Expansion Method. Adv. Chem. Phys. 1988, 74, 255–338.

(40) Murrell, J. N.; Mottram, R. E. Potential Energy Functions for Atomic Solids. Mol. Phys. 1990, 69, 571–585.

(31) Varandas, A. J. C.; Murrell, J. N. A ManyBody Expansion of Polyatomic Potential Energy Surfaces: Application to Hn Systems. Faraday Discuss. Chem. Soc. 1977, 62, 92–109.

(41) Balm, S. P.; Allaf, A. W.; Kroto, H. W.; Murrell, J. N. Potential-Energy Function for Large Carbon Clusters. J. Chem. Soc., Faraday Trans. 1991, 87, 803–806.

(32) Varandas, A. J. C.; Zhang, L. Dynamics of HO2+O3 Reaction using a Test DMBE Potential Energy Surface: Does It Occur Via Oxygen or Hydrogen Atom Abstraction? Chem. Phys. Lett. 2004, 385, 409–416.

(42) Eggen, B. R.; Johnston, R. L.; Murrell, J. N. Carbon Cluster Structures and Stabilities Predicted from Solid-State Potentials. J. Chem. Soc., Faraday Trans. 1994, 90, 3029–3037.

(33) Varandas, A. J. C.; Rocha, C. M. R. Cn (n=2-4): Current Status. Phil. Trans. R. Soc. A 2018, 376, 20170145.

(43) Tailor, P. M.; Wheatley, R. J.; Besley, N. A. An Empirical Force Field for the Simulation of the Vibrational Spectroscopy of Carbon Nanomaterials. Carbon 2017, 113, 299–308.

(34) Rocha, C. M. R.; Varandas, A. J. C. Energy-Switching Potential Energy Surface for Ground-State C3 . Chem. Phys. Lett. 2018, 700, 36–43.

(44) Lai, S. K.; Setiyawati, I.; Yen, T. W.; Tang, Y. H. Studying Lowest Energy Structures of Carbon Clusters by Bondorder Empirical Potentials. Theor. Chem. Acc. 2016, 136, 20–32.

(35) Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. A Second-Generation Reactive Empirical Bond Order (REBO) Potential Energy Expression for Hydrocarbons. J. Phys. Condens. Matter 2002, 14, 783– 802.

(45) Varandas, A. J. C.; Silva, J. D. Potential Model for Diatomic Molecules Including the United-Atom Limit and Its Use in a Multiproperty Fit for Argon. J. Chem. Soc., Faraday Trans. 1992, 88, 941–954.

(36) Brenner, D. W. Empirical Potential for Hydrocarbons for Use in Simulating the ACS Paragon Plus Environment

14

Page 15 of 27 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

The Journal of Physical Chemistry

(46) Varandas, A. J. C. Combined-HyperbolicInverse-Power-Representation of Potential Energy Surfaces: A Preliminary Assessment for H3 and HO2 . J. Chem. Phys. 2013, 138, 054120–054133.

Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087–1092. (55) Karton, A.; Tarnopolsky, A.; Martin, J. M. L. Atomization Energies of the Carbon Clusters Cn (n=2-10) Revisited by Means of W4 Theory as well as Density Functional, Gn , and CBS Methods. Mol. Phys. 2009, 107, 977–990.

(47) Galvão, B. R. L.; Mota, V. C.; Varandas, A. J. C. Modeling Cusps in Adiabatic Potential Energy Surfaces. J. Phys. Chem. A 2015, 119, 1415–1421. (48) Galvão, B. R. L.; Mota, V. C.; Varandas, A. J. C. Modeling Cusps in Adiabatic Potential Energy Surfaces Using a Generalized Jahn-Teller Coordinate. Chem. Phys. Lett. 2016, 660, 55–59.

(56) Gingerich, K. A.; Finkbeiner, H. C.; Schmude, R. W. Enthalpies of Formation of Small Linear Carbon Clusters. J. Am. Chem. Soc. 1994, 116, 3884–3888. (57) Gilbert, J.; Nocedal, J. Global Convergence Properties of Conjugate Gradient Methods for Optimization. SIAM J. Optimiz. 1992, 2, 21–42.

(49) Rocha, C. M. R.; Varandas, A. J. C. Multiple Conical Intersections in Small Linear Parameter Jahn-Teller Systems: the DMBE Potential Energy Surface of Ground-State C3 Revisited. Phys. Chem. Chem. Phys. 2018, 20, 10319–10331.

(58) Wilson Jr., E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations; McGraw Hill: New York, 1955.

(50) Rocha, C. M. R.; Varandas, A. J. C. Accurate Ab Initio-Based Double ManyDody Expansion Potential Energy Surface for the Adiabatic Ground-State of the C3 Radical Including Combined Jahn-Teller plus Pseudo-Jahn-Teller Interactions. J. Chem. Phys. 2015, 143, 074302–074318.

(59) Breier, A. A.; Büchling, T.; Schnierer, R.; Lutter, V.; Fuchs, G. W.; Yamada, K. M. T.; Mookerjea, B.; Stutzki, J.; Giesen, T. F. Lowest Bending Mode of 13 C-Substituted C3 and an Experimentally Derived Structure. J. Chem. Phys. 2016, 145, 234302–234311.

(51) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran 77: the Art of Scientific Computing. Second Edition; Cambridge University Press: New York, 1996; Vol. 1.

(60) Karton, A.; Rabinovich, E.; Martin, J. M. L.; Ruscic, B. W4 Theory for Computational Thermochemistry: In Pursuit of Confident Sub-kJ/mol Predictions. J. Chem. Phys. 2006, 125, 144108–144124. (61) Pitzer, K. S.; Clementi, E. Large Molecules in Carbon Vapor. J. Am. Chem. Soc. 1959, 81, 4477–4485.

(52) Nijenhuis, A.; Wil, H. S. Combinatorial Algorithms for Computers and Calculators. Second Edition; Academic Press: Cambridge, 1978; Vol. 1.

(62) Richard, R. M.; Lao, K. U.; Herbert, J. M. Aiming for Benchmark Accuracy with the Many-Body Expansion. Acc. Chem. Res. 2014, 47, 2828–2836.

(53) Corana, A.; Marchesi, M.; Martini, C.; Ridella, S. Minimizing Multimodal Functions of Continuous Variables with the Simulated Annealing Algorithm. ACM Trans. Math. Softw. 1987, 13, 262–280. (54) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E.

ACS Paragon Plus Environment

15

The Journal of Physical Chemistry 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

TOC Graphic

ACS Paragon Plus Environment

16

Page 16 of 27

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

The Journal of Physical Chemistry

VCk=VC(2)2+VC(3)3+VC(4)4

ACS Paragon Plus Environment

14

The Journal of Physical Chemistry

Page 18 of 27

12

Δ− R/10-2 a0

1 2 10 3 4 8 5 6 6 7 8 9 4 10 11 12 2 13 14 0 15 16 -2 17 18 19

CC/VQZ-expt. DMBE-expt. DMBE-CC/VQZ

2

3

ACS Paragon Plus Environment

4

5

6

7

cluster size (k)

8

9

10

0

Page 19 of 27

kC( P)

-1 -2 -3 -4

(2)

V

k

3 -1 EC /10 kJ mol

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

The Journal of Physical Chemistry 3

-5

V(2+3)

-6

V W4

(2+3+4)

expt. -7

2

3

4 ACS Paragon 5 Plus Environment 6 7

cluster size (k)

8

9

10

V 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

(2)

R1 2.483

C2

The Journal of Physical Chemistry (2+3)

V

V

(2+3+4)

R1 2.445

C3

R1 R2 2.524 2.435

C4

R1 R2 2.575 2.539

C5 R1 R2 R3 2.483 2.485 2.536

C6 R1 R2 R3 2.478 2.459 2.584

C7 C8 C9 C10

R1 R2 R3 R4 2.488 2.456 2.534 2.668

R1 R2 R3 R4 2.488 2.467 2.530 2.607

R1 R2 R3 R4 R5 2.484 2.465 2.549 2.596 2.554 ACS Paragon Plus Environment

Page 20 of 27

240

Page 21 of 27

40

# of n-body terms

# of n-body terms

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

The Journal of Physical Chemistry

180

120

2-body 3-body

30

4-body

20

10

0 2

4

5

6

7

cluster size (k)

60

0

3

2

3

4

5

6

ACS Paragon Plus Environment

7

cluster size (k)

8

9

10

14

∆− R/10-2 a0

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

The Journal of Physical Chemistry

Page 22 of 27

12 10 8 6

CC/VQZ-expt. DMBE-expt.

4

DMBE-CC/VQZ 2 0 -2

2

3

4

5

6

ACS Paragon Plus Environment

7

cluster size (k)

8

9

10

0.3

Page 23 of 27

The Journal of Physical Chemistry

3-body energy

0.2

Energy/Eh

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

0.1 0.0 -0.1

C3 - ref. (42)

TS

-0.2

C3 - ref. (43) 1

-0.3

Min

C3( A’) - this work

-0.4 -0.5 -0.6

20

40

ACS Paragon Plus Environment

60

80

100

∠CCC/deg

120

140

160

180

0.2

Energy/Eh

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

The Journal of Physical Chemistry

Page 24 of 27

2-body energy

0.1

0.0

-0.1

C2 - ref. (42) C2 - ref. (43) C2(a3Πu) - this work

-0.2

-0.3

C2(X1Σ+g) - ref. (19) 1

2

3

4

ACS Paragon Plus Environment

R/a0

5

6

7

8

2

Page 25 of 27

ΔEST/102 kJ mol-1

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

The Journal of Physical Chemistry

Triplet

1

k

ref.

k

ref.

1 3 5 7 9

(33) 2 (55) (33) 4 (20) (21) 6 (22) (23) 8 (6) (6) 10 (6)

0

-1

-2

Singlet -3

1

2

3ACS Paragon 4 Plus5Environment 6 7

cluster size (k)

8

9

10

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

R1

R2

1

3

12

8

1

9

6

15

10

16

11

R2

8 1

6

12 1

3 4

11 6

7

9

13

2 3 4 5

1

10

7

5

11

6

12

9

14

15

(e)

R1

(f) R1

R2

7 6

11 10

10 13

4

7 8

8

2 3

15

11

R2

5

13 14

5

(d)

R1

2

9

4 3

12 7

2

(c)

R1

R2

8

6

14

13

4

Page 26 of 27

(b)

R2

7

5 2

The Journal of Physical Chemistry R1

(a)

8 9

12

1 2

13

14

10

11 3 4

5

12 13 10

9

ACS Paragon Plus Environment

4

Page 27 of 27

The Journal of Physical Chemistry

DMBE-expt.

3

k

2 -1 ΔEC /10 kJ mol

1 2 3 2 4 5 6 1 7 8 9 10 0 11 12 13 -1 14 15 16 17 -2 18 19

W4-expt.

2

3

ACS Paragon Plus Environment

4

5

6

7

cluster size (k)

8

9

10