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.

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.

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.

FIG. 1.

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.

FIG. 1.

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.

Close modal

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

(1)

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

(2)

where

$s \equiv \frac{Q}{q}$
sQq (Q is the value of the point charge used to terminate the cluster) and
$r \equiv \frac{R}{a_0}$
rRa0
(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

(3)

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:

(4)

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

$\frac{1}{6\pi }\frac{1}{R_a^2}$
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, Wp 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., Ra = 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.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.

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:

(5)

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

(6)

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

(7)

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 Ucs is expanded as

(8)

where k2 and k4 are the spring constants and rcs 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 Å.

Table I.

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 Å6B−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−1r0 (Å)   
S shell–S shell 1.81 0.65 3.08   
Cd shell–S shell 1.74 1.88 2.34   
Shell k2 (eVÅ−2k4 (eVÅ−4Y(e  
Cd core–Cd shell 25.20   3.82   
S core–S shell 58.46 35 101 −3.14   
Buckingham A (eV) ρ (Å) C (eV Å6B−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−1r0 (Å)   
S shell–S shell 1.81 0.65 3.08   
Cd shell–S shell 1.74 1.88 2.34   
Shell k2 (eVÅ−2k4 (eVÅ−4Y(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.

Table II.

Material properties of wurtzite CdS determined experimentally and calculated using our interatomic potential model.

PropertyExperimentIP 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}$
ε110
 
8.53c 7.04 
$\epsilon _{33}^{0}$
ε330
 
9.62c 7.93 
$\epsilon _{11}^{\infty }$
ε11
 
5.23d 5.16 
$\epsilon _{33}^{\infty }$
ε33
 
5.29d 5.39 
PropertyExperimentIP 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}$
ε110
 
8.53c 7.04 
$\epsilon _{33}^{0}$
ε330
 
9.62c 7.93 
$\epsilon _{11}^{\infty }$
ε11
 
5.23d 5.16 
$\epsilon _{33}^{\infty }$
ε33
 
5.29d 5.39 
a

Reference 75.

b

Reference 76.

c

Reference 77.

d

Reference 78.

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.

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 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.

FIG. 2.

The calculated DOS as a function of energy relative to the Fermi level for a series of embedded clusters with different QM cluster sizes.

FIG. 2.

The calculated DOS as a function of energy relative to the Fermi level for a series of embedded clusters with different QM cluster sizes.

Close modal

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.

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 

FIG. 3.

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.

FIG. 3.

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.

Close modal

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.

FIG. 4.

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).

FIG. 4.

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).

Close modal

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 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

${\rm Cu}_{\mathrm{Cd}}^{\times }$
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.

FIG. 5.

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.

FIG. 5.

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.

Close modal

Forming an

${\rm In}_{\mathrm{Cd}}^{\times }$
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.

FIG. 6.

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.

FIG. 6.

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.

Close modal

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

${\rm Cu}_{\mathrm{Cd}}^{\times }$
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.

FIG. 7.

The calculated spin density associated with (a) a

${\rm Cu}_{\mathrm{Cd}}^{\times }$
Cu Cd × defect and (b) an
${\rm In}_{\mathrm{Cd}}^{\times }$
In Cd ×
defect in a CdS linear chain, determined using a plane-wave DFT approach. Larger yellow spheres represent S ions, smaller blue spheres represent Cd ions, and medium (blue and purple) spheres at the center represent the impurity ions (Cu and In, respectively).

FIG. 7.

The calculated spin density associated with (a) a

${\rm Cu}_{\mathrm{Cd}}^{\times }$
Cu Cd × defect and (b) an
${\rm In}_{\mathrm{Cd}}^{\times }$
In Cd ×
defect in a CdS linear chain, determined using a plane-wave DFT approach. Larger yellow spheres represent S ions, smaller blue spheres represent Cd ions, and medium (blue and purple) spheres at the center represent the impurity ions (Cu and In, respectively).

Close modal

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

${\rm Cu}_{\mathrm{Cd}}^{\times }$
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.

FIG. 8.

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.

FIG. 8.

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.

Close modal

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

${\rm In}_{\mathrm{Cd}}^{\times }$
In Cd × acceptor level (2.89 eV above the valence band maximum) could be compared with the 2 eV level reported for nanowires.84 Similar observations are made for Cu impurities,81 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 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

(A1)

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

(A2)

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

(A3)

Once induced, the interaction energy ΔWdip of this dipole with E is given by

(A4)

which implies that the total polarization energy ΔWpol of the line segment is given by

(A5)

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

(A6)

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

(A7)

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

(A8)

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.

1.
Y.-L.
Lee
and
Y.-S.
Lo
,
Adv. Funct. Mater.
19
,
604
(
2009
).
2.
S.
Rühle
,
M.
Shalom
, and
A.
Zaban
,
ChemPhysChem
11
,
2290
(
2010
).
3.
M.
Seol
,
H.
Kim
,
Y.
Tak
, and
K.
Yong
,
Chem. Commun.
46
,
5521
(
2010
).
4.
S.-M.
Liu
,
F.-Q.
Liu
,
H.-Q.
Guo
,
Z.-H.
Zhang
, and
Z.-G.
Wang
,
Solid State Commun.
115
,
615
(
2000
).
5.
H.
Zhao
,
E. P.
Douglas
,
B. S.
Harrison
, and
K. S.
Schanze
,
Langmuir
17
,
8428
(
2001
).
6.
W.-S.
Chae
,
J.-H.
Ko
,
I.-W.
Hwang
, and
Y.-R.
Kim
,
Chem. Phys. Lett.
365
,
49
(
2002
).
7.
J.
Zhang
,
L.
Sun
,
C.
Liao
, and
C.
Yan
,
Solid State Commun.
124
,
45
(
2002
).
8.
J. G.
Díaz
,
J.
Planelles
,
G. W.
Bryant
, and
J.
Aizpurua
,
J. Phys. Chem. B
108
,
17800
(
2004
).
9.
A.
Vorokh
and
A.
Rempel
,
Phys. Solid State
49
,
148
(
2007
).
10.
H.
Cao
,
G.
Wang
,
S.
Zhang
,
X.
Zhang
, and
D.
Rabinovich
,
Inorg. Chem.
45
,
5103
(
2006
).
11.
H. M.
Fan
,
X. F.
Fan
,
Z. H.
Ni
,
Z. X.
Shen
,
Y. P.
Feng
, and
B. S.
Zou
,
J. Phys. Chem. C
112
,
1865
(
2008
).
12.
H.
Li
,
X.
Wang
,
J.
Xu
,
Q.
Zhang
,
Y.
Bando
,
D.
Golberg
,
Y.
Ma
, and
T.
Zhai
,
Adv. Mater.
25
,
3017
(
2013
).
13.
Y.
Ye
,
Y.
Dai
,
L.
Dai
,
Z.
Shi
,
N.
Liu
,
F.
Wang
,
L.
Fu
,
R.
Peng
,
X.
Wen
,
Z.
Chen
,
Z.
Liu
, and
G.
Qin
,
ACS Appl. Mater. Interfaces
2
,
3406
(
2010
).
14.
C. J.
Barrelet
,
A. B.
Greytak
, and
C. M.
Lieber
,
Nano Lett.
4
,
1981
(
2004
).
15.
G.
Dai
,
B.
Zou
, and
Z.
Wang
,
J. Am. Chem. Soc.
132
,
12174
(
2010
).
16.
Y.
Huang
,
X.
Duan
, and
C.
Lieber
,
Small
1
,
142
(
2005
).
17.
R.-M.
Ma
,
L.
Dai
,
H.-B.
Huo
,
W.-J.
Xu
, and
G. G.
Qin
,
Nano Lett.
7
,
3300
(
2007
).
18.
X.
Duan
,
Y.
Huang
,
R.
Agarwal
, and
C. M.
Lieber
,
Nature (London)
421
,
241
(
2003
).
19.
A.
Pan
,
R.
Liu
,
Q.
Yang
,
Y.
Zhu
,
G.
Yang
,
B.
Zou
, and
K.
Chen
,
J. Phys. Chem. B
109
,
24268
(
2005
).
20.
H.
Pan
,
G.
Xing
,
Z.
Ni
,
W.
Ji
,
Y. P.
Feng
,
Z.
Tang
,
D. H. C.
Chua
,
J.
Lin
, and
Z.
Shen
,
Appl. Phys. Lett.
91
,
193105
(
2007
).
21.
T.
Zhai
,
X.
Fang
,
Y.
Bando
,
B.
Dierre
,
B.
Liu
,
H.
Zeng
,
X.
Xu
,
Y.
Huang
,
X.
Yuan
,
T.
Sekiguchi
, and
D.
Golberg
,
Adv. Funct. Mater.
19
,
2423
(
2009
).
22.
T.
Zhai
,
X.
Fang
,
Y.
Bando
,
Q.
Liao
,
X.
Xu
,
H.
Zeng
,
Y.
Ma
,
J.
Yao
, and
D.
Golberg
,
ACS Nano
3
,
949
(
2009
).
23.
L.
Li
,
P.
Wu
,
X.
Fang
,
T.
Zhai
,
L.
Dai
,
M.
Liao
,
Y.
Koide
,
H.
Wang
,
Y.
Bando
, and
D.
Golberg
,
Adv. Mater.
22
,
3161
(
2010
).
24.
J. S.
Jie
,
W. J.
Zhang
,
Y.
Jiang
, and
S. T.
Lee
,
Appl. Phys. Lett.
89
,
223117
(
2006
).
25.
J.
Zhai
,
L.
Wang
,
D.
Wang
,
H.
Li
,
Y.
Zhang
,
D. Q.
He
, and
T.
Xie
,
ACS Appl. Mater. Interfaces
3
,
2253
(
2011
).
26.
X.
Duan
,
C.
Niu
,
V.
Sahi
,
J.
Chen
,
J. W.
Parce
,
S.
Empedocles
, and
J. L.
Goldman
,
Nature (London)
425
,
274
(
2003
).
27.
Y.
Tak
,
S. J.
Hong
,
J. S.
Lee
, and
K.
Yong
,
J. Mater. Chem.
19
,
5945
(
2009
).
28.
J.
Li
and
L.-W.
Wang
,
Phys. Rev. B
72
,
125325
(
2005
).
29.
S.-P.
Huang
,
W.-D.
Cheng
,
D.-S.
Wu
,
J.-M.
Hu
,
J.
Shen
,
Z.
Xie
,
H.
Zhang
, and
Y.-J.
Gong
,
Appl. Phys. Lett.
90
,
031904
(
2007
).
30.
H.
Pan
and
Y. P.
Feng
,
ACS Nano
2
,
2410
(
2008
).
31.
W.
Sangthong
,
J.
Limtrakul
,
F.
Illas
, and
S. T.
Bromley
,
Nanoscale
2
,
72
(
2010
).
32.
S.
Karthikeyan
,
E.
Deepika
, and
P.
Murugan
,
J. Phys. Chem. C
116
,
5981
(
2012
).
33.
C. J.
Barrelet
,
Y.
Wu
,
D. C.
Bell
, and
C. M.
Lieber
,
J. Am. Chem. Soc.
125
,
11498
(
2003
).
35.
P.
Minary
,
J. A.
Morrone
,
D. A.
Yarne
,
M. E.
Tuckerman
, and
G. J.
Martyna
,
J. Chem. Phys.
121
,
11949
(
2004
).
36.
C. A.
Rozzi
,
D.
Varsano
,
A.
Marini
,
E. K. U.
Gross
, and
A.
Rubio
,
Phys. Rev. B
73
,
205119
(
2006
).
37.
A.
Warshel
and
M.
Levitt
,
J. Mol. Biol.
103
,
227
(
1976
).
38.
J. H.
Harding
,
A. H.
Harker
,
P. B.
Keegstra
,
R.
Pandey
,
J. M.
Vail
, and
C.
Woodward
,
Physica B
131
,
151
(
1985
).
39.
L.
Seijo
and
Z.
Barandiaran
,
J. Math. Chem.
10
,
41
(
1992
).
40.
M.
Sierka
and
J.
Sauer
,
Faraday Discuss.
106
,
41
(
1997
).
43.
T. A.
Wesolowski
and
A.
Warshel
,
J. Phys. Chem.
97
,
8050
(
1993
).
44.
ChemShell, a computational chemistry shell, 1999, see http://www.chemshell.org.
45.
P.
Sherwood
,
A. H.
de Vries
,
M. F.
Guest
,
G.
Schreckenbach
,
C. R. A.
Catlow
,
S. A.
French
,
A. A.
Sokol
,
S. T.
Bromley
,
W.
Thiel
,
A. J.
Turner
,
S.
Billeter
,
F.
Terstegen
,
S.
Thiel
,
J.
Kendrick
,
S. C.
Rogers
,
J.
Casci
,
M.
Watson
,
F.
King
,
E.
Karlsen
,
M.
Sjøvoll
,
A.
Fahmi
,
A.
Schäfer
, and
C.
Lennartz
,
J. Mol. Struct.: THEOCHEM
632
,
1
(
2003
).
46.
A number of alternative approaches to the derivation of the embedding potential, which might improve its description, have been reported and include wave function or charge density based techniques as already mentioned (See Refs. 37–43).
47.
M. A.
Nygren
,
L. G. M.
Pettersson
,
Z.
Barandiaran
, and
L.
Seijo
,
J. Chem. Phys.
100
,
2010
(
1994
).
48.
A. A.
Sokol
,
S. T.
Bromley
,
S. A.
French
,
C. R. A.
Catlow
, and
P.
Sherwood
,
Int. J. Quantum Chem.
99
,
695
(
2004
).
49.
A. A.
Sokol
,
S. A.
French
,
S. T.
Bromley
,
C. R. A.
Catlow
,
H. J. J.
van Dam
, and
P.
Sherwood
,
Faraday Discuss.
134
,
267
(
2007
).
50.
C. R. A.
Catlow
,
A. A.
Sokol
, and
A.
Walsh
,
Chem. Commun.
47
,
3386
(
2011
).
51.
A.
Walsh
,
C. R. A.
Catlow
,
M.
Miskufova
, and
A. A.
Sokol
,
J. Phys.: Condens. Matter
23
,
334217
(
2011
).
52.
G.
Dutta
,
A. A.
Sokol
,
C. R. A.
Catlow
,
T. W.
Keal
, and
P.
Sherwood
,
ChemPhysChem
13
,
3453
(
2012
).
53.
D. O.
Scanlon
,
C. W.
Dunnill
,
J.
Buckeridge
,
S. A.
Shevlin
,
A. J.
Logsdail
,
S. M.
Woodley
,
C. R. A.
Catlow
,
M. J.
Powell
,
R. G.
Palgrave
,
I. P.
Parkin
,
G. W.
Watson
,
T. W.
Keal
,
P.
Sherwood
,
A.
Walsh
, and
A. A.
Sokol
,
Nature Mater.
12
,
798
(
2013
).
54.
J. D.
Gale
and
A. L.
Rohl
,
Mol. Simul.
29
,
291
(
2003
).
55.
To automate the construction of CdS linear chain models with variable region sizes, we developed a shell script which reads in the total number of ions, the bond length, the size of the QM region, the size of the MM active and frozen regions, and the central ion type. A cluster is then cut with the regions divided into the appropriate sizes. In the MM region, due to the symmetry of the linear chain, atomic shells, which are used in modeling polarizable ions, are initially placed directly on cores (they may be displaced from the cores when a charged defect is present).
56.
We have developed a FORTRAN 90 program called TERMINATE to implement this procedure, which is available on request. The values of n, N, a0, and q are input, and Q and R are output.
57.
M. F.
Guest
,
I. J.
Bush
,
H. J. J.
Van Dam
,
P.
Sherwood
,
J. M. H.
Thomas
,
J. H.
Van Lenthe
,
R. W. A.
Havenith
, and
J.
Kendrick
,
Mol. Phys.
103
,
719
(
2005
).
58.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
59.
A.
Bergner
,
M.
Dolg
,
W.
Küchle
,
H.
Stoll
, and
H.
Preuss
,
Mol. Phys.
80
,
1431
(
1993
).
60.
K. A.
Peterson
,
J. Chem. Phys.
119
,
11099
(
2003
).
61.
M.
Ernzerhof
and
G. E.
Scuseria
,
J. Chem. Phys.
110
,
5029
(
1999
).
62.
C.
Adamo
and
V.
Barone
,
J. Chem. Phys.
110
,
6158
(
1999
).
63.
J. P.
Perdew
,
A.
Ruzsinszky
,
G. I.
Csonka
,
O. A.
Vydrov
,
G. E.
Scuseria
,
L. A.
Constantin
,
X.
Zhou
, and
K.
Burke
,
Phys. Rev. Lett.
100
,
136406
(
2008
).
64.
D.
Andrae
,
U.
Häussermann
,
M.
Dolg
,
H.
Stoll
, and
H.
Preuss
,
Theor. Chim. Acta
77
,
123
(
1990
).
65.
G.
Igel-Mann
, Doktorabeit (Institut für Theoretische Chemie, Stuttgart,
1987
).
66.
M.
Dolg
,
U.
Wedig
,
H.
Stoll
, and
H.
Preuss
,
J. Chem. Phys.
86
,
866
(
1987
).
67.
B.
Metz
,
H.
Stoll
, and
M.
Dolg
,
J. Chem. Phys.
113
,
2563
(
2000
).
68.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
558
(
1993
).
69.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
49
,
14251
(
1994
).
70.
G.
Kresse
and
J.
Furthmüller
,
Comput. Mater. Sci.
6
,
15
(
1996
).
71.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
72.
P. E.
Blöchl
,
Phys. Rev. B
50
,
17953
(
1994
).
73.
M.
Born
and
K.
Huang
,
Dynamical Theory of Crystal Lattices
(
Oxford University Press
,
Oxford
,
1956
).
74.
B. G.
Dick
and
A. W.
Overhauser
,
Phys. Rev.
112
,
90
(
1958
).
75.
A. W.
Stevenson
,
M.
Milanko
, and
Z.
Barnea
,
Acta Crystallogr., Sect. B
40
,
521
(
1984
).
76.
Numerical Data and Functional Relationships in Science and Technology
, edited by
O.
Madelung
,
M.
Schölz
, and
H.
Weiss
(
Springer
,
Berlin
,
1982
), Vol.
17
.
77.
I.
Kobiakov
,
Solid State Commun.
35
,
305
(
1980
).
78.
S.
Ninomiya
and
S.
Adachi
,
J. Appl. Phys.
78
,
1183
(
1995
).
79.
Alternatively, one can think in terms of the curvature of the orbitals, to which the kinetic energy is directly proportional. When the orbitals are out of phase, their curvature is greater in comparison to when they are in phase and overlap, meaning the kinetic energy is greater.
80.
T. Y.
Lui
,
J. A.
Zapien
,
H.
Tang
,
D. D. D.
Ma
,
Y. K.
Liu
,
C. S.
Lee
,
S. T.
Lee
,
S. L.
Shi
, and
S. J.
Xu
,
Nanotechnology
17
,
5935
(
2006
).
81.
V.
Ghiordanescu
,
M.
Sima
,
I.
Enculescu
,
M. N.
Grecu
,
L.
Mihut
,
M.
Secu
, and
R.
Neumann
,
Phys. Status Solidi A
202
,
449
(
2005
).
82.
K. S.
Rathore
,
D. D.
Patidar
,
N. S.
Saxena
, and
K. B.
Sharma
,
J. Ovonic Res.
5
,
175
(
2009
).
83.
I.
Shafiq
and
A.
Sharif
,
J. Phys. D: Appl. Phys.
42
,
035412
(
2009
).
84.
W.
Zhou
,
D.
Tang
, and
B.
Zou
,
J. Alloys Compd.
551
,
150
(
2013
).
85.
E.
Bertran
,
A.
Lousa
,
M.
Varela
,
M. V.
García-Cuenca
, and
J. L.
Morenza
,
Sol. Energy Mater.
17
,
55
(
1988
).
86.
R. M.
Ma
,
L.
Dai
,
H. B.
Huo
,
W. Q.
Yang
,
G. G.
Qin
,
P. H.
Tan
,
C. H.
Huang
, and
J.
Zheng
,
Appl. Phys. Lett.
89
,
203120
(
2006
).
87.
88.
D. G.
Seiler
,
D.
Heiman
,
R.
Feigenblatt
,
R. L.
Aggarwal
, and
B.
Lax
,
Phys. Rev. B
25
,
7666
(
1982
).
89.
S.
Yan
,
L.
Sun
,
P.
Qu
,
N.
Huang
,
Y.
Song
, and
Z.
Xiao
,
J. Solid State Chem.
182
,
2941
(
2009
).
90.
J. P.
Bosco
,
D. O.
Scanlon
,
G. W.
Watson
,
N. S.
Lewis
, and
H. A.
Atwater
,
J. Appl. Phys.
113
,
203705
(
2013
).
91.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Chem. Phys.
124
,
219906
(
2006
).
92.
D. M.
Eigler
and
E. K.
Schweizer
,
Nature (London)
344
,
524
(
1990
).
93.
P.
Caroff
,
K. A.
Dick
,
J.
Johansson
,
M. E.
Messing
,
K.
Deppert
, and
L.
Samuelson
,
Nat. Nanotechnol.
4
,
50
(
2009
).
94.
Computer Simulation of Solids
, edited by
C. R. A.
Catlow
and
W. C.
Mackrodt
, 1st ed. (
Springer
,
Berlin
,
1982
), Chap. 1.
95.
A.
Walsh
,
J.
Buckeridge
,
C. R. A.
Catlow
,
A. J.
Jackson
,
T. W.
Keal
,
M.
Miskufova
,
P.
Sherwood
,
S. A.
Shevlin
,
M. B.
Watkins
,
S. M.
Woodley
, and
A. A.
Sokol
,
Chem. Mater.
22
,
2924
(
2013
).