A computational study of photonic materials based on Ni bis(dithiolene) fused with benzene, possessing gigantic second hyperpolarizabilities

Abstract

We propose a family of photonic materials based on nickel bis(dithiolene) [NiBDT] and benzene (B). Linear and planar oligomers of 1D and 2D polymers, together with several rich π-electron linkers (e.g. C₁₄H₈, C₂₄H₈, and C₄O₂H₄), have been employed in order to demonstrate the superior performance of the proposed materials. A key property for the design of efficient photonic materials is the second hyperpolarizability. We have, thus, computed it together with several other properties, such as the dipole moment, polarizability, first excitation energy, etc. We have also used two indices to estimate the degree of the diradical character of the considered derivatives. The calculations have been performed by employing density functional theory. Considerable care has been taken to validate the computational procedure we have used. The main findings of the present work are as follows: (a) the 1D and 2D coordination polymers based on nickel bis(dithiolene) and benzene and in particular the linear ones, (B-NiBDT)_n, are very likely to be efficient photonic materials, since they have gigantic second hyperpolarizability. (b) The second hyperpolarizability greatly depends on the organic linker with which NiBDT is fused and on the structure; for example, the linear oligomers have much higher second hyperpolarizability than the planar ones. The present work together with a recent article (A. Avramopoulos et al., J. Phys. Chem. C, 2016, 120, 9419–9435) demonstrates a mechanism for designing efficient photonic materials.

---

I. Introduction

There is an ever-increasing need for novel photonic materials for use in a large number of applications (e.g. laser protection,¹ devices for all-optical control,² optical communication, optical signal processing³). Here we propose a novel family of photonic materials: 1D and 2D coordination polymers (nanosheets), based on nickel bis(dithiolene) and benzene. We recently reported that the oligomers (RA-NiBDT)_n, RA = radiaannulene, NiBDT = bis(ethylene-1,2-dithiolato)Ni, have giant second hyperpolarizabilities. We found in the literature that Kambe et al.^4–7 synthesized 1D and 2D (nanosheets) coordination polymers, based on NiBDT. It was a natural extension of our previous work to compute the (hyper)polarizabilities of some of their models, employing density functional theory. Specifically, we got some fragments of the 1D and 2D polymers and using these linear and planar oligomers, we performed a systematic study of their linear and nonlinear optical properties. We found that some of these oligomers have extremely large nonlinearities, which is one of the prerequisites for an efficient photonic material.

In order to demonstrate the superior non-linear optical (NLO) performance of the proposed materials, over related compounds, we have carried out extensive experimentation, which involved NiBDT, fused with a series of rich π-electron organic linkers:

(a) Benzene. This was fused with NiBDT to produce linear and planar oligomers (nanosheets). Extensive experimental work has been reported on the corresponding polymers.^4–7 The linear oligomer, symbolically, is denoted as (B-NiBDT)_n, where n = 1–6. The planar oligomers have been studied by employing 3 models (Fig. 1). This choice also allows studying the significant effect of the geometry on the properties of interest.

Fig. 1 The structures (1–3) of several planar nickel bis(dithiolene) complexes, computed with B3LYP/6-31G*.

(b) Linkers with extensive conjugation. We have used C₁₄H₈ (pyracyclene), C₂₄H₁₂ (dicyclopenta(cd,lm)perylene), and C₂₄H₈ (this is a radiaannulene, RA). These linkers allow studying the effect of the extensive π-electron network on the properties of interest.

(c) C ₄ O ₂ H ₄ , C ₄ S ₂ H ₄ , C ₄ OH ₄ , and C ₄ SH ₄ . We have used a series of small molecules to connect the NiBDT units, the effect of which has been compared with that of benzene. These are C₄O₂H₄ (1,4-dioxin), C₄S₂H₄ (1,4-dithiin), C₄OH₄ (furan), and C₄SH₄ (thiophene). Furan and thiophene are aromatic, 1,4-dioxin is non-aromatic and 1,4-dithiin is antiaromatic. It is known that thiophene is an important element of many photonic materials.⁸

Of major importance in this work are the (hyper)polarizabilities and in particular the second hyperpolarizability. For their rationalization we have used several other properties, such as E_HOMO (energy of the highest occupied molecular orbital), E_LUMO (energy of the lowest unoccupied molecular orbital), and the first excitation energy. A significant property of our studied compounds is their diradical character (DC). This is one of the major characteristics for understanding the intriguing properties of the NiBDT derivatives. Thus, we have used two indices which estimate, or they are related to, the strength of the DC. This work, in connection with our previous results, demonstrates a mechanism for designing photonic materials with exceptionally large nonlinearities. This mechanism includes a series of units consisting of NiBDT, π-electron rich linkers, Li substituents⁹ and the appropriate structure, for example a linear arrangement leads to a much larger second hyperpolarizability in comparison to the corresponding planar one.

This article is organized as follows: in Section II we briefly present the most relevant articles of the extensive literature, which is devoted to the various topics of this work; in Section III, we describe the methods we have used, in Section IV we analyze our results and finally in Section V, the concluding remarks are presented.

1. Literature survey

We briefly review some recent, relevant work, related to 2D nanosheets, polymers, NiBDT and thiophene.

Intensive research is currently being carried out on 2D nanosheets (e.g. graphene,¹⁰ BN,¹¹chalcogenides,¹² silicene,¹³ and germanane¹⁴), due to the possibility of tuning their properties (e.g. photoluminescence¹⁵ and photoconductivity¹⁶) in order to optimize their performance in many applications (e.g. photonics, biomedicine, electronics, and energy storage).^2,10,17,18 In particular various, π-conjugated (X = Ni, Pd, Pt, Co) bis(dithiolene) nanosheets have been studied.^19–21 High conductivity was observed for Ni bis(dithiolene) π-nanosheets.⁴ Tang and Zhou²² studied the binding of ethylene to nickel bis(dithiolene). A 2D nickel bis(dithiolene) sheet is an organic topological insulator.^4,16,23,24 Zhou found that this sheet behaves like a metal, when it interacts with graphene to form a heterobilayer material.²⁵ Shojaei and Kang found that between two layers of a π-conjugated Ni bis(dithiolene) complex (NiS₄C₄), covalent bonds are formed. The bilayer is potentially useful for the design of electromechanical and optoelectronic devices.²⁶ Shojaei et al.²³ investigated the electronic and mechanical properties of a 2D π-conjugated Ni bis(dithiolene) sheet (NiS₄C₄) and its Pd analogue (PdS₄C₄). Liu et al.²⁷ found, by using first principles calculations, that the π-conjugated Ni bis(dithiolene) complex with formula Ni₃S₁₂C₁₂ can be used to detect poisonous gases (e.g. CO, NO). A change in the properties (from semiconducting to metallic) of the Ni₃S₁₂C₁₂ nanosheet was induced by adsorbed molecules. Pal et al.²⁰ reported the synthesis of palladium bis(dithiolene) nanosheets. Clough et al.²¹ used 2D polymer nanosheets containing Co-bis(dithiolene) units for efficient hydrogen evolution from water. It has also been reported that a cobalt dithiolene metal–organic framework exhibits band-like metallic conductivity.²⁸ Other related materials have also been studied, for example, Ni₃(HITP)₂,²⁹ Cu₃(HITP)₂ [22]³⁰ and HTT-Pt,³¹ where HITP = 2,3,6,7,10,11-hexaiminotriphenylene and HTT = triphenylene hexathiol. Excellent NLO properties have been observed for low-dimensional materials, together with the possibility of developing new photonic devices based on them.^{1–3,32–35}

Polymers possessing delocalized π-electrons have been investigated extensively, as a potential source of photonic materials.^36–42

NiBDT. Many articles on Ni bis(dithiolene) related derivatives and several of their properties have been reported (e.g. HOMO and LUMO energy levels, third order NLO properties).^43–52 Deplano et al.⁵³ reviewed square-planar d⁸ metal mixed-ligand dithiolene complexes as second order nonlinear optical chromophores. One of their goals was to propose tools for optimizing the optical properties. Several markers for the prediction of the optical properties have been proposed (e.g. λ_max, ν(C=C)).⁵³ Garreau-de Bonneval et al.⁵⁴ reviewed the synthesis, electronic properties and applications of neutral d⁸ metal bis-dithiolene complexes. They noted that the NiBDT derivatives, due to their optical, electronic, magnetic and conductive properties, provide an excellent basis for the development of many devices (e.g. opto-electronic).⁵⁴

Thiophenes. Iyoda and Shimizu reviewed π-expanded oligothiophene macrocycles; they noted that their properties (e.g. optoelectronic) depend on the inserted π-systems. The authors discussed the remarkable optical and electronic properties of these systems.⁵⁵ Li et al. employed a hybrid computational scheme, including MD, QM/MM and the Frenkel exciton theory of vibronic coupling, to study the optical spectra of α-sexithiophene nanoparticles.⁵⁶ Gonzalez et al. discussed the self-assembly of fully conjugated π-expanded macrocyclic oligothiophenes, complexed with fullerenes. The role of non-covalent interactions in controlling their electronic properties was analysed.⁵⁷ Oligothiophene linkers were used to create a series of structural scaffolds. A series of molecules were studied and this allowed for fine-tuning of emissions in the range of 412–540 nm.⁵⁸

II. Methods

This section discusses the following: (i) a description and justification of the functionals we employed; (ii) definition of (hyper)polarizabilities; (iii) definition of the diradical character and the indices we used to compute it; (iv) description of the procedures we have employed to validate the methods we have used; (v) description of the methods we have employed to compute the molecular structures, (vi) calculation of the excitation spectrum; (vii) decomposition of the molecular polarizabilities of some selected systems into atomic contributions by using the fractional occupation Hirshfeld-iterative method and (viii) vibrational effects on the (hyper)polarizabilities.

(a) Functionals. The molecular properties were computed by employing a series of functionals which include the following:

(a) CAM-B3LYP.⁵⁹ This is the long range (LR) corrected version of the B3LYP functional; it has given satisfactory (hyper)polarizabilities for large and extended molecules.^9,60,61

(b) LC-BLYP.⁶² This long-range corrected functional (BLYP) has been shown to provide reasonable second hyperpolarizability values for singlet diradical molecules^63,64 and a series of fullerenes.⁶⁵

(c) ωB97X.⁶⁶ This LC functional gave satisfactory results for the second-order nonlinear optical responses of a series of ferrocene-tetrathiafulvalene hybrids⁶⁷ and several properties (e.g. affinities, dissociation energies, and dipole moments) of ytterbium doped silicon clusters.⁶⁸

