One-dimensional embedded cluster approach to modeling CdS nanowires

CdS quantum dots and wires supported on wide band gap semiconducting oxide substrates have recently emerged as promising materials for “third generation” solar cells.^1–3 Photoexcitation of charge carriers occurs when photons are incident on the CdS nanostructures. The induced current in the solar cell is driven by the injection of negative charge carriers into the conduction band of the oxide substrate (where they diffuse to the front electrode) and the recharging of the CdS nanostructures via a redox electrolyte connected to a counter electrode. Understanding this process in detail requires accurate predictions of the structures and optical properties of such CdS nanostructures and the nanostructure/substrate interface, and an analysis of the role of defects in the material.^4–12 

In addition to solar cell applications,^3,12,13 recent experimental work on one-dimensional (1D) CdS nanostructures has resulted in a range of other device applications such as waveguides,^14,15 light-emitting diodes,¹⁶ logic gates,¹⁷ lasers,^18–20 cathodoluminescence,²¹ field-emitters,^22,23 photodetectors,²³ gas sensors,^24,25 transistors,²⁶ and photoelectrochemistry.²⁷ Considering the importance of this material, both for its device potential and relevance to the understanding of the fundamental physics of anisotropic quantum confinement, there have been relatively few theoretical studies performed,^28–32 which have concentrated on nanowires of diameter ⩽5 nm and, though they have had success in modeling the electronic structure of such nanowires, they have not considered modeling defects, which are key to understanding their optoelectronic properties.³³ Each of these studies has used a plane-wave density functional theory (DFT) approach to model the material. As large supercells are required in order to reduce the interaction between neighboring periodically repeated nanowires, this approach can be quite computationally intensive. For charged defects, the size of supercell necessary to reduce the periodic image charge interaction sufficiently becomes intractable. There have been some attempts to tackle this problem in 1D systems (see, e.g., Refs. 34–36) which are, however, limited to treating thin wires, and fail to account for the polarizing effect of periodic image charges in the direction of the wire. Additionally, to model the electronic structure accurately the use of hybrid density functionals, where the exchange-correlation potential contains a proportion of Hartree-Fock (HF) exchange, is preferable. These methods typically increase the computational cost by approximately an order of magnitude. For these reasons it can be concluded that, for performing large-scale simulations and modeling defects in nanowire structures, plane-wave DFT is not suitable at present.

For bulk and two-dimensional (2D) systems, a useful method to determine accurate models of localized charged defects is to use the hybrid quantum mechanical/molecular mechanical (QM/MM) embedded cluster approach.^37–43 Our focus in the present work is on the implementation of this approach in the ChemShell package,^44,45 where the defect and immediately surrounding region (typically ∼100 atoms) are treated at a QM level of theory, using a Gaussian basis set and a hybrid functional for exchange and correlation to determine the electronic structure accurately. This cluster is embedded in a larger region (typically ∼10 000 atoms), which is treated at an MM level of theory using a polarizable shell model, to provide an accurate elastic and dielectric response of the infinite system to defect formation. The MM region has an “active” inner part, close to the QM region, which is allowed to relax explicitly, and a “frozen” outer part, which is polarized through long-range Coulomb interactions. An embedding potential is provided by large-core pseudopotentials on cationic sites in the interface between the QM and MM regions,⁴⁶ preventing spillage of electrons onto anionic sites in the MM region.⁴⁷ The entire system is surrounded by point charges, which are fitted to reproduce the Madelung potential of the infinite solid. The advantages of this approach are as follows:

• Artificial interactions between periodic images of supercells, which occur using plane-wave DFT and are particularly problematic in the case of charged defects, are avoided;

• There is no need for a charge compensating background when treating charged systems, as is the case using plane-wave DFT;

• An accurate account of polarization effects is given;

• Access to the vacuum level, which is not well defined in periodic methods using a plane-wave basis set, is provided, facilitating the calculation of ionization reactions with an absolute energy reference;

• Access to high levels of quantum chemistry theory to treat the QM region, e.g., coupled cluster, multi-reference or configuration interaction techniques, which are not implemented in standard plane-wave codes, is provided;

• The analysis of orbital contributions to each energy band is less complicated; in principle such analysis is easier using atom-centered local basis functions rather than plane waves, which is, however, realized in periodic codes with local basis sets.

