Bernardino Tirri

Gloria Mazzone

Alistar Ottochian

Jerôme Gomar

Umberto Raucci
email: umberto.raucci@unina.it

Carlo Adamo

Ilaria Ciofini
email: ilaria.ciofini@chimieparistech.psl.eu

A combined Monte Carlo/DFT approach to simulate UV-vis spectra of molecules and aggregates: merocyanine dyes as a case study

The combination of a Monte Carlo (MC) sampling of the configurational space with Time dependent Density Functional Theory (TD-DFT) to estimate vertical excitations energies has been applied to compute the absorption spectra of a family of merocyanine dyes in both their monomeric and dimeric forms. These results have been compared to those obtained using a static DFT/TD-DFT approach as well as to the available experimental spectra.

Though suffering of the limitations related to the use of DFT and TD-DFT for this type of systems, our data clearly show that the classical MC sampling provides a suitable alternative to classical molecular dynamics to explore the structural flexibility of these donor-acceptor (D-π-A) chromophores enabling a realistic description of the Potential Energy Surface of both their monomers and aggregates (here dimers) and thus of their spectra.

Overall, the combination of MC sampling with quantum mechanics (TD-DFT) calculations, carried out in implicit dioxane solvent on random snapshots, provides a workable compromise to solve the combined challenge of accuracy and time-consuming problem not only for merocyanines momers but also for their dimers, up to now less investigated. Indeed, the simulated absorption spectra fairly agree with the experimental ones, suggesting the general reliability of the method.

Introduction

Merocyanines are a class of polar dyes showing a variety of colors, with exceptional non linear optical responses, photo-dimerization, photo-isomerization, solvatochromism, and thermochromism properties, and unique self-assembly behavior arising from their large dipolar character, leading to have been employed extensively as spectral sensitizers for silver halide photography, as hybrid large band gap semi-conductor materials, in optical disks as recording media, in industrial paints, for trapping of solar energy, as laser materials, as photorefractive materials, as antitumor agents, as probes for biological systems, as sun filters. 1 This class of molecules is characterized by the presence of a donor and an acceptor group connected by an ethylene or polyethylene bridge (see Figure 1 for a schematic representation of these molecules). Thus, the so-called Bond-Length Alternation (BLA), accounting for the difference between single and double carbon-carbon bonds of the polyethylene bridge, represent a relevant structural parameter of merocyanines. Small (large) BLAs are an indication of an increased (reduced) delocalization of the  electrons and of an improved (worsen) coupling between the donor-acceptor moieties, thus leading to a red (blue) shift of the spectra. In terms of the electronic structure, a small BLA indicates a formally neutral structure, while a large BLA suggests an increasing weight of the charge-separation (zwitterionic) form. 2 Nonetheless, BLA is by itself a challenging property for computational approaches such as Density Functional Theory (DFT), since its correct reproduction requires the reduction (or even cancellation) of the Self-Interaction Error (SIE), which, inducing a large electronic delocalization, leads to shorter BLAs. 3,4 In the last decade, several efforts have been devoted to the theoretical description of the merocyanine family of dyes using DFT and TD-DFT, a class of numerical methods that allows for a reliable prediction of electronic absorption spectra not only in term of band position (i.e. max) but also in intensity and band shape in simple and complex media. [5][6][7] Several studies focus on the definition of the most suitable exchange-correlation functional for modelling the merocyanines electronic absorption spectra. Since the most intense transition in merocyanines is expected to show a CT character standard DFT approaches are expected to fail. 8 However, previous studies have shown that functionals normally curing (or mitigating) CT problems, including range separated hybrids and optimally-tuned functionals, do not significantly improve the results for merocyanines. 5,9,10 Actually, it has even been suggested that the insensibility of the computed transitions to the choice of the DFT approach is a fingerprint of the merocyanine. 11 Furthermore, the fine-structure of their absorption spectra is still matter of debate. Indeed, while from one side, theoretical model studies and fine analysis of experiment outcomes suggest that the broadening of the first absorption band is mainly related to the C=C vibration of the polymethine chain [11][12] , from another side, studies using complete modelling of the band shape of the absorption spectrum using a Franck-Condon (FC) approach indicate that the origin of this structured band can be ascribed to a number of different vibrational modes. 10,[13][14][15][16] In some cases, difficulties have been evidenced, so that approximate FC models are used. For instance, vibrational levels are considered as being equal for both the ground and excited states (adiabatic approximation) or only one degree of freedom is explicitly considered. 15,16 As rule of thumb, the vibronic coupling is expected to shift the max of about 0.3 eV. 15 Furthermore, simulating the merocyanines behavior in complex media is even more complicated considering that: i) merocyanines undergo to strong solvatochromism 17 , and the solvent effects cannot be in turn avoided in the modelling; 16,[18][19] ii) in solution merocyanines form H-and J-type aggregates (dimers), differing in the relative orientation of the dipole moment of the two constituting monomers (antiparallel for H-aggregates, parallel for the Jtype) [20][21] . All these issues are strongly entwined and, often, they have been investigated separately.

Here, focusing on both the monomeric and dimeric forms of the merocyanines 1-4 reported in Figure 1, we tried to disentangle the intricate relationships between the structural and electronic features of merocyanines exploring their complex phase space with a stochastic approach relying on Monte Carlo (MC) molecular simulations. The flexible structure of merocyanines adds another level of complexity in the modeling. Indeed, standard methods to obtain the vibrationally resolved optical spectra [22][23][24][25] cannot be applied in this situation due to the presence of several energetically accessible minima existing on the ground and/or excited state potential energy surfaces (PESs). [26][27][28] In such cases a single structure can hardly be representative of the real phase space and ensemble averages become mandatory. For this reason, we simulated the absorption spectrum as an average of those vertically computed for a statistically representative average over an ensemble of molecular structures, extracted from the ground state PES. In this respect, the widely applied strategy to build the statistical ensemble is based on Molecular Dynamics (MD) simulations, where the time evolution of the system is achieved by integrating the classical Newton's equations of motion, the corresponding potential energy being computed either at classical or quantum chemistry level of theory. [26][27][28][29][30] Although MD simulations are widely and successfully applied, some concerns still arise. Indeed, when large energy barriers (often of several kBT) to conformational changes (e.g. bond rotations in flexible molecules, geometrical rearrangements in aggregates) exist, it is easy to be trapped in local -low energyconformations, leading to a poor quality statistical sampling. Increasing the temperature of the simulation, in conjunction or not with simulated annealing approaches or replicate trajectories starting from different initial conditions are common procedures to improve the space-phase sampling. 31 Driven dynamics (metadynamics, umbrella sampling) are also viable alternatives, requiring the identification of the relevant reaction coordinate(s). 32,33 Nevertheless, spectra simulated at high temperature are clearly broadened with respect to the experimental ones, and more replicas of the trajectory or driven dynamics can be computationally difficult to achieve especially in the case of ab initio molecular dynamics.

For these reasons, we relied on an alternative philosophy to escape from local minima on the PES based on random sampling with the Monte Carlo method. In this approach, configurations in the phase space are generated by randomly changing the coordinates of a randomly chosen atom of the system, and accepting the new configuration according to the Metropolis sampling algorithm. [START_REF] Allen | Computer Simulation of Liquids[END_REF] The random moves and the possibility to accept higher energy configurations allow MC simulations to avoid being trapped on local minima reaching energetically distinct regions of the PES, which allows the simulation of the photo-physical properties of complex chromophores in realistic environments such as condensed phases (solution, crystal or amorphous phases), or in matrix (biological or polymeric).

Although this approach is very powerful, few examples are reported in literature coupling MC sampling with quantum chemistry calculation for the simulation of UV-Vis absorption spectra of molecules, mainly focused on the organization of the solvent molecules around the chromophore. [START_REF] Georg | [END_REF][36][37] The main aim of the present work is thus to show the potentiality of a combined Monte Carlo/Quantum Chemistry, here Time Dependent Density Functional Theory (TD-DFT), approach to simulate the electronic absorption spectra of the merocyanines family of flexible molecules and their aggregates (here dimers).

Computational details

