Exploring Chemical Space with Alchemical Derivatives: BN

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Exploring Chemical Space with Alchemical Derivatives: BNSimultaneous Substitution Patterns in C60 Robert Balawender,*,† Michał Lesiuk,‡ Frank De Proft,§ and Paul Geerlings*,§ †

Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland Faculty of Chemistry, University of Warsaw, Pasteura 1, PL-02-093 Warsaw, Poland § Algemene Chemie, Vrije Universiteit Brussel, Faculteit Wetenschappen, Pleinlaan 2, 1050 Brussels, Belgium ‡

ABSTRACT: With the idea of using alchemical derivatives to explore in an efficient, computer- and cost-effective way Chemical Space was launched several years ago. In the context of Conceptual DFT response functions, these energies vs nuclear charge derivatives permit the estimatation of the energy of transmutants of a given starting or reference molecule showing different nuclear compositions. After an explorative study on small and planar molecules (Balawender et al. J. Chem. Theory Comput. 2013, 9, 5327) by the present authors of this paper, the present study fully exploits the computational advantages of the alchemical derivatives in larger three-dimensional systems. Starting from a single reference calculation on C60, the complete BN substitution pattern, from single substituted C58BN via the belt (C20(BN)20 and the ball C12(BN)24 structures to the fully substituted (BN)30, is explored. Successive and simultaneous substitution strategies are followed and compared, indicating that both techniques yield identical results up to 13 substitutions but that for higher substitutions the simultaneous approach needs to be taken. Due to the cost-efficiency of the algorithm this path can indeed be followed as opposed to earlier work in the literature where for each step a full SCF calculation was at stake leading to prohibitively large computational demands for adopting the simultaneous approach. Previously formulated rules governing the substitution pattern by Kar and co-workers are scrutinized in this context and reformulated giving chemical insight in the gradual substitution process and the relative energies of the isomers. In its present form the method offers an interesting venue to study BN substitution patterns in higher fullerenes and graphene and in general paves the way for more efficient exploration of the Chemical Space.

I. INTRODUCTION Chemistry and chemists are exploring Chemical Compound Space (CCS), the space populated by all imaginable chemicals with natural nuclear charges and realistic interatomic distances for which chemical interactions exist.1−3 In their continuous efforts toward Molecular Design, chemists try to “detect” in this space stable molecules with interesting/optimal properties. Navigating through this space of unimaginable size is costly, especially for experimentalists in front of a myriad of synthetic problems but also for theoreticians. Indeed, exploring Chemical Space by even relatively moderate level ab initio or even semiempirical brute force techniques may already lead to an enormous increase of the number and complexity of calculations with increasing number and complexity of atoms characterizing the subspace investigated. This problem urges theoretical and computational chemists to look for efficient ways to drastically shorten the navigation time, as they are indeed supposed to take the lead in this kind of navigations due to the still lower cost of computational experiments as compared to synthetic work and so to assume their role as a guide for experimentalists. Ingenious alternatives to the brute force approach such as simulated annealing,4 genetic algorithms,5 linear combinations of atomic potentials,6 and the Inverse Molecular Design approach7,8 have been very promising. © XXXX American Chemical Society

These methods allow for the designing of rationally novel compounds with some desired characteristics.9−11 Another highly promising approach for a more efficient exploration of CCS, and closely related to the linear combinations of atomic potentials, is the alchemical coupling approach (AC).12−18 Within AC any two iso-electronic molecules in CCS can be coupled “alchemically” through the interpolation of their external potentials in the Hamiltonian. The energy derivatives can be obtained simply by virtue of the Hellmann−Feynman theorem. This approach has been used to study reaction energetics,15 to compute ab initio energy gradients, to predict orbital energies,14 and to estimate the change in covalent bond potentials in the most recent work.17 The alchemical derivatives,12,19,20 the underpinnings of the AC method, find a natural place in the context of Conceptual Density Functional Theory (CDFT)21,22 where the energy of a molecule is written as a functional of the number of electrons 5 and the external potential v(r) due to the nuclei. Upon perturbing the molecule in 5 and/or v(r) a perturbation expansion can be written in which in each term a derivative of the type ∂nδmE/ ∂5 nδv(r1)...δv(rm) acts as a coefficient of the perturbation, which Received: November 5, 2017 Published: January 4, 2018 A

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

substitution in this case is already 63, the gain in computing time vs the brute force approach is already 1 order of magnitude. As the number of BN substitutions is limited to 8 in pyrene as compared to 30 in C60, the expected increase in efficiency is definitely substantial. If successful, this approach could be used routinely in alchemical studies on BN substitution of carbon nanotubes and its graphene analogues.49−55 The choice was moreover made because systematic (timeconsuming) studies of the relative stability of the isomers of successive BN-substituted fullerenes C60−2n(BN)n appeared. Kar and co-workers continued their initial study56 for n = 1−7 in two ways: one with the 12 carbon atoms excluded from the substitution which leads finally to C12(BN)24 ball57 and another without imposing any constraints which results in C20(BN)20 belt.58 Upon analyzing the results, these authors were able to formulate a number of “rules“ governing the substitution pattern: (I) Carbon−Carbon Distance Rule. BN substitution always prefers double-bonded sites rather than single-bonded carbon pairs. Thus, the carbon pair of any hexagon−hexagon (doublebond type) junction of C60 is the first victim of BN substitution. The preference of forming B−N bonds instead of B−B or N−N bonds (in line with the bond energy values showing higher bond strength for B−N and C−C as compared to B−B and N−N59); (II) Hexagon-Filling Rule. BN units prefer to fill available, already partially mutated hexagons. This pattern increases the number of B−N bonds. Thus, once the first substitution takes place in a hexagon, the next two BN units replace the carbons of the same hexagon. Further substitution spreads to adjacent hexagons one by one; (III) N-Site Attachment Rule. The incoming BN group always prefers to attach to an existing BN unit via the BN−BN link. The other possibility NB−NB is less favorable; (IV) Continuity Rule. Joining existing BN chains spread over different hexagons significantly enhances the stability. This substitution pattern found in higher BN-fullerenes overrides the hexagonfilling rule. We will refer to these rules in the text as Rule I, II, etc. In view of the importance of a systematic evaluation of BNsubstitution patterns of fullerenes and nanotubes not only in the theoretical studies but also in the successful efforts in the synthesis of BN carbon nanotubes60−63 (note that the existence of Boron Nitride nanotubes was predicted in 1994 by Rubio et al.,64 only one year later followed by its synthesis by Chopra et al.65), the rules reported by Kar and co-workers56−58 based on successive substitution will be verified and scrutinized in the simultaneous substitution used in this work. The paper is organized as follows. In Section II, we present the theoretical and computational background. In Section III, the method of calculation is presented. Finally in Section IV the alchemical results are discussed. Naming conventions for the heteromonocycles and for the fused ring systems are in accordance with the IUPAC recommendations.66,67

can be identified as the response function of the system due to the perturbation, independent of the magnitude of the perturbation, and as such an intrinsic property of the system. Numerous studies appeared involving d5 and/or dv perturbations. Traditionally the former perturbation has been concerned primarily with the prediction of phenomena associated with electron transfer.21,22 The density response to the external potential change, the so-called linear response function, is a central example of the latter perturbation. Its physical and mathematical properties were studied,23−26 techniques for its numerical and analytical evaluation were developed,27,28 and its chemical importance was scrutinized,25,29−32 up to connections with molecular conductivity.33 The change in external potential is typically understood as a geometry change result,34−36 the case of the δv perturbation due to a nuclear charge change rarely being treated.12,13,20,37,38 Alternatively, the calculation involving the linear response function can be superseded by the calculation using Fukui-like functions, i.e. the derivatives of the density with respect to nuclear position or the nuclear charge. These derivatives can be calculated using the coupled perturbed (CP) scheme, which constitutes a powerful technique for the calculation of molecular response properties involving electric, magnetic, and geometric perturbations.39−41 This scheme was also successfully used in calculations of the relaxation effect due to the change in the electron number,42−47 the second and the third order derivatives with respect to nuclear charges.20,37 It is also applicable for systems where the effective core potential approximation is used.48 It was shown that on the basis of only a single SCF type calculation and the corresponding alchemical derivatives of the reference molecule, the CCS of first neighbors − implying changes in nuclear charge of ±1 − could be fully explored by simple arithmetic operations at negligible supplementary cost to a SCF calculation. As illustrative transmutations showing the potential of the CP method, some examples of increasing complexity were presented,37 starting with the deprotonation of small molecular systems with one and more hydrogen atoms, continuing with the transmutation of the nitrogen molecule, and ending with the substitution of isoelectronic (B,N) units for (C,C) units and N units for C−H units in carbocyclic systems. The overall trends observed for the alchemical deprotonation energy proved the usefulness of the alchemical indices as probe in chemical reactivity investigations. The results of calculations for the BN derivatives of benzene and pyrene show that this method has great potential for efficient and accurate scanning of chemical space. These successes prompted us to pursue our investigations on CC/BN substitutions, but now in the case of carbon, cage compounds for which C60, the archetypical fullerene, and by extension nanotube, will be considered in this paper. C60 shows a myriad of possible substitution patterns upon increasing n (the number of CC/BN substitutions) from 1 to 30, promoting it to an ideal benchmark system. Our aim is to see if the alchemical approach is able to predict qualitatively, starting from a single SCF calculation on the compound C60, the full transmutation route up to n = 30, using only a second order perturbation expansion involving only first and second order alchemical derivatives of the parent compound. The evaluation of the energy of all transmutation products involves only some simple, extremely fast arithmetic operations. In our previous study we mentioned that for azoborapyrene a ratio of 8 to 1 was found for the computing time of the alchemical derivatives vs the SCF part of the parent molecule. As the number of isomers for a single

II. THEORETICAL BACKGROUND II.A. Chemical Compound Space and Taylor Expansion in It. Every point of the chemical space is mathematically defined as a vector, (R,Z), in the M-dimensional conformational space (RA and ZA denote the A-th nuclear location and the corresponding nuclear charge, respectively) ⎧ ⎫ R = (R1, R 2, ..., R M), Z = (Z1, Z 2 , ..., ZM ) ⎪ ⎪ ⎬ *M ≡ ⎨(R, Z) ⎪ ⎪ ∀ A , B ∈ {1, 2, .., M }, A ≠ B : RA ≠ R B, ZA ≥ 0 ⎭ ⎩

(1) B

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

The transmutation energy can be decomposed also into the electronic and the nuclear parts:

which is associated with some electronic configuration characterized by the electron number 5 and some characteristics in the spin space, e.g. singlet, doublet, etc.37 To each point of the conformational space (R,Z) corresponds an external potential, defined as v(r) ≡ v(r, Z , R) = −∑ ZA /|r − RA| A

Dal [5, Z , R, d Z] = Dal,el [d Z] + Dal,nuc[d Z] ⎛ ⎞ 1 al,el = ⎜⎜∑ μAal,el dZA + ∑ ηAB dZAdZB⎟⎟ 2 A ,B ⎝ A ⎠

(2)

⎛ 1 + ⎜⎜∑ μAal,nuc dZA + 2 ⎝ A

and the nuclear−nuclear repulsion energy Vnn ≡ Vnn[Z , R] =

∑∑ A

B>A

ZAZB |RA − R B|

ρ

{

⎛ ∂W [5, Z , R] ⎞ μAal [5, Z , R] ≡ ⎜ = μAal,el + μAal,nuc ⎟ ∂ZA ⎠5, R ⎝ (8)

and it measures the transmutational tendency for the atom at position RA in the molecule. The first component, the electronic contribution to the electrostatic potential at site A (the nuclearelectron attraction energy per unit nuclear charge,70 also involved in the diamagnetic shielding of an atom in a molecule71), is due to the Hellmann−Feynman theorem

(4)

with Ẽ [ρ,v] ≡ F[ρ] + ∫ ρ(r)v(r)dr, the energy functional for a system moving in an external potential v(r). The minimizer, ρg, represents the ground state density of the system. F[ρ] is a universal functional, the sum of the kinetic energy and the electron−electron repulsion energy. The sum of the electronic energy E and the nuclear−nuclear repulsion energy Vnn (other external fields absent) is the total energy W of a system

⎛ ∂E ⎞ μAal,el [5, Z , R] ≡ ⎜ ⎟ ⎝ ∂ZA ⎠5, R

(5)

=

and consequently, the independent variables are the coordinates in the conformational space (R,Z) and the electron number, 5 . Moving from one point to another in the chemical space is characterized by the change in the electron number d5 , the change of the position vector dR = (dR1,dR2,...,dRM), and the change of the charge vector dZ = (dZ1,dZ2,...,dZM) − the socalled transmutation vector. So “the journey” in the chemical space can be expressed with the change vector (d5 , dZ, dR). Properties of the “end point journey” molecule can be efficiently estimated using Taylor series expansions (TSE), e.g. the total energy change, ΔW[d5 , dZ, dR]. In this paper, only the so-called alchemical transmutations characterized by the transmutation vector dZ at constant electron number (d5 = 0) and frozen geometry (dR = 0) will be considered. It was shown that the second order Taylor expansion based prediction preserves a realistic perspective from a chemical point of view, if the molecules at their equilibrium geometries are used as the basis for an expansion and the changes in the atomic charges are limited to dZi ∈ {−1,0,1}.20,37 The alchemical transmutation energy connected with the transmutation vector dZ is defined as Dal [5, Z , R , d Z] ≡

1 2

⎛

⎞ ⎛ ∂v(r) ⎞ ⎟ dr = − ⎜ ⎝ ∂ZA ⎠5, R 5

∫ ⎜⎝ δEδ[v5(r,)v] ⎟⎠

∫ |r ρ−(rR) | d r A

(9)

Here, the relation (δE[5 ,v]/δv(r))5 = ρ(r; 5 ,v) and the linear dependence of v(r) on the nuclear charge (∂v(r,ZM, RM)/∂ZA)RM = −1/|r − RA| are used. The second component, resulting from Vnn, is ⎛ ∂V ⎞ μAal,nuc [5, Z , R] ≡ ⎜ nn ⎟ = ⎝ ∂ZA ⎠5, R

∑ B≠A

ZB RAB

(10)

The second order derivative, the alchemical hardness12 (the effective “screened” potential) is ⎛ ∂ 2W ⎞ al al,el al,nuc al ηAB [5, Z , R] ≡ ⎜ = ηAB + ηAB = ηBA ⎟ ⎝ ∂ZA ∂ZB ⎠5, R (11)

It measures the transmutational resistivity for the atoms at the positions RA and RB in the molecule. The contribution from Vnn is simply ⎛ ∂Vnn ⎞ 1 al,nuc al,nuc ηAB = ηBA [Z , R] ≡ ⎜ ⎟ = (1 − δAB) RAB ⎝ ∂ZB∂ZA ⎠R

al dZAdZB = Dal, μ[d Z] + Dal, η[d Z] ∑ μAal dZA + ∑ ηAB A

(7)

The components in the last line will be abbreviated as D , Dal,η,el, Dal,μ,nuc, and Dal,η,nuc, respectively. If dZ = (0,...,dZA = 1,...,0) and there is a vacancy at the A-th position, ZA = 0, Dal is related to the proton affinity. II.B. First and Second Order Derivatives. The first order derivative, the alchemical potential of nucleus A at position RA,19,69 is

(3)

∫ ρ(r)d r = 5} = E[̃ ρg ; v]

W [5, Z , R] ≡ E[5, v[Z , R]] + Vnn[Z , R]

⎠

A ,B

al,μ,el

The external potential v(r) uniquely determines the electronic Hamiltonian in the second-quantization form, which allows, at least in principle, for solving the molecular Schrödinger equation. Thus, the potential function v(r) encodes all of the chemical information for a given electronic configuration. Due to the Hohenberg−Kohn variational principle,68 the ground state energy is given by E[5, v] = min E[̃ ρ ; v]

⎞

al,nuc dZAdZB⎟⎟ ∑ ηAB

A ,B

≈ ΔW [d 5 = 0, d Z , d R = 0]

(12) (6)

where RAB = |RA − RB|. The higher-order derivatives of Vnn are trivially equal to zero. The electronic contribution to the alchemical hardness is the effective “screened” potential defined by the linear response kernel

where Dal,μ and Dal,η are the contributions from the first and second derivatives, respectively. For simplicity, the dependence on (5 , Z, R) is suppressed on the right side of the equations. C

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation ⎛ ∂ 2E ⎞ al,el [5, Z , R] ≡ ⎜ = ηAB ⎟ ⎝ ∂ZB∂ZA ⎠5, R =

∫

∫

where v0(r) and v1(r) are the external potential of the reference and the target system, respectively. The coupling parameter λ goes from 0 to 1. The change of the nuclear−nuclear interaction energy in the AC approach is defined as

(∂ρ(r)/∂ZA)R dr |r − R B|

(∂ρ(r)/∂ZB)R al,el d r = ηBA |r − RA|

0 1 0 λ Vnn = V nn + λ(Vnn − V nn )

(13)

It can be rewritten as ⎛ ∂ E[5, v[Z , R]] ⎞ = ⎟ ⎜ ∂ZB∂ZA ⎠5, R ⎝ 2

=

⎛

⎞

∫ ∫ χ(r, r′)⎜⎝ ∂∂vZ(r) ⎟⎠ A

=

∫ ∫ ⎜⎝ δδρv((rr′)) ⎟⎠

⎛ ∂v(r′) ⎞ ⎜ ⎟ d rd r′ ⎝ ∂ZB ⎠5, R 5, R

∫ ∫ |r − Rχ(r||,r′r′−) R | d rd r′ A

Because the functional dependence of this energy component on the charge vector is different in the two compared methods, the physical interpretation of the Vnn derivatives is different. In the AC method, there exists only a first derivative over λ which is a simple difference between the Vnn energies at the “end points”, while in the method used here the first derivative over ZA is the nuclear contribution to the electrostatic potential at site A, eq 10; the second derivative, on the other hand, is the inverse of the distance between two atoms or zero, eq 12. It means through the scheme presented in this work, someone gets information about the influence of the environment and the distances (the first and second derivative, respectively) of the separate atoms involved in the alchemical change on the response due to perturbation. With fixed nuclear position RA (frozen geometry), the external potential of the AC method, eq 17, can be rewritten with the help of eq 2 as

⎞ ⎛ ∂v(r) ⎞ ⎛ ∂v(r′) ⎞ ⎟ d rd r′ ⎜ ⎟ ⎜ ⎝ ∂ZA ⎠5, R ⎝ ∂ZB ⎠5, R N

⎛

(14)

B

where χ(r,r′) = (δ E[5 ,v]/δv(r′)δv(r)))5 is the so-called linear response function mentioned in the Introduction. Its importance is easily recognized when written as δρ(r)/δv(r′)))5 , i.e., it gives the response of the density at a given point r when the system is perturbed in its external potential at another point r′, the situation at stake in any chemical reaction. It is well-known in its time dependent form in solid state physics. Its computation, visualization, and chemical interpretation29 were extensively explored in recent years (for a review see ref 25 and the Introduction). In this work, the calculation of the linear response function is superseded by the calculation of the Fukui-like function 2