Using this method, the defect properties of bulk ZnO^48–50 and GaN,^51,95 the adsorption of CO₂ on ZnO surfaces,⁵² and the band alignment of TiO₂ in the rutile and anatase phases⁵³ have been successfully modeled.

Due to the fact that periodic images of charged defects are avoided, the hybrid QM/MM embedded cluster approach is also ideal for modeling nanowires. To date, however, the model has not been applied to 1D systems. We therefore present in this paper an extension of the model to treat such systems. We use as our test case a 1D CdS linear chain, which has the advantage of being a simple system, while providing the features necessary to test 1D embedding, i.e., confinement in two orthogonal directions and periodicity in the third. We treat the pure CdS system, and the system containing an In or Cu impurity (in the appropriate charge states) substituting on a Cd site. We benchmark our model using plane-wave DFT calculations, finding excellent agreement between the calculated density of states (DOS) using both approaches (as the model system is simple, the plane wave DFT calculations are not excessively computationally expensive). We find good agreement in the description of the defect states (where applicable), and demonstrate the ability of our embedded cluster approach to determine ionization energies with an absolute energy reference. We also show that our results accord with available experimental data.

The rest of the paper is structured as follows: In Sec. II we describe our model and our calculations; in Sec. III we present our interatomic forcefield model; in Sec. IV we present the results obtained using our hybrid QM/MM model, comparing them with results we have obtained using a plane-wave DFT approach; and in Sec. V we summarize our study.

To treat a CdS 1D system using a hybrid QM/MM embedded cluster approach, our first requirement was an MM interatomic forcefield model, with polarizable shells on the ions, which accurately reproduces the bulk structural, elastic, and dielectric properties of CdS, and that stabilizes in the wurtzite phase at ambient conditions. We have fitted such an empirical forcefield model using the GULP code,⁵⁴ the parameters of which we present in Sec. III along with the resulting bulk properties of CdS.

Our QM/MM model consists of a line of alternating Cd and S ions, with a particular ion in the center (see Fig. 1). In this approach, the embedding potential was provided by large-core effective core potentials (ECPs) placed on cationic sites in the interface between the QM and MM regions, to facilitate which, we set our QM region to be anion-terminated, so that the interface region consisted of a single Cd ion either side. The MM region was divided into the active and frozen parts, the extent of which was allowed to vary.⁵⁵ We note here that the linear chain is in fact a saddle point for a 1D single chain of CdS. When the system is free to relax in all directions the chain “buckles”, which implies that, by only allowing relaxations along the direction of the chain, we are in fact modeling a strained system, where the “buckling” would be analogous to shear strain in bulk.

We next needed a means of terminating our clusters so that they accurately represent an infinite system, for which we placed point charges at either end of the entire cluster, beyond the frozen MM region and in line with the CdS chain. The value of charge and distance from the cluster that reproduces the Madelung potential of the infinite chain was determined analytically as follows. In this ideal infinite system, the Madelung field (V_{M, ideal}) acting on each ion is conventionally given by

where q is the charge of the ion and a₀ is the ion-ion separation (or bond length). In an embedded cluster system with 2n + 1 ions in the QM and active region, and 2N + 1 ions in the QM, active, and frozen regions (because the cluster has a central ion, it necessarily has an odd total number of ions), the Madelung field (V_{M, cluster, i}) acting on ion i in the active region is given by

where

$s≡Qq$ (Q is the value of the point charge used to terminate the cluster) and

$r \equiv \frac{R}{a_0}$

$r≡Ra0$ (R is the distance from the point charge to the center of the cluster). The values of Q and R were determined by minimizing the function

with respect to s and r.⁵⁶ 

Once the model, i.e., the atomic coordinates and charges, was set up, we employed the existing ChemShell architecture as an interface between the QM and MM drivers, allowing self-consistent field and geometry optimization calculations to be performed. We have used GULP⁵⁴ as the MM driver and GAMESS-UK⁵⁷ as the QM driver in our simulations. Our QM region was modeled using a valence triple-zeta plus polarization Gaussian basis set.^58–60 The Perdew-Burke-Ernzerhof fitting to the generalized gradient approximation (GGA), corrected for solids,^61–63 with 25% HF exchange (PBEsol0), hybrid functional, was used for XC. A 28 electron ECP was used to model core electrons in Cd.⁶⁴ The embedding potential was provided by a large core (46 electron) pseudopotential,⁶⁵ placed on the interface Cd ions. In modeling Cu and In impurities, we have used 10 and 28 electron ECPs, respectively.^66,67

