Group 16 substitution effects on the electronic and optical properties of 1-dimensional conjugated heterocycle oligomers: a relativistic self-interaction corrected real-space TDDFT study

Abstract

First-principles calculations employing real-space self-interaction corrected (time-dependent) density functional theory with incorporation of scalar and spin–orbit relativistic corrections have been used to elucidate the electronic characteristics and optical response of linear heterocycle-based oligomers as a function of Group 16 atom substitution. Analysis reveals intense substitution-dependent ultraviolet π–π* transitions in all oligomeric systems investigated. Upon p-type doping via single-hole introduction, these materials exhibit characteristic polaron absorptions in the near-infrared region. In the neutral species, absorption in the visible region is seen only for the heaviest Group 16 substitutions. However, upon oxidation all of the oligomers display significant visible light absorption. The systematic modulation of both optical and band gap properties through chalcogen substitution and charge-carrier doping provides a framework for the engineering of organic materials with applications in multiple technological domains.

---

1. Introduction

From the discovery by John Wesley Hyatt in 1869 of polymers to substitute the use of ivory, tortoiseshell, horn, and linen to the revolution of Wallace Carothers with the development of Nylon, the first synthetic fibre, the field of polymer science has expanded dramatically. This expansion has included developing organic materials with the intention of replacing current inorganic materials for their electronic and optical applications: photovoltaic cells,¹ transistors,² thermoelectric generators,³ light-emitting devices,⁴ and even sensors.^5,6

To meet the specific electronic and structural requirements for real-world applications, the development of organic materials has focus on several key characteristics. These materials should feature a narrow band gap to optimize light absorption and energy efficiency, as well as a low-lying highest occupied molecular orbital for stability. A specific example where these considerations are paramount is the case of organic materials for use in dye-sensitised solar cells. Here, the lowest unoccupied molecular orbital (LUMO) of the donor should be tuned to lie slightly above the LUMO of the acceptor material, reducing energy waste during charge transfer. These precise adjustments provide efficient charge transfer between the donor and acceptor while minimising energy losses.^7–9 Another key aspect is the ability to control the assembly and morphology of the organic materials, which is critical for efficient charge transport.^1,6,10

One research topic expanding with this focus is the development of organic conjugated polymers/oligomers. These pi-conjugated systems, distinct from sp hybridized bonds, are characterized by an extended framework of alternating sp² and/or sp-hybridized atoms wherein unhybridized p orbitals engage in lateral overlap to form a continuous π-electron cloud. This electronic architecture facilitates extensive delocalization of π-electrons throughout the molecular backbone. These materials show promising new avenues, as the length of the polymer can be modified and the individual moieties can be substituted depending on the specific characteristics required for the intended application.^3,6,11

Among the various modifications, substituting heteroatoms within the monomer units presents an effective strategy for tailoring the properties of the material. In the case of heterocyclic monomers, modifying the heteroatoms can alter the optical and electronic characteristics of the conjugated material. By replacing the oxygen heteroatom with a heavier chalcogen, such as sulfur, selenium, or tellurium, and moving down the periodic table, the changes in size, electronegativity, polarizability, and increased strength of intermolecular interactions also lead to a reduction in the HOMO–LUMO gap. This tunable band gap combined with the ease of processing organic materials make them attractive candidates for optoelectronic applications.^6,12,13

In our investigation of novel organic materials for photoelectric applications, we examined the natural product R-telomestatin (Fig. 1), produced by Streptomyces anulatus 3533-SV4, which inspired both our previous work¹⁴ and the current study. Previously, all-trans linear analogs of telomestatin were examined with chain lengths varying from one to ten oxazole units and evaluated under different charge conditions (neutral, electron-doped, and hole-doped).

Fig. 1 Structural comparison of telomestatin in its native cyclic conformation (top) and its linear all-trans derivative (bottom).

In the present study, we concentrate exclusively on heptameric azole chains, as this corresponds to the number of conjugated oxazole ring units in the natural product, whilst also being large enough to adequately display the electronic and optical properties characteristic of the longer chains studied previously. We begin with a summary of the theoretical approaches used and then move on to study the effect of substituting the oxygen in the linear oligomeric oxazole framework with heavier chalcogen elements from Group 16. This modification aims to preserve the inherent physical characteristics of the polymer while fine-tuning its photoelectric properties.^14–16

In order to investigate the effect of Group 16 substitution on ground-state properties of the polyazoles we study the evolution of their electronic density of states (DOS) using self-interaction corrected real-space Density Functional Theory (SIC-DFT) to gain insight into changes in the electronic structure on varying the heteroatom. Following this, the optical absorption properties of the substituted chains are studied using real-space time-propagation SIC-TDDFT to probe the excited-state electronic behaviour of the oligomers in their neutral and hole-doped (monocationic) states. Spin–orbit corrections are included throughout to provide a balanced description of the azole chains for all of the chalcogen-substituted azole systems.

In addition, we here study the effect of removal of an electron to produce hole-doped monocationic chains. The anionic electron-doped species were not considered as the previous work had shown that the effects of addition or removal of a single electron had qualitatively very similar effects on the resulting absorption spectra.

2. Theory

2.1. DFT in real space