Two computational protocols were considered for the calculation of UV-Vis spectra. The first one, referred to as static, consists of a full optimization of the molecular structures using the PBE0 functional 38 and the 6-31+G(d,p) basis set. 39 In the case of dimers, dispersive interactions were accounted through the original Grimme empirical correction (GD3). 40 Vertical absorptions were then computed at the TD-DFT level, using global (PBE0) and range separated (B97XD) 41 hybrid functionals in conjunction with the 6-31+G(d,p) basis set. The first functional, PBE0, is widely used for the computation of UV-Vis spectra involving valence transitions, 42 while the second, containing a larger quantity of HF exchange, is more suitable for Charge Transfer (CT) excitations and it was already applied to study model merocyanines. 5 Solvent effects were taken into account using the IEF-PCM model, considering dioxane as solvent as in the experimental protocol. 47 Due to the aprotic nature of this solvent, no specific solute-solvent interactions are expected, so that the IEF-PCM provides accurate results for the electronic transitions. 43 The second approach, hereafter referred to as stochastic, is based on the generation of ensemble configurations for both monomer and dimeric forms of merocyanines 1-4 (Figure 1) using MC Metropolis simulations. Our method is schematized in Figure 2 and works as it follows: i) the statistical ensemble is built through the MC sampling of the ground state PES, ii) vertical excitations on randomly selected configurations are carried out at quantum mechanics level, iii) the electronic absorption spectrum is obtained convoluting by a gaussian broadening the so computed vertical excitations.

The energy of the MC configurations, was computed at molecular mechanics level using standard Generalized Amber Force Field (GAFF). 44 Production MC simulations were carried out for 1*10 7 steps at a 300 K temperature. Then, 500 structures were randomly extracted to simulate the UV-Vis spectrum and their vertical electronic transitions were computed using the same TD-DFT computational level as for the static approach, including solvent effects.

Convergence of the sampling was tested for the monomer 1 molecule by using a sampling of 600 providing the same computed spectral properties.

Finally, the resulting absorption spectrum (A(E)) was simulated as:

𝐴(𝐸) = 1 N p ∑ [∑ 𝑓 IL (𝑹 𝒌 )𝑔(𝐸 -Δ𝐸 IL (𝑹 𝒌 ), 𝜹) N S L≠I ] N P 𝑘
where vertical excitation energies (ΔEIL(Rk)) between the ground electronic state (I) and excited state (L) and corresponding oscillator strengths (fIL(Rk)) were computed for each sampled configuration (Rk), Ns is the number of excited states considered, NP is the number of sampled points, and g represents a normalized gaussian function centered on ΔEIL of full width at half maximum () of 0.3 eV.

All DFT and TD-DFT calculations were carried out using the Gaussian program, 45 while MC computations were performed using the LAMMPS code. 46

Results and Discussion

UV-Visible spectra: static vs stochastic approach.

In Table 1 are collected the computed absorption energies maxima (max) for monomers and dimers of the merocyanines sketched in Figure 1, obtained with the two approaches described above (static and stochastic).

From literature data it's known that merocyanine dyes undergo vibronic coupling significantly affecting the position of the max computed by a static approach . The red-shift of the main peak, with respect to the vertical approximation, was estimated to be about 0.3 eV. 15 The comparison of data coming from static calculations with the experimental wavelengths, 47 included in Table 1, not surprisingly reveals that the vertical transitions energies computed at the PBE0 level for all the monomers are significantly overestimated (between 0.42 and 0.50 eV). The computed Mean Absolute Deviations (MADs) are quite large, but essentially in agreement with the errors reported in the literature for similar systems. 5,9,10 Furthermore, if the estimated 0-0 correction is subtracted (0.3 eV), the deviations are within the expected range of accuracy for TD-DFT methods (0.1-0.2 eV), always taking into account the above reported caveats. These discrepancies slightly increase if the range-separated B97XD functional is considered (+0.47 eV for 1), in agreement with data in literature. 5 For these reasons, B97XD will be not further considered for the other molecules and only PBE0 results will be discussed.

Despite these significant deviations between theory and experiments, the experimental relative trends are correctly reproduced by the PBE0 results. Indeed, the outcomes obtained (Table1)

show significant blue-shifts going from molecule 1 to 4 (between 0.09 and 0.56 eV), in the same range of the experimental observation obtained by Wortmann and co-workers (between 0.01 and 0.56 eV), 47 with the exception of molecule 2, for which a red-shift of 0.22 eV is found, in agreement with the experimental trend.

The topology of the Frontier Molecular Orbitals (FMOs) of the investigated molecules is depicted in Figure 3. According to the TD-DFT calculations, the most intense absorption corresponds to a 𝜋 → 𝜋 * transition, from the highest occupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO). These FMOs are delocalized over the whole core of the molecules, including the polymethine unit and the pyridinic residues (or the amino group for molecules 3). Therefore, any structural deformation modifying the electronic conjugation of the backbone impacts the computed absorption max. On the other hand the CT character appears limited.

The trend found for monomers are kept for the dimers. All the static calculations reported in 47 Comparing the dimers absorption spectra, the one of dimer 2 is characterized by a red-shift (0.30 eV) of the maximum absorption, with respect to 1, the same found in the experiments. 47 In contrast, the molecule 3b is slightly blue shifted at theoretical level (0.03 eV), while it results red-shifted (0.10 eV) with respect to the recorded spectrum. 47 Resuming, static TD-DFT calculations on the monomeric and H-type dimeric forms of merocyanines provide a picture in qualitative agreement with the experimental data, even if the obtained errors are quite large, though in line with previous literature data.

A sizable improvement, compared to the static approach, is obtained simulating the electronic absorption spectra as an ensemble average over multiple configurations. As matter of fact, the MAD is halved going from the static to the stochastic approach for both monomers and dimers.

In details, the MAD is 0.46 eV for the static monomers whereas it becomes 0.26 eV using a MC sampling. For a neat comparison the vertical transitions are computed at the same PBE0/6-31+G(d,p) level of theory in both the static and the MC spectra. As in the case of the static approach, the largest error (+0. The MC sampling leads to even better results in the case of dimers, for which the MAD decreases from 0.38 eV (static approach) to 0.17 eV (stochastic approach). In this case, the capability of MC to account for several configurations accessible is crucial to reduce the simulated/experimental error on max. Compound 3b represents the trickiest molecule to simulate (error of +0.22 eV), also in its dimeric form. Overall the MC results are in fairly agreement with the experimental spectra both in terms of absolute computed-experimental shifts than in relative shifts between the different merocyanines.

Structural effects underpinning spectra

After the discussion of the static and stochastic TD-DFT results on the band position, it is interesting to examine how the structural features of merocyanines and their dimeric species influence the max and how the two different approaches are able to describe them. As abovestated, the major limitation of the static approach essentially relies on the impossibility to consider all the energetically accessible minima existing on both ground and/or excited PESs.

The MC simulations, though cannot be able to provide the ZPE corrections due to their quantum nature, allows overcoming such a limitation taking into account several energetically accessible conformations, thus supplying the contribution to max coming from the distortions of flexible geometrical parameters.

As matter of facts, several degrees of freedom can drastically influence the electronic properties of this family of chromophores. Relevant parameters, to be taken into consideration, are: i) the bond length alternation (Figure 4a), ii) the relative orientation between the pyridone and the pyridine rings (amino group for 3), linked by the conjugated chain (defined by the angle , Figure 4b) and iii) for the dimer, the relative distance (and orientation) of the two monomers.

Small (large) BLAs are indication of increased (reduced) delocalization of the  electrons and of an improved (worsen) coupling between the two moieties, thus leading to a red (blue) shift of the spectra. A similar effect on the electronic structure could be expected by a coplanar rearrangement of the pyridone and the pyridine rings (or amino moiety for 3), where smaller  values correspond to smaller transition energies. Short intermolecular distances, represented as the distance between the centre of mass (CoM) of the two monomers, lead to low energy transitions while an opposite effect is expected for larger inter-dimer distances.

The main geometrical parameters computed for merocyanines 1-4 are collected in Table 2 whereas the structures of the monomers are depicted in Going from static to MC calculations, small variations of the angle  are observed, the largest ones being of about ± 4 degrees (for 2, 3a and 3b), thus suggesting that this parameter, at least for monomers, is not relevant.