To model charged defects, an additional energy term, W_p, due to the polarization energy induced by the defect, had to be accounted for. The first approximation to this term depends only on the charge and not the configuration of the defect, so it may be added a posteriori to our calculation. For 2D and 3D systems, Jost's formula is used to calculate the relevant term.⁴⁸ For 1D linear systems, a modified expression for this energy term must be used:

where Q_d is the charge of the defect, R_a is the distance from the defect to the boundary between the active and frozen regions, and ε is the dielectric permittivity of the 1D system. We present the derivation of this term in the Appendix, including details on how ε is determined. Equation (4) differs from Jost's formula by the factor

$16π1Ra2$⁠, which reflects the reduction of the system dimensionality: a single line cannot provide as much dielectric response as a system extended in more than 1D, such as homogeneous dielectrics or their interfaces. As a result of this factor, W_p constitutes a significantly smaller correction to the defect energies than in bulk. Even in the case of an infinitely polarizable medium (i.e. a metallic wire) and a relatively small active region, e.g., R_a = 10 Å, the polarization energy of the remainder of the system due to the presence of a singly charged defect will contribute only −0.1 meV to the total energy.

To benchmark the results using our embedded cluster model, we have performed plane-wave DFT calculations on a CdS linear chain using the VASP code.^68–71 The 4d states of Cd were treated explicitly as valence, as were the 3d and 4d states of Cu and In, respectively, when modeling defects. The interaction between the core electrons and the valence electrons was included by the standard frozen-core projector augmented-wave (PAW) potentials.⁷² We fixed the bond length at the same value as used in our ChemShell calculations and employed the PBEsol0 hybrid functional for exchange and correlation, in order to allow direct comparison between the calculated DOS from our embedded cluster approach and from the plane-wave DFT approach. We found that the total energy was converged within 10⁻⁵ eV for a plane-wave energy cut-off of 1200 eV, a 1 × 1 × 16k-point grid, and a vacuum region extending for 15 Å in the x and y directions for a chain lying along the z-axis. These were the settings used in all subsequent calculations, apart for DOS calculations where a 1 × 1 × 60 k-point grid was used to achieve an accurate integration of the Brillouin zone. For pure CdS we have used a 2-atom cell, and for defect modeling we have used 24-atom supercells, which is the maximum supercell size feasible on current supercomputers, given simulation run times and memory requirements. For the supercell calculations, a 1 × 1 × 2 k-point mesh was found to give total energies converged within 10⁻⁵ eV.

The interatomic potential (IP) model used in our study is based on the Born model of ionic solids.⁷³ We simulate ion-ion interactions using a sum of three terms: The first is a Coulomb sum:

where U_ij is the energy of interaction and r_ij is the separation between ions i and j, and q_i is the charge on ion i; the second is a two-body Buckingham potential, with a damped dispersion term, of the form

where the parameters A, ρ, C, and B depend on the particular interaction considered; and the third is a two-body Morse potential of the form

where D_e, a, and r₀ are the species-dependent parameters. In order to describe the dielectric and structural properties accurately, the polarizability of the ions must be taken into account. We do this using the shell model of Dick and Overhauser,⁷⁴ where each ion is separated into a core and shell. The massless shell (of charge Y) is connected to the core by a spring. The energy U_c−s is expanded as

where k₂ and k₄ are the spring constants and r_c−s is the distance between the core and shell. The free ion polarizability α is given by α = Y²/k₂. In a lattice, the ionic polarizability is affected by the local structure and electrostatic field distribution. By fitting Y, k₂, and k₄ (if required), we account for these effects in our IP model automatically. The sum of the core and shell charge is constrained to equal the formal charge of the ion. The pair potential parameters used in this study for CdS are presented in Table I. The range of the interactions in Eqs. (6) and (7) is 15 Å.

We give a summary of the bulk physical properties of CdS in the wurtzite phase predicted by our IP model compared with experimental results in Table II (where ε⁰ is the static dielectric tensor and ε^∞ is the high frequency dielectric tensor). The agreement between the model and experiment, in particular for ε^∞, is good, although the static dielectric response is underestimated. For the purposes of the current study the IP model is acceptable.