(d) LC-PBE.⁶⁹ This long-range corrected version of PBE was developed by Perdew et al. Lu used this functional to satisfactorily compute the dynamic quadratic polarizability of several organic molecules.⁷⁰ Jacquemin et al. computed the longitudinal dipole moments and static electronic first hyperpolarizabilities of increasingly long polymethineimine oligomers.⁷¹ These authors noted the efficiency of long-range DFT, even in very pathological cases.

(e) BHandHLYP.^72,73 This has been shown to provide satisfactory second order hyperpolarizability values for diradical molecules.^74,75

(f) B2-PLYP.⁷⁶ This functional provided reasonable second hyperpolarizabilities for NiBDT and several Ni-based derivatives, compared with those calculated with the more accurate CCSD(T) and CASPT2 methods.⁷⁴

(b) Definitions. When a molecule is set in a uniform static electric field, F, its energy, E(F), is given by the following expansion:

E(F) = E⁰ − μ_iF_i − (1/2)α_ijF_iF_j − (1/6)β_ijkF_iF_jF_k − (1/24)γ_ijklF_iF_jF_kF_l − …, (1)

where E⁰ is the field free energy; μ_i, α_ij, β_ijk and γ_ijkl are the dipole moment, polarizability, first hyperpolarizability and second hyperpolarizability components, respectively. A summation over repeated indices is implied. Analytic⁷⁷ and finite field approaches,⁷⁸ using eqn (1), have been used to compute the (hyper)polarizabilities. The property values have been calculated by using the Romberg technique, to safeguard the numerical stability for the (hyper)polarizability results.^79–81 Several field strengths, with magnitude 2^mF, m = 1–4, and F = 0.001 a.u., have been employed. The GAUSSIAN 09 software⁷⁷ has been used for the reported DFT calculations.

(c) Diradical character. Definitions of the diradical character may be found in the literature.^74,82,83

It has been noted that in a singlet diradical, there are two singly occupied molecular orbitals of equal energy; these are weakly antiferromagnetically coupled.⁸⁴ Singlet diradicals should be studied with a multi-configurational wavefunction, which can properly treat static correlation. However, such a method is computationally very demanding, in particular for molecules with the size of those considered in this work. It is known that the symmetry broken (BS) DFT can simulate the static correlation;⁸⁴ this method has adequately described singlet diradicals.^74,84–86

The diradical character may be estimated by using the following expressions:⁷⁴

(a) E₁ = E_RDFT(singlet) − E_UDFT(singlet)

E₁ gives the energy lowering obtained by using the DFT broken symmetry solution, in comparison with the RDFT value. A large DC is associated with a large E₁ value.

(b) E₃ = (2(E_T − E_(U)DFT)/(2 − 〈S²〉_BS)),

where E_T is the UDFT energy of the triplet, E_(U)DFT is the BS-DFT energy of the singlet and 〈S²〉 is the total spin expectation value.⁸⁶

(d) Validation. Most of our results have been calculated by employing the (U)CAM-B3LYP and (U)BHandHLYP functionals with the 6-31G* basis set for H, C, S, and O atoms.^87,88 For Ni the quasi-relativistic, effective core potential, ECP28MWB(SDD), was used.⁸⁹ Detailed studies on the performance of CAM-B3LYP and BHandHLYP, in connection with the 6-31G* basis set, for Ni derivatives, have been performed elsewhere.⁹ A comparison of the DFT results with the more rigorous CASSCF/CASPT2 [complete active space multiconfiguration SCF/CAS second order perturbation theory] and CCSD(T) [coupled-cluster with single, double and perturbative triple excitations (CCSD(T))] ones confirms the satisfactory performance of the employed DFT methods.⁷⁴ It should also be noted that the adequacy of the 6-31G* basis set, for the selected compounds, has been confirmed by employing a series of larger basis sets, including more diffuse and polarization functions (Table S1, ESI†).

(e) Molecular structures. All the molecular structures have been computed by employing the B3LYP/6-31G* method. It has been demonstrated that this approach gives satisfactory structures for NiBDT derivatives having conjugated linkers.⁹ Furthermore, B3LYP was shown to give satisfactory geometries for diradical systems.⁹⁰ It has been established that a broken symmetry method provides similar molecular geometries, compared with those computed with RDFT.⁸⁴ Thus, RDFT was used for the geometry optimization. Vibrational analysis was performed for all structures employed in this study, to verify that a stationary point has been found.

(f) Excitation spectrum. Time dependent density functional theory (TD-DFT) was used to compute the transition energies and the transition dipole moments,^91,92 with the CAM-B3LYP functional.⁵⁹ This gave satisfactory estimates (in comparison with the corresponding experimental values) of transition energies of some organic dyes.^93–96

(g) The fractional occupation Hirshfeld-iterative method. As in one of our recent studies⁹ on Ni-(bis)dithiolene derivatives, we decomposed the molecular polarizabilities of some selected systems of the current study into atomic contributions by means of the fractional occupation Hirshfeld-iterative method (FOHI).^97,98 Through this well tested partitioning scheme, we obtained pure atomic polarizability values (hereafter intrinsic polarizabilities), free of size and origin dependent charge transfer interactions.⁹⁹ In addition, we have computed and plotted contour maps of the mean intrinsic polarizability (MIP) densities in order to identify the most polarizable regions of some of the studied systems.⁹⁹ The aim of these treatments is to obtain useful size independent chemical insights into the polarizability evolution of the systems considered by identifying the most polarizable atoms in the considered molecules and study effects which are expected to arise from the various structural changes that have been applied during our designing strategy. More information and some demonstrative examples on the application of the FOHI method in various study cases can be found in ref. 100–102 and references therein. All electron density matrices as well as their derivatives with respect to an electric field were computed analytically as implemented in GAUSSIAN 09,⁷⁷ while atomic partitioning was carried out with the BRABO package,¹⁰³ the STOCK program¹⁰⁴ and homemade FORTRAN codes.

(h) Vibrational effects. The reliable calculation of the vibrational properties of molecules of the size investigated here is generally a rather difficult task. We computed the vibrational contributions of two relatively small models, related to the derivatives of interest. The perturbation theory, as pioneered by Bishop and Kirtman,^105–107 was used for these calculations. Only the static contribution to the so-called nuclear relaxation (nr) was computed, which takes into account that the properties of the optimum geometry in an electric field differ from those of the field-free optimized geometry. The nr contribution was computed by means of the Bishop–Hasan–Kirtman approach, which is based on geometry optimizations in an applied electric field.¹⁰⁸ The nr contribution is only the lowest order part to the total vibrational contribution. Further contributions, such as the zero-point average or the curvature contribution, are generally too costly to compute even for the smallest molecules considered here. This implies that no conclusions on the convergence behavior of the perturbational series can be drawn. The properties were computed by numerical differentiations from both energies and analytically computed dipole moments for several field strengths and applying the Romberg–Rutishauser differentiation method.^80,81 The range of the applied field strengths was either 0.0001–0.0064, or, in some difficult cases, 0.000025–0.0064 a.u.

The dipole moments and (hyper)polarizabilities are expressed in atomic units; E₁ and E₂ are given in kcal mole⁻¹, while the excitation energy is given in eV. Conversion factors to other systems of units are also provided.¹⁰⁹

III. Results and discussion

The results presented in this section will demonstrate the extremely large second hyperpolarizabilities of 1D coordination polymers, based on Ni-bis(dithiolenes). This will be shown by employing a model involving (B-NiBDT)_n, n = 1–6 (Table 1). Two important related questions that we deal with in this section are as follows: (a) how does the structure affect the properties of interest? In order to comment on this question we have employed several models, some of which involve planar fragments of the nanosheet Ni₃S₁₂C₁₂ (Table 1, (B-NiBDT)_n; Fig. 1). (b) We have found that an important element for the optimization of the L&NLO properties is the fused linkers with the NiBDT units. Thus, we have used several rich π-electron linkers to find their effect on the L&NLO properties (Fig. 2). It is well known that extended π-electron networks are of cardinal importance for the development of NLO materials and the rationalization of their properties.⁹ A molecular unit which is characterized as an important building block for the fabrication of photonic materials is thiophene and its related derivatives (Fig. 3). Thus, we consider its effect, as a linker, on the L&NLO properties of the resulting Ni-bis(dithiolene) compounds. This work mainly presents electronic contributions to (hyper)polarizabilities. However, for completeness, we will briefly comment on the vibrational properties, specifically on the nuclear relaxation contribution to (hyper)polarizabilities of two relatively small model derivatives, related to the compounds of interest.

Table 1 The dipole moment, (hyper)polarizabilities, E_HOMO/LUMO, E_HOMO − E_LUMO gap (|H–L|), indices E₁ and E₂ and the excitation energy (E_exc) of (B-NiBDT), n = 1–6. The properties (P) were computed at the B3LYP/6-31G* optimized geometry. All property values are given in a.u., except those of E_exc (eV) and E₁/E₂ (kcal mol⁻¹)

n	P
	μ x	α xx	α	β xxx × 10^−3	γ xxxx × 10^−5	(γxxxx × 10^−5)/n^a	(γxxxx × 10^−5)/N^b	E 1 E 2	E HOMO	E LUMO	|H–L|	E exc
a n = the number of benzene–NiBDT units*. b N: number of electrons. c E 1 = ER − EU, method: B3LYP/6-31G*. d E 2 = (2(ET − EU)/(2 − 〈S^2〉BS)), method: B3LYP/6-31G*. e Method: (U)CAM-B3LYP/6-31G*. f Method: (U)BHandHLYP/6-31G. g Method: (U)B2PLYP/6-31G*. h (U)LC-ωPBE. i Method: (U)LC-PBEPBE/6-31G*. j Method: (U)ωB97X/6-31G*. k Method: (U)CCSD/6-31G*. l Method: (U)CCSD(T)/6-31G*. m Method: (U)LC-BLYP/6-31G*.
1	−0.438^f	371^f	229^f	3.8^f	19^f	19.0^f	0.14^f	0.0	−0.264^f	−0.111^f	0.153^f	1.92
	−0.231^e	422^e	239^e	6.0^e	20^e	20.0^e	0.14^e	11.67	−0.272^e	−0.108^e	0.164^e
	−0.042^g	488^g		10.9^g	41^g	41.0^g	0.30^g		−0.262^g	−0.098^g	0.164
	−0.302^h	396^h		5.4^h	23^h				−0.334^h	−0.074^h	0.260^h
	−0.380^i	369^i	387^i	2.70^i	14^i				−0.293^i	−0.089^i	0.204^i
	−0.206^j	435^j	242^j	8.2^j	26^j
	−0.163^k	421^k		2.1^k	13^k
	−0.197^l	486^l		6.1^l	38^l