⎛ ∂ρ(r) ⎞ ⎜ ⎟= ⎝ ∂? ⎠

⎛

⎞ ⎛ ∂v(r′) ⎞ ⎜ ⎟d r′ , ⎝ ∂? ⎠

∫ ⎜⎝ δδρv((rr′)) ⎟⎠

vλ(r) = −∑ ZAλ /|r − RA| ≡ v(r, Zλ , R)

of the reference and target systems. Then, for the nuclear contribution we have (see eq 3) ⎛ ∂V λ ⎞ ⎜⎜ nn ⎟⎟ ⎝ ∂λ ⎠

? ∈ {5, Z1, ..., ZM , R1, ..., RM }

(15)

1 0 = Vnn − V nn λ= 0

In taking the derivative of the density in HF or DFT, it is sufficient to take the derivatives of the spin orbitals (SO) and the derivatives of the SO occupation numbers: ⎛

=

⎝

p

0

⎛ ∂E[5, vλ] ⎞ ⎜ ⎟ ⎝ ⎠ ∂λ

(16)

In the alchemical derivatives calculation, the contribution from the SO occupation numbers can be neglected in the majority of cases (the ordering of the SOs is assumed to be unchanged), and the remaining difficulty is the calculation of the wave functions relaxation. Assuming that all difficulties in the evaluation of the basis set derivatives and the weight factors derivatives (due to use of the linear combination of the atomic orbitals method and numerical integration for the exchange-correlation part) have been overcome, the remaining difficulty is the calculation of the derivatives of the SO coefficients. These calculations (the relaxation effects to chemical reactivity indices) can be done effectively in the framework of the coupled perturbed (CP) scheme (see ref 20 for details). II.C. AC Method vs Method Used in This Work. Now, we want to highlight the main difference between the AC method12−18 and the method presented here.37 In the former two molecules are coupled alchemically through interpolation of their external potentials in the Hamiltonian vλ(r) = v0(r) + λ(v1(r) − v0(r))