The equilibrium bond length of the linear system using our interatomic potential model has been determined (using the GULP code) to be a₀ = 2.37 Å. This value was used when cutting all clusters, and used in all plane wave DFT calculations.

We have tested the convergence in cluster size of our embedded cluster approach, performed calculations on the pure CdS linear chain system, and performed calculations on the CdS linear chain system containing an impurity at a cation site.

We have first performed tests on convergence in MM cluster size in our embedded cluster approach. We varied the sizes of the active and frozen MM regions for a system with a QM region consisting of 33 ions (with a S ion at the center) and a system with a QM region of 35 ions (with a Cd ion at the center). A schematic of our embedded cluster model is shown in Fig. 1. If we define the origin of coordinates to be at the center of the cluster, set the linear chain to lie along the x axis, and the interface between one region and the next to occur at the midpoint of the bond between the ions at their edges, the 33 ion QM cluster extends from −39.11 to 39.11 Å along x and the 35 ion QM cluster extends from −41.48 to 41.48 Å in x. The interface region between the QM and MM clusters consists of a single Cd ion, implying that the interface extends a distance of 2.37 Å beyond the QM region in the positive and negative x directions. Varying the extent of the active and frozen regions from 16 ions to 32 ions on either side of the QM region (i.e., from a distance of between 37.92 and 75.84 Å beyond the interface region in the positive and negative x directions), we find a difference in the calculated energy spectrum (defined as the maximum error in the energy levels) of less than one part in 10⁴, which indicates that an active and frozen MM region consisting of 16 ions each, on either side of the QM region, is well converged. We therefore set our MM regions to this size; the active MM region extends 37.92 Å beyond the interface region in both the positive and negative x directions, and the frozen MM region extends a further 37.92 Å beyond the active MM region in the positive and negative x directions. The entire cluster is terminated by point charges, whose charge and distance from the cluster are fitted (using TERMINATE⁵⁶) to reproduce the Madelung potential of the infinite linear chain.

To test convergence in QM cluster size we have performed a series of DOS calculations, varying the number of ions in the QM cluster from 11 ions (so that the cluster ranges from −13.04 to 13.04 Å in the x direction) to 43 ions (so that the cluster ranges from −50.96 to 50.96 Å in the x direction). We present the results in Fig. 2. For each calculation the active and frozen regions consist of 16 ions each, on either side of the QM cluster, with a single Cd ion acting as an interface between the QM and MM regions. A smearing value of 0.1 eV was used to calculate the DOS. From this plot it is evident that the DOS is well converged, particularly in terms of the band gap and the band widths, for QM clusters of 33 ions or greater.

In Fig. 2 there are several peaks in the valence band spectrum (highlighted by arrows) that decrease approximately linearly with QM cluster size, indicating that they are size-dependent artifacts. The charge density associated with the states corresponding to these peaks is largely confined to the QM cluster edge, which implies that they are related to the termination of the QM region. We therefore decided to vary the embedding potential, in order to remove or minimize these artifacts of the embedding model. The initial embedding potential was provided by placing a large core pseudopotential on the Cd ions in the interface region, with 46 electrons in the core.⁶⁵ This pseudopotential consisted of S, P, and D non-local projectors. To allow variations, we constructed a simpler local pseudopotential from the original non-local pseudopotential, where the local Gaussian functions consisted of the radial part of the P projectors. We first tested that this procedure gave the same results as using the original potential by performing a DOS calculation for an embedded cluster with 37 ions in the QM region, which gave identical results to a calculation using the original embedding pseudopotential. We next varied the prefactors in the Gaussian functions used and observed the effect on the peaks highlighted in Fig. 2. By decreasing the prefactors by a factor of six we found that these artifacts were reduced substantially. We therefore chose the resultant local pseudopotential as our embedding potential.