2	−0.670^f	1336^f	576^f	5.4^f	327^f	163.5^f	1.28^f	0.37	−0.254^f	−0.148^f	0.106^f	1.26
	−0.392^e	1684^e	766^e	17.1^e	570^e	285.0^e	2.23^e	7.05	−0.263^e	−0.139^e	0.124^e
	−0.282^g	2020^g		22.5^g	861^g	430.5^g	3.36^g		−0.249^g	−0.137^g	0.112^g
	−0.569^i	1347^i	1275^i	−16.1^i	640^i				−0.315^i	−0.117^i	0.198^i
	−0.347^j	1779^j	782^j	25.2^j	899^j				−0.286^j	−0.161^j	0.165^j
3	−0.784^f	2732^f	1076^f	−2.3^f	1693^f	564.3^f	4.50^f	1.14	−0.256^f	−0.157^f	0.099^f	1.13
	−0.497^e	3836^e	1617^e	13.7^e	5380^e	1793.3^e	14.30^e	5.72	−0.256^e	−0.148^e	0.118^e
	−0.455^g	4579^g		1.6^g	5670^g	1890.0^g	15.10^g		−0.251^g	−0.146^g	0.105^g
	−0.584^h	3463^h							−0.318^h	−0.126^i	0.192^i
	−0.670^i	2809^i	1525^i	−178.0^i	3810^i				−0.288^j	−0.129^j	0.159^j
	−0.441^j	4141^j	1667^j	23.3^j	9000^j
4	−0.846^f	4383^f	1667^f	−16.6^f	4650^f	1162.5^f	9.38^f	2.01	−0.257^f	−0.161^f	0.096^f	1.04
	−0.571^e	6743^e	2764^e	−29.6^e	23 600^e	5900.0^e	47.60^e	2.55	−0.266^e	−0.151^e	0.115^e
	−0.509^j	7407^j	2857^j	−40.2^j	42 000^j				−0.289^j	−0.133^j	0.156^j
5	−0.881^f	6173^f	2307^f	−31.8^f	9160^f	1832.0^f	14.87^f	2.88	−0.259^f	−0.165^f	0.094^f	0.98
	−0.624^e	10 212^e	4135^e	−117.5^e	66 000^e	13200.0^e	107.10^e	1.41	−0.269^e	−0.156^e	0.113^e
	−0.710^h	8800^h		−110.0^h					−0.321^h	−0.134^h	0.187^h
	−0.775^i	6579^i	2892^i	−36.7^i	10 200^i				−0.290^i	−0.137^i	0.153^i
	−0.56^j	11 372^j	4280^j	−181.8^j	120 000^j				−0.312^j	−0.126^j	0.186^j
	−0.829^m	6253^m	2390^n	−30.7^m	7280^m
6	−0.901^f	8033^f	2971^f	−44.7^f	14 800^f	2466.7^f	20.11^f	3.76	−0.259^f	−0.168^f	0.091^f	0.94
	−0.660^e	14 047^e	5606^e	−233.8^e	139 000^e	23166.7^e	188.86^e	1.01	−0.270^e	−0.158^e	0.112^e
	−0.743^h	11 890^h		−180.0^h	86 000^h				−0.291^j	−0.139^j	0.152^j
	−0.597^j	15 796^j	5858^j	−380.1^j	250 000^j

---

Fig. 2 Structures of NiBDT-based (right column) and their corresponding non-NiBDT-based (left column) derivatives, optimized with the B3LYP/6-31G*/Ni(SDD) method, in the gas phase.

Fig. 3 Structures of NiBDT-based derivatives, optimized with the B3LYP/6-31G* method, in the gas phase.

The basic parameters using which we study the above issues are, besides the (hyper)polarizabilities, the dipole moment, the E_HOMO, E_LUMO, |E_HOMO − E_LUMO|(|H–L|) gap and the distribution of spin density. Of particular importance is the diradical character (DC); for its quantification we have employed two indices. The variety of systems we have studied and the properties we have computed highlight the comparative character of the present work.

1. NiBDT fused with benzene

1.1 Linear models. The structures of (B-NiBDT)_n, where n = 1–6, were optimized by the B3LYP/6-31G* method (Table 1). The properties were computed by employing the following functionals: BHandHLYP, B2PLYP, CAM-B3LYP, LC-PBE, ωB97X and LC-BLYP. Pilot computations have also been performed with the CCSD and CCSD(T) methods. In connection with the above techniques the 6-31G* basis set was used.

E _HOMO gets a maximum at n = 2 for BHandHLYP and B2PLYP, while for CAM-B3LYP, E_HOMO gets a maximum at n = 3.

E _LUMO decreases with n for BHandHLYP and CAMB3LYP.

|H–L|. It decreases with n for BHandHLYP, CAM-B3LYP and ωB97X.

Effect of a uniform external electric field on |HLG|↑ and |HLG|↓. It has been shown that upon application of an electric field the spin states on the two edges of a ribbon may be affected.^110–112 We have considered the effect of a uniform electric field on |HLG| of (B-NiBDT)₆, for spin α(↑) and spin β(↓) orbitals, |HLG|↑ and |HLG|↓, respectively. The main points of Fig. 4 are:

Fig. 4 The |E_HOMO − E_LUMO| gap (HLG) for α spin (↑) and β spin (↓) orbitals of the (B-NiBDT)₆ derivative as a function of the field. All field values (F = F_x) were computed with the (U)B3LYP/6-31G* method at the B3LYP/6-31G* gas phase optimized geometry.

(a) |HLG|↑ = |HLG|↓, for F = 0.

(b) If F > 0, then |HLG|↑ ≠ |HLG|↓; for example for 14 × 10⁻⁴ < F < 48 × 10⁻⁴ a.u., then |HLG|↑ > |HLG|↓; for F = 48 × 10⁻⁴ a.u., we observe |HLG|↑ = |HLG|↓, and for 48 × 10⁻⁴ < F ≤ 64 × 10⁻⁴ a.u. |HLG|↑ < |HLG|↓.

(c) The difference |HLG|↑ − |HLG|↓ is small/negligible.

E _exc. This property decreases with increasing n and seems to reach a constant value (ca. 0.9 eV) (Table 1). In Fig. 5 the absorption spectrum of (B-NIBDT)_n, n = 1–6, is shown. A strong absorption peak was computed, ranging from 646 to 1323 nm. For higher energies (shorter wavelengths) the observed absorptions are associated with small oscillator strengths (<0.1).

Fig. 5 Absorption spectra of several (B-NiBDT)_n, n = 1–6, derivatives, computed with (U)CAM-B3LYP/6-31G* at the B3LYP/6-31G* gas phase optimized geometry.

Dipole moments. Significant differences are observed between the dipole moment values calculated by some of the employed functionals for n = 1 (Table 1). A comparison of the dipole moment values, computed by the CCSD and CCSD(T) methods, shows the significant effect of triples (T). The dipole moment data and n of [(B-NiBDT)_n] are fitted with the function: |μ_x| = −0.020n² + 0.227n, R² = 0.9874 (Fig. 6).

Fig. 6 The dependence of the dipole moment, polarizability and first hyperpolarizability on the number of (B-NiBDT) units (n). The values were computed with the (U)CAM-B3LYP/6-31G* method. The black line depicts the fitted curve.

Polarizabilities. A comparison of the CCSD and CCSD(T) α_xx values shows that the perturbative triples (T) have a significant effect (Table 1). The five employed functionals give values for α_xx in the range 369–488 a.u. for n = 1. There is a remarkable dependence of α_xx on the selected functional. For n = 1–3, the maximum α_xx value is computed by B2PLYP. There is reasonable agreement for n = 1–3 for α_xx, computed with CAM-B3LYP and B2PLYP. The polarizability data and n of [(B-NiBDT)_n] have been connected with the function α_xx = 429.83n^1.9674, R² = 0.9996 (Fig. 6).

The variation of α_xx(−ω;ω) for (B-NiBDT)₃ is shown in Fig. S1 (ESI†). This property gets a maximum value for λ = 1100 nm. The computed λ is in agreement with the corresponding value (1097 nm) of Table 1.

In Fig. 7 we present atomic intrinsic polarizabilities (INPs) of (B-NiBDT)₂ together with the computed FOHI atomic charges and the corresponding MIP densities including those of its constituent units. In this treatment, we used a triple ζ quality basis set, namely 6-311G*, for all atoms. At the respective level of theory, the computed total mean molecular polarizability of (B-NiBDT)₂, which amounts to 768.4 a.u., is decomposed into a dominant charge transfer part (CT) of 679.8 a.u. (∼77%) and a size independent intrinsic polarizability part (IN) of 106.6 a.u. (∼13%). As in all prolate shaped molecules, the CT part, which defines the magnitude of their dipole polarizability, is dominated by charge transfer interactions along the longitudinal axis of the molecule. On the other hand, the intrinsic polarizability of (B-NiBDT)₂ is dictated by the individual atomic INPs of Ni and S atoms. Their overall contribution accounts for 48% of the total molecular INP of (B-NiBDT)₂. By comparing the illustrated integrated values on each atom (Fig. 7), we see that for this molecule: INP(S) > INP(Ni) > INP(C) > INP(H). The respective ordering, also observed in the isolated NiBDT unit (shown at the bottom of Fig. 7), follows only partially the number of electrons (or size) of each atom, and thus their atomic polarizabilities, in their free forms. Hence, while Ni and S (Ni : 28e, S : 18e) are found to be more polarizable than all C and H atoms of this system, the INP of S and Ni appears to follow the opposite trend. The relative INPs loosely correlate with the charges on each of these atoms in the sense that the more positively (negatively) charged an atom is in the molecule the less (more) polarizable it becomes. This trend, though, does not apply to atoms of the same type. For instance, the negatively charged S atoms attached to the two carbon hexagons bear slightly smaller INPs than those linked to the terminal HCCH unit which remains intrinsically as polarizable as in free NiBDT. Nevertheless, the observed differences are small and are not expected to play an important role in the behavior of the molecule. On the other hand, the INP of Ni in (B-NiBDT)₂ increases considerably in comparison with free NiBDT. The opposite trend is observed in the case of S atoms which lose some of their polarizable character. Additional information about the character of (B-NiBDT)₂ in terms of polarizability can be retrieved by comparing its mean intrinsic polarizability densities (MIP) with those of its constituent units, namely, benzene and NiBDT. As explained elsewhere,⁹⁹ positive MIP densities imply easily polarizable regions of the molecules. On the other hand, negative densities stand for “nonpolarizable” zones in which the electron density is expected to reduce under an external field perturbation. The illustrated contour maps suggest that the most affected areas of (B-NiBDT)₂, with respect to the free molecules, are those in the vicinity of Ni atoms and on the central carbon hexagon. The latter unit is expected to bear a considerably reduced aromatic character as compared to the free molecule of benzene or aromatic hexagons of larger graphene like molecular systems.¹⁰²