A ,B

(1 − δAB) ((ZA0 + dZA)(ZB0 + dZB) − ZA0 ZB0) RAB

where dZ = Z − Z . It is easy to check that it equals D [dZ] = Dal,μ,nue[dZ] + Dal,η,nue[dZ]. The electronic energy in both methods is calculated in the same way, eq 4. The first derivative of it in the AC method is

∑⎜ p

∑

(20)

⎛ ∂Ψ*(r) ⎞ p ⎟ ⎟ ⎝ ∂? ⎠

⎛ ∂np ⎞ ⎟Ψp(r)Ψ*p(r) ⎝ ∂? ⎠

1 2

1

∑ np⎜⎜Ψp(r)⎜⎜

⎛ ∂Ψp(r) ⎞⎞ * ⎟⎟⎟ + + Ψ p(r)⎜⎜ ⎝ ∂? ⎠⎟⎠

(19)

A

where ZλA = Z0A + λ(Z1A − Z0A) with Z0 and Z1 denoting the charges

N

⎛ ∂ρ(r; {Ψp}, {np}) ⎞ ⎜⎜ ⎟⎟ = ∂? ⎝ ⎠

(18)

λ=0

al,nue

⎛ ∂E[5, v[Zλ , R]] ⎞ ⎟ =⎜ ∂λ ⎝ ⎠ = −∑ A

=

∫

λ=0

(Z1 − ZA0 ) dr ρλ = 0 (r) A |r − RA|

∑ (Z1A − ZA0 )μAal,el

(21)

A

al,μ,el

the same as D in eq 7 evaluated for dZA = second derivatives can be written formally as ⎛ ∂ 2E[5, v ] ⎞ λ ⎜ ⎟ 2 ∂ λ ⎝ ⎠

= λ= 0

Z1A

−

Z0A.

The

(Z1A − ZA0 )(ZB1 − ZB0) ∑ ∑ ηAal,el ,B A

B

(22) al,η,el

in eq 7 for dZA = Z1A − Z0A. In the method used in the

i.e., as 2D present paper, each row/column of the alchemical hardness matrix is calculated separately by means of the CP scheme (evaluations are needed only for symmetry unique atoms). In the AC method, if the CP scheme or the finite difference method is used to evaluate the second derivatives,17 such evaluation is to be performed separately for each alchemical “direction” Z1 − Z0 (a separate calculation ∂2E[5 ,vλ]/∂λ2 for each target system).

(17) D

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Table I. Alchemical Derivatives and Their Electronic and Nuclear Components for C60 at HF/STO-3Ga atom C atom j

first derivatives μali μal,el i −14.5312 −60.7323 second derivatives ηal1j ηal,el 1j

μal,nuc i 46.1920 ηal,nuc 1j

η1alj

connection type and used notation

al η1,9

1 9 2 7 8 3 25 23 20 24 22 21 41 40 38 39 37 36 58 54 59 53 52 60 a

−1.3019 0.3230 0.2602 0.1816 0.1671 0.1529 0.1487 0.1466 0.1437 0.1413 0.1392 0.1345 0.1303 0.1302 0.1266 0.1266 0.1245 0.1240 0.1222 0.1218 0.1214 0.1206 0.1200 0.1181

−1.3019 −0.0555 −0.1034 −0.0038 −0.0470 −0.0718 0.0010 0.0373 0.0007 0.0267 0.0107 0.0174 0.0431 0.0339 0.0250 0.0353 0.0266 0.0378 0.0433 0.0406 0.0453 0.0412 0.0441 0.0437

0.0000 0.3785 0.3636 0.1855 0.2141 0.2247 0.1478 0.1093 0.1430 0.1146 0.1284 0.1171 0.0871 0.0964 0.1016 0.0913 0.0979 0.0861 0.0789 0.0813 0.0760 0.0794 0.0759 0.0744

1.000 0.806 0.562 0.517 0.473 0.460 0.454 0.445 0.437 0.431 0.416 0.403 0.403 0.392 0.392 0.385 0.384 0.378 0.377 0.376 0.373 0.372 0.366

diagonal element of the η matrix hexagon−hexagon junction hexagon−pentagon junction para-position at the hexagon meta-position at the hexagon nonvicinal position at the pentagon

ηalCC ηalh‑h ηalh‑p ηalpara ηalmeta ηalnvp

Results ordered according to decreasing values of ηal1j, 1 > j. All values are in au.

Figure 1. Schlegel diagram of C60 and structures of 20-BN fullerene (belt structure) and 24-BN fullerene (ball structure). H and P represent hexagons and pentagons, respectively. Surfaces with corresponding primed and unprimed number are related by inversion with respect to the center of symmetry in C60.

E

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

It is clear that the contribution from the first term in eq 26 is constant for given n (ηAA = η11). After the examination of the values of the ηal matrix (diagonal elements are negative, offdiagonal elements are positive, see Table I), we can conclude that the first term is negative and constant for given n, the second term gives always a positive contribution to Dal (it increases the transmutation energy), whereas the last term (the heteroatoms contribution) has always a negative contribution to Dal (it decreases the transmutation energy). Based on this analysis, we can formulate the general alchemical BN substitution rule (shortly named the alchemical rule): at given n, maximize the contribution from the heteroatom part and minimize the contribution from the homoatom part. To verify the substitution rules formulated by Kar and coworkers,56−58 the calculation of all possible transmutation products should be performed. In the case of n = 1, it is easy, the 1770 distances (the number of distances possible out of the 60 carbons) comprise 23 different distance types, so there are 23 isomers. In the case of n = 2, there are 487 635 tetrahedrons composed of four different vertices. In each of these tetrahedrons vertices can be substituted in 6 ways, but because Dal[Nn, Bn] is a symmetric function, Dal[Nn,Bn] = Dal[Bn,Nn], the number of possible sequences reduce to 3. The use of the symmetry significantly reduces the number of possible isomers to 3574 isomers. Based on the results of the single BN substitution, that the hexagon−hexagon junction position is the most preferred site, the C56(BN)2 isomers should also be built by the stepwise substitution of C−C units at hexagon−hexagon junctions. As there are 30 such junctions, the number of considered tetrahedrons reduces to 29. Based on the C60 symmetry, they comprise only 8 different tetrahedrons, and the number of potentially stable isomers is 24. The number of possible isomers of the successive values of n increases very sharply: the calculation of all possible isomers − the force method − is impractical, and its number is nontrivial to predict on the basis of symmetry rules. To simplify calculations, we proposed to replace the force method by the k-memory method, in which for n substitutions, no more than the k energetically lowest isomers from the preceding, n − 1 substitutions, were used as the starting point. Let us explain this method through the example for k = 200 and n = 4. There are 23 isomers of C58BN. For each isomer of the n = 2 substitution, all possible isomers of the following substitution n = 3 are calculated. After this step, no more than the k energetically lowest isomers of C54(BN)3 are saved for the next substitution, n = 4. To check the influence of the k value on the obtained structure of the most stable isomer for a given n, several values of k were used. For k = 200 and k = 500, the result was the same, so we can assume that the results from the 200memory method are “memory f ree“ results. To verify the influence of the geometry relaxation and to confront our alchemical results with the previous Kar and coworkers investigations,56−58 the geometry of n-BN fullerenes was fully optimized at the B3LYP/3-21G level (used in ref 58). As Dal in our case is a symmetric function with respect to the replacement of Nn with Bn, Dal[Nn,Bn] = Dal[Bn, Nn], two structures are generated for each isomer, the original structure for dZn = {Nn,Bn} and its mirror structure for dZn = {Bn,Nn}. To distinguish them, the second term in eq 26 was split into the nitrogen−nitrogen interaction and the boron−boron interaction terms

Their number may be huge if the number of transmuted atoms is significant, e.g. in the case of C60, there is only one symmetry unique atom, while there are 23 possible C58BN isomers (the target systems). Note that after completing the calculation for the alchemical derivatives for atoms, one is able to make the prediction for any target system (with the same accuracy as the AC method) and the analysis of the contributions due to the change of the nuclear charge at particular atomic sites (not possible in the AC method).

III. METHOD OF CALCULATION The main idea of this work was to show the usefulness of the alchemical derivatives as a tool for exploring chemical space. This is illustrated by a qualitative prediction of the most stable isomers of BN-substituted fullerenes C60−2n(BN)n, where n = 1−30, based on a single SCF calculation of the C60 molecule. The experimental geometry of the prototype fullerene from the gasphase electron diffraction was used.72 The Cartesian coordinates for each C atom in the cage were obtained analytically with the procedure of ref 73, and all calculations were carried out using the locally modified version of the GAMESS program package74 at the HF/STO-3G level of theory. If the structure has icosahedral symmetry, the geometry of the molecule is defined by only two parameters, which for convenience we chose as the bond length in the pentagons, r5, and the length of bond connecting pentagons, r6. C60 at Ih symmetry is the minimal arrangement such that no two pentagons share an edge (isolated pentagons): the edges of each pentagon join only hexagons, and the edges of each hexagon alternately join pentagons and hexagons. All carbons are equivalent by the symmetry, each carbon (vertex) is double bonded to one carbon (hexagon−hexagon junction, h-h) and sigma bonded to two other carbons (hexagon−pentagon junction, h-p). The alchemical information (alchemical derivatives and their electronic and nuclear components) is presented in Table I. The carbon atoms are numbered in accordance with the IUPAC recommendations75 as it is illustrated via the Schlegel diagram in Figure 1. In case of the (C−C)n → (B−N)n transmutation, the set of the transmutation vectors is ‐ 2n(BN)n + CC60 60

⎧ = ⎨d Z = (dZ1, ..., dZ60)|dZA ∈ {− 1, 0, 1}, ⎩ ⎪

⎪

60

⎫

60

⎪

∑ dZA = 0, ∑ |dZA| = 2n⎬ ⎭

⎪

A=1

A=1

(23)