Fig. 3 shows the calculated energy band diagram for the pure CdS linear chain (with a bond length of 2.37 Å), determined using our plane-wave DFT approach. The bands display a cosine dispersion, as would be expected from a simple Kronig-Penney model. Interestingly, we find that the band gap (2.38 eV) is indirect. The top of the valence band is formed by S 3p states, which are split into a doublet and singlet, with the doublet being the higher energy state. The doublet consists of p orbitals perpendicular to the direction of the chain, and is lower in energy (at k = 0 in the Brillouin zone) when alternating orbitals are in phase, due to the fact that, when in phase, the electron distributions sum together where they overlap, causing their densities to be more diffuse, and hence, the electrons have lower momentum and therefore, lower kinetic energy. At the edge of the Brillouin zone (k = π/2a₀), where the orbitals are out of phase, the electron distributions cancel where they overlap, and hence the densities occupy lower volumes, resulting in the electrons having greater momentum and hence higher kinetic energies; hence, the energy at the edge of the Brillouin zone is higher than at the Γ-point, and we observe an indirect gap.⁷⁹ 

Next we compared the DOS of a CdS linear chain obtained with our hybrid QM/MM approach with a 51 ion QM cluster using the reparameterized embedding pseudopotential, and using a plane-wave DFT approach. In both cases, we set the Cd–S bond length equal to 2.37 Å, and use the PBEsol0 hybrid XC functional. Our results are shown in Fig. 4. Both calculated DOS have been aligned so that their Fermi levels are equal. The agreement between the two calculations, in particular at the top of the valence band (below 0 eV) and bottom of the conduction band (above 1.8 eV), is striking. There is some disagreement in the higher conduction bands (at ∼7 eV), which is probably due to misrepresentation of very high unoccupied bands by the Gaussian basis set used in the embedded cluster approach. There is also a very slight disagreement in the deep “valence” band at ∼−12 eV; the embedded cluster approach determines this band to be shifted upwards by 0.1 eV relative to that determined by the plane-wave DFT approach. By calculating the energy levels of a sulfur atom using the two approaches, we confirm that this offset originates from the different approaches to the description of the interaction between the core and valence electrons employed by the two methods: In our plane-wave DFT approach, we use the PAW method to treat the core-valence interaction,⁷² while in our localized orbitals approach to treat the QM region in our embedded cluster model the core-valence interaction is facilitated by the energy-consistent Stuttgart-Dresden small-core ECP.⁶⁴ However, the agreement between the two calculations in describing the upper valence bands, the band gap, and the lower conduction bands is excellent. There are some disagreements in the intensities of the peaks, which are related to how the DOS are normalized. Overall, we see that our embedded cluster approach gives a very similar description of the most important features of the DOS in determining properties of interest of the material to that of the plane-wave DFT description.

On inspection of the energy levels in the different valence bands and comparing to molecular calculations, we can attribute the lowest plotted band (at ∼−12 eV) to S 3s states. The band lying above (at ∼−10 eV) can be attributed to Cd 4d states, with the next two lowest bands (one at ∼−4 eV and the other just below 0 eV) attributed to S 3p states. There is a splitting in this band due to the symmetry of the system, with S 3p states along the linear chain (i.e., σ-bonds) being lower in energy than the 3p states perpendicular to the linear chain (i.e., π-bonds).

The addition of impurity ions in 1D CdS nanostructures has led to some interesting properties. Doping of CdS and Cd_1−xZn_xS nanoribbons with Cu has been found to have consequences for photoactivity and photoconductivity, with evidence of de-localized and localized donor and acceptor states in the band gap.⁸⁰ Photoemission from states in the band gap has also been observed when Cu is added to Mo-doped CdS nanowires,⁸¹ and increases in conductivity have been observed in Cu-doped CdS nanoparticles.⁸² Similar results for photoactivity and photoconductivity, showing evidence of deep defect levels, have been observed in In-doped CdS nanoribbons⁸³ and nanowires.⁸⁴ It has been found, however, from optical transmission measurements on CdS thin films⁸⁵ and conductivity and electroluminescence measurements on CdS nanobelt-based field-effect transistors,⁸⁶ that In acts as a shallow donor, with no deep levels in the gap. Theoretical and computational studies on such systems are required in order to elucidate the exact role of the substitutional defect, as well as the stability of charge carriers with respect to the formation of charge-compensating intrinsic point defects. The key challenge in the description of defect states in extended systems is an accurate treatment of defect levels, or ionization energies required to transfer a charge carrier from its defect trap to a conduction state.

In this section we study the effect of adding a Cu and an In impurity on a Cd site in our linear chain model, to form Cu_Cd and In_Cd substitutional defects. We perform relaxations using our QM/MM approach (with a 51 ion QM region) and plane-wave DFT approach (with a 24 ion supercell), comparing the two. Our goal here is to demonstrate the ability of our QM/MM model to treat such defects, as the application of our QM/MM model to more realistic systems is more feasible in contrast to the plane-wave DFT approach.