Fig. 7 FOHI charges, intrinsic atomic polarizabilities (INPs) and MIP density contour maps of (B-NiBDT)₂, C₆H₆ and NiBDT computed at the (U)CAM-B3LYP/6-311G* level.

First hyperpolarizabilities. A comparison of CCSD(T) with the CCSD β_xxx values shows that the triples (T) have a significant effect (Table 1). The five considered functionals give values for β_xxx in the range 2.1–10.9 × 10³ a.u. There is a remarkable change of sign for β_xxx. Curve fitting (Fig. 6) has disclosed a rather complex function connecting β_xxx with n[(B-NiBDT)_n].

We have employed a two-state model (TSM)¹¹³ in order to rationalize the first hyperpolarizability results of (B-NiBDT)_n. The geometries were computed with the B3LYP/6-31G* method, while the β_xxx values were calculated with the CAM-B3LYP/6-31G* approach (Table 2). For this model we have computed the dipole moment of the ground and excited states, μ_g and μ_e, respectively, the transition dipole moment, μ_ge, and the excitation energy, ΔE. For n = 1, we have employed several excited states; for n = 2, we have used 2 excited states and for the rest of the oligomers, one. Of course, as one would expect the value of β_xxx may greatly depend on the selected excited state. Both the numerically computed and the TSM β_xxx values follow the same trend; however, the first gets a maximum value for n = 2, while the second one approaches a maximum for n = 3. ΔE and (μ_g)_x decrease with n, while |(μ_ge)_x| increases; (μ_e)_x and (μ_g)_x − (μ_e)_x get a maximum value for n = 2, that is they follow a similar trend to β_xxx (Table 2).

Table 2 The transition dipole moment (μ_ge), excitation energy (ΔE) and dipole moments of the ground (μ_g) and the excited state (μ_e)^a of (B-NiBDT)_n, n = 1–6. For comparison the corresponding data for hydrogen fluoride are also reported. The property values were computed at the B3LYP/6-31G* optimized geometry with the (U)CAM-B3LYP/6-31G* method. The description of states is given in Table S4 (ESI). All property values are given in a.u.

P	N						HF
	1	2	3	4	5	6
a This state is associated with the first allowed electronic transition. b State 5. c State 42. d State 49. e State 53. f State 11. g β xxx = [3(μe − μg)(μge)^2]/(ΔE)^2, ΔE = Eexc − Eg. h Numerical value.
ΔE	0.070^a 0.053^b 0.169^c 0.188^d 0.197^e	0.046^a 0.050^f	0.041^a	0.038^a	0.036^a	0.034^a	0.558^a
|(μge)x|	1.491^a 1.279^b 0.955^c 2.206^d 1.152^e	3.924^a 2.366^f	7.047^a	9.629^a	11.911^a	13.936^a	0.935^a
(μe)x	−0.907^a 2.137^b 0.432^c −0.090^d 0.050^e	3.028^a −4.032^f	0.983^a	−0.092^a	−0.650^a	−0.965^a	0.555^a
(μg)x	−0.231	−0.392	−0.497	−0.571	−0.624	−0.661	−0.747
(μe)x − (μg)x	−0.676^a 2.368^b 0.663^c 0.141^d 0.281^e	3.420^a −3.640^f	1.480^a	0.479^a	−0.026^a	−0.304^a	1.302^a
β xxx	−906^a 4190^b 63^c 58^d 28^e	73 320^a −24 222^f	127 580^a	91 021^a	−8540^a	−149 540^a	10.2^a
β xxx	6047	17 100	13 700	−29 560	−117 460	−233 850	16.9

---

For comparison we also report the first hyperpolarizability value of HF (hydrogen fluoride). The |β_xxx| value of (B-NiBDT)₆ is 13 837 times larger than that of HF. How can we interpret this large difference using the TSM? It appears that the properties that could be used to rationalize the much larger |β_xxx| of the oligomer are ΔE and |(μ_ge)_x|; the first(denominator) is much smaller and the second(numerator) is much larger than the corresponding properties of HF (Table 2).

Second hyperpolarizabilities. Triples (T) have a significant effect (Table 1). The maximum γ_xxxx value for n = 1–3 is computed with B2PLYP. The value of γ_xxxx depends, remarkably, on the functional. For most of the considered cases the minimum γ_xxxx is given by BHandHLYP and the maximum by ωB97X. The ratio γ_xxxx (ωB97X)/γ_xxxx (CAM-B3LYP) equals 1.3 for n = 1, and, stabilizes at 1.8 for n ≥ 5. Both ωB97X and CAM-B3LYP are long-range corrected functionals. We observe that the ratio γ_xxxx (n + 1)/γ_xxxx (n) decreases as n increases. So, n₂/n₁ = 28.9, but n₆/n₅ = 2.1 (CAM-B3LYP).

We have considered two indices in order to quantify the change of γ_xxxx with the size of the derivative: (a) the first index is I₁ = γ_xxxx/n, where n is the number of benzene (B)–NiBDT units. I₁ demonstrates that CAM-B3LYP leads to much more polarizable units than BHandHLYP. (b) The second index is I₂ = γ_xxxx/N, where N is the number of electrons of the derivative. Both I₁ and I₂ clearly demonstrate that the considered structures (i) are extremely polarizable and (ii) there is a rapid, nonlinear, increase of I₁ and I₂ with the size of the structure. One of the main points, emerging from Table 1, is that some of the oligomers (B-NiBDT)_n have giant second hyperpolarizabilities. The hyperpolarizability data and n[(B-NiBDT)_n] have been connected with the function γ_xxxx = 19.641n^5.0248, R² = 0.999 (Fig. 8).

Fig. 8 The dependence of the second hyperpolarizability on (a) the number of (B-NiBDT) units (n), (b) the excitation energy (E_exc), and (c) E₂. The values were computed with the (U)CAM-B3LYP/6-31G* method, except those of E₂, where the (U)B3LYP/6-31G* method was employed.

Effect of benzene on the (hyper)polarizabilities. We have found that α_xx and γ_xxxx of NiBDT are 221.9 and 558 × 10³ a.u., respectively.⁷⁴ The corresponding properties of NiBDT-B ((U)HBandHLYP) are 371 and 19 × 10⁵ a.u., respectively. Thus, benzene has a very large effect on both α_xx and γ_xxxx. Furthermore, one may consider two derivatives of NiBDT in which (i) one H has been substituted with CH₃ and (ii) 2 hydrogen atoms have been substituted with the ligand CS₃. These derivatives have the following γ_xxxx values: 739 × 10³ a.u. and 1535 × 10³ a.u., respectively. All these properties have been computed with the (U)BHandHLYP/6-31G*.⁷⁴ The significant effect of benzene, as a substituent, is observed.

Diradical character (DC). We used two indices, E₁ and E₂, to estimate the DC of the considered derivatives (Table 1). These indices show that DC increases with increasing n. Moreover E₂ is associated with the exchange interaction, J, (E₂ = 2J), which can be experimentally measured by electron spin resonance (ESR).⁸⁶ It has been found that E₁ and E₂, for NiBDT, are 0.05 and 9.91 kcal mol⁻¹, respectively.⁷⁴ Comparing E₂ of NiBDT with that of B-NiBDT, we infer that the fusion of benzene with NiBDT leads to a small decrease of DC. However, for 2 ≤ n, all the oligomers have larger DC than NiBDT.

Effect of a uniform external electric field on DC. Su et al.¹¹⁴ reported that they controlled the singlet–triplet energy gap of a diradical in the solid state by temperature, demonstrating that an external stimulus can modify the DC. Here we shall comment on the effect of a uniform external electric field on DC. Most considered field values lead to an increase of ΔE(F) = E(F)_(RDFT) − E(F)_(UDFT), that is an increase of F leads to an increase of DC, in comparison to the value that is associated with F = 0 (Fig. 9).

Fig. 9 The dependence of ΔE = E(RDFT) − E(UDFT) on the field (F = F_x) for the (B-NiBDT)₆ derivative. All values were computed with the (U)B3LYP/6-31G* method at the B3LYP/6-31G* gas phase optimized geometry.

Spin density distribution (triplet state). It was suggested that, in a singlet diradical, the weakly coupled electrons, and in a triplet radical, the corresponding electrons are located at the same position of the molecule.⁸⁴ We shall adopt the same hypothesis and comment on the Mulliken spin densities of the lowest triplet state of (B-NiBDT)_n, where n = 1–3. For n = 1, we observe (Fig. 10) a large spin density at nickel and sulfur atoms; at the carbon atoms the spin density is smaller or negligible. For comparison, we have also computed the spin densities of NiBDT. We found that the spin density at Ni is very small, while at the sulfur atoms the spin density is, approximately, twice that of carbon atoms. For n = 2, the total spin densities at sulfur and nickel atoms are 1.004 and 0.68, respectively; one of the two nickel atoms has negligible spin density. At the carbon atoms the spin density is negligible. For n = 3, the total spin densities at sulfur and nickel are 1.046 and, 0.597, respectively. The spin density is concentrated at the middle Ni atom. At the carbon atoms the spin density is very small or negligible. Thus, these computations show that the unpaired electrons are localized at the sulfur and nickel atoms. However, our VB computations for NiBDT suggest that the unpaired electrons are localized at the carbon atoms.⁷⁴ In Fig. S2 (ESI†) the spin density distribution for (B-NiBDT)_n, n = 2, 6, is also depicted.

Fig. 10 Computed Mulliken triplet spin densities with the UB3LYP/6-31G* method.

Substitution of Ni with Pd. The results of Table 3 present some properties of (B-MBDT)_n, where M = Ni, Pd and n = 2, 3. An increase of the atomic number of the metal (going from Ni to Pd) decreases the values of μ_x, α_xx, β_xxx, and γ_xxxx (Table 3). However, the Pd derivatives are associated with much larger DC. This finding confirms our previous computations related to the DC and (hyper)polarizabilities of A-[(SC)₂-M-(SC)₂]-A, where A = SC₃C₂ and M = Ni, Pd.¹¹⁵ The results of Table 1 show a different trend, that is increase of γ_xxxx is connected with an increase of DC.