In the current work the OCTOPUS code was used for all optical response and density of states (DOS) calculations. OCTOPUS employs regular spatial grids within a pre-defined volume of space and therefore avoids the use of a conventional (e.g. Gaussian-type) basis-set expansion of the Kohn–Sham wavefunctions.^17–20 This has the advantage of making convergence of properties such as energies dependent only on the grid spacing and simulation volume, and avoiding errors due to basis-set incompleteness.

Density functional theory (DFT),^21,22 with non-hybrid functionals of the electronic charge density is well suited to discretization on real-space grids due to the (semi-)local form of the functionals.²³ In addition to ease of convergence, a further advantage of this approach is that real-space calculations scale very favourably on (massively) parallel computing platforms and can therefore provide a level of accuracy for extremely large systems that can only be obtained for relatively small molecules with conventional quantum chemical DFT implementations.²⁴

The (time-independent) electronic ground-state of the system being studied is obtained through solution of a Kohn–Sham system of N coupled one-particle Schrödinger equations

h_KSϕ_m(r) = ε_mϕ_m(r) (m = 1,…,N) (1)

where r = (x,y,z) is the vector of coordinates in 3-dimensional space. The Kohn–Sham effective potential v_KS(r) is

v_KS(r) = v_ext(r) + v_H[n(r)] + v_xc[n(r)] (3)

This allows us to compute the electronic energy E_elec[n(r)] for our system by combining the energies arising from v_KS(r) acting on the total charge density n(r) with the sum of the individual kinetic energies ε^kin_m for all of the occupied states²³

Here, the external potential has been neglected for simplicity and the kinetic energy is obtained as

2.1.1. Collinear spin-DFT. In the collinear approximation within spin density functional theory (SDFT), it is assumed that the electron spins are aligned (anti)parallel with a reference axis for the system (usually taken to be the z-axis). Thus, the total charge density, n(r), can be expressed as a sum of the two separate spin densities

n(r) = n_↑(r) + n_↓(r) (6)

with

2.1.2. Non-collinear spin-DFT. Generalisation of the above leads to non-collinear SDFT, where the spin angular momentum vector s is no longer constrained to take one of two values (↑ or ↓). Instead, the electron's magnetisation can point in any direction and the wavefunctions, ϕ_m, are replaced with two-component spinors ϕ_mwhose components, ϕ_±, are complex functions, the combination of which permits full flexibility in the description of the system's magnetisation.

The energy of the system now depends on both the charge density, n(r), and the magnetisation density, m(r)

where σ = (σ_x,σ_y,σ_z) is the vector of 2 × 2 Pauli spin-matrices.

The additional freedom permitted by this formalism means that non-collinear calculations are able to describe in a natural way the response of electronic spin to applied magnetic fields. Similarly, non-collinear SDFT provides a route to non-perturbative (self-consistent) description of spin–orbit coupling effects in relativistic calculations.

The need for accurate descriptions of magnetic behaviour in advanced quantum materials has led to a considerable amount of recent work on theoretical methods and computational implementations of non-collinear (TD)DFT methods.^25–29 In addition to the OCTOPUS implementation used here, other codes such as GPAW and INQ permit self-consistent non-collinear SOC within the real-space (TD)DFT formalism for the accurate treatment of heavy-atom effects in materials science.^30,31

2.2. Self-interaction correction

In practical calculations with approximate exchange–correlation functionals, such as the local spin density approximation (LSDA) used in the present work, it is important to be aware of the existence of the so-called self-interaction error (SIE). SIE arises due to the interaction of an electron with itself through the total-density dependence of the electrostatic Hartree term in the KS potential (see eqn (3)).

The Hartree energy is given by

and clearly shows that in the limiting example of a one-electron system, with the electron at position r the energy still contains a contribution from the same electron's density at position r′.

Whilst the exact XC functional would fully cancel this spurious self-interaction, approximate functionals such as LSDA do not, leading to the result that for a given density functional approximation (DFA)

E^DFA_xc[n^1e(r)] + E_H[n^1e(r)] ≠ 0 (12)

where n^1e(r) in eqn (12) is a one-electron charge density. Thus, the presence of the Hartree self-interaction combined with its incomplete cancellation by the approximate XC functional leads to the remaining net self-interaction error found with common DFAs.

Different schemes have been developed over the years in an attempt to develop a method for self-interaction correction (SIC). Whilst exact SIC formulations do exist, such as the famous Perdew–Zunger method (PZSIC),³² these are numerically complicated and can be difficult to use in practice. Instead of performing the SIC on an orbital-by-orbital basis, as is done in PZSIC, an attractive alternative is to make use of the average SIE of the electrons in the system to correct the error in the total density.³³

For a given DFA in the collinear SDFT approximation, this average density self-interaction correction (ADSIC) is given as

where both the energy and the SIC are functionals of the up and down spin densities n_↑ and n_↓.

Whilst the ADSIC method has been available for some time in spin-unpolarised and spin-polarised forms, the latter has been restricted to the collinear SDFT approximation. This has meant that systems displaying non-trivial magnetic behaviour or requiring spin–orbit coupling corrections could not easily be corrected to remove the SIE.

Both PZSIC and ADSIC were, however, recently reformulated to make them applicable to the non-collinear SDFT situation.³⁴ Whilst PZSIC remains troublesome to use in many cases, the non-collinear ADSIC method is as simple to use and as stable as the more familiar collinear implementation.