Replacing a Cd ion with a Cu impurity in the neutral charge state (i.e., forming a

$Cu Cd ×$ defect, where we use the standard Kröger-Vink notation) leads to a decrease in the distance to the nearest-neighbor S ions from 2.37 Å to 2.19 Å (2.24 Å) using our QM/MM (plane-wave DFT) approach. Although the distances calculated by the different approaches are in reasonable agreement, their difference can be largely attributed to the fact that the relatively small supercell (24 ions) used in the periodic plane-wave DFT approach effectively reduces the relaxations of the ions. In Fig. 5, we show the resulting spin DOS compared with the pure CdS system, calculated using both approaches. The descriptions of the upper valence bands (from ∼−5 to 0 eV) by the two approaches are in good agreement. There are several peaks in the deeper “valence” states (at −6.5, −10.6, and −13.7 eV) in the QM/MM result (Fig. 5(b)) which are not reproduced in the plane-wave DFT result (Fig. 5(a)); these peaks are associated with Cu 3s and 3p states, which are treated explicitly in the QM/MM approach but modeled as core states in the plane-wave DFT approach. There is a small discrepancy in the description of the bottom of the conduction band by the two approaches. Both show evidence of an additional state just below the bottom of the conduction band, but the QM/MM

${\rm Cu}_{\mathrm{Cd}}^{\times }$

$Cu Cd ×$ defect calculation results in a blue shift of the lowest unoccupied states compared with the pure CdS calculation, which is not observed in the plane-wave DFT approach. This is probably related to the larger relaxation of the ions in the QM/MM approach.

Forming an

$In Cd ×$ defect leads to an increase in the distance to the nearest-neighbor S ions from 2.37 Å to 2.46 Å (2.45 Å) using our QM/MM (plane-wave DFT) approach. As with the

${\rm Cu}_{\mathrm{Cd}}^{\times }$

$Cu Cd ×$ defect, the two approaches are in good agreement, with the difference largely arising from the small supercell (24 ion) used in the plane-wave DFT approach. We show the resulting spin DOS, compared with the pure CdS DOS, from both approaches in Fig. 6. There is good agreement between the two approaches in describing the valence bands, as in the

${\rm Cu}_{\mathrm{Cd}}^{\times }$

$Cu Cd ×$ defect case, apart from several peaks (at −1.5, −6.5, −10.6, and −13.7 eV) in the QM/MM result (Fig. 6(b)) which are absent in the plane-wave DFT result (Fig. 6(a)). These peaks are associated with In 4s and 4p states, which are explicitly treated in the QM/MM approach but are modeled as core states in the plane-wave DFT approach. In contrast to the

${\rm Cu}_{\mathrm{Cd}}^{\times }$

$Cu Cd ×$ defect case, the description of the bottom of the conduction band by the two approaches is consistent, reflecting the close agreement in the relaxation of the ions in both calculations. Both approaches show evidence of an extra state just above the top of the valence band and one just below the bottom of the conduction band, associated with the In impurity.

Further, to study the nature of the states introduced by the

$Cu Cd ×$ and

${\rm In}_{\mathrm{Cd}}^{\times }$

$In Cd ×$ defects, we have plotted the spin density of both in Fig. 7 (determined using plane-wave DFT). In the case of

${\rm Cu}_{\mathrm{Cd}}^{\times }$

$Cu Cd ×$⁠, we observe a strong localization of the hole state on the Cu

$3d_{z^2}$

$3dz2$ orbital, accompanied by a pronounced spin polarization of the nearest-neighbor S ions, from which we can infer that, in the case of high-level Cu-doping, the system would exhibit ferromagnetic ordering. For the

${\rm In}_{\mathrm{Cd}}^{\times }$

$In Cd ×$ defect, the spin density is significantly more delocalized over In 5p and Cd 5p orbitals, which constitute the π-bonding orbital system that forms the lower conduction band.