Table 3 The dipole moment, (hyper)polarizabilities, E_HOMO/LUMO, E_HOMO − E_LUMO gap (|H–L|) of (B-MBDT)_n, n = 2, 3 and M = Ni, Pd. The properties (P) were computed at the B3LYP/6-31G* optimized geometry. All property values are given in a.u.

n/M	P
	μ x	α xx	β xxx × 10^−3	γ xxxx × 10^−5	E 1	E HOMO	E LUMO	|H–L|
a Method: (U)CAM-B3LYP/6-31G*. b E 1 (kcal mol^−1) = E(RDFT) − E(UDFT). Method: B3LYP/6-31G*.
2
Ni Pd	−0.392 −0.536	1684 1496	17.1 7.1	570 332	0.37 1.63	−0.263 −0.264	−0.139 −0.140	0.124 0.124
3
Ni Pd	−0.497 −0.641	3836 3162	13.7 3.7	5380 1946	1.14 3.26	−0.256 −0.266	−0.148 −0.146	0.118 0.120

---

From the results of Table 3 it is also observed that (γ_xxxx(n = 3)/γ_xxxx(n = 2)) [Ni] = 9.43, (γ_xxxx(n = 3)/γ_xxxx(n = 2)) [Pd] = 5.86. Similarly for E₁ it was found that (E₁(n = 3)/E₁(n = 2)) [Ni] = 3.08, (E₁(n = 3)/E₁(n = 2)) [Pd] = 2. The larger increase of γ_xxxx [Ni], compared with that of γ_xxxx [Pd], with the number of (B-MBDT) units, is associated with the larger increase of E₁[Ni] compared with that of E₁[Pd].

1.2 Planar models. The γ_xxxx values computed with the considered functionals are in qualitative agreement; smaller differences have been found for the α_xx values (Table 4). The structure of 1 (Fig. 1) is a model of the Ni₃S₁₂C₁₂ nanosheet.⁵ Derivatives 1 and 2 have 12 and 6 NiBDT units, respectively (Fig. 1). Thus, it is not surprising that γ_xxxx(1)/γ_xxxx(2) = 23.3 (ωB97X). Although both 2 and 3 have 6 units of benzene and NiBDT, their structure is very different (Fig. 1). This structural difference greatly affects the second hyperpolarizability, that is γ_xxxx(3)/γ_xxxx(2) = 23.9; a smaller difference is observed for the α_xx values; the corresponding ratio is 2.4 (ωB97X). It is observed that 3 has a remarkably greater DC than 2 (CAM-B3-LYP). A comparison of the second hyperpolarizabilities of the linear (B-NiBDT)₆ (Table 1) and the corresponding planar derivative (2; Table 3 and Fig. 4) shows that γ_xxxx(linear) = 483 γ_xxxx(planar); method: UCAM-B3LYP/6-31G*. The linear compound has larger E₁ (larger DC) and smaller E_exc, in comparison with the planar one (Tables 1 and 4), highlighting the great effect of the structure on the second hyperpolarizability.

Table 4 Components of the (hyper)polarizabilities (a.u.), E₁ (kcal mol⁻¹) and the excitation energy (E_exc/eV) of several planar nickel bis(dithiolene) complexes. All the reported properties were computed using the B3LYP/6-31G* optimized geometry

Method/basis set	Property
	α xx	α yy	β xxx (×10^−3)	γ xxxx (×10^−5)^a	E 1	E exc
a Value computed by taking the first derivative of βxxx/yyy with respect to the field. b E 1 (kcal mol^−1) = E(RDFT) − E(UDFT), method: B3LYP/6-31G*. c Fig. 4. d β yyy component. e γ yyyy component.
Compound: 1^c
(U)BHandHLYP/6-31G* (U)ωB97X/6-31G* (U)CAM-B3LYP/6-31G*	5895 7350	6455 7322	—	3860 8005	0.12	1.16
Compound: 2^c
(U)BHandHLYP/6-31G* (U)CAM-B3LYP/6-31G* (U)ωB97X/6-31G*	2550 2837 2858	2550 2837 2858	—	204 288 343	0.63	1.57
Compound: 3^c
(U)CAM-B3LYP/6-31G* (U)ωB97X/6-31G*	6327 6761	3959 2477	99.3 (−8.2)^d 163.2 (−17.3)^d	5014 (−3500)^e 8200 (−1480)^e	1.41	1.19

---

2. NiBDT fused with linkers having extensive conjugation

We have selected 4 different conjugated organic linkers (L): C₆H₆, C₁₄H₈, C₂₄H₁₂, and C₂₄H₈ (Table 5). We have then designed 8 oligomers involving, first, one of the above linkers, fused with two NiBDT units and, second, 3 NiBDT units fused with two linkers (Fig. 2). These computational experiments have been performed in order to show the effect of the linker, and the conjugation pattern in particular, on the (hyper)polarizabilities. For (NiBDT)₂-L, we observe that γ_xxxx = 270 ± 46 × 10⁵ a.u. (UCAM-B3LYP/6-31G*), while for (NiBDT)₃-L₂, three of the linkers (C₁₄H₈, C₂₄H₁₂, and C₂₄H₈) give rather similar γ_xxxx values; a significantly different and much larger second hyperpolarizability is associated with C₆H₆ (Table 4). An approximate estimate of the interaction induced γ_xxxx for NiBDT-C₂H₂-NiBDT is 25.4 × 10⁶ a.u. (UCAM-B3LYP/6-31G*). The γ_xxxx values for NiBDT and C₆H₆ are 77 200 a.u. and 1946 a.u., respectively. The largest α_xx values for both (NiBDT)₂-L and (NiBDT)₃-L₂ are associated with C₂₆H₁₂. An increase of n of the oligomer (NiBDT)_m-L_n leads to an increase of DC and a decrease of the excitation energy (Table 5).

Table 5 The (hyper)polarizabilities (a.u.), E₁ (kcal mol¹), E_exc (eV) of NiBDT fused with various conjugated linkers. ^a The properties of the linkers are also reported. The values were computed at the B3LYP/6-31G* optimized geometry

Derivative	Property^b					Derivative	Property^b
	α xx	γ xxxx × 10^−5	(γxxxx × 10^−5)/N^c	E 1	E exc		α xx	γ xxxx × 10^−5	(γxxxx × 10^−5)/N^c	E exc
a Fig. 2. b Method: (U)CAM-B3LYP/6-31G*. c N: total number of electrons. d E 1 (kcal mol^−1) = EDFT(singlet) − EUDFT(singlet); method: B3LYP/6-31G*. e RA: radiaannulene, NiBDT: bis(ethylene-1,2-dithiolato)nickel. f Values were taken from J. Phys. Chem. C, 2016, 120, 9419–9435. g z-Component.
NiBDT-C2H2-NiBDT^e NiBDT-C2H4-NiBDT^e	1333 704	256 45	1.11 0.41	0.50 0.19	1.28 1.97	C6H6	71	0.019	0.00046	7.47
(NiBDT)3-(C2H2)2 (NiBDT)3-(C2H4)2	3318 1215	3480 125	9.94 0.35	1.34	1.12 1.91	C10H8	150	0.180	0.00265	6.29
NiBDT-C10H4-NiBDT	1297	280	1.00	1.26	1.99	C14H8	179	0.354	0.00385	7.35
(NiBDT)3-(C10H4)2	2755	1425	3.17	2.63	1.66	C26H12	459	8.143	0.04847	6.40
NiBDT-C20H8-NiBDT	1775	316	0.92	0.63	2.04	C24H12	465	6.567	0.04209	6.60
(NiBDT)3-(C20H8)2	3834	1506	2.61	1.39	1.88	C46H20	1364	168.438	0.56904	2.34
NiBDT-RA-NiBDT^e	1592^f	224^f	0.66	0.50^f	1.52^f	C24H8	339 (712)^g	0.952 (15.440)^g	0.00626	2.70
(NiBDT)3-(RA)2	3535^f	1517^f	2.66	0.63^f	1.43^f	C46H12	845 (1212)^g	19.137 (5.749)^g	0.06645	2.83

---

In order to show the effect of the delocalized π electrons – in particular those connected with the linkers – on the properties of interest, we fused NiBDT with C₆H₈ (Fig. S3, ESI†). This linker allows the disconnection of the π-systems of NiBDT units. The results show that there is a very large effect on γ_xxxx, for example:

γ_xxxx[(NiBDT)₃–(C₂H₂)₂]/γ_xxxx[(NiBDT)₃–(C₂H₄)₂] = 28.

This effect increases with n. The other properties (e.g. α_xx) are associated with a significant, but less pronounced, effect.

In Fig. 11 we contrast the intrinsic atomic polarizabilities and the +MIP densities of NiANi and NiBNi (Fig. 11). In these planar polyaromatic systems, the CAM-B3LYP/6-311G* CT part of the dipole polarizabilities accounts for 84% (714.27 a.u.) and 82% (522.23 a.u.) of their total mean molecular polarizabilities, respectively. The stronger CT contribution obtained for NiANi is in accord with its larger size in comparison with NiBNi in both its molecular directions. As in (B-NiBDT)₂, Ni and S are intrinsically the most polarizable atoms in both NiANi and NiBNi, with INP(S) > INP(Ni). In either molecule, out of the eight S atoms those which are bonded to the polyaromatic linkers appear to be the most polarizable ones. In terms of element contributions to the global molecular intrinsic polarizability of NiANi/NiBNi, FOHI analysis suggests that Ni and S atoms contribute 53%/62% of the total molecular intrinsic polarizability. In this case the contribution of the latter elements follows the inverse of the corresponding molecular sizes. In contrast to what we saw in (B-NiBDT)₂, when NiBDT becomes a part of either DCPP or pyracyclene, there is a notable decrease in the intrinsic polarizability of the Ni atoms as compared to free NiBDT. On the other hand, while the inner S atoms show an increase in their intrinsic polarizabilities, the side S atoms retain most of their character in terms of INP. Considering the main fragments of either system, our results suggest that the global intrinsic polarizabilities of the carbon frames of the linkers in NiANi/NiBNi, calculated as the sum of all atomic INPs, slightly decrease after they are functionalized with NiBDT. The most affected C atoms are those connected to the S atoms which lose about 60% of their initial INPs. Turning our attention to Fig. 11, where the computed MIP density distributions are illustrated, we see that both linkers have similar effects on the attached NiBDT. In these plots we also see that the MIPs of the linkers exhibit small variations in NiANi and NiANi with respect to their free forms. The most affected areas are spotted in the vicinity of the terminal CC bonds as is also suggested by the INPs discussed above. Finally it would be interesting to note that when these densities are compared to those of the (B-NiBDT)₂ compound we see that there is a change in the distribution of the positive and negative areas which should be related to the decreased Ni INPs in the NiANi and NiBNi compounds in comparison with (B-NiBDT)₂.