Rewriting eqn (13) in terms of the local charge and magnetization densities leads to the simple form shown below

However, it was shown by the authors that this form as it stands does not produce the correct collinear limit, which requires that the total energy and SIC depend only on n_↑ and n_↓ (see eqn (13) in ref. 34).

The approach taken by the authors in their OCTOPUS implementation avoids this issue by first diagonalising the spin-density matrix to obtain the spin-densities n_↑ and n_↓. These are then averaged, normalising by their integrals to define N_↑ and N_↓ in the frame defined by the local magnetisation. Thereafter, the potential is computed independently for n_↑ and n_↓ in this local frame. Finally, the potential is rotated back to the global frame using the total magnetisation.

3. Methods

3.1. Geometry optimisation

Geometry optimizations were conducted at the density functional theory (DFT) level employing the Orca computational chemistry package (version 5.0.4).³⁵ The re-regularized SCAN (r²SCAN) meta-GGA functional was selected based on its demonstrated capacity to yield geometries of equivalent or superior quality compared to more computationally intensive correlated wavefunction methods.^36,37

To address the inherent limitation of (semi-)local functionals in describing long-range dispersion forces, the density-dependent D4 dispersion correction was incorporated with r²SCAN throughout all geometry optimization procedures.^38–41 The r²SCAN-D4 methodology was used in conjunction with the flexible Def2-TZVP basis set.^42,43 The resulting geometries are available in the SI for this article.

As demonstrated in our previous research, this computational protocol provides sufficient accuracy for the present study with only very small alteration in the geometries on going to the much larger Def2-QZVP basis set.¹⁴

Re-optimisation of the hole-doped cationic species was not performed, as in our previous work on the polyoxazole chains very little difference had been found when separate geometry optimisations were performed (approximately 0.001 Å in bond lengths and 0.05° in angles). A test of this using the tellurazole chain found that the same was true of the heavier derivatives. This can be explained by extensive delocalisation of the hole, and the associated positive charge, throughout the conjugated π-system of the azole rings (Fig. S1) which effectively dilutes the geometry effects that might be expected in cases displaying greater localisation of the charge.

3.2. Time-propagation TD-ADSIC-SO

Real-space grid discretization of the (TD)DFT equations was performed using OCTOPUS version 14.1, with exchange and correlation terms approximated using the modified Perdew–Zunger local spin-density approximation (LSDA) functional^32,44,45 as provided by the libxc functional library version 5.1.0.^46,47

Real-space methods have the advantage of scaling efficiently, both with the size of the electronic system, and with the number of compute processes (thus making them particularly efficient for large calculations in distributed computing environments). Furthermore, the time-propagation approach to excited state electron dynamics yields the absorption spectrum of the system up to arbitrary frequencies without specifying the number of excitations required, and using only the occupied electronic states, unlike linear response methods such as solution of the Casida equation. Convergence of the properties calculated in this way depend only on the grid spacing and simulation volume for ground state DFT. Convergence of excitation energies also depends on the length of the time-propagation, but for transitions of interest in UV-Vis spectroscopy, this convergence is rapid. Nonlinearity and complex time-dependent processes such as collective electron dynamics are also naturally described by this approach.

In line with our previous work,¹⁴ the LSDA calculations were corrected using the average-density self-interaction correction (ADSIC) to eliminate spurious SIE and restore the physically correct −1/r asymptotic charge density decay.⁴⁸ The use of the ADSIC correction has been shown to produce results of similar quality to much more expensive methods such as CAS-SCF and EOM-CCSD, with a large part of the improvement over standard TD-L(S)DA being attributed to the improved description of long-range effects due to the correct asymptotic decay of the density tail.⁴⁹ Also, in our previous work on telomestatin derivatives TD-ADSIC was found to lead to excellent reproduction of the experimental absorption spectrum of the natural product.¹⁴

Spin–orbit coupling effects were incorporated self-consistently via the non-collinear spinor SDFT representation, and within the ADSIC-corrected LSDA theoretical framework outlined above, yielding the TD-ADSIC-SO method used throughout the current work. The method employs relativistic Kohn–Sham reference states subsequently modified to account for self-interaction errors anticipated in conventional density functional approaches. The resulting framework adequately addresses non-collinear spin structures, which prove essential for accurate spectroscopic predictions involving heavy-element-containing systems. In contrast to our previous methodology, which was applied only to relatively light elements, spin–orbit coupling effects were necessary in the present study to obtain physically correct descriptions of the transitions leading to the absorption spectra.³⁴

Fully-relativistic Hartwigsen–Goedecker–Hutter norm-conserving LDA pseudopotentials were used for all (TD)DFT calculations.⁵⁰ Based on prior convergence testing with linear six-membered oxazole oligomers,¹⁴ we employed a grid spacing of 0.2 Å and an atom-centered sphere radius of 6 Å, yielding Fermi energy convergence better than 0.01 eV in the calculation of the electronic ground states. This was checked for the tellurazole oligomer and the Fermi energy convergence was found to be the same, indicating that these parameters were suitable for even the heaviest of the substituted chains. Time-propagation TDDFT calculations utilized a time step of 0.003 ħ eV⁻¹ (0.00197 fs) for stable propagation over 20 000 steps, corresponding to total propagation times of 60 ħ eV⁻¹ (39.5 fs). The use of longer propagation times than in the previous work yielded better resolution of the individual peaks (0.1 eV), which was expected to be important in the study of potential peak splitting due to spin–orbit effects as the shorter propagation time used in the previous work did not produce full resolution of some of the peaks.¹⁴