Finally, we have calculated the vertical ionization energies (for electrons and holes) of the pure CdS system and the system containing the Cu and In impurity, using our QM/MM approach (see Fig. 8). Such calculations are not possible with the plane-wave DFT approach, as one does not have access to the vacuum level. We first determined the electron ionization energy of the pure system, by taking the energy difference between the system in the neutral and positive charge states, with all electronic degrees of freedom relaxed, which placed the top of the valence band 6.97 eV below the vacuum level. Next we determined the electron affinity of the pure system (or alternatively the hole ionization energy), by taking the energy difference between the system in the neutral and negative charge states, with all electronic degrees of freedom relaxed, which placed the bottom of the conduction band 3.88 eV below the vacuum level. This gave us a band gap of 3.09 eV, substantially larger than that determined from the Kohn-Sham levels (2.38 eV), but interestingly in reasonable agreement with the direct band gap determined using plane-wave DFT (3.17 eV, see Fig. 3), which may be related to the nature of the self-consistent field solver used by molecular QM methods (based on atomic orbitals that favor localized states), and its inherent difficulty in modeling the hole state at the zone edge. Turning to the Cu impurity, we find that the

$Cu Cd ×$ defect is electronically stable, requiring 3.30 eV to excite an electron to the edge of the conduction band. We determine that the energy required to excite a hole to the top of the valence band (or, equivalently, to trap an electron from the top of the valence band) is 3.22 eV. In both cases the ionization energy exceeds the band gap, indicating that the defect is resonant. Relaxing the atomic structure, we observe that, if an electron is trapped by this defect, forming

${\rm Cu}_{\mathrm{Cd}}^{\prime }$

$Cu Cd ′$⁠, it is a deep donor, requiring 0.22 eV for the electron to be excited to the conduction band. Turning now to the In impurity, as in the case of the Cu impurity we find that the neutral

${\rm In}_{\mathrm{Cd}}^{\times }$

$In Cd ×$ defect is electronically stable, requiring 3.14 eV (2.89 eV) to excite an electron (hole). This places the defect in the valence band with respect to electron ionization, and in the band gap with respect to hole ionization, indicating that it acts as a very deep acceptor. When a hole is trapped by the defect, forming

${\rm In}_{\mathrm{Cd}}^{\bullet }$

$In Cd •$⁠, the state relaxes until it lies below the top of the valence band, indicating that it will not be stable, as ionizing the hole gains 0.89 eV. When an electron is trapped, forming

${\rm In}_{\mathrm{Cd}}^{\prime }$

$In Cd ′$⁠, it is a deep donor, lying 0.91 eV below the conduction band minimum.

In presenting the above results we have concentrated on methodological aspects: We have demonstrated the capabilities of the hybrid QM/MM approach to tackle the main challenges in the defect characterisation in 1D, with a focus on the electronic structure. It is of interest to compare our current results with those known experimentally for the bulk material to assess the extent to which confinement in 2D affects the fundamental electronic properties of a typical wide-gap semiconductor such as CdS. As one would expect, the localized valence states are affected least, with the ionization potential of the linear chain (6.97 eV) being only 4% lower than the experimentally observed value for bulk (7.26 eV).⁸⁷ Both the calculated optical gap from the hybrid QM/MM approach (3.09 eV) and the direct gap from the plane-wave DFT approach (3.17 eV) could be compared with the band gap of 2.58 eV for the bulk material at 1.8 K.⁸⁸ The pronounced blue shift we predict can be attributed to the quantum confinement, and indeed has been observed in real CdS nanowires (see, e.g., Ref. 89). The width of the upper valence band, ΔV_B = 4 eV, (see Fig. 4) is in turn lower by ∼1 eV than the experimentally observed bulk value.⁹⁰ We determined ΔV_B = 4.8 eV for the bulk system using our plane-wave DFT approach, whereas Bosco et al.⁹⁰ reported a value ΔV_B = 4.5 eV, determined using the HSE06 hybrid functional⁹¹ with the fraction of Hartree-Fock exchange fitted to reproduce the bulk band gap. The calculated defect levels of Cu and In impurities can be broadly correlated with those observed, e.g., by photoluminescence in real systems,^81–84 but their exact position is of course strongly affected by the idealized atomic structure of the linear chain. For example, the

$In Cd ×$ acceptor level (2.89 eV above the valence band maximum) could be compared with the 2 eV level reported for nanowires.⁸⁴ Similar observations are made for Cu impurities,⁸¹ but further analysis will require calculations on more realistic systems.