Fig. 11 FOHI charges, intrinsic atomic polarizabilities (INPs) and MIP density contour maps of NiBDT-C₂₀H₈-NiBDT (NiANi), NiBDT-C₁₀H₄-NiBDT (NiBNi), DCPP and pyracyclene computed at the (U)CAM-B3LYP/6-311G* level.

3. H₄C₄O₂, H₄C₄S₂, H₄C₄O and H₄C₄S fused with NiBDT

We have experimented with a further series of linkers: H₄C₄O₂ (I), H₄C₄S₂ (II), H₄C₄O (III) and H₄C₄S (IV). H₄C₄O₂ and H₄C₄S₂ were selected because we would like to find the effect of replacing CH (in benzene) with oxygen and sulfur (Table 6 and Fig. 3). H₄C₄O and H₄C₄S were selected in order to study the effect of replacing a six-membered ring with a five-membered one. All the molecules are planar, except those of II, which has a boat-shape (Fig. S4, ESI†). These derivatives have their dipole moment oriented along the x axis. For computation of all the molecular properties the CAM-B3LYP/6-31G* method has been employed. In addition, for the properties of IV(a,b), the BHandHLYP/6-31G* functional has been used. These functionals gave similar property values except for β_zzz.

Table 6 The dipole moment (a.u.), (hyper)polarizabilities (a.u.), E₁ (kcal mol⁻¹) and E_exc (eV) of NiBDT units, fused with various linkers.^a All property values were computed at the B3LYP/6-31G*/SDD(Ni) optimized geometry

Method/basis set	Property
	μ x	α zz	α xx	β xxx	β zzz	γ zzzz × 10^−5	γ xxxx × 10^−5	E 1	E exc
a Fig. 3. b The dipole moment is oriented along the x-axis (Fig. 3). c E 1 (kcal mol^−1) = E(RDFT) − E(UDFT), method: B3LYP/6-31G*.
(B-NiBDT)3
(U)CAM-B3LYP/6-31G*	0.497	804	3836	−13 700	1	−1100.0	5380.0	1.14	1.13
C14H4Ni3S12O6 (I)
(U)CAM-B3LYP/6-31G*	1.242	463	1685	29 880	0	−0.4	645.8	5.65	1.47
C14H4Ni3S12S6 (II)
(U)CAM-B3LYP/6-31G*	1.185	1259	1642	26 079	18 882	285.9	549.8	3.14	1.25
C14H4Ni3O3S12 (III)
(U)CAM-B3LYP/6-31G*	1.004	1724	940	14 935	−1010	574.9	74.2	3.14	1.06
C14H4Ni3S15 (IVa)
(U)BHandHLYP/6-31G* (U)CAM-B3LYP/6-31G*	1.081 1.124	1665 1912	864 923	8794 6021	−1260 −8811	591.0 889.0	43.4 59.5	1.25	1.88
C14H4Ni3S15 (IVb)
(U)BHandHLYP/6-31G* (U)CAM-B3LYP/6-31G*	0.664 0.661	572 577	2009 2366	18 300 22 910	453 240	−0.9 −1.2	1013.0 1100.0	1.57	1.82

---

I and II have similar properties along the x-direction, but present significant differences in their properties along the z-direction (α_zz, β_zzz, γ_zzzz). Their E₁ values are also remarkably different. The above pattern of changes could be attributed (i) to the different linkers, that is H₄C₄O₂ is the linker in I and H₂C₄S₂ is the linker in II and (ii) the different structures of I and II. It was found that I has a linear structure, while II has a boat-shaped one (Fig. S4, ESI†). Benzene (C₆H₆) as a linker leads to much larger γ_xxxx values in comparison with 1,4-dioxin (H₂C₄O₂) and 1,4-dithiin (H₂C₄S₂). The same trend is observed for α_xx.

In Fig. 12 we show the intrinsic atomic polarizabilities of compounds I and II (Fig. 3). Contrary to what was observed in the previous set of models the Ni atoms do not feature the largest intrinsic polarizability. Depending on their position in the framework of these chains, Ni loses large portions of its NPs. In terms of element contribution for either II/I, more than 50% of the molecular intrinsic mean polarizability originate from S atoms. In both cases, the atomic polarizabilities of Ni and S atoms notably vary depending on their position. The main trend holding for both systems is that the Ni and S atoms located at the outer chain sides are always more polarizable, and thus more sensitive to external field perturbations than those located close to the center of each chain. Also, there is a clear difference in the character of the various S atoms since FOHI analysis suggests that out of the 18 S atoms present in I, those which are directly bonded to Ni appear to be more polarizable than the S atoms of the embedded dithiin units. Finally, by comparing the average INP of O and S in free (1,4)-dithiine and in (1,4)-dioxine with the corresponding values in compounds I and II, it is revealed that there is a decrease of the average INP in all atoms of both fragments when they are present in I and II linear chains.

Fig. 12 FOHI intrinsic atomic polarizabilities and MIP density contour maps of C₁₄H₄Ni₃S₁₂S₆ (II; Fig. 3) and C₁₄H₄Ni₃S₁₂O₆ (I; Fig. 3) compounds computed at the UCAM-B3LYP/6-311G* level. Numbers in parentheses correspond to the INPs of (1,4)-dithiine and (1,4)-dioxine computed at the RCAM-B3LYP/6-311G* level.

It is remarkable that II and III have, approximately, the same DC, as this is computed by E₁, but several of their other properties have significant differences (e.g. β_zzz, γ_xxxx). A similar trend has also been observed in other cases (Table 5; e.g. (NiBDT)₂RA/(NiBDT)₃(RA)₂). These observations may be associated with the fact that the (hyper)polarizabilities are more sensitive probes, than the DC, for changes taking place in the structure of the molecule. A very small γ_xxxx value for III and IVa is observed. The (hyper)polarizabilities of IVa and IVb are quite different; these oligomers have a similar DC.

It is also noted that a comparison of the pairs I/II and III/IVa shows that the linkers with oxygen lead to oligomers with higher DC, in comparison to those with sulfur.

The employed linkers highlight the superiority of benzene for the production of oligomers with much larger second hyperpolarizability (γ_xxxx) in comparison with the other ones.

4. Vibrational contributions

We have chosen two relatively small models (Fig. S5, ESI†) to study both the electronic and vibrational contributions to α_xx and γ_xxxx. The computed static nr contributions to the non-vanishing longitudinal diagonal component of the electrical properties computed for the spin-broken unrestricted DFT wavefunction are shown in Table 7, together with the purely electronic contributions. They were computed with the same basis set/pseudopotential combination as the electronic properties in the rest of this paper. To be consistent with the majority of the electronic results, the CAM-B3LYP functional was first tried for both molecules. It turned out, however, that the spin-broken solution of the C₈O₄NiS₄H₄ molecule (Fig. S4, ESI†) obtained with this functional showed no convergence for the property values obtained from the Romberg method. In order to have a full set of comparable values, the B3LYP functional was then applied for both molecules, where the problems observed with CAM-B3LYP did not occur. As shown in Table 7 for molecule C₈S₄NiS₄H₄ (Fig. S4, ESI†), the differences in the properties calculated with the two functionals are not very large (up to a factor of ∼2). While the electronic contributions to α and γ are of similar magnitude, the nr contributions are much larger for molecule C₈S₄NiS₄H₄ than for C₈O₄NiS₄H₄, especially for γ, which is about four orders of magnitude larger for C₈S₄NiS₄H₄ than for C₈O₄NiS₄H₄.

Table 7 Electronic (el) and nuclear relaxation (nr) contribution to the non-vanishing longitudinal diagonal components of the electrical properties of C₈X₄NiS₄H₄, X = S, O, computed with the spin-broken unrestricted DFT wavefunction (U) as well as with the restricted wavefunction (R) using the B3LYP and CAM-B3LYP functionals. All property values are given in a.u.

Molecule	C8S4NiS4H4^a		C8O4NiS4H4^a
Method	(U)B3LYP		(U)B3LYP		RB3LYP
Property	el	nr	el	nr	el	nr
a Fig. S1 (ESI). b Method: (U)CAM-B3LYP. c Method: RCAM-B3LYP.
α xx	812	527	430	42.6	599	123
	659^b	533^b
	804^c	399^c
γ xxxx (×10^6)	2.5	317	1.65	0.13	−0.16	−0.69
	4.2^b	696^b
	0.92^c	41^c

---

In addition to the nr contribution to the electrical properties for the spin-broken unrestricted DFT wavefunction, we also computed the corresponding properties for the restricted DFT wavefunction, with CAM-B3LYP for molecule C₈S₄NiS₄H₄ and with B3LYP for molecule C₈O₄NiS₄H₄. It turns out that the two sets of properties differ substantially, both for the electronic as well as for the vibrational (nr) contributions (Table 7). While the polarizabilities are of similar order of magnitude, the second hyperpolarizabilities are about one order of magnitude larger with the unrestricted wavefunction.

IV. Conclusions

It is well known that there is a great need for novel photonic materials for use in a large number of applications (e.g. optical communications), which require an ever increasing improvement of their performance.

In this work we propose a novel family of such materials: 1D and 2D coordination polymers, based on nickel bis(dithiolene) and benzene. Factors of cardinal importance for the photonic materials are their NLO properties; these are computed by employing density functional theory.

Extensive computational experimentation has been performed, by employing linear and planar oligomers of the 1D and 2D polymers, in order to document and demonstrate the superior NLO performance of the proposed materials. For comparison, we have also used several rich π-electron linkers (e.g. C₆H₆,C₁₄H₁₈, C₂₄H₁₂, and C₄O₂H₄).

A property of key importance for this work is the second hyperpolarizability (γ). However, we computed several other properties (e.g. polarizabilities, E_HOMO, E_LUMO, and excitation energies), which have been correlated with γ and contributed to its rationalization. The diradical character, which was characterized by employing two indices, was also employed for the understanding of the properties of interest. This work mainly deals with the electronic contribution to (hyper)polarizabilities. However, we have also commented on the nuclear relaxation contribution to α_xx and γ_xxxx of some relatively small model compounds. While the nr contribution for the O-containing molecule was modest, quite large contributions were found for the S-containing molecule. The electronic contributions for both molecules, on the other hand, were of similar magnitude. In addition, it was found that the spin-broken wavefunction yields quite different vibrational and electronic properties than the restricted wavefunction. Considerable care has been taken for the validation of the employed computational procedure. This involved several functionals and state-of-the-art techniques (e.g. UCCSD(T)).