Concerning BLAs, the small lengths of the polymethine chain should prevent large errors at DFT level, generally observed in longer oligomeric chains 48,49 and, indeed, PBE0 and ωB97XD functionals give very close values for molecule 1, even if their errors become significantly different when larger systems are considered. Larger variations are, instead, observed in going from static to MC calculations, with different behavior for the considered systems. In particular, the BLAs increase for 2, 3b and 4 and decrease for 1 and 3a, although in some cases the variations are quite small. Overall, a small energy difference (about 3 kcal/mol) between the PBE0 optimized and the MC average structures is found for all the molecules. This energy difference is not only related to the BLA since the two structures differ also in other geometrical parameters (i.e., lateral chains conformation), as can be inferred from the structures reported in Figure 4c-f. More interestingly, it cannot be traced from our results a direct relationship between any single selected structural parameter, including the BLAs or  angles, and the red-shifts observed in going from static to MC values. Such a shift results from a global deformation rather than a single dominant degree of freedom. In other words, the hypothesis of a dominant effect of the C=C double bond variation on the merocyanine spectra cannot be confirmed by our calculations. Furthermore, it should be recalled that structures obtained via the stochastic approach are produced at classical molecular mechanics level thus the observed effect on the spectra is the sum of small structural effect due to the difference between the full DFT and the MM based MC approach and the dominant effect related to the need of a more exhaustive sampling allowed by the MC approach. Of course favorable error cancellation in this procedure cannot be excluded.

Therefore, the net effect of stochastic calculations is the red-shift of the max providing better agreement with the experimental spectra. It should be also pointed out that the resulting MAD is of the same order of magnitude as the E 0-0 estimated for similar systems (0.3 eV). 15 In an attempt to provide a comparison even with the vibrationally resolved absorption spectrum, the structure of the lowest electronic excited state (S1) of molecule 1 was optimized at the PBE0/6-31G level (see Supporting Information Figure SI.1 for their representation). However, because of the flexibility of the molecule, the obtained S1 minimum energy structure is strongly distorted with respect to the ground state one ( close to 90°). As a consequence, the FC approximation is no longer valid and it is not possible, at least with usual approaches, to compute the vibronic contribution.

Once defined the global molecular flexibility as the main factor affecting the absorption energies of the merocyanine monomers, it is interesting to establish how these structural features can impact the aggregate formation.

All the molecules considered are push-pull systems possessing a very high dipole moment on the ground electronic state (20.3 Debye for 1). Previous studies have shown that dipole-dipole interactions are indeed the driving force for their self-assembly in dilute solution leading to dimers showing an antiparallel orientation of the dipole moment of the constituting monomers and thus an overall null dipole moment. Furthermore, the presence of a null dipole moment of dimers avoids their further clustering and thus the formation of larger aggregates. 21 In agreement with the experimental observation, the MC sampling on the merocyanine dimers indicates that most of the stable dimeric forms are cofacial and antiparallel aligned (H aggregates) while only a few of them have a parallel arrangement of the two units (J dimer). In particular static calculations, performed at PBE0/6-61++G(2d,2p) level of theory, indicate that the H aggregates are most stable than the J ones of 18.1 kcal/mol and this is confirmed by the MC calculations showing less than 1% of parallel dimers. Furthermore, due to the relative high difference in energy between the two forms it is clear that it is difficult to go from one to the other arrangement.

Due to the statistical sampling which massively favor the H dimers, in agreement with the experimentally observed blue-shift, only these latter will be considered in the following. These dimers appear even more flexible than the isolated monomers, since they undergo large structural modifications, as illustrated by the data reported in Table 2 andFigure The intermolecular distance was computed as the distance between the centers of mass (CoM) of the two moieties. Shorter distances correspond to a more compact structure, whereas the longer distances arise from a relative shift of the two units to reduce the steric hindrance. Thus, the intermolecular distance strongly depends on the steric repulsion in the packed structures, so that a large panel of values is found, ranging from 3.9 (molecule 3b, the most compact) to 5.6 Å (molecule 2). It is interesting to note that the MC structures are relatively higher in energy compared to the optimized ones. The difference in total energies between the two sets is just twice than the value observed for the monomer (7-9 kcal/mol, see Table 2), the exception being the dimer 2 where very compact rearrangement (CoM = 4.3 Å) lead to significantly higher energies (E = 14.2 kcal/mol).

Going from the monomers to the dimers, all the transitions are blue-shifted, so it seems that the mutual polarization of the two monomers in the aggregate plays the major role. Therefore, we have simulated the electronic spectra of the monomer at the geometry in the corresponding dimer, for both optimized static and MC structures. As it appears from the results collected in Table 2, all the transitions are close (within a few thousandths of eV) to those computed for the corresponding isolated monomers, regardless of the approach used to obtain such a structure (static or MC). As for the isolated monomer, it is not possible to draw a direct relationship between the BLA and the  parameters, the main factor affecting the blue-shifts of the dimers seems to be their formation, accounted by the distance of the CoMs. Nevertheless, there is not a correlation between the blue-shifts in the transitions and the energy of formation of the dimers.

Similarly, no correlations were found between the static and MC transitions for the dimers and their corresponding CoMs distances (or their differences). However, the numerical improvement in computed absorption energies appears evident looking at the computed MAD (0.17 eV for the MC dimers), well below the expected value for such systems. Figure 6, where the computed spectra for the dimer 1 are reported, also exemplifies this improved agreement

Conclusions

In conclusion, the overall agreement between MC results and experimental data suggests that the main structural effects on the electronic spectra are captured by the MC approach combined with vertical TD-DFT calculations. Although classical MC simulations do not allow to reconstruct the vibronic progression, they take into account the flexible nature of the molecules and its effect on the max which is overall shifted leading to a better agreement with the experimental spectra. Our results reveal that despite its simplicity a combined MC/TD-DFT approach is reliable to describe flexible systems that cannot be properly described in a FC static framework. Indeed, in the case of the merocyanines dyes the characterization of the FC minimum leads to an error bar (0.3-0.5 eV) for the spectra of the monomers obtained with standard functionals, including global or range-separated hybrids, in perfect agreement with literature data. 22,50 Finally, to the best of our knowledge, a detailed FC analysis of merocyanines dimers has never been performed, the only study concerning two units bound by aliphatic chains which significantly damp intermolecular soft modes. This is not surprising, in view of the difficulties that can be expected in the optimization of the dimers on the excited electronic states. In this respect, the MC stochastic sampling here proposed is a reliable alternative to sample the dimer conformational space. Promising applications can be expected in the future, where the MC approach could be used to generate initial conditions for ab initio molecular dynamics.

  3 eV) is obtained for the compound 3b, suggesting the subtle effects played by the length of the polymethine bridge and the BLA. The relative shifts computed for the different merocyanines in their monomeric forms also fairly agree with the static approach and the experimental results. The character of the brightest electronic transition corresponds to a 𝜋 → 𝜋 * transition, where HOMO and LUMO orbitals are delocalized over the whole core of the molecule (FMOs of a representative geometry are reported in Supporting Information see Figure SI.1 and SI.2).

Figure 4 .

 4 In the case of MC simulations average values are reported (see Supporting Information Figure SI.2 for their representation).

  5. Indeed, significant variations are observed for both  and BLAs. Starting from the first parameter, the static structures of dimers 2, 3a and 3b are characterized by large torsional angles ( around 50 degrees), while dimers 1 and 4 still retain an almost planar arrangement of their moieties (around 10 degrees). As above-mentioned, their corresponding isolated monomers are all almost planar. The weak non-covalent interactions between the pyridine rings force the planar structure of the monomers in molecules 1 and 4, while the flexible amino terminal chain in dimers 3 prevents such arrangement. In molecule 2, instead, the head-tail rearrangement is characterized by a steric repulsion between the ethyl and butyl chains substituting group of the two aromatic rings, which induces a shearing of the two moieties. As a consequence, the leucoline moieties can rotate around a C-C bond of the polymethine chain. These distorted structures break the extended  conjugation and induce a significant increase of the BLAs in the dimers exhibiting the largest values of  (for molecules 1, 2, 4). Therefore, the zwitterionic form is favored over the neutral one.

  a) computed as the average value of the two monomers at their static geometry in the dimer; b) computed as the average value of the two monomers at their MC geometry in the dimer; c) from reference 47; d) Mean Absolute Deviation

  0 a) refer to Figure 4b. b) Refer to Figure 4a , c) refer to text for further explanation.

Figure captions Figure 1 .Figure 2 .

 captions12 Figure captions Figure 1. Sketch of merocyanine dyes considered in the present paper. Bu=butyl, Hex=hexyl Figure 2. Schematic illustration of the stochastic approach.

Figure 3 .

 3 Figure 3. Illustration of the HOMO/LUMO computed for merocyanines monomers. Blue/ yellow regions indicate the negative/positive isosurfaces.