For the hole-doped species, a single electron was removed and the electronic ground state recalculated. All of the monocationic azole chains displayed doublet ground states with total magnetisation m = 1. The hole density was found to be fully delocalised along the length of the azole chains (Fig. S1).

In the calculation of the Density of States (DOS) for the azole chains, unoccupied states totalling 50% of the number of occupied states were included so that a reasonable coverage of the conduction band region was obtained. In order to provide better resolution of individual features, a broadening of 0.004 Hartree was applied when creating the DOS plots. The decay of the DOS in the gap region does not reach zero, however this is an artifact of the Lorentzian fitting (even though the broadening factor used was half of the default). In the gap regions between the valence and conduction band peaks, no extra states were found.

4. Results and discussion

4.1. Density of states

The ADSIC-SO density of states for the neutral chains in which the Group 16 heteroatoms have been substituted by moving down the periodic table from oxygen to tellurium is presented in Fig. 2. To ensure clarity and consistency in data interpretation, the energy axis has been normalized by subtracting the Fermi energy, E_Fermi, from the state energies.

Fig. 2 ADSIC-SO density of state (DOS) plots for neutral oligomeric chains. Energy scale normalised by shifting E_Fermi to zero for each oligomer. The non-zero value of the plot in the gap region is a result of the Lorentzian fit used to produce the plots but no states are present in this region.

Fig. 2 illustrates a decrease in the band gap as the chalcogen series transitions from oxygen to sulfur, selenium, and tellurium. This behaviour can be attributed to a systematic increase in atomic radii within Group 16; as one moves from oxygen to tellurium, atomic and covalent radii increase. This expansion creates greater electron-nucleus distances in heavier chalcogens, where the nuclear charge experiences significant screening from inner electron shells. The reduced effective nuclear charge combined with increased distance substantially weakens the coulombic attraction between the nucleus and valence electrons. These electrons experience diminished electrostatic binding forces, occupy higher energy states, and become increasingly polarizable. The progressive shielding effect and attenuated coulombic interactions directly manifest in the declining electronegativity values observed down the group, from 3.44 for oxygen to 2.10 for tellurium, reflecting the decreasing ability of these atoms to attract and retain their valence electrons.⁵¹

In contrast, as shown in Fig. 3, the contraction of the band gap is less marked in the hole-doped species, indicating that positive charge conditions lead to an altered electronic structure. With one fewer electron, the reduction in electron–electron repulsion allows the remaining valence electrons to relax into slightly lower-energy orbitals. Furthermore, the highest occupied molecular orbital (HOMO) in the neutral species is transformed into a singly occupied molecular orbital in the cation, with its energy decreasing due to the diminished exchange and Coulomb repulsions. The holes introduced are positioned within delocalized valence-band states (Fig. S1), effectively generating vacant states at the top of the valence band that can facilitate subsequent low energy electron transitions. This behaviour provides a rationale for the observed increase in valence band density in the hole-doped species relative to the neutral species of all the oligomers tested.⁵²

Fig. 3 ADSIC-SO density of state (DOS) plots for hole-doped oligomeric chains. Energy scale normalised by shifting E_Fermi to zero for each oligomer. The non-zero value of the plot in the gap region is a result of the Lorentzian fit used to produce the plots but no states are present in this region.

In direct comparison to telomestatin, the natural product serving as the initial structure for this investigation, the oxazole chain (both neutral and positively charged) exhibits an expansion of the band gap. Conversely, the thiazole, selenazole, and tellurazole neutral chains exhibit a reduction of the band gap, supporting the direction of increasing electronic delocalization down Group 16. The thiazole-based positively charged species maintains a band gap similar to telomestatin, contrasting the selenazole and tellurazole positively charged chains exhibiting further decrease of the band gap.

4.2. Optical absorption spectra

The optical absorption spectra calculated at the TD-ADSIC-SO level for the neutral oligomers are shown in Fig. 4. Table 1 contains energies in eV of the most important excitations observed in the spectra (and discussed below).

Fig. 4 Isotropic TD-ADSIC-SO absorption spectra of neutral oligomeric azole chains. Main peaks for each derivative discussed in text are labelled with asterisks and vertical dashed lines indicate visible region. Spectra shifted by +0.45 eV (see section 4.2).

Table 1 Major optical absorption feature energies in neutral chalcogen substituted polyazole chains. The features referred to here as the main peaks are highlighted in Fig. 4

Chalcogen	Absorption onset/eV	Main peak/eV
O	3.8	4.7
S	3.3	4.2
Se	3.0	4.0
Te	2.8	3.6

---

As in our previous work on polyoxazole chains of varying lengths,¹⁴ calculation of the absorption spectrum for the neutral form of the natural product telomestatin was performed in order to evaluate any shift of the calculated spectrum relative the experimental UV-Vis (SI of Doi et al., 2006⁵³). When the TD-ADSIC-SO telomestatin spectrum was compared with experiment, the correct peak distribution was obtained, but a shift of +0.45 eV was required to overlay the absorption maxima (Fig. S2). The previous calculations, which employed different pseudopotentials and neglected spin–orbit effects, required a smaller shift of +0.25 eV. However, it was felt that even this slightly larger shift was within the deviation to be expected from TDDFT when using semi-local functional approximations, and so this was applied to all subsequent absorption spectra obtained in the current work.