The main findings of this work are the following:

(a) The 1D and 2D coordination polymers, and in particular the linear ones, (B-NiBDT)_n, are likely to form the basis of a new class of efficient photonic materials with superior performance. There is a steep increase of γ with n and even for relatively small n, gigantic second hyperpolarizability values have been computed for (B-NiBDT)_n.

(b) We have found a great dependence of γ_xxxx on (i) the linker with which the NiBDT units are fused (benzene, as a linker, leads to the maximum observed effect) and (ii) the structure, for example γ(linear; n = 6)/γ(3; planar) = 28 (Tables 1 and 3).

(c) This work, in connection with the results of our previous article,⁹ demonstrates a mechanism for designing photonic materials with exceptionally large nonlinearities. This mechanism involves a series of units consisting of NiBDT, π-electron rich linkers, Li substituents⁹ and the appropriate structure, for example a 1D arrangement leads to a much larger second hyperpolarizability in comparison to the corresponding 2D one.

Conflicts of interest

There are no conflicts to declare.

References

1. Y. Chen, T. Bai, N. N. Dong, F. Fan, S. F. Zhang, X. D. Zhuang, J. Sun, B. Zhang, X. Y. Zhang, J. Wang and W. J. Blau, Prog. Mater. Sci., 2016, 84, 118–157.
2. X. F. Liu, Q. B. Guo and J. R. Qiu, Adv. Mater., 2017, 29, 1605886.
3. M. Feng, H. B. Zhan and Y. Chen, Appl. Phys. Lett., 2010, 96, 033107.
4. T. Kambe, R. Sakamoto, T. Kusamoto, T. Pal, N. Fukui, K. Hoshiko, T. Shimojima, Z. F. Wang, T. Hirahara, K. Ishizaka, S. Hasegawa, F. Liu and H. Nishihara, J. Am. Chem. Soc., 2014, 136, 14357–14360.
5. T. Kambe, R. Sakamoto, K. Hoshiko, K. Takada, M. Miyachi, J. H. Ryu, S. Sasaki, J. Kim, K. Nakazato, M. Takata and H. Nishihara, J. Am. Chem. Soc., 2013, 135, 2462–2465.
6. K. Hoshiko, T. Kambe, R. Sakamoto, K. Takada and H. Nishihara, Chem. Lett., 2014, 43, 252–253.
7. R. Matsuoka, R. Sakamoto, T. Kambe, K. Takada, T. Kusamoto and H. Nishihara, Chem. Commun., 2014, 50, 8137–8139.
8. Handbook of Thiophene-Based Materials: Applications in Organic Electronics and Photonics, ed. I. F. Perepichka and D. F. Perepichka, John Wiley & Sons, Ltd, Chichester, UK, 2009.
9. A. Avramopoulos, H. Reis, N. Otero, P. Karamanis, C. Pouchan and M. G. Papadopoulos, J. Phys. Chem. C, 2016, 120, 9419–9435.
10. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666–669.
11. K. Watanabe, T. Taniguchi and H. Kanda, Nat. Mater., 2004, 3, 404–409.
12. J. N. Coleman, M. Lotya, A. O'Neill, S. D. Bergin, P. J. King, U. Khan, K. Young, A. Gaucher, S. De, R. J. Smith, I. V. Shvets, S. K. Arora, G. Stanton, H. Y. Kim, K. Lee, G. T. Kim, G. S. Duesberg, T. Hallam, J. J. Boland, J. J. Wang, J. F. Donegan, J. C. Grunlan, G. Moriarty, A. Shmeliov, R. J. Nicholls, J. M. Perkins, E. M. Grieveson, K. Theuwissen, D. W. McComb, P. D. Nellist and V. Nicolosi, Science, 2011, 331, 568–571.
13. B. Lalmi, H. Oughaddou, H. Enriquez, A. Kara, S. Vizzini, B. Ealet and B. Aufray, Appl. Phys. Lett., 2010, 97, 223109.
14. E. Bianco, S. Butler, S. S. Jiang, O. D. Restrepo, W. Windl and J. E. Goldberger, ACS Nano, 2013, 7, 4414–4421.
15. A. Splendiani, L. Sun, Y. B. Zhang, T. S. Li, J. Kim, C. Y. Chim, G. Galli and F. Wang, Nano Lett., 2010, 10, 1271–1275.
16. K. F. Mak, C. Lee, J. Hone, J. Shan and T. F. Heinz, Phys. Rev. Lett., 2010, 105, 136805.
17. Y. M. Lin, C. Dimitrakopoulos, K. A. Jenkins, D. B. Farmer, H. Y. Chiu, A. Grill and P. Avouris, Science, 2010, 327, 662.
18. F. Bonaccorso, Z. Sun, T. Hasan and A. C. Ferrari, Nat. Photonics, 2010, 4, 611–622.
19. T. Y. Lu, H. Feng, S. W. Yang and J. C. Zheng, Comput. Mater. Sci., 2017, 126, 170–175.
20. T. Pal, T. Kambe, T. Kusamoto, M. L. Foo, R. Matsuoka, R. Sakamoto and H. Nishihara, ChemPlusChem, 2015, 80, 1255–1258.
21. A. J. Clough, J. W. Yoo, M. H. Mecklenburg and S. C. Marinescu, J. Am. Chem. Soc., 2015, 137, 118–121.
22. Q. Tang and Z. Zhou, J. Phys. Chem. C, 2013, 117, 14125–14129.
23. F. Shojaei, J. R. Hahn and H. S. Kang, Chem. Mater., 2014, 26, 2967–2974.
24. Z. F. Wang, N. H. Su and F. Liu, Nano Lett., 2013, 13, 2842–2845.
25. J. Zhou, RSC Adv., 2014, 4, 13361–13366.
26. F. Shojaei and H. S. Kang, J. Phys. Chem. C, 2014, 118, 25626–25632.
27. H. M. Liu, X. L. Li, L. Chen, X. L. Wang, H. Z. Pan, X. M. Zhang and M. W. Zhao, J. Phys. Chem. C, 2016, 120, 3846–3852.
28. A. J. Clough, J. M. Skelton, C. A. Downes, A. A. de la Rosa, J. W. Yoo, A. Walsh, B. C. Melot and S. C. Marinescu, J. Am. Chem. Soc., 2017, 139, 10863–10867.
29. D. Sheberla, L. Sun, M. A. Blood-Forsythe, S. Er, C. R. Wade, C. K. Brozek, A. Aspuru-Guzik and M. Dinca, J. Am. Chem. Soc., 2014, 136, 8859–8862.
30. M. G. Campbell, D. Sheberla, S. F. Liu, T. M. Swager and M. Dinca, Angew. Chem., Int. Ed., 2015, 54, 4349–4352.
31. J. S. Cui and Z. T. Xu, Chem. Commun., 2014, 50, 3986–3988.
32. S. Biswas, C. S. Tiwary, S. Vinod, A. K. Kole, U. Chatterjee, P. Kumbhakar and P. M. Ajayan, J. Phys. Chem. C, 2017, 121, 8060–8069.
33. P. Kumbhakar, A. K. Kole, C. S. Tiwary, S. Biswas, S. Vinod, J. Taha-Tijerina, U. Chatterjee and P. M. Ajayan, Adv. Opt. Mater., 2015, 3, 828–835.
34. K. P. Wang, J. Wang, J. T. Fan, M. Lotya, A. O'Neill, D. Fox, Y. Y. Feng, X. Y. Zhang, B. X. Jiang, Q. Z. Zhao, H. Z. Zhang, J. N. Coleman, L. Zhang and W. J. Blau, ACS Nano, 2013, 7, 9260–9267.
35. N. Liaros, A. B. Bourlinos, R. Zboril and S. Couris, Opt. Express, 2013, 21, 21027–21038.
36. L. R. Dalton, A. W. Harper, B. Wu, R. Ghosn, J. Laquindanum, Z. Y. Liang, A. Hubbel and C. Z. Xu, Adv. Mater., 1995, 7, 519–540.
37. S. Ellinger, K. R. Graham, P. J. Shi, R. T. Farley, T. T. Steckler, R. N. Brookins, P. Taranekar, J. G. Mei, L. A. Padilha, T. R. Ensley, H. H. Hu, S. Webster, D. J. Hagan, E. W. Van Stryland, K. S. Schanze and J. R. Reynolds, Chem. Mater., 2011, 23, 3805–3817.
38. P. K. Hegde, A. V. Adhikari, M. G. Manjunatha, C. S. S. Sandeep and R. Philip, Adv. Polym. Technol., 2011, 30, 312–321.
39. X. W. Zhan, Y. Q. Liu, D. B. Zhu, X. C. Liu, G. Xu and P. X. Ye, Chem. Phys. Lett., 2001, 343, 493–498.
40. M. Liu, H. S. Quah, S. C. Wen, J. Y. Wang, P. S. Kumar, G. Eda, J. J. Vittal and W. Ji, J. Mater. Chem. C, 2017, 5, 2936–2941.
41. Y. Matsuzaki, M. Nakano, K. Yamaguchi, K. Tanaka and T. Yamabe, Chem. Phys. Lett., 1996, 263, 119–125.
42. X. W. Zhan, Y. Q. Liu, D. B. Zhu, X. C. Liu, G. Xu and P. X. Ye, Chem. Phys. Lett., 2002, 362, 165–169.
43. C. M. Amb, C. L. Heth, S. J. Evenson, K. I. Pokhodnya and S. C. Rasmussen, Inorg. Chem., 2016, 55, 10978–10989.
44. L. Pilia, M. Pizzotti, F. Tessore and N. Robertson, Inorg. Chem., 2014, 53, 4517–4526.
45. T. B. Faust, P. M. Usov, D. M. D'Alessandro and C. J. Kepert, Chem. Commun., 2014, 50, 12772–12774.
46. X. B. Yang, L. Zhou, L. B. Huang, J. J. Xu, Y. Zhou, S. T. Han, Z. X. Xu, V. C. Y. Lau, M. H. W. Lam, W. Y. Wong and V. A. L. Roy, RSC Adv., 2013, 3, 12075–12079.
47. L. Hu, J. Qin, N. Zhou, Y. F. Meng, Y. Xu, J. L. Zuo and X. Z. You, Dyes Pigm., 2012, 92, 1223–1230.
48. T. T. Bui, O. Thiebaut, E. Grelet, M. F. Achard, B. Garreau-deBonneval and K. I. M. C. Ching, Eur. J. Inorg. Chem., 2011, 2663–2676,  DOI:10.1002/ejic.201001288.
49. X. Q. Wang, Q. Ren, J. W. Chen, W. T. Yu, H. L. Fang, T. B. Li, H. J. Cong, X. T. Liu, L. Y. Zhu, G. H. Zhang and D. Xu, Solid State Sci., 2011, 13, 896–903.
50. S. Dalgleish, J. G. Labram, Z. Li, J. P. Wang, C. R. McNeill, T. D. Anthopoulos, N. C. Greenham and N. Robertson, J. Mater. Chem., 2011, 21, 15422–15430.
51. P. Basu, A. Nigam, B. Mogesa, S. Denti and V. N. Nemykin, Inorg. Chim. Acta, 2010, 363, 2857–2864.
52. X. Q. Wang, Q. Ren, H. L. Fan, J. W. Chen, Z. H. Sun, T. B. Li, X. T. Liu, G. H. Zhang, D. Xu and W. L. Liu, J. Cryst. Growth, 2010, 312, 2206–2214.
53. P. Deplano, L. Pilia, D. Espa, M. L. Mercuri and A. Serpe, Coord. Chem. Rev., 2010, 254, 1434–1447.
54. B. Garreau-deBonneval, K. I. M. C. Ching, F. Alary, T. T. Bui and L. Valade, Coord. Chem. Rev., 2010, 254, 1457–1467.
55. M. Iyoda and H. Shimizu, Chem. Soc. Rev., 2015, 44, 6411–6424.
56. W. Q. Li, Q. Peng, H. L. Ma, J. Wen, J. Ma, L. A. Peteanu and Z. G. Shuai, Chem. Mater., 2017, 29, 2513–2520.
57. J. D. C. Gonzalez, M. Iyoda and J. P. Rabe, Nat. Commun., 2017, 8.
58. B. Holzer, J. Bintinger, D. Lumpi, C. Choi, Y. Kim, B. Stoger, C. Hametner, M. Marchetti-Deschmann, F. Plasser, E. Horkel, I. Kymissis and J. Frohlich, ChemPhysChem, 2017, 18, 549–563.
59. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
60. B. Kirtman, V. Lacivita, R. Dovesi and H. Reis, J. Chem. Phys., 2011, 135, 154101.
61. T. Miletic, A. Fermi, I. Orfanos, A. Avramopoulos, F. De Leo, N. Demitri, G. Bergamini, P. Ceroni, M. G. Papadopoulos, S. Couris and D. Bonifazi, Chem. – Eur. J., 2017, 23, 2363–2378.
62. Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai and K. Hirao, J. Chem. Phys., 2004, 120, 8425–8433.
63. S. Bonness, H. Fukui, K. Yoneda, R. Kishi, B. Champagne, E. Botek and M. Nakano, Chem. Phys. Lett., 2010, 493, 195–199.
64. R. Kishi, S. Bonness, K. Yoneda, H. Takahashi, M. Nakano, E. Botek, B. Champagne, T. Kubo, K. Kamada, K. Ohta and T. Tsuneda, J. Chem. Phys., 2010, 132, 094107.
65. O. Loboda, R. Zalesny, A. Avramopoulos, J. M. Luis, B. Kirtman, N. Tagmatarchis, H. Reis and M. G. Papadopoulos, J. Phys. Chem. A, 2009, 113, 1159–1170.
66. J. D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128.
67. C. G. Liu, M. L. Gao and Z. J. Wu, RSC Adv., 2014, 4, 38300–38309.
68. X. H. Xie, D. S. Hao and J. C. Yang, Chem. Phys., 2015, 461, 11–19.
69. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
70. S. I. Lu, Theor. Chem. Acc., 2015, 134, 1589.
71. D. Jacquemin, E. A. Perpete, M. Medved, G. Scalmani, M. J. Frisch, R. Kobayashi and C. Adamo, J. Chem. Phys., 2007, 126, 144105.
72. A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377.
73. C. T. Lee, W. T. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
74. L. Serrano-Andres, A. Avramopoulos, J. B. Li, P. Labeguerie, D. Begue, V. Kello and M. G. Papadopoulos, J. Chem. Phys., 2009, 131, 134312.
75. H. Fukui, R. Kishi, T. Minami, H. Nagai, H. Takahashi, T. Kubo, K. Kamada, K. Ohta, B. Champagne, E. Botek and M. Nakano, J. Phys. Chem. A, 2008, 112, 8423–8429.
76. S. Grimme, J. Chem. Phys., 2006, 124, 034108.
77. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci and G. A. Petersson, et al., Gaussian 09, Revision A.02, Gaussian, Inc., Wallingford, CT, Gaussian 09, Revision D.1, Gaussian, Inc., Wallingford, CT, 2009.
78. H. D. Cohen and C. C. Roothaan, J. Chem. Phys., 1965, 43, S034–S039.
79. P. J. Davis and P. Rabinowitz, Numerical Intergration, Blaisdell, London, 1967.
80. W. Romberg, Kgl. Norske Vid. Selsk. Forsk, 1955, 28, 30–36.
81. H. Rutishauser, Num. Math., 1963, 5, 48–54.
82. L. Salem and C. Rowland, Angew. Chem., Int. Ed., 1972, 11, 92.
83. J. Wirz, Pure Appl. Chem., 1984, 56, 1289–1300.
84. V. Bachler, G. Olbrich, F. Neese and K. Wieghardt, Inorg. Chem., 2002, 41, 4179–4193.
85. L. Serrano-Andres, D. J. Klein, P. Schleyer and J. M. Oliva, J. Chem. Theory Comput., 2008, 4, 1338–1347.
86. M. Abe, Chem. Rev., 2013, 113, 7011–7088.
87. P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, 213–222.
88. M. M. Francl, W. J. Pietro, W. J. Hehre, M. S. Gordon, D. J. DeFrees and J. A. Pople, J. Chem. Phys., 1982, 77, 3654–3665.
89. M. Dolg, U. Wedig, H. Stoll and H. Preuss, J. Chem. Phys., 1987, 86, 866–872.
90. T. Kubo, A. Shimizu, M. Uruichi, K. Yakushi, M. Nakano, D. Shiomi, K. Sato, T. Takui, Y. Morita and K. Nakasuji, Org. Lett., 2007, 9, 81–84.
91. M. E. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem., 2012, 63, 287–323.
92. A. Dreuw and M. Head-Gordon, Chem. Rev., 2005, 105, 4009–4037.
93. D. Jacquemin, E. A. Perpete, G. E. Scuseria, I. Ciofini and C. Adamo, J. Chem. Theory Comput., 2008, 4, 123–135.
94. D. Jacquemin, A. Planchat, C. Adamo and B. Mennucci, J. Chem. Theory Comput., 2012, 8, 2359–2372.
95. D. Jacquemin, V. Wathelet, E. A. Perpete and C. Adamo, J. Chem. Theory Comput., 2009, 5, 2420–2435.
96. C. A. Guido, S. Knecht, J. Kongsted and B. Mennucci, J. Chem. Theory Comput., 2013, 9, 2209–2220.
97. D. Geldof, A. Krishtal, F. Blockhuys and C. Van Alsenoy, J. Chem. Theory Comput., 2011, 7, 1328–1335.
98. P. Bultinck, C. Van Alsenoy, P. W. Ayers and R. Carbo-Dorca, J. Chem. Phys., 2014, 140, 144104.
99. N. Otero, C. Van Alsenoy, C. Pouchan and P. Karamanis, J. Comput. Chem., 2015, 36, 1831–1843.
100. N. Otero, C. Van Alsenoy, P. Karamanis and C. Pouchan, Int. J. Comput. Theor. Chem., 2013, 1021, 114–123.
101. P. Karamanis, N. Otero, C. Pouchan, J. J. Torres, W. Tiznado, A. Avramopoulos and M. G. Papadopoulos, J. Comput. Chem., 2014, 35, 829–838.
102. P. Karamanis, N. Otero and C. Pouchan, J. Phys. Chem. C, 2015, 119, 11872–11885.
103. C. van Alsenoy and A. Peeters, J. Mol. Struct. THEOCHEM, 1993, 286, 19–34.
104. B. Rousseau, A. Peeters and C. Van Alsenoy, Chem. Phys. Lett., 2000, 324, 189–194.
105. D. M. Bishop and B. Kirtman, J. Chem. Phys., 1991, 95, 2646–2658.
106. D. M. Bishop and B. Kirtman, J. Chem. Phys., 1992, 97, 5255–5256.
107. D. M. Bishop, J. M. Luis and B. Kirtman, J. Chem. Phys., 1998, 108, 10013–10017.
108. D. M. Bishop, M. Hasan and B. Kirtman, J. Chem. Phys., 1995, 103, 4157–4159.
109. 1 a.u. of energy ≈ 27.21070 eV, ≈ 627.51 kcal mol⁻¹, 1 a.u. of electric field ≈ 1.71514 × 10⁷ esu, ≈ 5.14192 × 10¹¹ V m⁻¹, 1 a.u. dipole moment ≈ 2.54174 D, ≈ 8.47831 × 10⁻³⁰ C m, 1 a.u. of polarizability ≈ 0.148176 × 10⁻²⁴ esu, ≈ 0.164867 × 10⁻⁴⁰ C² m² J⁻¹, 1 a.u. of first hyperpolarizability ≈ 0.863993 × 10⁻³² esu, ≈ 0.320662 × 10⁻⁵² C³ m³ J⁻², 1 a.u. of second hyperpolarizability ≈ 0.503717 × 10⁻³⁹ esu, ≈ 0.623597 × 10⁻⁶⁴ C⁴ m⁴ J⁻³.
110. Y. W. Son, M. L. Cohen and S. G. Louie, Nature, 2006, 444, 347–349.
111. E. Rudberg, P. Salek and Y. Luo, Nano Lett., 2007, 7, 2211–2213.
112. O. Hod and G. E. Scuseria, ACS Nano, 2008, 2, 2243–2249.
113. J. L. Oudar and D. S. Chemla, J. Chem. Phys., 1977, 66, 2664–2668.
114. Y. Su, X. Wang, L. Wang, Z. Zhang, X. Wang, Y. Song and P. P. Power, Chem. Sci., 2016, 7, 6514–6518.
115. A. Avramopoulos, H. Reis, G. A. Mousdis and M. G. Papadopoulos, Eur. J. Inorg. Chem., 2013, 4839–4850.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7tc05047j

---