We present an embedded cluster model to treat one-dimensional nanostructures, using a hybrid quantum mechanical/molecular mechanical (QM/MM) approach. A segment of the nanowire (circa 50 atoms) is treated at a QM level of theory, using density functional theory (DFT) with a hybrid exchange-correlation functional. This segment is then embedded in a further length of wire, treated at an MM level of theory. The interaction between the QM and MM regions is provided by an embedding potential located at the interface. Point charges are placed beyond the ends of the wire segment in order to reproduce the Madelung potential of the infinite system. We test our model on the ideal system of a CdS linear chain, benchmarking our results against calculations performed on a periodic system using a plane-wave DFT approach, with electron exchange and correlation treated at the same level of approximation in both methods. We perform our tests on pure CdS and, importantly, the system containing a single In or Cu impurity. We find excellent agreement in the determined electronic structure using the two approaches, validating our embedded cluster model. As the hybrid QM/MM model avoids spurious interactions between charged defects, it will be of benefit to the analysis of the role of defects in nanowire materials, which is currently a major challenge using a plane-wave DFT approach. Other advantages of the hybrid QM/MM approach over plane-wave DFT include the ability to calculate ionization energies with an absolute reference and access to high levels of theory for the QM region which are not incorporated in most plane-wave codes. Our results concur with available experimental data.
I. INTRODUCTION
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,16 logic gates,17 lasers,18–20 cathodoluminescence,21 field-emitters,22,23 photodetectors,23 gas sensors,24,25 transistors,26 and photoelectrochemistry.27 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.33 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,46 preventing spillage of electrons onto anionic sites in the MM region.47 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 ZnO48–50 and GaN,51,95 the adsorption of CO2 on ZnO surfaces,52 and the band alignment of TiO2 in the rutile and anatase phases53 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.
II. METHODOLOGY
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,54 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.55 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.
Schematic of the embedded cluster model of the CdS linear chain. Horizontal dotted lines represent repeated patterns of ions. The model is repeated in the negative x direction.
Schematic of the embedded cluster model of the CdS linear chain. Horizontal dotted lines represent repeated patterns of ions. The model is repeated in the negative x direction.
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 (VM, ideal) acting on each ion is conventionally given by
where q is the charge of the ion and a0 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 (VM, cluster, i) acting on ion i in the active region is given by
where
with respect to s and r.56
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 GULP54 as the MM driver and GAMESS-UK57 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.64 The embedding potential was provided by a large core (46 electron) pseudopotential,65 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, Wp, 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.48 For 1D linear systems, a modified expression for this energy term must be used:
where Qd is the charge of the defect, Ra 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
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.72 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−5 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−5 eV.
III. INTERATOMIC POTENTIAL MODEL
The interatomic potential (IP) model used in our study is based on the Born model of ionic solids.73 We simulate ion-ion interactions using a sum of three terms: The first is a Coulomb sum:
where Uij is the energy of interaction and rij is the separation between ions i and j, and qi 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 De, a, and r0 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,74 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 Uc−s is expanded as
where k2 and k4 are the spring constants and rc−s is the distance between the core and shell. The free ion polarizability α is given by α = Y2/k2. In a lattice, the ionic polarizability is affected by the local structure and electrostatic field distribution. By fitting Y, k2, and k4 (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 Å.
Interatomic pair potential parameters for bulk CdS, including shell polarizations on Cd and S ions (e is the electronic charge).
Buckingham | A (eV) | ρ (Å) | C (eV Å6) | B (Å−1) |
Cd shell–S shell | 123.53 | 0.85 | 71.16 | 26.95 |
Cd core–Cd shell | 21.56 | 40.09 | ||
S shell–S shell | 277.04 | 15.66 | ||
Morse | De (eV) | a (Å−1) | r0 (Å) | |
S shell–S shell | 1.81 | 0.65 | 3.08 | |
Cd shell–S shell | 1.74 | 1.88 | 2.34 | |
Shell | k2 (eVÅ−2) | k4 (eVÅ−4) | Y(e) | |
Cd core–Cd shell | 25.20 | 3.82 | ||
S core–S shell | 58.46 | 35 101 | −3.14 |
Buckingham | A (eV) | ρ (Å) | C (eV Å6) | B (Å−1) |
Cd shell–S shell | 123.53 | 0.85 | 71.16 | 26.95 |
Cd core–Cd shell | 21.56 | 40.09 | ||
S shell–S shell | 277.04 | 15.66 | ||
Morse | De (eV) | a (Å−1) | r0 (Å) | |
S shell–S shell | 1.81 | 0.65 | 3.08 | |
Cd shell–S shell | 1.74 | 1.88 | 2.34 | |
Shell | k2 (eVÅ−2) | k4 (eVÅ−4) | Y(e) | |
Cd core–Cd shell | 25.20 | 3.82 | ||
S core–S shell | 58.46 | 35 101 | −3.14 |
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 ε0 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.
Material properties of wurtzite CdS determined experimentally and calculated using our interatomic potential model.
Property . | Experiment . | IP model . |
---|---|---|
Lattice parameter a, b (Å) | 4.13a | 4.14 |
Lattice parameter c (Å) | 6.71a | 6.68 |
Wurtzite parameter u | 0.377a | 0.379 |
Bulk modulus (GPa) | 62b | 74 |
$\epsilon _{11}^{0}$ | 8.53c | 7.04 |
$\epsilon _{33}^{0}$ | 9.62c | 7.93 |
$\epsilon _{11}^{\infty }$ | 5.23d | 5.16 |
$\epsilon _{33}^{\infty }$ | 5.29d | 5.39 |
Property . | Experiment . | IP model . |
---|---|---|
Lattice parameter a, b (Å) | 4.13a | 4.14 |
Lattice parameter c (Å) | 6.71a | 6.68 |
Wurtzite parameter u | 0.377a | 0.379 |
Bulk modulus (GPa) | 62b | 74 |
$\epsilon _{11}^{0}$ | 8.53c | 7.04 |
$\epsilon _{33}^{0}$ | 9.62c | 7.93 |
$\epsilon _{11}^{\infty }$ | 5.23d | 5.16 |
$\epsilon _{33}^{\infty }$ | 5.29d | 5.39 |
The equilibrium bond length of the linear system using our interatomic potential model has been determined (using the GULP code) to be a0 = 2.37 Å. This value was used when cutting all clusters, and used in all plane wave DFT calculations.
IV. RESULTS
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.
A. Convergence tests for embedded cluster approach
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 104, 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 TERMINATE56) 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.
The calculated DOS as a function of energy relative to the Fermi level for a series of embedded clusters with different QM cluster sizes.
The calculated DOS as a function of energy relative to the Fermi level for a series of embedded clusters with different QM cluster sizes.
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.65 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.
B. Pure CdS linear chain
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 = π/2a0), 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.79
Band dispersion for the pure CdS linear chain, determined using plane-wave DFT. Red and blue lines indicate conduction and valence bands respectively. The corresponding DOS is also shown.
Band dispersion for the pure CdS linear chain, determined using plane-wave DFT. Red and blue lines indicate conduction and valence bands respectively. The corresponding DOS is also shown.
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,72 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.64 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.
The calculated DOS as a function of energy relative to the Fermi level for an embedded cluster with a 51 ion QM region (red line) and that of a linear chain calculated using a plane-wave DFT approach (black line).
The calculated DOS as a function of energy relative to the Fermi level for an embedded cluster with a 51 ion QM region (red line) and that of a linear chain calculated using a plane-wave DFT approach (black line).
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).
C. CdS linear chain containing an impurity
The addition of impurity ions in 1D CdS nanostructures has led to some interesting properties. Doping of CdS and Cd1−xZnxS 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.80 Photoemission from states in the band gap has also been observed when Cu is added to Mo-doped CdS nanowires,81 and increases in conductivity have been observed in Cu-doped CdS nanoparticles.82 Similar results for photoactivity and photoconductivity, showing evidence of deep defect levels, have been observed in In-doped CdS nanoribbons83 and nanowires.84 It has been found, however, from optical transmission measurements on CdS thin films85 and conductivity and electroluminescence measurements on CdS nanobelt-based field-effect transistors,86 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 CuCd and InCd 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
The calculated spin DOS as a function of energy relative to the Fermi level for a pure CdS linear chain (black line), and that of a linear chain containing a single CuCd defect (red and green lines). The DOS were determined using (a) a plane-wave DFT approach and (b) a hybrid QM/MM approach with a 51 atom QM region. Negative spin DOS indicates spin-down. Pure CdS is spin degenerate, therefore only spin-up DOS is shown.
The calculated spin DOS as a function of energy relative to the Fermi level for a pure CdS linear chain (black line), and that of a linear chain containing a single CuCd defect (red and green lines). The DOS were determined using (a) a plane-wave DFT approach and (b) a hybrid QM/MM approach with a 51 atom QM region. Negative spin DOS indicates spin-down. Pure CdS is spin degenerate, therefore only spin-up DOS is shown.
Forming an
The calculated spin DOS as a function of energy relative to the Fermi level for a pure CdS linear chain (black line), and that of a linear chain containing a single InCd defect (red and green lines). The DOS were determined using (a) a plane-wave DFT approach and (b) a hybrid QM/MM approach with a 51 atom QM region. Negative spin DOS indicates spin-down. Pure CdS is spin degenerate, therefore only spin-up DOS is shown.
The calculated spin DOS as a function of energy relative to the Fermi level for a pure CdS linear chain (black line), and that of a linear chain containing a single InCd defect (red and green lines). The DOS were determined using (a) a plane-wave DFT approach and (b) a hybrid QM/MM approach with a 51 atom QM region. Negative spin DOS indicates spin-down. Pure CdS is spin degenerate, therefore only spin-up DOS is shown.
Further, to study the nature of the states introduced by the
The calculated spin density associated with (a) a
The calculated spin density associated with (a) a
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
The calculated defect ionization energy levels associated with a Cu and In impurity substituting on a Cd site, in various relevant charge states, shown relative to the vacuum level. The conduction band (CB) is represented by red shading, while the valence band (VB) is represented by blue shading. Black levels represent electron ionizations and blue levels represent hole ionizations. Neutral defects are denoted by a cross, positive defects by a dot, and negative defects by a prime.
The calculated defect ionization energy levels associated with a Cu and In impurity substituting on a Cd site, in various relevant charge states, shown relative to the vacuum level. The conduction band (CB) is represented by red shading, while the valence band (VB) is represented by blue shading. Black levels represent electron ionizations and blue levels represent hole ionizations. Neutral defects are denoted by a cross, positive defects by a dot, and negative defects by a prime.
D. Relevance to experiment
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).87 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.88 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, ΔVB = 4 eV, (see Fig. 4) is in turn lower by ∼1 eV than the experimentally observed bulk value.90 We determined ΔVB = 4.8 eV for the bulk system using our plane-wave DFT approach, whereas Bosco et al.90 reported a value ΔVB = 4.5 eV, determined using the HSE06 hybrid functional91 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
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.
V. CONCLUSIONS
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.
ACKNOWLEDGMENTS
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.
APPENDIX: Polarization energy in 1D
To model the dielectric response of the 1D linear chain due to the presence of a point charge Qd 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 (ΔWind) is
Once induced, the interaction energy ΔWdip of this dipole with E is given by
which implies that the total polarization energy ΔWpol 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 = Ra to x = ∞, where Ra is the distance from the point charge to the edge of the active region, and multiplying by 2, we determine the total polarization energy Wp in the frozen region to be
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 HCart is the matrix of second derivatives of the potential energy with respect to particle coordinates.