Consistent with the trend observed in Fig. 2, a red-shift in the absorption is seen as one moves down the periodic table, correlating to the systematic shrinking of the band gap. The oxazole chain displays the onset of absorption at 3.8 eV, which shifts to 3.3 eV for the thiazole chain. The absorption onset for the selenazole chain is similarly redshifted to a value of 3.0 eV, just inside the visible part of the spectrum. The tellurazole chain follows this trend, but its absorption onset is shifted to 2.8 eV, well inside the visible region. In addition, the tellurazole chain displays a second absorption peak in the visible at just over 3 eV, suggesting that it should be the most highly coloured of all the neutral systems studied here.

All of the systems studied display very strong absorption in the UV region, with the O, S, and Se compounds producing clear single peaks at 4.7, 4.2 and 4.0 eV, respectively. Our previous work on polyoxazole chains of varying length suggested that these peaks could be attributed to plasmonic excitations, due to the length-dependence of their intensity, and also the fact that they are polarised along the length of the chains and have negligible intensity when the excitation ‘kick’ is applied in a direction transverse to the chain (Fig. S4).¹⁴ In the case of the polytellurazole chain, however, the major peak at 3.6 eV is greatly diminished in intensity. Whilst the peak at 3.4 eV appears to be a shoulder similar to those that can be seen in similar positions for the other compounds, this reduction in the main peak height, along with the appearance of the shouldered peak at 3.8 eV, suggests involvement of the pronounced spin–orbit coupling associated with tellurium. An alternative explanation for this observation might be the broad conduction band characteristic of the neutral tellurazole chains, which can accommodate a wider range of excitation wavelengths.

In order to investigate this, the absorption spectrum of the neutral tellurazole compound was calculated using the same non-collinear SDFT approach including ADSIC, but neglecting the spin–orbit correction (Fig. S3). It can be seen that only a small part of the peak splitting described above can be attributed to spin–orbit coupling, meaning that the majority of the disruption of the expected single collective plasmon-like excitation must originate from a combination of scalar relativistic effects and shell-structure differences between Te and the lighter chalcogens.

In the hole-doped chains (Fig. 5), the large absorption peaks in the UV are shifted to lower values than those of the neutral chains, consistent with the reduction in the band gaps observed in Fig. 3. The intensity of these peaks undergoes distinct changes upon doping; the oxazole and selenazole chains experience a decrease in intensity, whereas the thiazole and tellurazole chains show an increase. Notably, the tellurazole chain, which initially exhibited three closely-spaced peaks in the neutral state, now only displays two peaks. Given the increase in intensity of the tellurazole peak at ∼3.8 eV, it is feasible that the missing peak has redshifted and combined with the central peak of the triad.

Fig. 5 Isotropic TD-ADSIC-SO absorption spectra of hole-doped oligomeric azole chains. Vertical dashed lines indicate visible region. Spectra shifted by +0.45 eV (see section 4.2).

The hole-doped species demonstrate a significant absorption feature in the infrared (IR) region, a behaviour that has been previously reported and attributed to a low energy polaron-like excitation (Fig. 6).^14,54,55 The existence of this new peak suggests that these compounds should display redox-switchable near-IR absorption, and that it is sensitive to the Group 16 substitution employed in the material.

Fig. 6 Schematic representation of semiconductor electronic structure depicting the fundamental band gap between the occupied valence band (VB) and the unoccupied conduction band (CB), with an emergent spin polaron (SP) state induced by hole doping. The empty circle denotes the excess hole. Band gap modifications due to lattice relaxation effects are not depicted.

In addition, on hole doping all of the compounds studied were found to display absorption in the visible. Furthermore, excitations across the whole visible region are seen, making it possible to optimise the visible absorption characteristics of a material based on these building blocks in a tunable way.

The computational findings are consistent with both experimental and theoretical studies reported in the literature.^56–58 Specifically, as lighter chalcogens are replaced with heavier ones, there is a systematic shift in the optical absorption edge toward longer wavelengths (lower energies). For example, in dye-sensitized TiO₂ complexes, substituting the chalcogen in a donor molecule from oxygen to sulfur, selenium, and ultimately tellurium results in a progressive red shift of the absorption maximum from 2.37 eV to 2.07 eV.⁵⁹ Similarly, in conjugated polymers, replacing thiophene with selenophene or tellurophene reduces the π–π* transition energy, leading to a red-shift in the absorption spectrum.^60–62

Inspection of the frontier charge distribution in the polyazoles studied here indicates that, in addition to general changes in the potential influencing the near-gap states, increasing involvement of the heavy chalcogen valence p-states in the delocalised π-orbitals can be expected for Se- and Te-containing conjugated systems (as can be seen in Fig. S1).

5. Conclusion