Although the current model is idealized, with advances in microscopy and related techniques the fabrication and manipulation of ultrathin (one atom thick) chains of atoms has become a reality.^92,93 The preparation of a CdS linear chain may therefore be possible in the near future, and the present study offers insight into the electronic properties of such a system.

We have developed an embedded cluster approach to treat 1D systems, based on the hybrid QM/MM approach implemented in ChemShell. We have taken a CdS linear chain as our test case, and compared calculated DOS with those obtained using plane-wave DFT calculations implemented in VASP, using the same hybrid functional (PBEsol0) for exchange and correlation and the same Cd–S bond length. We find excellent agreement between the two approaches, indicating that our embedded cluster approach is valid. We have also used our approach to treat Cu and In impurities in the system, comparing our results with those obtained using a plane-wave DFT approach. We find good agreement between the calculated defect levels using both approaches, and demonstrate the power of the embedded cluster approach in determining ionization energies with an absolute energy reference, which is not possible using a plane-wave DFT approach. As the hybrid QM/MM method avoids spurious interactions between periodic images of charged defects, it will be better suited to treating defects in more realistic 1D systems, which are key to understanding their optoelectronic properties. Another advantage of the QM/MM method is that access to high levels of theory to treat the QM region, such as meta-GGA or the coupled cluster or configuration interaction techniques, is available, in contrast to the majority of plane-wave codes.

We stress that, although we have used a CdS linear chain as our test case, the method will be applicable to any material with a band gap. The method will also be applicable to extended systems with 1D features, such as dislocations in solids, steps on surfaces, and supported nanowires and nanotubes. When treating such systems, the interaction with the environment will determine the details of the model set up, such as extensions of the QM region to the atoms at the interface and the level of the MM treatment of the environment. For example, the interaction between typical wide-gap semiconductor wires and polymeric or metallic supports will mostly be of a mechanical nature, with a possible addition of image forces. On the other hand, the modeling of a step on a CdS surface would necessarily involve the extension of the QM/MM boundary to the surface itself, while maintaining the one-dimensional periodicity.

The work has been performed under the HPC-EUROPA2 project (Project No. 228398) with the support of the European Commission–Capacities Area–Research Infrastructures. J.B., A.W., C.R.A.C., and A.A.S. acknowledge funding from EPSRC Grant No. EP/IO1330X/1. S.T.B. acknowledges support from the Spanish government (Grant Nos. FIS2008-02238 and MAT2012-30924) and the Generalitat de Catalunya (Grant Nos. 2009SGR1041 and XRQTC). The authors also acknowledge the use of the UCL Legion High Performance Computing Facility (Legion@UCL) and associated support services, the IRIDIS cluster provided by the EPSRC funded Centre for Innovation (EP/K000144/1 and EP/K000136/1), and the HECToR supercomputer through membership of the UK's HPC Materials Chemistry Consortium, which is funded by EPSRC Grant No. EP/F067496. We are also grateful for funding obtained from the Royal Society under their International Joint Projects 2009/R3 scheme that helped establish the collaborative project between the co-authors based in London and Barcelona.

To model the dielectric response of the 1D linear chain due to the presence of a point charge Q_d at x = 0, we introduce a continuous medium approximation in full analogy with the 3D case. We recall that the polarization density P of a linear dielectric is

where χ is the susceptibility and E is the electric field, and the electric displacement D is

where ε is the dielectric constant (throughout this derivation we use atomic units). Whereas for the 3D case the susceptibility is a property related to a unit volume, in the present 1D case we relate it to a unit length. Considering a very small segment of the linear chain, Δx, over which the external electric field can be assumed to be constant, the energy required to induce a dipole (ΔW^ind) is

Once induced, the interaction energy ΔW^dip of this dipole with E is given by

which implies that the total polarization energy ΔW^pol of the line segment is given by

The electric field of the point charge, screened by the dipoles along the line, is given by

Integrating Eq. (A5) from x = R_a to x = ∞, where R_a is the distance from the point charge to the edge of the active region, and multiplying by 2, we determine the total polarization energy W_p in the frozen region to be

which is Eq. (4) given in Sec. II of the main text.

In the context of atomistic molecular mechanical simulations, the dielectric constant can be obtained from 1D periodic calculations using interatomic potentials.^54,94 For the 1D linear chain, we have

where a is the unit cell length, q is the vector of point charges in the unit cell, and H_Cart is the matrix of second derivatives of the potential energy with respect to particle coordinates.