Figure 4 .Figure 5 .Figure 6 .Figure 6 .

 4566 Figure 4. Optimized static (left) and Monte Carlo averaged structures (right) for molecule 1 (a and b) and 4 (c and d) Figure 5. Optimized static (left) and Monte Carlo averaged structures (right) for dimers ofmolecules 1 (a and b), 2 (c and d) and 3b (e and f) Figure 6. Absorption spectra of dimer 1 in 1,4 dioxane computed with the static (blue dashed line) and stochastic (green solid line) approaches. The experimental spectrum is also reported for comparison (red line)

Table 1

 1 refer to co-facial antiparallel dimers (H-aggregate, see next paragraph for further explanation on structures), representing the global minimum on the ground state PES. A significant overestimation of the max is found at the PBE0 level, even if the deviations, falling within 0.32 eV (1 and 2) and 0.45 eV (3b), are slightly smaller than the corresponding monomers. Indeed, the final MAD for dimeric species is 0.38 eV, to be compared with 0.46 eV for the monomeric ones. As experimentally observed, all vertical transitions are blue-shifted in going from the monomer to the dimer, although this variation is slightly overestimated with respect to experiments.

Table 1 .

 1 Computed max (in eV) considering the dioxane as solvent and using a static and a MC approach All the electronic transitions were computed at the PBE0/IEF-PCM level of theory, except for the values in parenthesis, which were evaluated using the B97XD/IEF-PCM model. Mean Absolute Deviation (MAD in eV) with respect to the experimental values is also reported.

	Monomer

Table 2 .

 2 Computed main geometrical parameters corresponding to the static approach or averaged MC simulations. The energy differences (E) between optimized PBE0 and MC structures are also reported. Values in parenthesis refer to structures optimized with the B97XD functionals.

	Monomer

Acknowledgements I.C. gratefully acknowledge support from the European Research Council (ERC) for funding under the European Union's Horizon 2020 research and innovation programme (grant agreement No 648558, STRIGES CoG grant).

Data availability statement. The authors confirm that the data supporting the findings of this study are available within the article and its supplementary materials. Raw data are available from the corresponding author, upon request.


Davide Luise

Liam Wilbraham

Frédéric Labat
email: frederic.labat@chimie-paristech.fr

Ilaria Ciofini
email: ilaria.ciofini@chimie-paristech.fr

Modeling UV-Vis spectra of low dimensional materials using electrostatic embedding: the case of CdSe

come

Modeling UV-Vis spectra of low dimensional materials using electrostatic embedding: the case of CdSe

Introduction

A large effort has recently been devoted to devise new photovoltaic technologies which could potentially outperform silicon-based solar cells in terms of efficiency and fabrication costs [1]. In particular, photosensitization of wide band gap n-type semiconductors (SC) such as TiO2 has paved the way to the development of various sensitized solar cell (SSC) technologies. Different classes of sensitizers have been tested, such as organic or inorganic dye molecules [2;3], quantum-dots (QD) [4] or perovskites [5;6]. The resulting photo-electrochemical devices usually use a monolayer of photosensitive material adsorbed on a thin layer of nanocrystalline oxide to harvest light, while a photoinduced electron transfer successively takes place from the sensitizer to the SC. In the case of QD-sensitized solar cells (QDSC), inorganic II-VI semiconductor nanoparticles (such as CdSe, CdS, CdTe, PbS and PbSe) have been used as sensitizers as their size-dependent optical properties can be leveraged to target specific applications. The nanoscopic nature of these particles induces a quantum confinement of the exciton wave function, increasing the nanoparticle's bandgap with decreasing size [7][8]. Furthermore, the nanocrystallites can generate multiple electron-hole pairs upon absorption of a single photon, which, in principle, allows one to exceed the Shockley-Queisser limit. Among the possible inorganic II-VI semiconductors for QDSC application, CdSe has been a material of choice.

CdSe nanocrystallites can be readily obtained via colloidal synthesis both in aqueous and organic media [8;9], with various shapes and dimensions depending on the surfactant molecules used during synthesis, leading especially to zero dimensional (0D) quantum dots (QD) [7][8], 1D nanorods and 2D nanoplatelets (NPL) [10], with 1D and 2D objects having dimensions outside of the nanoscale range. In particular, QD are 3D confined spherical nanoparticles with typical diameters in the range of 2-10 nm, whereas CdSe NPL are 2D extended nanosheets confined due to their finite thickness along the normal to a well-controlled exposed plane at the experimental level [10]. Although both systems therefore present photophysical properties which can be tuned according to their dimensions and to the application targeted, CdSe NPL present some advantages over QD, including sharper photoluminescence, higher quantum yield and larger absorption cross section. [11;12;13] They can also be easily synthesized by atomic deposition methods from monodisperse colloidal solutions of nanosheets, with a finely controlled thickness, layer by layer with atomic precision along the normal to the exposed (100) plane [10;14]. In addition, their excitonic and charge carrier properties can also be easily tuned by varying the passivating ligands during the synthesis. [15] Due to their importance for photovoltaic applications, a better understanding of the light absorption properties of such objects is fundamental. However, the resulting CdSe sensitizers are either finite but composed of thousands of atoms as in the case of the QD species or extended 2D periodic objects such as NPL systems, preventing the use of the most advanced electronic structure methods to simulate and better understand their unique optical properties.

From a computational point of view, the study of excited state properties and phenomena such as the simulation of UV-Vis. spectra remains significantly challenging and computationally demanding when treating either large finite systems or extended periodic objects. For example, for the latter case, formulations have been proposed only recently, making use of either Gaussian-type Orbitals (GTO) [16;17;18] or projector augmented wave [19] basis sets. In practice however, most studies addressing the calculation of the photophysical properties of such materials rely on extracting one part of the system, or cluster, and treating it at the time-dependent density functional theory (TD-DFT) level, which is well-established in the molecular field [see 20 for instance]. Further refinement of the cluster model can be performed to treat larger molecular aggregates by partitioning this molecular cluster into different layers to which different levels of theory are applied, either belonging to the quantummechanics (QM) or molecular mechanics (MM) fields, resulting in the so-called QM/MM [21] or QM/QM' [22; 23;24;25;26;27;28] techniques. In most cases however, the effect of the crystalline environment surrounding the extracted cluster is either absent or only partially addressed.

To simultaneously overcome the barrier of cluster size and retain a high-fidelity computational model, we recently proposed the SC-Ewald approach [29] to simulate the photophysical properties of 3D periodic materials, inspired by the original work of Derenzo et al. [30] which focused on the study of localized structural defects in ionic crystals. SC-Ewald is based on the concept of simulating the environment produced by a periodic system by an array of point charges, which are fitted to reproduce the exact electrostatic potential (at the ground or excited state) of an infinite periodic system. The photophysical properties of an infinite periodic solid can thus be computed at moderate computational cost by carrying out calculations on a molecular cluster treated at the QM level, while its surroundings are approximated by point charges to reproduce the electrostatic interaction of the extended system with the cluster. In addition, since the array of point charges can be self-consistently equilibrated with a given excited state of the molecular cluster, not only absorption but also emission spectra can be simulated. So far, this approach has been successfully applied to different phenomena occurring in 3D periodic systems such as excited-state intramolecular proton-transfer [31], aggregation induced emission [32] and mechanochromism [33]. A similar method has been lately developed and applied in the Crespo-Otero group (see for instance ref. 34). Here, we present an extension of the SC-Ewald approach to systems with reduced periodicity, that is systems that are periodic in one (polymers, helices, nanotubes) or two directions (surfaces). As a case study, we consider the application of this method to investigate both quantum-confined QD and NPL of CdSe, along with bulk CdSe. We note that, to the best of our knowledge, only one similar approach based on electrostatic embedding and applicable to periodic systems in one, two or three directions has already been proposed: the Periodic Electrostatic Embedded Cluster Method (PEECM) [35], although no applications to excited-state properties calculations have been reported so far.

The paper is structured as follows: in Section 2, the SC-Ewald approach is briefly recalled and its extension to systems with reduced periodicity is presented. Computational details are given in Section 3, while in Section 4, results obtained on CdSe bulk, QD and NPL are presented and discussed.

Finally, in Section 5, some conclusions and perspectives of this work are drawn.

The methods