The first order terms in the Taylor expansion of the transmutation energy cancel, eq 7 (due to the symmetry, μalA is independent of A), and the transmutation energy is Dal [5, Z , R, d Z] =

1 2

al dZAdZB ∑ ηAB

(24)

A ,B

To simplify the analysis of the above equation, the dZ will be replaced by two vectors N n = {(N1, ..., Nn)|Ni < Nj , dZNi = 1}, Bn = {(B1, ..., Bn)|Bi < Bj , dZBi = −1}

(25)

so eq 24 can be rewritten as n

Dal [N n, Bn] = nη11al +

n

n

n

∑ ∑ (ηNaliNj + ηBaliBj) − ∑ ∑ ηNaliBj i=1 j>i

i=1 j=1

(26) F

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Table II. Transmutation Energies of BN-Fullerenes from the Alchemical Calculation at the HF/STO-3G Level and after Geometry Optimization at the B3LYP/3-21G Levela alchemical prediction HF/STO-3G n

ΔDball n

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

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 −28.24 0.00 −20.08 0.00 0.00 0.00 0.00 −43.30 −105.42 −126.76 −198.92 −195.78 −141.82 −69.03 0.00 0.00 0.00

i

ΔDball n

1 1 1 1 1 1 1 1 1 1 1 1 3 1 14 1 1 1 1 7 23 43 71 90 91 78 1 1 1

0.00 0.00 0.00 0.00 0.00 0.00 0.00 −58.73 −59.49 −117.72 −60.18 −112.57 −69.78 −36.83 −2.76 −69.53 −77.50 −91.24 −114.46 −73.23 −51.27 0.00 0.00

B3LYP/3-21G i

ΔD̃ ball n

ΔEsimult n

type

i

1 1 1 1 1 1 1 3 11 54 20 86 25 4 2 22 11 8 8 14 4 1 1

0.00 0.00 0.00 0.00 0.00 0.00 0.00 −15.06 −23.53 27.23 −47.88 −64.51 −145.39 −222.39 −274.47 −185.05 −95.63 −78.63 −59.61 11.42 1.19 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 −0.25 −1.44 0.00 −1.26 0.00 −2.26 0.00 0.00 0.00 −2.38 0.00 0.00 0.00 0.00 −0.38 −0.94 0.00 −2.64 −2.76

eq eq N eq N eq N eq N eq N N eq N eq eq eq eq eq N eq eq eq eq eq eq N eq eq

1 1 1 1 1 1 1 1 1 1 6 2 1 2 1 3 1 1 1 3 1 1 1 1 5 6 1 3 2

a belt ΔDbelt n and ΔDn are the difference of transmutation energies between the simultaneous and successive substitution methods in the belt case and − Dsucc,belt and Dsimult,ball − Dsucc,ball , respectively. The index i means the i-th position of the considered isomer in the decreasing the ball case, Dsimult n n n n stability series based on the simultaneous alchemical prediction. ΔD̃ ball n is the difference of transmutation energies between the successive substitution sequence obtained in this work from the alchemical prediction and the one proposed by Kar and co-workers for the ball case. The differences of = Ein − E1n. Type means the most stable isomer pattern after the geometry transmutation energies at the B3LYP/3-21G level are defined as ΔEsimult n optimization. All values are in kcal/mol.

n

Dnal, NN =

n

i=1 j>i

Dal,NN n

n

n

As it was discussed above, a simple brute force technique is not applicable for such problems in view of the fact that computational timings become prohibitively large. In this work, the 200-memory method is assumed as a memory free method (the results obtained in this way are denoted by superscript “simult”). The successive substitution used by Kar and co-workers is equivalent to the 1-memory method introduced in this work (denoted as “succ belt”). The ability of these methods to qualitatively predict the most stable structural arrangement of BN-fullerenes is presented in Table II (see column labeled by ΔDbelt = Dsimult − Dsucc,belt ; a negative sign n n n indicates that our simultaneous method yields the most stable isomer) and in Figure 2. The crucial conclusion is that the most stable isomer is not always simply a combination of the incoming BN unit with the preceding most stable isomer. Comparing the alchemical prediction at the HF/STO-3G level, the successive substitution method which leads to the belt structure is inconsistent with the simultaneous substitution method for n = 14, 16. Passing beyond the belt system (n > 20) the successive substitution method fails almost systematically. Moreover, it is unable to reproduce the ball structure for n = 24 as the most stable structure. The

∑ ∑ ηNaliNj , Dnal, BB = ∑ ∑ ηBaliBj i=1 j>i

(27)

Dal,BB n ,

If = two structures are equivalent by symmetry, only one geometry optimization was done, and the structure was labeled “eq”. If Dal,NN < Dal,BB (Dal,NN > Dal,BB n n n n ), the structure was labeled “N” (“B”), where its mirror structure was labeled “B” (“N”). The geometry optimization was done for both of them.

IV. RESULTS AND DISCUSSION IV.A. Overall Picture. As it was mentioned, the main objective of this work is to show the qualitative accuracy of the alchemical method using the BN substitution of fullerenes as an illustrative example. The methodology used in this work assumes that all substitutions are made at the same time, simultaneously. We will refer to this method as simultaneous substitution. The results so obtained will be confronted with the results reported by Kar and co-workers56−58 As mentioned above, their initial study56 for n = 1 − 7 was continued in two ways: one with 12 carbon atoms excluded from the substitution which leads finally to C12(BN)24 ball57 and the second one without imposing any constraints which results in the C20(BN)20 belt.58 G

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 2. Difference of the transmutation energies (in kcal/mol) for the most stable isomer from simultaneous and successive substitution methods in the belt and ball cases.

positions. The 1,9-azaborafullerene is more stable than other BNfullerenes in which the boron and nitrogen atom are located on one and the same face of C60 (hexagonal or pentagonal). The most unstable isomer is the 1,60-azaborafullerene: the N and B atoms are located at the positions which are related by the inversion with respect to the center of symmetry in C60. It should be noted that the stability of C58BN does not regularly increase with the decrease of the distance between the heteroatoms, e.g. al,nuc al,nuc al,nuc ηalh‑h > ηalh‑p > ηalpara > ηalmeta > ηalnvp, but ηal,nuc h‑h > ηh‑p > ηnvp > ηmeta > al,nuc ηpara (see notation in Table I). The single substitution results lead to the conclusion that the most stable isomer of C56(BN)2 is to be built by the stepwise substitution of a C−C unit at a hexagon−hexagon junction. In the case n = 2, the last term in eq 26 is a combination of four components, and the second term is a sum of the nitrogen− nitrogen alchemical hardness and the boron−boron alchemical hardness. To maximize the contribution from the heteroatom part to Dal, the two other largest elements in the sequence of Table I should be considered for locating the second BN entity: the hardness for the h-p junction and the para-position at the hexagon. Therefore, the next BN unit should also be located in the same hexagon, e.g. 1,6,5,9-diazadiborafullerene (Figure 3, 2-1 structure). The transmutation energy is