Using real-space time-propagation TDDFT with explicit self-consistent correction for both self-interaction error and the influence of spin–orbit coupling, this investigation has shown that the electronic properties of oligomers containing Group 16 atoms correlate systematically with the identity of the incorporated chalcogen. The results demonstrate a progressive decrease in band gap on going down the chalcogen group from oxygen to tellurium. This contraction of the band gap is less pronounced in hole-doped species, indicating that the presence of the positive charge leads to alterations in the potential determining the nature of the valence and conduction bands leading to electronic structures that partially mitigate the chalcogen-induced effects.

Comparative analysis with telomestatin reveals distinct patterns based on which chalcogen is incorporated into the linear structures. The oxazole chain, in both neutral and positively charged states, exhibits band gap expansion relative to the reference compound. Conversely, thiazole, selenazole, and tellurazole neutral chains demonstrate progressive band gap reduction, supporting the trend of increasing electronic delocalization down Group 16.

The observed correlation between electronic and optical properties and the nature of the chalcogen substitution reflects fundamental atomic characteristics of the chosen Group 16 element. As heavier chalcogens replace oxygen in the molecular framework, the optical absorption edge undergoes systematic red-shifts along with successive reductions in the band gap. The effect reaches maximum expression in the tellurazole system which displays spectral features distinct from those of the lighter chalcogens, including broadened and split absorbance peaks extending into the visible region. Although strong spin–orbit effects can be expected in Te-containing compounds, the absorption spectra studied here did not change significantly when these effects were neglected, suggesting that the influence of Te in this case is more likely to be due to shell-structure and scalar-relativistic effects. Spin–orbit interactions are likely to be much more important in other properties of the tellurazole chain where inter-state coupling is involved, such as fluorescence.

Charge-carrier doping effects include increase of the band gap, alteration of absorbance peak intensities, and introduction of new transitions, particularly in the infrared region, consistent with polaronic behaviour. The observed simplification of peak splitting in tellurazole upon doping suggests suppression of spin–orbit effects, providing an additional dimension of electronic tunability.

The findings of the current study complement previously obtained data regarding oligomer chain length effects, where extended conjugation promotes increased delocalization, leading to additional red-shifts in optical absorption behaviour and further band gap narrowing.¹⁴ The combination of chain-length modulation with chalcogen substitution and controlled doping therefore represents a powerful strategy for systematic tuning of optoelectronic properties in oligomeric heterocycle-based systems.

Future work will include evaluation of solvent/microenvironmental effects. These are expected to be most pronounced in the doped species due to charge redistribution, but even in the case of the neutral chains it will be interesting to see how their inclusion shifts the absorption spectra of the azole chains (and how much the empirical shift used in this work can be attributed to solvation). In addition, levels of theory beyond adiabatic semi-local TDDFT should be applied in order to obtain a more accurate descriptions of complex phenomena like the charge-induced polaronic excitations suggested by the present results and our previous work.¹⁴

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5dt02361k.

Acknowledgements

The authors acknowledge financial support from the Industrial Biotechnology Innovation Centre (IBioIC, doctoral training grant BB/W059899/1). This work was supported by UK Research and Innovation (UKRI) under the UK Government's Horizon Europe funding guarantee grant No. IFS 10057167 (University of Aberdeen) for EU project BlueRemediomics, which is funded by the European Union under the Horizon Europe Programme, Grant Agreement No. 101082304. The authors would also like to acknowledge the support of the Maxwell Compute Cluster funded by the University of Aberdeen. All 2D chemical structures in this work were created using Marvin version 23.1, ChemAxon (https://www.chemaxon.com).

References