We first briefly recall the SC-Ewald approach originally developed for 3D periodic systems, before introducing modifications necessary to deal with 1D or 2D periodic systems. For 3D periodic systems, a large supercell obtained by replicating (2Nx x 2Ny x 2Nz) crystallographic unit cells in the positive and negative directions about the origin of the crystallographic axis system is built, with atomic charges at the positions of the crystallographic atoms. This is further divided into three distinct zones:

(i) an inner first zone (Z1) containing a cluster of interest with 𝑁 1 atoms to be treated explicitly at the QM level; (ii) an intermediate second spherical zone (Z2) with 𝑁 2 sites (𝑁 2 > 𝑁 1 ), containing point charges held constant at values that are equal to the partial Mulliken atomic charges provided by periodic calculations, surrounding the first zone; (iii) an outer third zone (Z3) with 𝑁 3 sites (𝑁 3 ≫ 𝑁 2 ) between the boundary of Z2 and the entire array formed by the (2Nx x 2Ny x 2Nz) unit cells, containing variable point charges. The intermediate second zone with point charges held constant is introduced in order to avoid artificial over-polarization of the electron density in Z1 by the nearby point charges close to the QM boundary. In order to define the optimal array of charges, the following procedure is next applied:

1. The electrostatic potential of the infinite crystal lattice VE is computed at all sites of Z1 and Z2, using Ewald lattice formula [36]. For a 3D periodic system with 𝑁 atoms per unit cell with charges 𝑞 𝑛 (𝑛 ∈ [1; 𝑁]), a unit cell volume 𝑉 and a set of direct lattice vectors 𝑻, the electrostatic potential at site with coordinates 𝒙 𝑚 is given by:

𝑉 𝐸 3𝐷 (𝐱 𝑚 ) = ∑ ∑ 𝑞 𝑛 erfc(𝜉|𝐱 𝑚 -𝐱 𝑛 + 𝐓|) |𝐱 𝑚 -𝐱 𝑛 + 𝐓| 𝑁 ′ 𝑛=1 𝑻 + + 4𝜋 𝑉 ∑ 𝑒 - 𝑘 2 4𝜉 2 𝑘 2 ∑ 𝑞 𝑛 𝑒 -𝑖𝐤⋅(𝐱 𝑚 -𝐱 𝑛 ) 𝑁 𝑛=1 𝒌≠𝟎 - 2𝜉 √𝜋 𝑞 𝑚 (1)
where 𝑁 ′ indicates that the term (𝑛 = 𝑚, 𝑻 = 𝟎) is excluded from the real space sum, 𝜉 is a parameter controlling the relative decay of the real and reciprocal space sums and 𝒌 form a discrete set of vectors in the reciprocal space.

2. The values of the point charges in Z3 are determined by solving a set of N3 linear simultaneous homogeneous equations to reproduce VE at the crystallographic sites in Z1 and Z2:

∑ 𝑞 𝑖 + ∆𝑞 𝑖 |𝒓 𝒌 -𝒓 𝒊 | 𝑁 𝑇 𝑖≠𝑘 = 𝑉 𝐸 (𝒓 𝒌 ) (2)
where 𝑘 ∈ [1; 𝑁 1 + 𝑁 2 ] and 𝑁 𝑇 = 𝑁 1 + 𝑁 2 + 𝑁 3 the total number of point charges in the three zones. Since N1+N2<N3, this is an underdetermined system for which we seek the minimum-norm solution to a real linear least squares problem:

min ||𝐴𝑋 -𝐵||
This solution selects the smallest possible variations ∆𝑞 𝑖 of the charges in Z3. Additionally, the fitting procedure is constrained to produce a neutral charge array with a dipole moment of zero:

{ ∑ ∆𝑞 𝑖 𝑁 𝑇 𝑖>𝑁 1 +𝑁 2 = 0 ∑ ∆𝑞 𝑖 𝒓 𝑖 𝑁 𝑇 𝑖>𝑁 1 +𝑁 2 = 0 (3)
3. The quality of the charge-fitting is finally estimated by: (i) computing the root-mean square (RMS) error between the direct Coulomb sum potential and the exact Ewald potential for NR sites randomly chosen in Z1 and Z2. Typically, NR=1000; (ii) evaluating the range of variation of charges in Z3. A wide variation indicates a failure at some point in the procedure.

Furthermore we have tested that largest deviations of charges in Z3 occurs at its border, that is the farthest from the QC cluster region (zone 1).

When generalizing this approach to 1D or 2D periodic systems, some modifications need to be introduced. Although no application to 1D systems is presented in this work, formulas and modifications required in this case are given below for completeness. In particular, supercells obtained by replicating (2Nx) and (2Nx x 2Ny) crystallographic unit cells are built for 1D and 2D systems, respectively, and the site partitioning into three distinct zones mentioned above is carried out.

To determine 𝑉 𝐸 (𝒓 𝒌 ) at all sites in Z1 and Z2, Ewald formulas in reduced periodicity are used [37].

More precisely, for a system periodic in the 𝑥 direction with 𝑁 atoms with charges 𝑞 𝑛 (𝑛 ∈ [1; 𝑁]) per unit cell, a lattice parameter 𝑎 and a set of direct lattice vectors 𝐓, the electrostatic potential at site with coordinates 𝐱 𝑚 is given by:

𝑉 𝐸 1𝐷 (𝐱 𝑚 ) = ∑ ∑ 𝑞 𝑛 erfc(𝜉|𝐱 𝑚 -𝐱 𝑛 + 𝐓|) |𝐱 𝑚 -𝐱 𝑛 + 𝐓| 𝑁 ′ 𝑛=1 𝑻 + 1 𝑎 ∑ ∑ 𝑞 𝑛 𝑒 -𝑖𝑘 1 (𝑥 𝑚 -𝑥 𝑛 ) 𝐾 0 ( 𝑘 1 2 4𝜉 2 , 𝜌 𝑚𝑛 2 𝜉 2 ) 𝑁 𝑛=1 𝑘 1 ≠0 - 1 𝑎 ∑ 𝑞 𝑛 {𝛾 + log(𝜌 𝑚𝑛 2 𝜉 2 ) + 𝐸 1 (𝜌 𝑚𝑛 2 𝜉 2 )} 𝑁 𝑛=1 - 2𝜉 √𝜋 𝑞 𝑚 ( 4
)
where 𝐾 0 (𝑢, 𝑣) is an incomplete modified Bessel function of the second kind defined as:

𝐾 0 (𝑢, 𝑣) = ∫ 𝑑𝑡 𝑡 𝑒 -𝑢𝑡-𝑣 𝑡 ∞ 1 (5)
and evaluated using an adaptive Romberg method, 𝐸 1 (𝑣) is the exponential integral defined as:

𝐸 1 (𝑣) = ∫ 1 𝑡 𝑒 -𝑣𝑡 𝑑𝑡 ∞ 1 = ∫ 1 𝑡 𝑒 -𝑡 𝑑𝑡 ∞ 𝑣 , (6)
𝛾 is the Euler Mascheroni constant, 𝜌 𝑚𝑛 = √(𝑧 𝑚 -𝑧 𝑛 ) 2 + (𝑦 𝑚 -𝑦 𝑛 ) 2 and 𝑘 1 forms a discrete set of 𝐤 vectors in the reciprocal space.

Similarly, for a slab periodic in the 𝑥 and 𝑦 directions with unit cell area 𝐴 and a set of direct lattice vectors 𝑻:

𝑉 𝐸 2𝐷 (𝐱 𝑚 ) = ∑ ∑ 𝑞 𝑛 erfc(𝜉|𝐱 𝑚 -𝐱 𝑛 + 𝐓|) |𝐱 𝑚 -𝐱 𝑛 + 𝐓| 𝑁 ′ 𝑛=1 𝑻 + 𝜋 𝐴 ∑ 𝑞 𝑛 ∑ 𝑒 -𝑖𝐤⋅(𝐫 𝑚 -𝐫 𝑛 ) 1 𝑘 𝑔(𝑘, 𝑧 𝑚 -𝑧 𝑛 , 𝜉) 𝒌≠0 𝑁 𝑛=1 - 2√𝜋 𝐴 ∑ 𝑞 𝑛 ( 1 𝜉 𝑒 -𝜉 2 (𝑧 𝑚 -𝑧 𝑛 ) 2 + √𝜋 (𝑧 𝑚 -𝑧 𝑛 ) erf(𝜉(𝑧 𝑚 -𝑧 𝑛 ))) 𝑁 𝑛=1 - 2𝜉
√𝜋 𝑞 𝑚 (7) where:

𝑔(𝑘, 𝑧, 𝜉) = 𝑒 𝑘𝑧 erfc ( 𝑘 2𝜉 + 𝜉𝑧) + 𝑒 -𝑘𝑧 erfc ( 𝑘 2𝜉 -𝜉𝑧) (8)
and 𝒌 form a discrete set of vectors in the reciprocal space.

Finally, while the charge neutrality condition in Z3 (second equation in Eq. ( 3)) is still maintained, the dipole equation then becomes

∑ ∆𝑞 𝑖 𝑥 𝑖 𝑁 𝑇 𝑖>𝑁 1 +𝑁 2 = 0 and ∑ ∆𝑞 𝑖 𝑥 𝑖 + 𝑁 𝑇 𝑖>𝑁 1 +𝑁 2 ∑ ∆𝑞 𝑖 𝑦 𝑖 𝑁 𝑇 𝑖>𝑁 1 +𝑁 2 = 0
for 1D and 2D periodic systems, respectively.

Computational details

Ground state periodic calculations:

All ground state periodic calculations were carried out at the DFT level with the Crystal17 code [38][START_REF] Dovesi | CRYSTAL17 User's Manual[END_REF], which solves self-consistently both Hartree-Fock (HF) and Kohn-Sham equations using atomcentered Gaussian orbitals as basis set. Based on a previous benchmark on the performances of different exchange-correlation functionals for CdSe bulk [START_REF] Szemjonov | [END_REF] performed by some of us, the global hybrid B3PW91 functional has been selected for all periodic calculations. Additional semi-empirical London-type dispersion corrections were employed for the optimization of CdSe bulk, considering the DFT-D3 scheme [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements H-Pu[END_REF].

A double-zeta basis set augmented by one d polarization function, combined to the SBKJC pseudopotentials, have been used both for Cd and Se [42;43]. As a consequence, the following electrons were treated explicitly: 4s 2 4p 6 4d 10 5s 2 for Cd and 4s 2 4p 4 for Se. All ligands considered in the adsorption on the CdSe (100) surface were treated with all-electron basis sets, with the following contractions schemes: (9s, 3p, 1d) → [3s, 2p, 1d] for C, (14s, 6p, 1d) → [4s, 3p, 1d] for O, (7s, 1p) → [3s, 1p] for H, and (14s, 8p, 1d) → [5s, 4p, 1d] for S. [44] The following tolerance criteria were set for the accuracy of Coulomb and exchange series: 10 -7 as the overlap threshold of Coulomb integrals, 10 -7 as penetration threshold of Coulomb integrals, 10 -7

as overlap threshold for HF exchange integrals, 10 -9 and 10 -20 as pseudo-overlap for HF exchange

series. An extra-large DFT integration grid was employed throughout, with 75 radial and 974 angular points. As a result of a comprehensive test carried out on each system on the convergence of the total energy with the increase of the k-points number, shrinking factors of 6, 16 and 14 were selected for the NPL systems, Zinc Blende bulk and Wurtzite bulk respectively, thereby corresponding to 20, 145 and 192 k-points in the irreducible Brillouin zones of these systems. Full geometry optimizations (atomic coordinates and lattice parameters) were performed for all the periodic systems involved in the investigation, using default values [45] for convergence criteria.

CdSe NPLwere modeled as a zinc blende 2D periodic slab consisting of a film with a finite thickness exposing two free surfaces and composed of alternating layers of Cd and Se atoms parallel to the experimentally-exposed polar (100) crystallographic plane [46]. Since the experimental anionic ligands are known to link to Cd-rich surfaces, this surface was terminated by Cd layers on both top and bottom sides of the slab. Inorganic (SH -and OH -) and organic (HCOO -) anions were chosen as passivating ligands and the bidentate adsorption mode was considered following the results previously obtained in Ref [47]. HCOO was chosen to model the oleic acid since the long alkyl chain does not influence the electronic and optical properties of the ligand/nanoplatelet interface [47].

Embedded excited-state calculations:

All excited states calculations were performed at the TD-DFT level using a development version of the Gaussian code [48], considering the long-range corrected PBE (LC-PBE) functional [49;50]. All clusters used for excited-state calculations were extracted from the optimized structures obtained with periodic boundary conditions, either in three or two directions, for the bulk and QD models and for the slab models, respectively. Point charges used for the electrostatic embedding were determined according to the procedure described in Section 2. All spectra presented below were obtained by convoluting computed vertical transitions with Gaussian functions with half-width of 0.01 eV.

In addition, in order to more quantitatively estimate the spatial extent of a photo-induced charge transfer, the DCT index has been used. Developed by some of us [51], it quantifies the length of the hole-electron separation due to a charge transfer excitation computed as the distance of the hole and electron charge density distributions baricenters (R+ and R-, respectively).

𝐷 𝐶𝑇 = |𝑅 + -𝑅 -|
More details and examples of application can be found in ref [52;53].

Results and Discussions

To the best of our knowledge, there have been no previous attempts to model the excited state properties of CdSe bulk. On the other hand, several studies can be found in the literature of modeling the photophysical properties of CdSe quantum dots [54][55][56][57][58][59][60][61][62][63]. These works usually apply TD-DFT on bare and/or ligand-saturated nano-clusters composed by Cd and Se atoms [54,55,56]. This protocol has been used to model in particular the effect of the various capping molecules (e.g. OPMe3, NH2Me)

and of the solvent on the optical properties (such as UV-Vis. Absorption and emission spectra) of

CdSe nanoparticles [54]. However, this approach is computationally limited by the size of the cluster, and has been applied to model CdSe nanoclusters composed by up to tens of atoms. Recently, Nadler

and Sanz [57] simulated the optical spectra of CdSe nanoclusters both naked and capped with ligands composed of up to 68 atoms, using a real time TD-DFT approach (RT-TDDFT) [57]. They showed that RT-TDDFT can successfully reproduce the absorption spectrum of CdSe quantum dots at low computational cost. Consequently, they indicate the RT-TDDFT approach is more efficient for modeling medium-to large-sized systems than conventional linear-response TD-DFT implementations. However, this method's main drawback is that it does not allow for the assignment of the spectral features. As an alternative to first principles calculations, tight-binding models have also been used to model the size-dependence of the optical band-gap of CdSe QD [58, 59 , 60, 61].

More recently, many-body perturbation theory calculations have also proved to perform well when applied to investigate the excitonic effects in CdSe nanocrystals, predicting optical gaps in good agreement with experiment [62]. On the other hand, concerning CdSe NPL, few theoretical works have been reported so far [62,63], all of them using a tight-binding approach [62], eventually combined [63] with Configuration Interaction to simulate the excitonic properties of NPLs.

In the present work, we propose a computational protocol to model the excited state properties of both extended systems (such as CdSe bulk or 2D CdSe NPLs) and finite size nano-clusters, relying on the same theoretical approach based on DFT and TD-DFT. By using the SC-Ewald methodology we were able to calculate vertical excitation energies and reproduce the UV-Vis. absorption spectra of the systems under investigation at low computational cost.

In this section, we will first discuss the results obtained on CdSe bulk, considering two phases: zinc blende and wurtzite. Next, CdSe systems with lower dimensionalities will be described, starting with QD of increasing sizes followed by NPL. For this latter, in agreement with a previous study, a welldefined (100) polar orientation stabilized by three passivating anionic ligands will be considered.

CdSe bulk

CdSe is a II-VI semiconductor which can be obtained at atmospheric pressure both in the hexagonal wurtzite (P63mc) and cubic zinc blende (F-43m) phases, although CdSe single crystals are experimentally known to preferentially grow in wurtzite form. [64]. In both cases, Cd and Se are tetrahedrally coordinated and the two phases differ only in their stacking sequences: zinc blende shows an ABCABC sequence along the [111] direction while wurtzite has an ABABAB one along the c axis so that only the third nearest neighbor configuration is different. For this reason, their electronic properties are very similar. [START_REF] Szemjonov | [END_REF] Several studies can be found concerning the electronic and structural properties of bulk CdSe. [65][66][67][68] In particular, a previous study by some of us [START_REF] Szemjonov | [END_REF] led to the development of a computational procedure, based on the use of global hybrid functionals, that serves as a good description of the structural and electronic features of CdSe bulk structures. The same protocol is applied here with the aim of choosing the best DFT approach. However, while in ref. As mentioned above, to simulate the UV-Vis. spectra of the two CdSe bulk phases considered, the SC-Ewald approach was employed. In the case of Zinc Blende, its primitive unit cell with 2 atoms was chosen as the quantum cluster (QC) for all SC-Ewald calculations. A (6x6x6) supercell was built around the QC, accounting for a total of 3460 point charges, with 2937 point charges in Z3 and 523 in Z2 (see Figure SI.2). The resulting charges produce a RMS deviation from the Ewald potential of the order of 0.1 μV at all sites of Z1 and Z2 while the range of variation of charges in Z3 is between -0.475 |e| and + 0.466 |e|. Since the Zinc Blende polymorph is metastable at room temperature and pressure [64], we were unable to find suitable experimental UV-Vis. absorption spectra to compare with our computed results. However, our calculations reveal that Zinc Blende should not absorb in the visible region (see Figure 1a). Indeed, the computed absorption spectrum presents the first intense electronic transition (f=0.36 a.u) at 320 nm. This transition corresponds mainly to a one electron excitation from the HOMO-2 to the LUMO involving Cd s and Se p orbitals with an associated σ→σ* character. The associated DCT index (2.955 Å) confirms its weak charge transfer character. Computed UV-Vis. absorption spectra and experimental absorption band maxima positions (red vertical sticks, Ref. [START_REF] Parsons | [END_REF]) of CdSe wurtzite bulk. The molecular orbitals involved in the electronic transitions of each band are also plotted. Se and Cd atoms are depicted in yellow and cyan, respectively.

Following an analogous procedure, the Wurtzite phase UV-Vis. Absorption spectrum was also simulated. In this case, the primitive unit cell composed of Se atoms at the top of the VB and at the bottom of the CB. On the other hand, the second band computed at 512 nm arises from an electronic transition mainly corresponding to a one electron excitation from the HOMO-2 to the LUMO showing a n→σ* character and less intense than the first one (f=0.08 a.u.). This transition shows a more pronounced charge transfer character than the first one, as quantified by the associated CT distances of 2.972 Å for the 647 nm transition and 5.145 A for the 512 nm one. Indeed, in this latter case the hole-orbital is localized on one of the two selenium atoms, whereas the electron is delocalized over the other selenium and the nearby Cd atom present in the unit cell.

Overall, for CdSe bulk, a good qualitative and quantitative agreement of computed data with available experimental data has been obtained. In particular, for the UV-Vis. spectra simulation, the electrostatic embedding approach chosen here allowed to recover both short-and long-range electrostatic interactions at low computational cost, by considering a finite set of well-chosen point charges to reproduce the electrostatic effect of the infinite crystalline environment.

CdSe QD

To evaluate the performances of the SC-Ewald approach when applied to quantum-confined systems such as QD, the size-dependent UV-Vis. absorption spectra of CdSe QD was investigated, considering both zinc blende and wurtzite QDs. Since experimental QDs are known to be spherical with large sizes (in the 2-10 nm range) here, a 3D electrostatic embedding will be added around the cluster considered. Indeed, since the experimental QDs investigated here contain thousands of bulklike ordered atoms, their optical properties are simulated as being affected by the electrostatics produced by a bulk-like system. Additionally, since these two phases crystallize into the F-43m and P63mc space groups belonging to the cubic and hexagonal crystal systems, respectively, we also explored the impact of the asymmetry of the central QC on computed UV-Vis. spectra. This is here investigated by constructing QC of increasing sizes along a given lattice direction and verifying that the computed spectra are identical for symmetry-equivalent directions. SI.2) is reported, along with the associated largest QC dimension.