weakness and the limitations of the successive substitution method are unequivocal, particularly when some carbons are excluded from replacement. The differences between the transmutation energies predicted by both methods in such a simult,ball case are presented in Figure 2 and in Table II as ΔDball n = Dn succ,ball − Dn . The successive substitution works correctly only for n ≤ 8 and for the ball structure for n = 24 (and of course its preceding structure). Trustworthy atomic arrangements for n > 24 are impossible to obtain with the alchemical method, but it provides the possibility to point out the small set of potential candidates for the most stable structure, as it is clearly shown in the last column of Table II (after the geometry optimization). After geometry optimization of the six most stable isomers predicted by the alchemical method at the B3LYP/3-21G level (the results presented in Table II (B3LYP/3-21G part), the is less than −3 kcal/mol. In maximal energy difference ΔEsimult n most cases, the alchemical predictions were consistent with the results obtained after geometry optimization (see last column in Table II). The results presented in the B3LYP/3-21G part of Table II validate our conclusion from previous work37 on the usefulness of HF/STO-3G alchemical predictions as directive results. Of course, the absolute magnitudes of the energies or other properties strongly depend on the calculation level, but their trends and qualitative relations are similar; but the main object of this work is to show the qualitative correctness of the alchemical method. In the remainder part of this paper, the results obtained from the alchemical methods will be analyzed in more detail, and the successive substitution rules will be verified in the case of the simultaneous substitutions. IV.B. Verification of Rules I and II. Based on the discussion above eq 26, Rule I is very easy to validate. Based on the alchemical derivatives the transmutation energy, eq 26, is Dal[(1), (j)] = ηal11 − ηal1j for the single BN-substituted C60 molecule, where the nitrogen atom substitutes the carbon at C1 and the boron atom substitutes the carbon at the jth position (N1 = (1) and Bn = (j)). As the off-diagonal elements of the ηal matrix are positive in the C60 case, the transmutation energy of the substitution at the hexagon−hexagon junction is lower than those of the other

Figure 3. Structures of 2-BN and 3-BN fullerenes. H

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

This result is in line with the alchemical rule (maximization of the contribution from the heteroatom part). Also the hexagonalfilling rule can be justified based on these two cases. In case of the 3-1 structure, the heteroatom part contains ηalh‑h, ηalh‑p, and ηalpara elements, each of them three times. In the 3-2 structure, one ηalh‑p and one ηalpara element are replaced by ηal1.20 and ηal1,23, and the chain structure (3-2) is less stable then the cyclic pattern (3−1). Two others isomers, 3-3 and 3-4, represent a branch structure. These structures are very similar, in both the nitrogen/boron atom is surrounded by three boron/nitrogen atoms. The homoatoms parts are equal for both isomers, 2ηalnvp = 4ηalmeta. The main difference is that in 3-3, the surrounded BN pair is located at the h-p junction, and in 3-4, the surrounded BN pair is located at the h-h junction. This results in 3-3 containing two h-h and three h-p junctions and 3-4 containing one h-h and four h-p junctions. This difference causes 3-4 to be less stable than 3-3. The part above presented evidence for rehashing two rules through the alchemical derivative based method: the hexagon− hexagon junction (double-bond-type) as the preferable position for the replacement by BN pairs (Rule I) and the hexagon-filing rule (Rule II). Further analysis will be less detailed, due to the increase in the number of mutual links between atoms and for clarity of discussion. The four most stable isomers for n = 4−8 after geometry optimization are presented in Figure 4. Rule II, the hexagon-filing rule, is clearly manifested in the case of n = 5, 7; the structures with completely filled hexagons are the most stable, e.g. 5-1, 7-1. The former rule is necessary but not sufficient. We will illustrate this on the example of the 6-BN fullerenes. Starting from the most stable 5-BN fullerene (5-1, the

al al al al al al al Dal [(1, 6), (5, 9)] = 2η1,1 + (η1,6 + η5,9 ) − (η1,5 + η1,9 + η5,6 + η6,9 ) al al al = 2ηCC + 2ηmeta − (2ηhal‐ h + ηhal‐ p + ηpara ) ≡ D2al‐ 1

(28)

This result is an analog of the most stable diazaborines.37 The very similar transmutant, 2-2 (1,8,5,9-diazadiborafullerene) is less stable than 2-1, and its transmutation energy is al al al al al al al Dal [(1, 8), (5, 9)] = 2η1,1 + (η1,8 + η5,9 ) − (η1,5 + η1,9 + η5,8 + η8,9 ) al al al = 2ηCC + 2ηmeta − (ηhal‐ h + 2ηhal‐ p + ηpara ) ≡ D2al‐ 2

(29)

The difference in energy between the two structures is − Dal2‑1 = ηalh‑h − ηalh‑p = 0.0628 au − it is due to the replacement of one h-h interaction in 2-1 by one h-p interaction in 2-2. Other isomers where the second BN group is attached to the first one and located in any of the h-p junctions (Figure 3, 2-3 structure), e.g. the 8-25 junction, are less stable Dal2‑2

al al al Dal [(1, 8), (9, 25)] = 2ηCC + (ηmeta + ηnvp ) al − (ηhal‐ h + 2ηhal‐ p + η1,25 ) ≡ D2al‐ 3

(30)

The other example is 1,3,2,4-diazadiborafullerene in which all substituted vertices are located on the pentagon face (Figure 3, 24 structure). The transmutation energy is al al Dal [(1, 3), (2, 4)] = 2ηCC + ηnvp − 3ηhal‐ p ≡ D2al‐ 4

(31)

To see the importance of the h-h junction substitutions, we introduce the h-h junction index, which is defined as n

n

w[N , B ] =

n

n

∑∑ i=1 j=1

ηNali Bj ηhal‐ h

n

−

⎛ η al

n

∑ ∑ ⎜⎜ i=1 j>i

⎝

NiNj ηhal‐ h

+

ηBali Bj ⎞ ⎟ ≡ w lab ηhal‐ h ⎟⎠ (32)

ηali,j/ηalh‑h

(the values are listed in Table I, “lab” means a structure index in the figures), and the transmutation energy can be rewritten in terms of ηalCC and ηalh‑h as ⎛

al

n

n

n

n

D [N , B ] =

−

∑∑ i=1 j=1

al nηCC

−

⎛ η al ηBali Bj ⎞ NiNj ⎜ ∑ ∑ ⎜ al + al ⎟⎟ ηh ‐ h ⎠ ⎝ i = 1 j > i ⎝ ηh ‐ h

ηhal‐ h⎜ ⎜

n

n

ηNali Bj ⎞ ⎟ = nη al − η al w[Nn, Bn] CC h‐h al ⎟ ηh ‐ h ⎠

(33)

Comparing this index for the structures discussed above, we have w 2 ‐ 1 = 2.33,

w 2 ‐ 2 = 2.14,

w 2 ‐ 3 = 2.08,

w 2 ‐ 4 = 1.94

(34)

and for the 1,9-azaborafullerene, w[(1),(9)] = 1. These numbers indicate that the number of the h-h junction substitutions is crucial for the isomer stability. The decreasing number of the h-h junction substitutions, from 2 in 2-1 to 0 in 2-4, lowers the value of the w index, which means decreasing stability thus confirming Rule I. Rule II can be justified in similar manner. Four possible isomers of 3-BN fullerene are presented in Figure 3. In case of the 3-1 and 3-2 structures, the number of the h-h junction substitutions is maximal (three). The most stable isomer places all six heteroatoms in the same ring. In this case, the number of the h-p junction substitutions is three compared to two in 3-2; this changes the w index from 4.00 for isomer 3-1 to 3.70 for 3-2.

Figure 4. Structures of the 4-, 5-, 6-, 7-, and 8-BN fullerenes. The second number in the labels indicates the position in the decreasing stability series after geometry optimization. I

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation H1 and H2 rings are filled), the replacement of the carbon pair in 22-23 and 4-18 (h-h junction) is more favorable than in 22-21 and 4-3 (h-p junction). The next BN pair can be attached to four possible heteroatoms: B5, N6, B7, and N8, (symmetry equivalent to B2, N12, B11, N10) to form 6-4, 6-1, 6-3, and 6-2, respectively. It should be noted that 6-3 and 6-4 are the mirror structure (Bstructure) of 6-1 and 6-2 (N-structure), respectively. Short analysis of these structures shows that the difference in the h-h junction index is very small, w6‑2 − w6‑1 = 0.02. Based on the h-h junction rule, we are not able to justify why 6-1 is more stable than 6-2. The number of saturated nitrogen or boron atoms (surrounded by three boron or nitrogen atoms) is the same in both structures, resulting in the same number of C−C, C−B, C− N, and B−N bonds. This means that the most stable structure cannot be determined by a simple bond counting rule as it was proposed in ref 59. Recalling our discussion for 2-BN fullerenes, we can notice that the incoming BN pair creates the 2-1- and 2-3like structures with remaining heteroatoms in both cases, but in the case of 6-2, the 2-4-like structure (the pentagon with only one carbon atom) appears. The difference in energy between the two structures is only

Figure 5. Isomers of 9-BN fullerene. 9-1 and 9-2 are the first two most stable isomers by the simultaneous method. The last two isomers 9-3 and 9-4 are the most stable isomer by the successive method in the ball case, from this work and proposed by Kar and co-workers, respectively.

denoted as 9-4 at Figure 5 was found as the most stable isomer by the successive method in the ball case. In this work, the most stable 9-BN fullerene based on the alchemical method was found to be 9-3, the product of adding a BN pair to B13 of 8-1 (the replacement of the carbon pair at positions 30-31). Further successive substitutions definitely yield different structures from those reported by Kar and co-workers until n = 23 (see Figure 6,

al al al al al D6al‐ 2 − D6al‐ 1 = −ηnvp + 2η1,25 − η1,24 + η1,21 − 2η1,38 al + η1,39 = 0.0628 au

(35)

Note the significant contribution from the alchemical hardness between not bonded atoms (η1j, j > 9). This means that Rule II needs some reformulation; the incoming BN prefers to maximize the number of consecutive filled hexagons while minimizing the number of pentagons with only one carbon atom. We will denote this rule as Rule IIa. This extension will be deeply discussed for the case of n = 12-16 substitutions. Note that 7-1 can be considered as a product of the successive substitution of the BN pair to the 6-BN fullerenes presented in Figure 4. As in the case of the 5-BN fullerenes, the most stable structure is a structure with maximum number of completely filled hexagons. The four most stable structures of the 8-BN fullerene are combinations of 7-1 and the incoming BN pair. Their ordering can be explained in the same way as in the case of the 6-BN fullerenes. The number of filled hexagons is 3 for 8-1 and 8-2, but the number of pentagons with one carbon is 1 in 8-1 and 2 in 8-2. IV.C. Verification of the N-Site Attachment Rule and the Continuity Rule. Based on the results for n = 1−7,56 Kar and co-workers found a preference for N-site attachment (NBN) of the incoming BN group to the existing BN units (Rule III). Keeping this substitution pattern, they proceed to replace with some restriction carbon pairs to form C12(BN)24, the BNC ball,57 or without imposing any constraints to form C20(BN)20, the BNC belt.58 As it was mentioned, if the successive substitution has to yield the ball structure at n = 24, the 12 concerned carbon atoms have to be excluded from the substitution. Without this restriction, the successive substitution process leads to the belt structure at n = 20. The 8-1 fullerene is at a crossroads for the successive substitution path. The incoming BN pairs can continue to fill a fourth hexagon to form 9-1 (see Figure 5), a starting structure of the successive substitution process which leads to the belt structure. With the restriction of fixing the position of some carbons pairs and Rule III, the next substitution can take place at positions 41-42 (see Figure 5, the incoming group spreads to the adjacent hexagon) or replacement of the carbon pair attached to N12 or N16 of 8-1. The structure

Figure 6. Schlegel diagram of 24-fullerene (C12(BN)24, the ball structure). The green and red numbers indicate the order of BN substitution in the sequence proposed by Kar and co-workers and obtained from alchemical calculation. The first eight substitutions are marked by H3, H4, and H5 labels for fully substituted hexagons and by the black number 8. No number is assigned to the 12 carbon atoms that are blocked for replacement.

where green and red numbers do not coincide after n = 8). In Table II, the difference of the transmutation energies between the most stable transmutant obtained in this work and in ref 57 is presented, ΔD̃ ball n (a negative sign indicates the structure found in this work is more stable than proposed by Kar and co-workers). It is worth paying attention to ΔD̃ ball n . For n = 11, 21, 22, this value is positive; for others it is negative or zero, with the biggest difference at n = 16. This observation confirms the validation of Rule IIa. The most stable 11-BN fullerene from the successive J

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation method in this work is 11-2 with three filled hexagons (see Figure 7), which is less stable than 11-3 with four completely filled

substitutions). BN chains which are spread over different hexagons are not observed in the simultaneous substitutions, as was proposed by Kar and co-workers in Rule IV. In its place, the following rule is proposed: the most stable structure of the higher BN fullerenes is always some combination of phenanthrene-, triphenylene-, or chrysene-like structures (7-1, 9-2, and 9-1, respectively) which are connected with each other directly or by one BN unit. IV.D. Beyond n = 24. For n > 24, further substitution which follows Rule 1 creates two pentagons and two filled hexagons (Rule IV is no more applicable). The ball structure contains two triphenylene-like structures, and Rule IV does not seem to yield any additional information. The pentagons without carbons are energetically not favorable, and they have at least one N−N or B−B bond. To formulate some rules for the path from C12(BN)24 to (BN)30, the transmutation energy can be split into the contribution from the individual pentagons J, and eq 24 is rewritten as

Figure 7. Most stable isomers of 11-BN and 16-BN fullerenes by the simultaneous substitution method (with label 1) and by the successive substitution method (with labels 2 and 3) for the ball case. The isomers with label 2 are from this work; those with label 3 were proposed by Kar and co-workers.