1. A. Moliton and J.-M. Nunzi, Polym. Int., 2006, 55, 583–600.
2. C. D. Sheraw, L. Zhou, J. R. Huang, D. J. Gundlach, T. N. Jackson, M. G. Kane, I. G. Hill, M. S. Hammond, J. Campi, B. K. Greening, J. Francl and J. West, Appl. Phys. Lett., 2002, 80, 1088–1090.
3. I. H. Eryilmaz, Y.-F. Chen, G. Mattana and E. Orgiu, Chem. Commun., 2023, 59, 3160–3174.
4. J. Zhang, Y. Wang, S. Liu, H. Yu, L. Zhang and W. Xie, Appl. Phys. Lett., 2022, 120, 171102.
5. K. Aliqab, J. Wekalao, M. Alsharari, A. Armghan, D. Agravat and S. K. Patel, Biosensors, 2023, 13, 759.
6. M. Jeffries-El, B. M. Kobilka and B. J. Hale, Macromolecules, 2014, 47, 7253–7271.
7. A. K. Edward, K. Hosseini, K. L. Perry and D. S. Seferos, Polym. Chem., 2025, 16, 3187–3210.
8. H. He, Z. Zhong, P. Fan, W. Zhao and D. Yuan, Small, 2025, 21, 2405156.
9. M. Yang, J. Tan, Z. Wang, D. Li, W. Xu, H. Wang, L. Yuan, C. Yang and Z. Meng, J. Am. Chem. Soc., 2025, 147, 24152–24161.
10. X. Gong, M. Tong, F. G. Brunetti, J. Seo, Y. Sun, D. Moses, F. Wudl and A. J. Heeger, Adv. Mater., 2011, 23, 2272–2277.
11. C. K. Chiang, C. R. Fincher, Y. W. Park, A. J. Heeger, H. Shirakawa, E. J. Louis, S. C. Gau and A. G. MacDiarmid, Phys. Rev. Lett., 1977, 39, 1098–1101.
12. G. C. Hoover and D. S. Seferos, Chem. Sci., 2019, 10, 9182–9188.
13. S. Ye, V. Lotocki, H. Xu and D. S. Seferos, Chem. Soc. Rev., 2022, 51, 6442–6474.
14. J. Mérgola-Greef and B. F. Milne, Phys. Chem. Chem. Phys., 2023, 25, 12744–12753.
15. K. Shin-ya, K. Wierzba, K.-i. Matsuo, T. Ohtani, Y. Yamada, K. Furihata, Y. Hayakawa and H. Seto, J. Am. Chem. Soc., 2001, 123, 1262–1263.
16. M.-Y. Kim, H. Vankayalapati, K. Shin-ya, K. Wierzba and L. H. Hurley, J. Am. Chem. Soc., 2002, 124, 2098–2099.
17. N. Tancogne-Dejean, M. J. T. Oliveira, X. Andrade, H. Appel, C. H. Borca, G. Le Breton, F. Buchholz, A. Castro, S. Corni, A. A. Correa, U. De Giovannini, A. Delgado, F. G. Eich, J. Flick, G. Gil, A. Gomez, N. Helbig, H. Hübener, R. Jestädt, J. Jornet-Somoza, A. H. Larsen, I. V. Lebedeva, M. Lüders, M. A. L. Marques, S. T. Ohlmann, S. Pipolo, M. Rampp, C. A. Rozzi, D. A. Strubbe, S. A. Sato, C. Schäfer, I. Theophilou, A. Welden and A. Rubio, J. Chem. Phys., 2020, 152, 124119.
18. X. Andrade, D. Strubbe, U. De Giovannini, A. H. Larsen, M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas, I. Theophilou, N. Helbig, M. J. Verstraete, L. Stella, F. Nogueira, A. Aspuru-Guzik, A. Castro, M. A. L. Marques and A. Rubio, Phys. Chem. Chem. Phys., 2015, 17, 31371–31396.
19. J. Alberdi-Rodriguez, M. J. T. Oliveira, P. García-Risueño, F. Nogueira, J. Muguerza, A. Arruabarrena and A. Rubio, Computational Science and Its Applications – ICCSA 2014, Cham, 2014, pp. 607–622.
20. X. Andrade, J. Alberdi-Rodriguez, D. A. Strubbe, M. J. T. Oliveira, F. Nogueira, A. Castro, J. Muguerza, A. Arruabarrena, S. G. Louie, A. Aspuru-Guzik, A. Rubio and M. A. L. Marques, J. Phys.: Condens. Matter, 2012, 24, 233202.
21. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, 864–871.
22. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133.
23. F. Nogueira, A. Castro and M. A. L. Marques, in A Tutorial on Density Functional Theory, ed. C. Fiolhais, F. Nogueira and M. A. L. Marques, Springer Berlin Heidelberg, Berlin, Heidelberg, 2003, pp. 218–256.
24. J. Jornet-Somoza, J. Alberdi-Rodriguez, B. F. Milne, X. Andrade, M. A. Marques, F. Nogueira, M. J. Oliveira, J. J. Stewart and A. Rubio, Phys. Chem. Chem. Phys., 2015, 17, 26599–26606.
25. A. Grubišić-Čabo, M. H. Guimarães, D. Afanasiev, J. H. G. Aguilar, I. Aguilera, M. N. Ali, S. Bhattacharyya, Y. M. Blanter, R. Bosma, Z. Cheng, Z. Dan, S. P. Dash, J. M. Dueñas, J. Fernandez-Rossier, M. Gibertini, S. Grytsiuk, M. J. Houmes, A. Isaeva, C. Knekna, A. H. Kole, S. Kurdi, J. L. Lado, S. Mañas-Valero, J. M. J. Lopes, D. Marian, M. Na, F. Pabst, S. B. Pierantoni, M. Regout, R. Reho, M. Rösner, D. Sanz, T. van der Sar, J. Sławiñska, M. J. Verstraete, M. Waseem, H. S. van der Zant, Z. Zanolli and D. Soriano, 2D Mater., 2025, 12, 031501.
26. V. Blum, R. Asahi, J. Autschbach, C. Bannwarth, G. Bihlmayer, S. Blügel, L. A. Burns, T. D. Crawford, W. Dawson, W. A. de Jong, C. Draxl, C. Filippi, L. Genovese, P. Giannozzi, N. Govind, S. Hammes-Schiffer, J. R. Hammond, B. Hourahine, A. Jain, Y. Kanai, P. R. Kent, A. H. Larsen, S. Lehtola, X. Li, R. Lindh, S. Maeda, N. Makri, J. Moussa, T. Nakajima, J. A. Nash, M. J. Oliveira, P. D. Patel, G. Pizzi, G. Pourtois, B. P. Pritchard, E. Rabani, M. Reiher, L. Reining, X. Ren, M. Rossi, H. B. Schlegel, N. Seriani, L. V. Slipchenko, A. Thom, E. F. Valeev, B. V. Troeye, L. Visscher, V. Vlček, H. J. Werner, D. B. Williams-Young and T. Windus, Electron. Struct., 2024, 6, 042501.
27. Z. Pu, H. Li, N. Zhang, H. Jiang, Y. Gao, Y. Xiao, Q. Sun, Y. Zhang and S. Shao, Phys. Rev. Res., 2023, 5, 013036.
28. H. Li, Z. Pu, Y. Q. Gao and Y. Xiao, J. Chem. Theor. Comput., 2024, 20, 10477–10490.
29. T. Wang, H. Li, Z. Pu, Y. Q. Gao and Y. Xiao, J. Chem. Phys., 2025, 162, 214104.
30. J. J. Mortensen, A. H. Larsen, M. Kuisma, A. V. Ivanov, A. Taghizadeh, A. Peterson, A. Haldar, A. O. Dohn, C. Schäfer, E. Örn Jónsson, E. D. Hermes, F. A. Nilsson, G. Kastlunger, G. Levi, H. Jónsson, H. Häkkinen, J. Fojt, J. Kangsabanik, J. Sødequist, J. Lehtomäki, J. Heske, J. Enkovaara, K. T. Winther, M. Dulak, M. M. Melander, M. Ovesen, M. Louhivuori, M. Walter, M. Gjerding, O. Lopez-Acevedo, P. Erhart, R. Warmbier, R. Würdemann, S. Kaappa, S. Latini, T. M. Boland, T. Bligaard, T. Skovhus, T. Susi, T. Maxson, T. Rossi, X. Chen, Y. L. A. Schmerwitz, J. Schiøtz, T. Olsen, K. W. Jacobsen and K. S. Thygesen, J. Chem. Phys., 2024, 160, 92503.
31. J. Simoni, X. Andrade, W. Fang, A. C. Grieder, A. A. Correa, T. Ogitsu and Y. Ping, arXiv, 2025, preprint, arXiv:2506.21908v2, 1–36,  DOI:10.48550/arXiv:2506.21908v2.
32. J. P. Perdew and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 5048–5079.
33. C. Legrand, E. Suraud and P. G. Reinhard, J. Phys. B: At., Mol. Opt. Phys., 2002, 35, 1115.
34. N. Tancogne-Dejean, M. Lüders and C. A. Ullrich, J. Chem. Phys., 2023, 159, 224110.
35. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
36. J. Sun, A. Ruzsinszky and J. Perdew, Phys. Rev. Lett., 2015, 115, 036402.
37. J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew and J. Sun, J. Phys. Chem. Lett., 2020, 11, 8208–8215.
38. E. Caldeweyher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2017, 147, 034112.
39. E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2019, 150, 154122.
40. E. Caldeweyher, J.-M. Mewes, S. Ehlert and S. Grimme, Phys. Chem. Chem. Phys., 2020, 22, 8499–8512.
41. S. Ehlert, U. Huniar, J. Ning, J. W. Furness, J. Sun, A. D. Kaplan, J. P. Perdew and J. G. Brandenburg, J. Chem. Phys., 2021, 154, 061101.
42. A. Schäfer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571–2577.
43. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297.
44. P. A. M. Dirac, Math. Proc. Cambridge Philos. Soc., 1930, 26, 376–385.
45. F. Bloch, Z. Phys., 1929, 57, 545–555.
46. M. A. Marques, M. J. Oliveira and T. Burnus, Comput. Phys. Commun., 2012, 183, 2272–2281.
47. S. Lehtola, C. Steigemann, M. J. Oliveira and M. A. Marques, SoftwareX, 2018, 7, 1–5.
48. C. Legrand, E. Suraud and P.-G. Reinhard, J. Phys. B: At., Mol. Opt. Phys., 2002, 35, 1115–1128.
49. M. Oliveira, B. Mignolet, T. Kus, T. Papadopoulos, F. Remacle and M. Verstraete, J. Chem. Theory Comput., 2015, 11, 2221–2233.
50. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641.
51. L. C. Allen, J. Am. Chem. Soc., 1989, 111, 9003–9014.
52. M. Z. Maialle, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1967–1974.
53. T. Doi, M. Yoshida, K. Shin-Ya and T. Takahashi, Org. Lett., 2006, 8, 4165–4167.
54. J. L. Brédas, J. C. Scott, K. Yakushi and G. B. Street, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 30, 1023–1025.
55. J. L. Bredas and G. B. Street, Acc. Chem. Res., 1985, 18, 309–315.
56. M. Planells, B. C. Schroeder and I. McCulloch, Macromolecules, 2014, 47, 5889–5894.
57. M. Al-Hashimi, Y. Han, J. Smith, H. S. Bazzi, S. Y. A. Alqaradawi, S. E. Watkins, T. D. Anthopoulos and M. Heeney, Chem. Sci., 2016, 7, 1093–1099.
58. A. Bourouina, M. Rekhis and M. Trari, Polyhedron, 2017, 127, 217–224.
59. G. Deogratias, O. S. Al-Qurashi and N. Wazzan, J. Mol. Model., 2023, 29, 86.
60. B. Amna, H. M. Siddiqi, A. Hassan and T. Ozturk, RSC Adv., 2020, 10, 4322–4396.
61. M. Kawakami, K. H. G. Schulz, A. J. Varni, C. F. Tormena, R. R. Gil and K. J. T. Noonan, Polym. Chem., 2022, 13, 5316–5324.
62. T. Oyama, Y. S. Yang, K. Matsuo and T. Yasuda, Chem. Commun., 2017, 53, 3814–3817.

---