Experimentally, a redshift of the first absorption peak is observed when increasing the QD dimension. This is qualitatively reproduced by the computed data, although a direct comparison between the computed and experimental data is difficult since the sizes considered are not identical in both cases.

In particular, the QD sizes indicated in Figure 3 refer to the QD diameter in the case of the experimental data, while they correspond to the end-to-end length in the case of the modeled QC, which can then be considered as only a portion of the real quantum dots. Nevertheless, the two absorption peaks positions obtained experimentally for the smallest and the largest QD at 410 nm and 503 nm are well-reproduced by the computed data at 410 and 501 nm, although calculations are underestimating the red shift with increasing size. However, we note that both the shape of the whole spectrum and the relative intensities of the different bands computed are in line with the experimental data, constituting encouraging results for the application of the SC-Ewald approach to the modeling of the photophysical properties of difficult systems such as real QD made of thousands of atoms.

On the other hand, looking at the spectra obtained on naked QD models (i.e. without the application of the SC-Ewald embedding), reported in This kind of transition is actually affected by the size of the QC considered, decreasing in energy with the increase of the particle size. In particular by analyzing the DCT values reported in Table 1, we can note an increase of the charge transfer extent associated to the first bright transition as the size of the QC model considered increases, thereby resulting in the experimentally observed redshift of the UV-Vis. spectra.

Table 1. DCT values, excitation wavelength and oscillator strength computed for the first bright transition in the case of CdSe QC models of different sizes. From the computed UV-Vis. spectra reported in Figure 4, it is clear that, as expected, the spectra obtained for the a-and b-QC are identical to the ones previously presented in the zinc blende case (see Figure SI.4-SI.6). On the other hand, c-QC data is significantly different which can be related to the 6-membered CdSe rings in chair-like conformation stacked one upon the other along this direction, thereby assuming a layered arrangement in this case, which is different from the zinc blende c-QC case (see Figure 3). Another notable difference can be observed upon increasing the QC size since new peaks appear at increasingly higher wavelengths. To better clarify this point, the largest c-QC with 18 atoms was chosen for analysis of the transitions involved. In particular, the analysis of the lowest-energy band shows that the band maximum located at 627 nm (1.98 eV) is ascribed to a single electronic transition (S0→S5) of charge transfer character with an associated oscillator strength of 0. Overall, the electrostatic embedding procedure considered here allowed to successfully recover important environmental effects impacting the photophysical properties of large systems such as QD, at low computational cost.

Number

In particular, we were able to recover the size-dependent optical properties observed experimentally.

Indeed, the bathochromic shift of the first absorption peak of CdSe QD of increasing sizes could be well-reproduced. In addition, for the CdSe zinc blende c-QC considered here, the hole-electron distance increases with the number of layers of CdSe rings, thereby indicating that peaks appearing at higher wavelengths originate from long-range electrostatic interactions, which weaken with the increase of the separation between the charge carriers. Finally, we note that when the environmental effects obtained by the SC-Ewald method are not included in the calculations, computed UV-Vis.

CdSe nanoplatelets (NPL)

In this section, we explore the application of the 2D electrostatic embedding scheme to compute the UV-Vis. absorption spectra of polar CdSe NPL, passivated with different ligands. We also evaluate the effect of the quantum confinement by considering slab models with three different thicknesses along the non-periodic direction. Indeed, CdSe NPL are known to preferentially grow in the zinc blende phase, exposing the (100) polar orientation at the experimental level. The slab models considered are composed of alternating charged layers of Cd 2+ and Se 2-ions, so that an N-layer slab Compared to CdSe.SH, the absorption onset is blue-shifted for both CdSe.OH and CdSe.HCOO. As a result, the absorption spectra of CdSe NPL passivated with OH -and HCOO -ligands cover a very tiny portion of the visible region, thereby suggesting a decrease in the light absorption properties of the CdSe NPL when these two latter ligands are used for passivation instead of SH -, with possible important consequences for photovoltaic application of such systems. In addition, it is interesting to note that the S0→S1 electronic transitions computed at 403 and 393 nm for OH -and HCOO -, respectively, are characterized by weak intensity (low oscillator strengths) and are mainly localized on the core of the NPL (Table SI.9-SI.10). This is in clear contrast with the CdSe.SH case, where passivating SH-ligands were involved in the electronic transition, thereby indicating that SH -anions are better candidates than HCOOH -or OH -to passivate CdSe NPL sensitizers.

Conclusions

We have applied the SC-Ewald method to model the excited state properties of extended periodic objects by replicating their crystalline potential by means of an electrostatic embedding approach. As a result of the application of the SC-Ewald method, the UV-Vis. absorption spectra of the CdSe bulk (3D periodic), CdSe quantum dots (0D periodic) and CdSe NPL (2D periodic) were computed and showed very good agreement compared with available experimental UV-Vis. absorption spectra.

Calculations showed that the choice of the quantum cluster (QC) was an important step and clearly affecting the results obtained. The results obtained with minimal cluster for the bulk phases allowed to recover the experimentally observed spectra at minimal computational cost. Our approach is also able to reproduce the main features of QD and, in particular, the size-dependent optical properties observed experimentally (ie the bathochromic shift of the first absorption peak of CdSe QD of increasing sizes). Of note this is not the case when using naked cluster, that is when excluding the charge embedding. Finally, complex 2D structures, such as CdSe nanoplatelets passivated with different ligands could be studied and the effect of the thickness of the platelets on the absorption spectra correctly reproduced. Overall, these results are encouraging for the use of embedding techniques when predicting photophysical properties of low dimensional materials at low computational cost.

  combined to the double-zeta LANL2DZ basis set along with associated effective core potentials for all atoms (Cd, Se, C, H, O, S). LC-PBE belongs to the class of range-separated functionals, which have been developed to overcome some of the TD-DFT limitations originating from standard DFT when dealing with long range interactions. In particular, standard global hybrid functionals fail to correctly reproduce through-space particle-hole charge transfer interactions (arising upon electron excitation) due to the incorrect description of the asymptotic tail of the exchange potential at long range. On the other hand, long-range corrected functionals, including a variable percentage of HF like non local exchange depending on electron-electron distance, recover the correct asymptotic behavior of the exchange potential thereby allowing to better describe electron excitations of charge transfer character[49;50].

  40 a D2 dispersion correction scheme[69] was adopted in combination with each functional tested, we employed here the D3[70] version implemented in Crystal17 thereby providing slightly different numerical results but the same qualitative description. As reported in TableSI.1 both optimized lattice parameters and computed band gaps obtained at the B3PW91-D3 level for both phases are in an overall very good agreement with the experimental data and previous literature data[START_REF] Szemjonov | [END_REF]. In particular, the maximal errors on the lattice parameters are of -0.016 and -0.005 Å for the wurtzite and zinc blende phases, respectively, while the corresponding errors on the band gaps are of +0.278 and +0.137 eV, therefore validating the chosen DFT model for periodic calculations. As expected, similar Densities of States (DOS) are obtained for the two phases as reported in Figure SI.1: a sharp contribution of Se4sp orbitals around -13 eV, followed by a broader one around -8 eV mainly due to Cd4d orbitals with a small contribution from the Se4sp orbitals, indicating a p-d coupling characteristic of tetrahedrally-coordinated compounds. Indeed, in binary crystalline solids having Td site symmetry, the metal's d orbitals have the same representation at Γ than the anion's p orbitals, due to the crystal field splitting imposing a break of the metal's d states into two groups: the lower energy eg states and higher energy t2g state (coupling with the chalcogenide p orbitals). The top of the valence band (VB) is mainly due to Se4sp contributions, with a non-negligible overlap with Cd5sp orbitals, indicating a partial covalent character of the Cd-Se bond in these bulk materials. This is further confirmed by the analysis of the computed Mulliken atomic charges, with values of ±0.381 and ±0.383 |e| for the two species in the wurtzite and zinc blende phases, respectively. On the other hand, the bottom of the conduction band (CB) is mainly dominated by Cd5sp orbitals, with a small overlap with Se4sp contributions.

Figure 1

 1 Figure 1 (a) Computed UV-Vis. absorption spectra of CdSe Zinc Blende bulk. The molecular orbitals involved in the electronic transitions associated to the first absorption band are also plotted. (b)

Figures SI. 4 -SI. 6 .Figure 2 :Figure 3 , 2 .Figure 3 :

 462323 Figures SI.4-SI.6. It is clear that, for a given QC size, the computed spectra are independent from the expansion direction, nicely recovering the cubic symmetry of the zinc blende phase and highlighting the isotropic surroundings of the QC in such a material.

Figure SI. 7 ,

 7 broad low-lying (582-630 nm) first absorption bands whose intensities decrease with the QC size are obtained. More importantly, it is evident that the calculated spectra do not replicate the bathochromic trend with the QD size observed experimentally (see Figure3b) and overall, their shape and peak position do not resemble the experimental ones.To better understand the optical properties of QD nanoparticles, the nature of the electronic transitions contributing to the first absorption band have been further analyzed, selecting the 26 atoms QC as an example. As reported in Supporting Information (Figure SI.8) for this cluster the lowest-energy band located at 501 nm (2.47 eV) is dominated by one single bright transition (S0→S1) corresponding to a HOMO→LUMO excitation of charge transfer character as indicated by a DCT value of 4.980 Å. The HOMO is mainly composed by Se 4p orbitals while the LUMO is a combination of Se 4p and Cd 4s orbitals. Overall, for all the QC models (except for the 6 atoms QC being too small), the first bright transition which is responsible for the absorption onset implies a charge transfer from Se 4p to Cd 4s as it is evident from the analysis of the electron density difference between the first (in the case of 26 atoms, 14 atoms, 10 atoms and 6 atoms QC) or the second (in the case of 22 atoms and 18 atoms QC) excited state and the ground state (Figure SI.9).

Figure 4 :

 4 Figure 4: Computed UV-Vis. Absorption spectra of Wurtzite with QC composed of (a) 10/12 atoms and (b) 14/18 atoms, and built by expansion along the a (solid black), b (dashed red) and c (solid blue) lattice directions.

  19 a.u. This transition stems from an electronic excitation involving the HOMO-4, composed mainly by Se 4p orbitals, and the LUMO made up of Cd 4s orbitals. The charge transfer character of the electronic transition was quantified by evaluating the associated DCT as done for the Zinc Blende case (see Figure SI.10). The 7.256 Å value obtained further confirms the marked charge transfer character associated to this electronic transition. We note that a similar investigation carried out on the 12 atoms c-QC led to similar results which are detailed in Figure SI.11 of the supporting material. Indeed, the first band at 525 nm (2.36 eV) is associated to a single electronic transition (S0→S1) of HOMO→LUMO type and charge transfer character, with an associated DCT value of 4.013 Å.

spectra (see Figure SI. 7 )

 7 are unable to reproduce the experimental trend observed for the first absorption peak position, especially not showing the size-dependent optical shift typical of quantum confined nanoparticles.

5 .Figure 5 :

 55 Figure 5: 7-layers CdSe zinc blende (100) slab model passivated with SH -(a), OH -(b) and HCOO - (c)ligands with 4 layers of Cd 2+ (in violet) alternating with 3 layers of Se 2-(in green). White, yellow, brown and red balls correspond to H,S, C and O atoms respectively. A 2D periodic geometry optimization of the resulting CdSe.SH, CdSe.OH and CdSe.HCOO slabs was first performed at the B3PW91-D3 level, by fully relaxing both lattice parameters and atoms, for slab models having 5, 7 and 9 layers that we refer to as 5L, 7L and 9L in the following. As a result of the geometrical optimization, as can be evidenced in Figure 5, in the case of the CdSe.SH slab model, sulfide ions occupy similar positions as the Se atoms in the inner layers, thereby respecting the tetrahedral arrangement of the CdSe structure. Indeed, the computed S1-Cd1-Se1 valence angle (110°, see Figure 5 for labeling) is very close to the Cd1-Se1-Cd2 one (111°). Furthermore, Cd-S distances are very similar (2.598 Å) to Cd-Se (2.671 Å) ones. Interestingly, the H-S-Cd-Se dihedral angles in

Figure 6 :Figure 7 :

 67 Figure 6: Computed (solid lines) and experimental (dashed lines) UV-Vis. absorption spectra of CdSe.SH NPL of different thickness (5 layers in black, 7 layers in green and 9 layers in orange). Experimental data from Ref. 74. Experimentally, the absorption coefficient presents a step-like structure which can be interpreted by resorting to the model of the infinite one-dimensional quantum well. The first absorption band

Acknowledgments

D.L, I.C. and F.L. acknowledge the use of computational resources at Très Grand Centre de Calcul (TGCC) provided by GENCI through project A0050810135. I.C. gratefully acknowledges support from the European Research Council (ERC) for funding under the European Union's Horizon 2020 research and innovation programme (grant agreement No 648558, STRIGES CoG grant). I.C., F.L. and D.L. acknowledge support from the Marie Sklodowska Curie Innovative Training Network MSCA-ITN-2017 (grant number 765739 COSINE).

Data availability statement. The authors confirm that the data supporting the findings of this study are available within the article and its supplementary materials.