Dal [5, Z , R, d Z] =

hexagons, reported by Kar and co-workers, and ΔD̃ ball n is positive; but both structures are less stable than 11-1 with five filled fused rings, the most stable 11-BN fullerene from the simultaneous substitution. The biggest difference in energy between the Kar and coworkers prediction and ours is observed in case of n = 16. 16-1 with seven rings (a phenathrene-like structure connected directly with a triphenylene-like structure) is the most stable 16-BN fullerene. The most stable isomer for the successive substitutions method for the ball case from this work is 16-2 with six rings (two phenathrene-like structures connected by two BN units). The isomer 16-3, reported by Kar and co-workers as the most stable isomers, has only four rings. Returning to the 9-BN fullerenes, the four isomers of the 9-BN fullerenes in decreasing stability order are presented in Figure 5. Structures 9-1 and 9-2 are the first two most stable isomers by the simultaneous method. The last two isomers 9-3 and 9-4 are the most stable isomer from this work and proposed by Kar and coworkers,57 respectively (both obtained by the successive method). Based on Rule IIa, 9-1 is more stable than 9-2; both have the same number of filled hexagons, but 9-1 has fewer pentagons with one carbon than 9-2. The 9-3 structure was found in our work as the most stable 9-BN fullerene based on the alchemical method. The geometry optimization at B3LYP/321G levels confirmed that it is more stable than 9-4. The results obtained in this work evidently indicate that the N-site attachment rule, Rule III in the form as it was proposed by Kar and co-workers, is incorrect and yields a wrong substitution pattern (see ΔD̃ ball n values in Table II). The detailed analysis of the successive substitution patterns where some restrictions (positions of six carbon pairs are kept fixed during the substitution process) have been imposed to reach the CBN ball shows that the continuity rule (joining existing BN chains spread over different hexagons significantly enhances the stability) is not valid. For substitutions with ΔD̃ ball n ≤ 0, the obtained structures were always some combination of phenanthrene-, triphenylene-, or chrysene-like structures (7-1, 92, and 9-1, respectively) connected with each other directly, sometimes by one BN unit. The BN chains were observed only at n = 11, 21, 22, where ΔD̃ ball n was positive. Concluding this part, in the case when Dal,NN ≠ Dal,BB n n , the most stable structure is the N type structure (this rule reformulates Rule III: the N-site attachment rule is not valid for simultaneous

∑ DJal[5, Z, R, d ZJ] J∈7

+

1 2

∑ ∑

al [5, Z , R, d ZJ , d ZK ] DJK

J∈7 K∈7

(36)

where 7 is a set of pentagons with label P1,...,P6, P′1,...,P′6. As is shown in Figure 1 (green labels), DalJ is the contribution from the atoms located on the same pentagon (A,B ∈ J), and DalJk is the contribution from the interaction between atoms located on two pentagons (A ∈ J, B ∈ K): DJal [5, Z , R, d ZJ ] =

1 2

∑

al ηAB dZAdZB ,

A ,B∈J

al DJK [5, Z , R, d ZJ , d ZK ] =

al dZAdZB ∑ ∑ ηAB A∈J B∈J

(37)

In the case n = 30, the values of DalJ are between −3.569 au for dZJ = {1,−1,1,−1,1} or dZJ = {−1,1,−1,1,−1} (see P1 or P2 pentagons in Figure 9, respectively) and −1.189 au for the homoatomic pentagons. The DalJK values range from −3.834 au up to 3.834 au. The extremes occur for two connected pentagons, the minimum when one contains only nitrogen atoms, while the second one only boron atoms, the maximum when both pentagons contain the same atoms. In the case of homoatomic pentagons (6 with only nitrogens and 6 with only borons), the most stable structure is when pentagons with the same type of atoms form a chain, e.g. the carbons in P2, P1, P5, P′3, P′1, P′6 are replaced by nitrogens and other carbons by borons. The most unstable structure is when one pentagon is surrounded by the same type of pentagons (with the same kind of atoms), e.g. all atoms in P1,...,P6 are nitrogens. The same analysis can be done for the case when in the pentagon four atoms are nitrogen and one is boron. This analysis shows that the minimalization of the first sum in eq 36 is more important than the minimalization of the second sum (interaction between pentagons) in the process of finding the most stable structure, the lowest Dal. The minimum value for DalJ occurs when there is only one homoatom bond in the pentagon. The number of all possible dZJ for such pentagons is 10 (two vectors with homoatoms at two initial positions, dZJ = {±1,±1,∓1,±1,∓1} and their cyclic permutations from 1 to 4 places to the right). For all of them DalJ = (5ηalCC − 3ηalh‑p + ηalnvp)/2 = −3.569 au, leading to the first sum in eq 36 to be equal to −42.823 au. In the second summation in eq 36 it is found that K

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

the directly connected pentagons. Four pairs of such pentagons with the lowest values of DalJK (from −0.385 au to −0.301 au) are presented in the bottom part of Figure 9. These structures indicate that pentagons in the most stable (BN)30 have to be connected by a B−N bond; it is in line with Rule I. We can have at most only 6 structures A shown in Figure 9. The other 34 possible pairs cannot be A type, which makes the answer to the question “which structure of (BN)30 is most stable” not trivial, but we can formulate a procedure which can verify the (BN)30 results obtained by the 200-memory method used in this work. This procedure is based on two rules: two pentagons are connected by a B−N bond (Rule I), and the completely substituted pentagons have only one homoatoms bond. This last one rule will be denoted Rule V. All possible isomers which fulfill these rules can be constructed in the following way: (i) let us assume that pentagon P1 after replacement has a nitrogen at positions 1,2,4 and boron at positions 3,5, dZP1 = {−1,1,−1,1,− 1}; (ii) replace all remaining carbons at the h-h junction connected to the P1 pentagon according to Rule I; (iii) take the next pentagon with carbon atoms which is connected to completely substituted pentagons and replace carbons by nitrogen and boron according to Rule V (only one homoatoms bond); (iv) repeat step (ii) for this pentagon; (v) repeat steps (iii) and (iv) for all remaining pentagons. After step (iii) and (iv), we have 5 possible structures for P1 and P2 pentagons connected by a B−N bond, all of them only with one homoatomic bond per pentagon. Repeating step (iii) and (iv) for the P3 pentagon yields 17 possible structures for P1, P2, and P3 pentagons. Finally, 1600 structures are found which fulfill Rule I and Rule V (note that in a very naive approach the transmutation energy is calculated for all possible 30 unordered outcomes from 60 possible positions, i.e. for 1.18 × 1018 structures; using Rule I, the h-h junction should be a B−N bond, which reduces this problem to 230 ≈ 109 structures). The obtained results are identical to those obtained with the 200-memory method. Further examination of the results after the geometry optimization confirms our opinion that the alchemical method is able to point out the small set of potential candidates for the most stable structure also in the case of extremely substituted fullerenes. In Figure 9, the most stable 30BN fullerenes obtained from the alchemical prediction and after geometry optimization are presented. The structure predicted as the most stable by the alchemical method, 30-1 in Figure 9, has six hexagons with a N−N−B−B chain (white hexagons, structure A), 4 structures B (e.g., P3 and P4), and two type C and type D structures (e.g., P3 and P6′ ). After the geometry optimization, the most stable structure is 30-2. It has also six structures A, but only 3 structures B, and no structures of types C and D. The most visible difference between these structures is the type of pair on the pentagons faces which are related by the symmetry center in C60. In 30-1 these pentagons are related by the symmetry center, e.g. P2 and P′2, and in 30-2, the P2 pentagon has a B−B bond, whereas the P′2 pentagon has a N−N bond shifted by one place to left. Concluding this part, the alchemical method for n > 24 can also correctly point out a small set of isomers among which the most stable structure is present (see Table II, final column) IV.E. Qualitative Correctness of the Alchemical Method at the HF/STO-3G Level. From the overall results both for n ≤ 24 and beyond n = 24, the results presented in Table II (especially the B3LYP/3-21G part) lead to the conclusion that the HF/STO-3G alchemical predictions are powerful as directive preliminary results. After geometry optimization, the natural weakness of the alchemical methods (if the position of nitrogen

Figure 8. Most stable n-BN fullerenes after geometry optimization with label n-1 for n = 12−16. 12-2 is the most stable isomer obtained from the simultaneous alchemical prediction. Two types of pentagons with one carbon are highlighted: the first one fused with three filled hexagons (marked by green circle) and the second one fused with two filled hexagons (marked by green empty circle) are presented. In the text these patterns are referred to as the first and second type pentagon, respectively.

Figure 9. Most stable 30-BN fullerenes obtained from the alchemical prediction (left) and after geometry optimization (right). The pentagons with a N−N bond are in blue; those with a B−B bond are in pink. The hexagons without homoatoms bonds are in gray (upper part). Pairs of directly connected pentagons with the lowest values of DalJK (−0.385, −0.310, −0.3077, and −0.301 au for A, B, C, and D, respectively) are shown.

there are 66 pairs of pentagons: 40 pairs with two pentagons connected directly, e.g. P1 and P2 pentagons; 6 pairs with two pentagons related by inversion with respect to the center of symmetry in C60, e.g. P1 and P1′ ; and 20 pairs with two pentagons not directly connected and not related by the symmetry, e.g. P1 and P′2. For a given type of pairs, the lowest value of DalJK is −0.3856 au for the first type of pairs listed above; −0.1295 au for the second type; and −0.1389 au for the last type. The main finding from our analysis is as follows: minimize the DalJK value for L

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

pentagons, and beyond n = 24 the completely substituted pentagons have only one homoatom bond (Rule V).

atoms and the boron atoms are interchanged, eq 26 yields exactly the same transmutation energy at frozen geometry) used here can be overcome. The alchemical method cannot distinguish between structures labeled by “N” and those labeled by “B” (see for example Figure 4: 4-1 and 4-2 are equivalent in the alchemical method). From geometry optimization of the structures of the simultaneous approach the following general feature emerges: the most stable structure is the N type structure (see eq 27 and the text below for definitions), e.g. 6-1 and 6-2 are N type structures, when 6-3 and 6-4 are B type structures. This feature was at the basis of the N-site attachment rule, but in the simultaneous substitution method where the most stable isomer is not always a descendant of the preceding most stable isomer, such an interpretation is misleading. Taking 12-1 in Figure 8 as an illustrative example for the N type structures, we can notice that the BN units which connect two phenanthrene-like structures are with boron atoms displaced inward and nitrogen atoms moved outward similar to 13-1 and 15-1 where the nitrogen bulges are slightly outward. Based on the data from Table II, two deficiencies of the alchemical method in the prediction of the most stable isomer can be categorized as follows: those up to 24 substitutions and those above this number. The first one is related to the appearance of a structure with two bonded carbons which are separated from the remaining carbons. The most stable isomers after geometry optimization for n = 12-16 are shown in Figure 8, labeled as n-1. In addition, 12-2 is the most stable n = 12 isomer obtained from the alchemical prediction. The analysis of the Taylor expansion for 12-1 and 12-2 structures shows the same interaction between the atoms located on the same face (the same contribution from the six largest hardness matrix elements listed in Table II). Two types of pentagons with one carbon are responsible for the difference in stability of these structures: the first one fused with three filled hexagons (marked by a green circle at Figure 8) and the second one fused with two filled hexagons (marked by an empty green circle at Figure 8). These patterns will be called first and second type pentagon, respectively. It may be noted that 12-2 has five completely filled hexagons, while 12-1 has only four; but 12-2 has three first type pentagons, whereas 12-1 has only two second type pentagons. In the successive substitution, the incoming BN units fill the next adjacent hexagons, to form 15-BN-fullerene with seven hexagons and five first type pentagons; but 15-1 is the most stable 15-BNfullerene, as it has six hexagons but only two pentagons of the first and two pentagons of the second type. Similar situations are encountered in the case of 17-BN-fullerene and 21-BN-fullerene; the most stable structure after geometry optimization is the thirdposition structure from the alchemical method. In both cases, the unfavorable two BN units chains appeared in the most stable forms predicted by the alchemical method. After geometry optimization, the structure with the lowest number of first type pentagons was found as the most stable BN fullerene. Note that the energy difference between the most stable structure after the geometry optimization and the one predicted by the alchemical method is less than 3 kcal/mol, even if the number of substitutions is beyond 24 pairs of carbon (see ΔEsimult n in Table II). Concluding this part, the substitution pattern emerging from our investigation reaffirms the need to rewrite Rule II as it was proposed in Section IV.B: the most stable BN fullerene for n ≤ 24 is characterized by the maximum number of filled hexagons (not necessarily consecutive) with the minimum number of first type

V. CONCLUSION In this paper, the usefulness of the alchemical derivatives in the prediction of the substitution pattern of carbon pairs of C60 fullerene by iso-electronic BN moieties was tested. Two methodologies of substitution were used: all substitutions made at the same time or substitutions made step by step: the simultaneous method and the successive method, respectively. The successive method was realized in two ways: one with 12 carbon atoms excluded from the substitution finally leading to the C12(BN)24 ball and the second without imposing any constraints resulting in the C20 (BN)20 belt. In all cases, the HF/ STO-3G alchemical predictions are qualitatively trustworthy and can be used as directive results for higher level studies. A comparison of the alchemical derivative based results with those obtained after geometry optimization at the B3LYP/3-21G level was more than satisfactory. The absolute magnitudes of the energies or other properties obviously depend strongly on the calculation level, but their trends and qualitative relations are similar. Thanks to the great computational efficiency of the alchemical method which follows from the fact that the evaluation of the energy of all transmutation products demands only some simple, extremely fast arithmetic operations, we were able to tackle the very large chemical space of the BN fullerenes and examine their stability. Starting from a very simple alchemical rule: maximize the contribution from the heteroatom part and minimize the contribution from the homoatom part, the substitution rules reported by Kar and co-workers were revisited and reformulated. A substitution pattern in the simultaneous method emerging from this work includes the following features: (I) the simultaneous BN substitutions prefer to replace the carbon−carbon pairs located at the h-h junction. The substitutions of the h-h junction occur in a way that eliminates the possibility of unstable homonuclear bonds and maximizes the number of BN bonds (this is the analog of rule I); (II) the most stable BN fullerene is characterized by the maximum number of filled hexagons (not necessarily consecutive) with the minimum number of first type pentagons (this rule includes rule II); (III) in the case when Dal,NN ≠ Dal,BB n n , the most stable structure is the N type structure (this rule reformulates Rule III: the N-site attachment rule); (IV) the most stable structure of the higher BN fullerenes is always some combination of the phenanthrene-, triphenylene-, or chrysene-like structures which are connected with each other directly or by one BN unit (this rule replaced Rule IV formulated by Kar and co-workers). It was shown that the procedure based on two rules: two pentagons are connected by a B−N bond (Rule I) and the completely substituted pentagons have only one homoatoms bond (Rule V) can be effectively used in finding of the most stable structures for more than 24 substitutions. As a result, the presented method can be effectively applied to study BN substitution and stability patterns in higher fullerenes, graphene, BN planar structures, and similar systems, exploring sometimes already large associated Chemical Spaces using a trustworthy yet a cost-effective procedure.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. M

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation ORCID

(18) Baben, M. T.; Achenbach, J. O.; von Lilienfeld, O. A. Guiding ab initio calculations by alchemical derivatives. J. Chem. Phys. 2016, 144, 104103. (19) von Lilienfeld, O. A.; Lins, R. D.; Rothlisberger, U. Variational particle number approach for rational compound design. Phys. Rev. Lett. 2005, 95, 153002. (20) Lesiuk, M.; Balawender, R.; Zachara, J. Higher order alchemical derivatives from coupled perturbed self-consistent field theory. J. Chem. Phys. 2012, 136, 034104. (21) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: 1989. (22) Geerlings, P.; De Proft, F.; Langenaeker, W. Conceptual density functional theory. Chem. Rev. 2003, 103, 1793−1873. (23) Liu, S.; Li, T.; Ayers, P. W. Potentialphilicity and potentialphobicity: Reactivity indicators for external potential changes from density functional reactivity theory. J. Chem. Phys. 2009, 131, 114106. (24) Geerlings, P.; Boisdenghien, Z.; Proft, F. D.; Fias, S. The E = E[N, v] functional and the linear response function: a conceptual DFT viewpoint. Theor. Chem. Acc. 2016, 135, 135. (25) Geerlings, P.; Fias, S.; Boisdenghien, Z.; De Proft, F. Conceptual DFT: chemistry from the linear response function. Chem. Soc. Rev. 2014, 43, 4989−5008. (26) Geerlings, P.; De Proft, F.; Fias, S. The Mathemical Properties of Linear Response Functions and their Implication for Chemistry. Acta Phys.-Chim. Sin. 2018, 34, xx. (27) Sablon, N.; De Proft, F.; Ayers, P. W.; Geerlings, P. Computing Fukui functions without differentiating with respect to electron number. II. Calculation of condensed molecular Fukui functions. J. Chem. Phys. 2007, 126, 224108. (28) Yang, W.; Cohen, A. J.; De Proft, F.; Geerlings, P. Analytical evaluation of Fukui functions and real-space linear response function. J. Chem. Phys. 2012, 136, 144110. (29) Sablon, N.; De Proft, F.; Geerlings, P. The Linear Response Kernel: Inductive and Resonance Effects Quantified. J. Phys. Chem. Lett. 2010, 1, 1228−1234. (30) Fias, S.; Geerlings, P.; Ayers, P.; De Proft, F. sigma, pi aromaticity and anti-aromaticity as retrieved by the linear response kernel. Phys. Chem. Chem. Phys. 2013, 15, 2882−2889. (31) Boisdenghien, Z.; Fias, S.; Da Pieve, F.; De Proft, F.; Geerlings, P. The polarisability of atoms and molecules: a comparison between a conceptual density functional theory approach and time-dependent density functional theory. Mol. Phys. 2015, 113, 1890−1898. (32) Boisdenghien, Z.; Van Alsenoy, C.; De Proft, F.; Geerlings, P. Evaluating and Interpreting the Chemical Relevance of the Linear Response Kernel for Atoms. J. Chem. Theory Comput. 2013, 9, 1007− 1015. (33) Stuyver, T.; Fias, S.; De Proft, F.; Fowler, P. W.; Geerlings, P. Conduction of molecular electronic devices: qualitative insights through atom-atom polarizabilities. J. Chem. Phys. 2015, 142, 094103. (34) Cohen, M. H.; Ganduglia-Pirovano, M. V.; Kudrnovsky, J. Electronic and Nuclear-Chemical Reactivity. J. Chem. Phys. 1994, 101, 8988−8997. (35) Cohen, M. H.; Ganduglia-Pirovano, M. V.; Kudrnovsky, J. Reactivity Kernels, the Normal-Modes of Chemical-Reactivity, and the Hardness and Softness Spectra. J. Chem. Phys. 1995, 103, 3543−3551. (36) Baekelandt, B. G. The nuclear Fukui function and Berlin’s binding function in density functional theory. J. Chem. Phys. 1996, 105, 4664− 4667. (37) Balawender, R.; Welearegay, M. A.; Lesiuk, M.; De Proft, F.; Geerlings, P. Exploring Chemical Space with the Alchemical Derivatives. J. Chem. Theory Comput. 2013, 9, 5327−5340. (38) Balawender, R.; Holas, A.; De Proft, F.; Van Alsenoy, C.; Geerlings, P. In Many-Body Approaches at different Scales: a Tribute to N. H. March on the Occasion of his 90th Birthday; Angilella, G. G. N., Amovilli, C., Eds.; Springer: 2017; p xxx. (39) Yamaguchi, Y.; Osamura, Y.; Goddard, J. D.; Schaefer, H. F. A New Dimension to Quantum Chemistry; Oxford University Press: New York, 1994.

Robert Balawender: 0000-0002-4216-5255 Funding

The authors acknowledge financial support by the VUB (Vrije Universiteit Brussel) under the form of a Strategic Research Program (SRP) (P.G. and F.D.P.) and the Interdisciplinary Centre for Mathematical and Computational Modelling computational grant (R.B.). F.D.P. also acknowledges the Francqui foundation for a position as a Francqui research professor. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The authors are grateful to A. Holas for helpful discussions. REFERENCES

(1) Kirkpatrick, P.; Ellis, C. Chemical space. Nature 2004, 432, 823− 823. (2) Dobson, C. M. Chemical space and biology. Nature 2004, 432, 824−828. (3) von Lilienfeld, O. A. First principles view on chemical compound space: Gaining rigorous atomistic control of molecular properties. Int. J. Quantum Chem. 2013, 113, 1676−1689. (4) Franceschetti, A.; Zunger, A. The inverse hand-structure problem of finding an atomic configuration with given electronic properties. Nature 1999, 402, 60−63. (5) Johannesson, G. H.; Bligaard, T.; Ruban, A. V.; Skriver, H. L.; Jacobsen, K. W.; Norskov, J. K. Combined electronic structure and evolutionary search approach to materials design. Phys. Rev. Lett. 2002, 88, 255506. (6) Wang, M.; Hu, X.; Beratan, D. N.; Yang, W. Designing molecules by optimizing potentials. J. Am. Chem. Soc. 2006, 128, 3228−3232. (7) Kuhn, C.; Beratan, D. N. Inverse strategies for molecular design. J. Phys. Chem. 1996, 100, 10595−10599. (8) Balamurugan, D.; Yang, W.; Beratan, D. N. Exploring chemical space with discrete, gradient, and hybrid optimization methods. J. Chem. Phys. 2008, 129, 174105. (9) De Vleeschouwer, F.; Chankisjijev, A.; Geerlings, P.; De Proft, F. Designing Stable Radicals with Highly Electrophilic or Nucleophilic Character: Thiadiazinyl as a Case Study. Eur. J. Org. Chem. 2015, 506− 513. (10) De Vleeschouwer, F.; Geerlings, P.; De Proft, F. Molecular Property Optimizations with Boundary Conditions through the Best First Search Scheme. ChemPhysChem 2016, 17, 1414−1424. (11) Teunissen, J. L.; De Proft, F.; De Vleeschouwer, F. Tuning the HOMO-LUMO Energy Gap of Small Diamondoids Using Inverse Molecular Design. J. Chem. Theory Comput. 2017, 13, 1351−1365. (12) von Lilienfeld, O. A.; Tuckerman, M. E. Molecular grandcanonical ensemble density functional theory and exploration of chemical space. J. Chem. Phys. 2006, 125, 154104. (13) von Lilienfeld, O. A.; Tuckerman, M. E. Alchemical Variations of Intermolecular Energies According to Molecular Grand-Canonical Ensemble Density Functional Theory. J. Chem. Theory Comput. 2007, 3, 1083−1090. (14) von Lilienfeld, A. O. Accurate ab initio energy gradients in chemical compound space. J. Chem. Phys. 2009, 131, 164102. (15) Sheppard, D.; Henkelman, G.; von Lilienfeld, O. A. Alchemical derivatives of reaction energetics. J. Chem. Phys. 2010, 133, 084104. (16) Chang, K. Y.; von Lilienfeld, O. A. Quantum mechanical treatment of variable molecular composition: from ’alchemical’ changes of state functions to rational compound design. Chimia 2014, 68, 602− 608. (17) Chang, K. Y.; Fias, S.; Ramakrishnan, R.; von Lilienfeld, O. A. Fast and accurate predictions of covalent bonds in chemical space. J. Chem. Phys. 2016, 144, 174110. N

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation (40) Baker, J.; Wolinski, K.; Malagoli, M.; Pulay, P. Parallel implementation of Hartree-Fock and density functional theory analytical second derivatives. Mol. Phys. 2004, 102, 2475−2484. (41) Bast, R.; Ekstrom, U.; Gao, B.; Helgaker, T.; Ruud, K.; Thorvaldsen, A. J. The ab initio calculation of molecular electric, magnetic and geometric properties. Phys. Chem. Chem. Phys. 2011, 13, 2627−2651. (42) Balawender, R.; Komorowski, L. Atomic Fukui function indices and local softness ab initio. J. Chem. Phys. 1998, 109, 5203−5211. (43) Balawender, R.; Komorowski, L.; De Proft, F.; Geerlings, P. Derivatives of molecular valence as a measure of aromaticity. J. Phys. Chem. A 1998, 102, 9912−9917. (44) Balawender, R.; Safi, B.; Geerlings, P. Solvent effect on the global and atomic DFT-based reactivity descriptors using the effective fragment potential model. Solvation of ammonia. J. Phys. Chem. A 2001, 105, 6703−6710. (45) Balawender, R.; De Proft, F.; Geerlings, P. Nuclear Fukui function and Berlin’s binding function: Prediction of the Jahn-Teller distortion. J. Chem. Phys. 2001, 114, 4441−4449. (46) Balawender, R.; Geerlings, P. Nuclear Fukui function from coupled perturbed Hartree-Fock equations. J. Chem. Phys. 2001, 114, 682−691. (47) Geerlings, P.; Proft, F. D.; Balawender, R. Reviews of Modern Quantum Chemistry 2002, 1053−1070. (48) Lesiuk, M.; Zachara, J. Molecular electrostatic potential at the atomic sites in the effective core potential approximation. J. Chem. Phys. 2013, 138, 074107. (49) Nag, A.; Raidongia, K.; Hembram, K. P.; Datta, R.; Waghmare, U. V.; Rao, C. N. Graphene analogues of BN: novel synthesis and properties. ACS Nano 2010, 4, 1539−1544. (50) Panchakarla, L. S.; Govindaraj, A.; Rao, C. N. R. Boron- and nitrogen-doped carbon nanotubes and graphene. Inorg. Chim. Acta 2010, 363, 4163−4174. (51) Raidongia, K.; Gomathi, A.; Rao, C. N. R. Synthesis and Characterization of Nanoparticles, Nanotubes, Nanopans, and Graphene-like Structures of Boron Nitride. Isr. J. Chem. 2010, 50, 399−404. (52) Raidongia, K.; Nag, A.; Hembram, K. P.; Waghmare, U. V.; Datta, R.; Rao, C. N. BCN: a graphene analogue with remarkable adsorptive properties. Chem. - Eur. J. 2010, 16, 149−157. (53) Kumar, N.; Moses, K.; Pramoda, K.; Shirodkar, S. N.; Mishra, A. K.; Waghmare, U. V.; Sundaresan, A.; Rao, C. N. R. Borocarbonitrides, BxCyNz. J. Mater. Chem. A 2013, 1, 5806−5821. (54) Govindaraj, A.; Rao, C. N. R. Doping of Graphene by Nitrogen, Boron, and Other Elements. In Functionalization of Graphene; 2014; DOI: 10.1002/9783527672790.ch10. (55) Moses, K.; Shirodkar, S. N.; Waghmare, U. V.; Rao, C. N. R. Composition-dependent photoluminescence and electronic structure of 2-dimensional borocarbonitrides, BCXN (x= 1, 5). Mater. Res. Express 2014, 1, 025603. (56) Pattanayak, J.; Kar, T.; Scheiner, S. Boron-nitrogen (BN) substitution patterns in C/BN hybrid fulterenes: C60−2x(BN)(x) (x = 1−7). J. Phys. Chem. A 2001, 105, 8376−8384. (57) Pattanayak, J.; Kar, T.; Scheiner, S. Boron-nitrogen (BN) substitution of fullerenes: C-60 to C12B24N24CBN ball. J. Phys. Chem. A 2002, 106, 2970−2978. (58) Kar, T.; Pattanayak, J.; Scheiner, S. Rules for BN-substitution in BCN-fullerenes. Separation of BN and C domains. J. Phys. Chem. A 2003, 107, 8630−8637. (59) Fan, X. F.; Zhu, Z. X.; Shen, Z. X.; Kuo, J. L. On the use of bondcounting rules in predicting the stability of C12B6N6 fullerene. J. Phys. Chem. C 2008, 112, 15691−15696. (60) Golberg, D.; Bando, Y.; Tang, C. C.; Zhi, C. Y. Boron Nitride Nanotubes. Adv. Mater. 2007, 19, 2413−2432. (61) Terrones, M.; Romo-Herrera, J. M.; Cruz-Silva, E.; Lopez-Urias, F.; Munoz-Sandoval, E.; Velazquez-Salazar, J. J.; Terrones, H.; Bando, Y.; Golberg, D. Pure and doped boron nitride nanotubes. Mater. Today 2007, 10, 30−38.

(62) Golberg, D.; Bando, Y.; Huang, Y.; Terao, T.; Mitome, M.; Tang, C.; Zhi, C. Boron nitride nanotubes and nanosheets. ACS Nano 2010, 4, 2979−2993. (63) Zhi, C.; Bando, Y.; Tang, C.; Golberg, D. Boron nitride nanotubes. Mater. Sci. Eng., R 2010, 70, 92−111. (64) Rubio, A.; Corkill, J. L.; Cohen, M. L. Theory of graphitic boron nitride nanotubes. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 5081−5084. (65) Chopra, N. G.; Luyken, R. J.; Cherrey, K.; Crespi, V. H.; Cohen, M. L.; Louie, S. G.; Zettl, A. Boron nitride nanotubes. Science 1995, 269, 966−967. (66) Moss, G. P. Nomenclature of fused and bridged fused ring systems - (IUPAC Recommendations 1998). Pure Appl. Chem. 1998, 70, 143− 216. (67) Powell, W. H. Revision of the Extended Hantzsch-Widman System of Nomenclature for Heteromonocycles. Pure Appl. Chem. 1983, 55, 409−416. (68) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864−B871. (69) von Lilienfeld, O. A.; Tuckerman, M. E. Molecular grandcanonical ensemble density functional theory and exploration of chemical space. J. Chem. Phys. 2006, 125, 154104. (70) Nalewajski, R. F.; Parr, R. G. Legendre Transforms and Maxwell Relations in Density Functional Theory. J. Chem. Phys. 1982, 77, 399− 407. (71) Ray, N. K.; Parr, R. G. Diamagnetic Shieldings of Atoms in Molecules and Their Relation to Electronegativity. J. Chem. Phys. 1980, 73, 1334−1339. (72) Hedberg, K.; Hedberg, L.; Bethune, D. S.; Brown, C. A.; Dorn, H. C.; Johnson, R. D.; De Vries, M. Bond lengths in free molecules of buckminsterfullerene, c60, from gas-phase electron diffraction. Science 1991, 254, 410−412. (73) Senn, P. Computation of the Cartesian Coordinates of Buckminsterfullerene. J. Chem. Educ. 1995, 72, 302−303. (74) 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. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. General Atomic and Molecular Electronic-Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (75) Cozzi, F.; Powell, W. H.; Thilgen, C. Numbering of Fullerenes (IUPAC Recommendations 2004). Pure Appl. Chem. 2005, 77, 843.

O

DOI: 10.1021/acs.jctc.7b01114 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX