DFT-Derived Reactive Potentials for the Simulation of Activated

Article pubs.acs.org/JPCB

DFT-Derived Reactive Potentials for the Simulation of Activated Processes: the Case of CdTe and CdTe:S Xiao Liang Hu,† Riccardo Ciaglia,‡,§ Fabio Pietrucci,†,‡ Grégoire A. Gallet,†,‡ and Wanda Andreoni*,†,‡ †

Institut de Théorie des Phénomènes Physiques, Ecole Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland Centre Européen de Calcul Atomique et Moléculaire, Ecole Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland

‡

S Supporting Information *

ABSTRACT: We introduce a new ab initio derived reactive potential for the simulation of CdTe within density functional theory (DFT) and apply it to calculate both static and dynamical properties of a number of systems (bulk solid, defective structures, liquid, surfaces) at finite temperature. In particular, we also consider cases with low sulfur concentration (CdTe:S). The analysis of DFT and classical molecular dynamics (MD) simulations performed with the same protocol leads to stringent performance tests and to a detailed comparison of the two schemes. Metadynamics techniques are used to empower both Car−Parrinello and classical molecular dynamics for the simulation of activated processes. For the latter, we consider surface reconstruction and sulfur diffusion in the bulk. The same procedures are applied using previously proposed force fields for CdTe and CdTeS materials, thus allowing for a detailed comparison of the various schemes.

■

INTRODUCTION The tremendous progress in nanotechnology and the urgency for new clean energy solutions have recently led to the emergence of new compounds, architectures, and structures, thus opening unprecedented and intriguing questions about their physicochemical behavior at diverse size and time scales. In particular, the need of simulations has become apparent, enabling us to understand and predict the functioning of real materials. In spite of the recent remarkable advancement in computational materials science, the application of atomistic models to simulate physicochemical processes on the nanoscale continues to face several difficulties. Indeed, a reliable representation of a real system is highly demanding: it requires not only an adequate size of the model but also the possibility to follow its evolution over relevant times and an accurate description of interatomic interactions at relevant scales. The “time problem” has in part benefitted from enhancement of the configuration sampling for which several techniques have been developed. In this way, one can obtain an understanding of “rare events” that unbiased molecular dynamics (MD) techniques cannot capture (for a review, see ref 1). The inadequacy of classical potentials for the description of systems with either covalent or metallic bonding has long been a critical problem that is alleviated by the use of molecular dynamics based on density functional theory (DFT) (Car−Parrinello method).2 Ultimately, however, the accuracy of DFT-based MD depends on that of the DFT functional selected. Moreover, in order to overcome severe problems related to size and time scales, computations of material properties often have to rely on classical schemes using model potentials. Therefore, the construction of interatomic potentials, their validation, and in particular the control of their © XXXX American Chemical Society

transferability are still necessary steps of the simulation approach to real materials. Here we consider the case of CdTe, a semiconductor of wide application in optoelectronic devices3 and of increasing interest as a key component of thin-film-based solar cells (see, e.g., refs 4 and 5) and nanostructured photovoltaics.6,7 The performance of solar cells based on CdTe thin films, which are in general polycrystalline, is known to depend on the morphology of the interfaces, especially the one with CdS in the p−n junction (for a review, see ref 8), and also to be correlated with the nature and distribution of grain boundaries.9 However, the microscopic mechanisms leading to the observed (and often uncontrolled) effects on the functioning of a solar cell elude experimental investigation. On the other hand, the complexity and size scale of the systems involved is not affordable yet with ab initio methods and requires, at least at the initial stage, simulations entirely based on classical schemes. A few force fields have been published over the years for CdTe.10−12 In particular, two recent proposals will be considered in this paper, namely, a bond order potential (BOP)11 and a Stillinger−Weber-type (SW) potential.12 In the former, the parametrization was adjusted to reproduce a number of characteristics of bulk and cluster structures of Cd, Te, and CdTe calculated within the Perdew−Burke−Ernzerhof (PBE) DFT functional13 including empirical corrections for van der Waals interactions,14 whereas the latter had as reference Special Issue: William C. Swope Festschrift Received: December 31, 2013 Revised: March 7, 2014

A

dx.doi.org/10.1021/jp412808m | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

empirical data for solid phases of Cd, Te, S, and their compounds. In this paper, we introduce an augmented Tersoff (AT) potential15 describing the interatomic interactions in CdTe also in the presence of sulfur impurities at low concentration. We follow the AT scheme developed in ref 15 with application to SiONH compounds, and slightly modify it to allow for an increased flexibility. The augmentation scheme15 is explicitly intended to improve on the original coordination-dependent Tersoff potential16,17 in the description of mixed systems with coordination defects. The parametrization is obtained by simultaneously mapping forces and energies onto the results of DFT calculations using the PBE exchange-correlation functional.13 This latter choice is motivated by the pervasive use and global performance of PBE for the representation of condensed matter systems, like solids, surfaces, interfaces, and nanostructures. For the fitting, we follow the same procedure as in ref 15. In contrast with all previous derivations, we check the validity of our potentials not only on ground-state geometries, structural energy differences, and vibrational spectra but also on the characterization of activated processes which are of great importance for the study of complex binding situations that form, e.g., during material processing, deposition, growth, and light absorption. In particular, using metadynamics18 as a procedure for the enhancement of configurational sampling, we simulate structural transformations of the CdTe Te-terminated (001) surface and sulfur diffusion in bulk CdTe, both within DFTbased19,20 and classical MD, and compare free-energy barriers and pathways thus obtained. All calculations are repeated with the available force fields,11,12 thus allowing for a detailed comparison.

where ζijIJ =

with IJK eijk = exp[(μIJ rij − μIK rik)mI ]

∑ Vij + ∑ i≠j

+

I

∑

I tijk =1+

fijIJ

(

nI −1/2n I

)

)

dI2 + (hiI − cos(θjik))2

(7)

j≠i

(8)

Note that, in analogy with ref 15, we fix the values of mI in eq 6 and of χIJ in eq 4: mI = 1 and χIJ = 1. The cutoff radii in eq 3 were first determined approximately from pair correlation functions of the solid phases at 400 K; then, the full fitting of the other parameters was performed for a few such values. They were then fixed on the basis of the performance of the global fitting of the parametrization. Indeed, the range of physically reasonable cutoff radii is limited. In eq 1, E0I is a “core” energy meant to establish different energy scales for different elements and Eci represents the energy penalty associated with deviation from the “normal” coordination z0i : Eic = cI ,1Δzi + cI ,2Δzi 2

(9)

where Δzi =

(1)

zi =

zi − zI0 |zi − zI0|

∑ fijIJ bijIJ j≠i

fs (|zi − zI0|)

(10)

(11)

⎧0 if |z| < z T − zB ⎪ ⎪1⎡ ⎛ |z| − z T ⎞⎤ ⎪ ⎢1 + sin⎜π ⎟⎥ if zT − zB |fs (z)| = int(|z|) + ⎨ 2 ⎢⎣ zB ⎠⎥⎦ ⎝ < |z| ⎪ < z T + zB ⎪ ⎪ if zT − zB < |z| ⎩1

(2)

(12)

where zT and zB are independent from the specific atom. Fitting: the Reference Data. As reference data, we have used the results of Car−Parrinello MD at finite temperature21 (see the Supporting Information). For each system, a limited set of configurations was selected from a typical total of 10,000−20,000 (corresponding to a run of 1−2 ps). In each case, we verified that the selected configurations for the fitting were sufficiently different so as to span the relevant regions of the potential energy landscape. The reference model systems and the number of frames taken for the fitting are reported in Table 1.

(3)

with R and S denoting the cutoff radii. The damping factors bij of the two-body attractive interaction terms in eq 2 include three-body effects and depend on an “effective coordination number” βIζIJij : bijIJ = χIJ (1 + βI ζijIJ

dI2

cI2

J

where the function f IJij defines their action cutoff as ⎧1 if rij < RIJ ⎪ ⎪ ⎡ ⎛ rij − RIJ ⎞⎤ ⎪1 ⎟⎟⎥ if RIJ < rij < SIJ = ⎨ ⎢1 + cos⎜⎜π S R − ⎪ 2 ⎢⎣ ⎝ IJ IJ ⎠⎥ ⎦ ⎪ ⎪0 if SIJ < rij ⎩

−

hiI = hI ,1 + hI ,2∑ ∑ fijIJ

where I runs on the different elements and i on the NI atoms. The Vij’s are explicit functions only of the interatomic distance rij and have the form of (generalized) Morse potentials: Vij = fijIJ [AIJ e−λIJ rij − bijIJ BIJ e−μIJ rij]

cI2

The term tIijk incorporates the effect of the bending angle θjik between atoms j and k around atom i. In ref 15, hIi is a constant for each atom of the same element. Here we modify it to further differentiate atoms of the same element having a different environment:

Eic

i

(6)

and

METHODS Augmented Tersoff Potential. The expression of the augmented Tersoff potentials we use is from Billeter et al.15 and includes only one small modification that we will specify below. The potential energy E is NI EI0

(5)

k≠i,j

■

1 E= 2

∑ fikIK eijkIJK tijkI

(4) B



Article

We also recall that the fitting protocol implies that force mismatch and energy error are optimized simultaneously using a weighted sum

Table 1. Systems Considered in the Fitting of the Potential Parameters and Number of Reference Configurations Selected from the Trajectories of the DFT-PBE Molecular Dynamics structure

no. of atoms

no. of MD frames

label

bulk (400 K) vac(Cd); vac(Te) (400 K) interstitial Te (400 K) (001)-(2 × 1) (400 K) (001)-c(2 × 2) (400 K) (110) (700 K) bulk (300 K)

216 214 217 192 192 192 217

12 26 16 40 24 5 16

a b c d e f g

system CdTe

CdTe:S

E W = wEE E + (1 − wE)rBE F

where rB is the Bohr radius and wE is the weight attributed to the energy mismatch. In the fitting presented here, we found a good compromise when using wE = 0.1. MD and Metadynamics Simulations. Metadynamics simulations were used to test the AT potential (see the Supporting Information). Different choices were made for the collective variables, namely, the Te coordination number for the change of reconstruction of the (001) CdTe surface (see ref 19) and path CVs for sulfur diffusion.20,22 The use of path CVs implies a reference path as the starting approach, from which, however, the subsequent evolution of the system can (and does indeed) deviate significantly. We adopted the same reference path for all our simulations, namely, for those based on Car− Parrinello MD as well as those using the three force fields mentioned above (AT-PBE, BOP, and SW).

Bulk CdTe was considered also in the presence of defects (keeping stoichiometry) as well as slabs representing Teterminated (001) and (110) surfaces. (We shall denote the former as (001)-Te from now on.) Our models are periodically repeated cubic supercells of edge 19.92 Å for the bulk and orthorhombic slabs for the surfaces. Twelve-layer-slab models were used to represent the truncated solid with periodic boundary conditions on the (x, y) plane (4 × 4) supercell and a vacuum of 12 or 24 Å. More details of the PBE calculations can be found in ref 19 and in the Supporting Information. With Cd−Cd, Cd−Te, and Te−Te potentials fixed to CdTe, we moved on to the fitting of the interactions of a CdTe matrix with sulfur atoms incorporated at low concentration (∼1%). Our reference data were thus MD simulations of this mixed system and not of CdS. Moreover, with no S−S pairs being present, only Cd−S and Te−S interaction parameters were optimized. Following ref 15, we estimate the quality of the fitting in terms of (i) the energy error EE,Cd (renormalized to the number of Cd atoms NCd,f in each configuration f) (eq 13), (ii) the force mismatch DF (eq 14), and (iii) the correlation coefficient CF between the forces Fpot,f and the DFT values FDFT,f: E E,Cd

2

EF2 =

1 = Nfr

⎛ Epot(x f ) − E DFT(x f ) ⎞2 ⎟⎟ ∑ ⎜⎜ NCd, f ⎠ f ⎝

■

RESULTS AND DISCUSSION All parameters of our AT potentials are reported in Tables 2, 3, and 4. Note that some parameters are apparently missing for Table 2. Parameters of the AT-PBE Potential, Part I: the Diagonal Termsa

Nfr

(13)

1 N 3 ∑ f fr

a

Nat, f 3

(14) N

CF =

∑ f fr ⟨Fpot, f , FDFT, f ⟩ N

parameter

cadmium

tellurium

sulfur

AI (eV) BI (eV) λI (Å−1) μI (Å−1) RI (Å) SI (Å) βI nI mI cI dI hI,1 hI,2

1093.64 174.905 2.66580 2.49760 3.5 4.0 0.011 × 10−4 0.38039 1 3.18118 × 108 1173.48 −0.98305 0.04883

5720.46 425.826 3.00421 1.90800 3.5 4.0 0.026 × 10−4 1.60070 1 3066.41 1.88724 −0.10578 0.00996

3.0 3.5 1.135 × 10−4 1.82946 1 274.829 1.75938 −0.37490 −0.05183

As specified in the text, the term χIJ in eq 4 was fixed to 1, as in ref 15.

Table 3. Parameters of the AT-PBE Potential, Part II: the Mixed Terms

⎛ ∂ ⎞2 ∂ ⎜ E DFT(x f )⎟⎟ ∑ ∑ ∑ ⎜ Epot(x f ) − ∂xi , d ⎠ f i d ⎝ ∂xi , d Nfr Nat, f

(17)

N

[∑ f fr ⟨Fpot, f , Fpot, f ⟩]1/2 [∑ f fr ⟨FDFT, f , FDFT, f ⟩]1/2

parameter

Cd−Te

Cd−S

Te−S

AIJ (eV) BIJ (eV) λIJ (Å−1) μIJ (Å−1)

3260.26 172.748 2.84030 1.55110

2398.64 238.497 2.93200 1.76874

9208.74 243.820 3.77077 2.10216

(15)

⎡ ∂ ⎤ Fpot, f = ⎢ − Epot(x f )⎥ , ⎢⎣ ∂xi , d ⎥⎦ i,d ⎡ ∂ ⎤ FDFT, f = ⎢ − E DFT(x f )⎥ ⎢⎣ ∂xi , d ⎥⎦ i,d

sulfur: this is only because the S−S interaction was irrelevant in these systems. Here we compare results obtained within this AT scheme and within DFT-PBE for systems not included in the fitting set. This work is meant to be a severe test of the force-field transferability, which is necessary to ensure before any further application. A comparison with our own calculations employing previously proposed potentials for CdTe and CdSTe is also drawn in most cases.

(16)

where Nfr is the number of frames. The results are quoted in the Supporting Information (Table 1). C



Article

Defects: Cd Vacancies. Cd vacancies are common defects in CdTe. The set used for the fitting of energies and forces included some frames of a system (see Table 1) with one Cd vacancy and one far-away Te vacancy. We test here the performance of the potential on the optimized structure (Figure 2a) and the geometry and relative energy of the most

Table 4. Parameters of the AT-PBE Potential, Part III: the Augmentation Terma

a

parameter

cadmium

tellurium

sulfur

E0I (eV) z0I cI,1 (eV) cI,2 (eV)

−1252.37 2.83636 0.05473 1.17577

−221.309 4.28011 0.27175 0.10414

−276.844 4.23381 0.24918 0.03498

zT = 0.66865 and zB = 0.37443 for all elements.

Bulk CdTe. Unlike most approaches to the construction of force fields, we do not include either experimental or DFT equilibrium geometries in the fitting set. The calculated lattice constants and bulk moduli within all approaches here considered are reported in the Supporting Information (Table 2). In particular, the experimental value for the lattice constant (at room temperature) is 6.48 Å,23 whereas the PBE value from geometry optimization is 6.64 Å. It is not surprising that, like PBE,24,25 the AT-PBE potential predicts a lower equilibrium density (lattice constant = 6.68 Å) and also underestimates the bulk modulus. We have checked on the influence that this discrepancy might have on the prediction of structural energy differences by performing additional PBE calculations with the experimental density. The results are reported in the following sections. The vibrational spectrum constitutes a more stringent test. Figure 1 shows the experimental spectrum, measured with

Figure 2. CdTe: The environment of a Cd vacancy (black spot) in the lowest-energy structure (a) and with Te in the anti-site position ((b) and (c) in PBE and AT-PBE, respectively). Green circles denote Te atoms.

common defect associated with it: Te in the anti-site position (Figures 2b and c) (see also the Supporting Information (Figure 2) for BOP and SW). The results are reported in Table 5. The representation of these defects is a very severe test for a classical potential, given the electronic redistribution following the loss of coordination. The most rewarding result is the correct prediction of trends in the changes of bond lengths and bond angles and especially in the structural energy difference. The performance of the AT-PBE is remarkable. Calculations using the SW predict a strong effect on the energy differences quoted above on passing from the PBE to the experimental density (decrease of the lattice constant). We checked that DFT-PBE would instead predict minor effects of the change of density (see the Supporting Information, Table 3), i.e., negligible changes for both structural features and relaxation energy, and only 0.3 eV decrease in the structural energy difference. Surface Properties. CdTe surfaces are examples of systems where severe coordination changes take place. The fitting set contains a number of frames from the MD trajectories taken at 400 K for the (001)-Te surface and only very few for the (110) surface at 700 K. We can then expect good results for the 0 K geometries of the former but not necessarily for the latter, and especially no prediction can be made at this point on the performance of our force field for the simulation of structural transformations. Therefore, the results we present below constitute an important step of the evaluation. All energy changes refer to our model supercells. (110) Surface. This surface is nonpolar and as such does not tend to reconstruct. On the other hand, a strong relaxation of the topmost layers takes place. Results are in Table 6. DFTPBE predicts a raise of the Te layer 0.23 Å above its ideal position and depression of the Cd layer by 0.57 Å. This behavior is well captured by our potential and also the energy gain associated with it. On the contrary, previous potentials fail to do it, having both Cd and Te layers relaxing by the same amount and in the same direction. We have to specify again that our model has a 2D unit cell with the edge being a multiple of the PBE value at equilibrium. Therefore, for the SW potential, a comparison should also be made on a supercell of the experimental unit cell, to which ref 12 refers. We did this test and found that in this case no relaxation is predicted in the SW scheme for the atoms of the top layers relative to the ideal

Figure 1. CdTe bulk (300 K): Vibrational spectra from MD simulations, compared with experiment.

neutron scattering,26 and the power spectrum obtained from the velocity autocorrelation function calculated on MD trajectories, sampled at room temperature, with our potential and also with the force fields in refs 11 and 12. The PBE functional is known to reproduce well the vibrational properties of covalent and ionic materials apart from a (rather) global softening. This is also the case for CdTe (Figure 1), a characteristic mirrored in the AT-derived spectrum. As illustrated in Figure 1 of the Supporting Information, our PBE calculations at 0 and 300 K show that the main temperature effect is also a red-shift of the whole spectrum (by 20−40 cm−1). The BOP potential gives an even lower density than PBE (lattice constant = 6.83 Å) whereas SW is fitted to reproduce the experimental value. Both BOP and SW spectra exhibit a different behavior than DFT-PBE, and especially tend to harden the stretching modes (Figure 1). D



Article

Table 5. CdTe: The Environment of a Cd Vacancy. (a) Is the Lowest-Energy Structure (Figure 2a) with All Four underCoordinated Anions in the Lattice Sites (Te*); (b) Has One Anion in the Anti-Site (Te**) (Figure 2b and c)a method

Te*−Cd

ΔE(rel)

Te**−Te

α(Te**)

ΔE(a−b)

PBE AT-PBE BOP SW

2.81 2.86 2.91 (2.95) 2.86 (2.90)

−1.82 −1.12 −0.63 (−0.01) −0.29 (−0.82)

2.97 2.86 3.09 (3.09) 3.16 (3.31)

119 102 99 (114) 107 (102)

−2.28 −1.56 −0.76 (−0.50) −0.26 (−1.32)

a

All calculations refer to the PBE lattice constant. For BOP and SW, the values with the corresponding lattice constants are also given in parentheses. ΔE(rel) denotes the gain in energy associated with the structural relaxation relative to the ideal geometry. ΔE(a−b) gives the energy difference between the a and b configurations. Distances are in Å, angles in deg, energies in eV.

of coexistence of the two structural domains.27,28 Two freeenergy barriers characterize the transition, which, within the accuracy of our calculations in all schemes, are the same. Both AT-PBE and SW are able to reproduce this type of transition, while in the BOP scheme the c(2 × 2) reconstructed surface transforms toward a disordered state. While the periodicity of the top layer changes, the bonding environment in the underlying layers is also modified. This variation is clearly shown in Figure 3 where the bond-angle

Table 6. CdTe-(110) Surface: Relaxation Energy ΔE and Layer Displacement Δz Relative to the Truncated Bulk in the PBE Supercella method

ΔE (eV)

ΔzCd (Å)

ΔzTe (Å)

PBE AT-PBE BOP SW

−3.50 −2.89 −3.92 (−0.01) −2.37 (0.00)

−0.57 −0.30 +1.36 (0.00) −1.03 (0.00)

+0.23 +0.49 +1.34 (−0.02) −1.03 (0.00)

a

The values in parentheses refer to calculations using the lattice constant of the corresponding interatomic potential (see text).

truncated surface. We have also verified that, if we used the experimental lattice constant, PBE would have predicted a slightly larger structural relaxation and a significantly higher relaxation energy (see the Supporting Information, Table 4). (001)-Te Surface: Reconstruction. This surface exhibits two dominant reconstructions: the (2 × 1) (ground state) and the c(2 × 2). Transitions from one to the other take place at temperatures of 400−600 K, and domains with one or the other are often observed to coexist.27,28 Within the PBE approximation, we have studied these transformations using metadynamics in ref 19. We chose a coordination number (CN) of the top-layer Te atoms as a collective variable, which accounts for the formation and cleavage of Te−Te bonds. The change in the pattern of the latter on the top layer determines the change of periodicity. Geometry optimization shows that the bond length in the dimers decreases slightly from the c(2 × 2) to the (2 × 1) modification: from 2.85 to 2.83 Å in DFTPBE and from 2.90 to 2.85 Å in the AT-PBE scheme. The energy difference is reduced from 1.1 to 0.8 eV. Still, a better comparison is difficult to attain with a classical potential. Comparison also with the previously proposed potentials is reported in Table 7. The transformation from one reconstruction to the other (see ref 19 and Figure 3 in the Supporting Information) goes through an intermediate (metastable) state having the characteristics of both. This is consistent with the observations

Figure 3. CdTe (001)-Te surface. Bond-angle distribution (on the top three layers) averaged over the trajectories of MD simulations at 400 K.

(Te−Cd−Te) distribution on the top three layers is plotted in both reconstructed surfaces at 400 K. Only the AT-PBE is able to reproduce the dramatic change from one to the other that accompanies the change of periodicity. Indeed, some configurations were part of the fitting set for our potential, which certainly guarantees a certain degree of accuracy from the start, but the structural properties of the surfaces were not fitted. Liquid CdTe. As a final test of the AT-PBE potential for CdTe outside the range of densities and temperatures employed in the fitting, we considered the liquid phase of CdTe. Our model consisted of 1728 atoms and had an experimental density of 0.027 atoms/Å3.29 After heating the crystalline sample from 0 K up to 1390 K in 1 ns, we ran a MD simulation for 1 ns at 1390 K and obtained the radial distribution function g(r) illustrated in Figure 4. Comparison with the one extracted from neutron scattering data at the same temperature29,30 shows excellent agreement in the range of the first shell: for the position both of the first peak of g(r) (at 2.87 Å vs 2.85 Å30) and of its first minimum (at 3.53 Å vs 3.60 Å30). Also, for the integral of g(r) under the first peak, we obtain 4.1 using the data from ref 30 and 4.0 for our MD simulation (note that ref 30 reports a value of 4.27: this is a detail that depends

Table 7. CdTe-(001)-Te Surface: Dimer Bond Length for the Two Reconstructed Surfacesa c(2 × 2)

2×1

method

Te−Te (Å)

Te−Te (Å)

ΔE (eV)

ΔF1 (eV)

ΔF2 (eV)

PBE AT-PBE BOP SW

2.83 2.85 2.91 3.21

2.85 2.90 2.98 3.23

1.11 0.82 1.11 0.58

0.5 0.6

0.5 0.6

1.0

1.0

a Energy difference ΔE = E(c(2 × 2)) − E(2 × 1) and free energy barriers for the c(2 × 2) → c(2 × 2) + (2 × 1) transition (ΔF1) and for the c(2 × 2) + (2 × 1) → (2 × 1) transition ΔF2.

E



Article

Figure 5. CdTe:S - Sulfur in (a) interstitial and (b) substitutional positions in the presence of Te. Green and yellow circles denote Te and S, respectively. Figure 4. CdTe-liquid phase. Radial distribution function calculated at 1390 K with the AT-PBE force field compared with experiment.

on the integration procedure). The relevant result is that our potential correctly predicts that in the liquid phase CdTe retains the same coordination as in the solid. Discrepancies beyond the first coordination shell must be expected, given the localized range of the AT-PBE potential. Sulfur Impurity in CdTe. Sulfur diffusion into CdTe from CdS is an important and unavoidable process in the p−n junction of solar cells. This motivates the study of the effects induced by the presence of sulfur on the properties of the host matrix. We start here by considering a bulk system with low sulfur concentration (a single additional S over about 100 Te atoms). The fitting set contains some configurations from the DFT-MD trajectories taken over 1.8 ps of the equilibrated system at 300 K (see Table 8). First, we consider the results of geometry optimization. Two local energy minima are found within DFT-PBE, one corresponding to an interstitial position in which S forms a pair with a Te atom (Figure 5a) and the other where S replaces a Te in the lattice (apart from a certain displacement due to its smaller size) and a Te−Te pair forms (Figure 5b). As reported in Table 8, the AT-PBE potential correctly reproduces the fact that configurations with interstitial S atoms forming S−Te pairs are thermodynamically favored over those with substitutional S. This is not the case for the SW potential. Moreover, we have verified20 that the (a) structure is favored in the DFT approach also when using the hybrid exchange-correlation functional PBE0.31 Path metadynamics22 was applied to simulate the diffusion of sulfur from one equilibrium position (a) to another equivalent one. Figure 6 illustrates the path resulting from simulations using our force field at room temperature, and using the same reference path as DFT-PBE. The effect of temperature is seen in the “spreading” of the initial state especially in an ensemble of configurations. The final state appears more narrow; however, this is only due to the fact that we did not continue

Figure 6. CdTe:S - Path of the S diffusion from metadynamics using the AT-PBE potential (see text). Green and yellow circles denote positions of Te and S atoms, respectively.

the simulation for a long enough time after the landing of the impurity in this basin. We were mainly interested in the approach to it and especially in estimating the free-energy barrier of this process. The latter turns out to be 0.6 eV in DFT-PBE, 0.3 eV in AT-PBE, and 0.4 eV in the SW scheme. A good agreement can thus be claimed, given that these values have a common uncertainty of 0.2 eV.

■

CONCLUSION In this paper, we have derived a new force field for CdTe and CdTe:S systems (low S concentration), which is grounded on the Tersoff paradigm for reactive, coordination-dependent potentials and the augmentation scheme introduced in ref 15. We were especially concerned about its verification on the prediction of diverse structural patterns, in which under- and overcoordinated atoms are present, and of the dynamics of microscopic processes that are relevant to the study of

Table 8. CdTe:S - Sulfur in (a) the Interstitial and (b) the Substitutional Position in the Presence of Tea (a)

(b)

method

d(S−Te)

d(S−Cd)

α(S)

d(Te−Te)

d(Te−Cd)

α(Te)

ΔE

PBE AT-BPE SW

2.45 2.51 2.53

2.60 2.62 2.60

100 98 109

2.78 2.67 3.19

2.66 2.57 2.67

98 103 89

−0.37 −0.26 +0.08

a

Anion−anion distances refer to the pairs, as shown in Figure 5; d(Te−Cd) (α(Te)) is an average over the four bond lengths (bond angles) of the atoms in the dimer. ΔE is the energy difference between the (a) and (b) structures. Distances are in Å, angles in deg, energies in eV. F



■

interfaces between heterogeneous materials, doping and growth. Our tests consisted of the calculation and analysis of vibrational spectra, structural changes induced by the formation of defects (Cd vacancies and S impurities), energy differences between different configurations in bulk and surface systems as well as free-energy differences for structural transformation and S diffusion, and the radial distribution function in the liquid phase. We have shown that the AT-PBE force field provides indeed a description of all these very different bonding situations in good global agreement with the one used as reference. Our next step will consist of extending our model to describe also situations with higher sulfur concentration and disordered alloys. Some words of caution must be added. We remark that, as any other fitting, also ours is still limited to the type of bonding present in the reference systems. For example, the Cd−Cd potential is repulsive in the whole range we considered, which is consistent with the positive charge of Cd in II−VI compounds. Therefore, it cannot be transferred to represent the Cd metal or, e.g., the formation of Cd clusters and islands in the composite materials. Also, the AT-PBE Te−Te potential cannot be expected to represent the complex polymorphism of solid tellurium. The best way to make use of our DFT-derived potentials is in preliminary long-time runs of large-scale systems to be followed by DFT simulations of reduced-size models and in hybrid quantum-classical (QM/MM) approaches.

■

ASSOCIATED CONTENT

Here we report on the computational scheme, including assessment of the quality of the fitting, and on additional results mentioned in the main text (bulk lattice constant, temperature effects on the bulk vibrational spectra, effects of change of lattice constant on the values calculated for geometry characteristics, relaxation energies, and structural energy differences). This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION

Corresponding Author

*E-mail: wanda.andreoni@epfl.ch. Present Address §

(R.C.) Luxottica Group, Sedico Belluno, Italy.

Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Dellago, C.; Bolhuis, P. Transition Path Sampling and Other Advanced Simulation Techniques for Rare Events. Adv. Polym. Sci. 2009, 167, 221−233. (2) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471−2474. (3) Triboulet, R.; Siffert, P. CdTe and Related Compounds; Physics, Defects, Hetero- and nano-structures, Crystal Growth, Surfaces and Applications; European Materials Research Society Monographs; Elsevier: Amsterdam, The Netherlands, 2009. (4) Green, M. A.; Emery, K.; Hishikawa, Y.; Warta, W.; Dunlop, E. D. Solar Cell Efficiency Tables (version 43). Prog. Photovoltaics 2014, 22, 1−9. (5) Kranz, L.; et al. Doping of Polycrystalline CdTe for Highefficiency Solar Cells on Flexible Metal Foil. Nat. Commun. 2013, 4, 2306−2312. (6) Zhang, M.; Wang, Y.-N.; Moulin, E.; Grützmacher, D.; Chien, C.J.; Chang, P.-C.; Gao, X.; Cariusb, R.; Lu, J. G. Core-shell CdTe-TiO2 Nanostructured Solar Cell. J. Mater. Chem. 2012, 22, 10441−10443. (7) Tong, S. W.; Mishra, N.; Su, C. L.; Nalla, V.; Wu, W.; Ji, W.; Zhang, J.; Chan, Y.; Loh, K. P. High-Performance Hybrid Solar Cell Made from CdSe/CdTe Nanocrystals Supported on Reduced Graphene Oxide and PCDTBT. Adv. Funct. Mater. 2014, DOI: 10.1002/adfm.201303010. (8) Kumar, S. G.; Rao, K. S. R. K. Physics and Chemistry of CdTe/ CdS Thin Film Heterojunction Photovoltaic Devices: Fundamental and Critical Aspects. Energy Environ. Sci. 2014, 7, 45−102. (9) Taylor, A. A.; Mendis, B. G. In Transmission Electron Microscopy Characterization of Nanomaterials; Kumar, C. S., Ed.; Springer: Berlin, Heidelberg, 2014; Chapter Electron Microscopy of Thin Film Inorganic and Organic Photovoltaic Materials. (10) Wang, Z. Q.; Stroud, D.; Markworth, A. J. Monte Carlo Study of the Liquid CdTe Surface. Phys. Rev. B 1989, 40, 3129−3132. (11) Ward, D. K.; Zhou, X. W.; Wong, B. M.; Doty, F. P.; Zimmerman, J. A. Analytical Bond-order Potential for the Cadmium Telluride Binary System. Phys. Rev. B 2012, 85, 115206−115224. (12) Zhou, X. W.; Ward, D. K.; Martin, J. E.; van Swol, F. B.; CruzCampa, J. L.; Zubia, D. Stillinger-Weber potential for the II-VI elements Zn-Cd-Hg-S-Se-Te. Phys. Rev. B 2013, 88, 085309−085322. (13) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (14) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787−1799. (15) (a) Billeter, S.; Curioni, A.; Fischer, D.; Andreoni, W. Ab Initio Derived Augmented Tersoff Potential for Silicon Oxynitride Compounds and Their Interfaces with Silicon. Phys. Rev. B 2006, 73, 155329−155343;(b) Phys. Rev. B 2009, 79, 169904 (E). (16) Tersoff, J. Empirical Interatomic Potential for Silicon with Improved Elastic Properties. Phys. Rev. B 1988, 38, 9902−9905. (17) Tersoff, J. Modeling Solid-state Chemistry: Interatomic Potentials for Multicomponent Systems. Phys. Rev. B 1989, 39, 5566−5568. (18) Laio, A.; Parrinello, M. Escaping Free-energy Minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562−12566. (19) Pietrucci, F.; Gerra, G.; Andreoni, W. CdTe Surfaces: Characterizing Dynamical Processes with First-principles Metadynamics. Appl. Phys. Lett. 2010, 97, 141914−141917. (20) Gallet, G. A.; Pietrucci, F.; Andreoni, W. To be published. (21) CPMD, Copyright IBM Corp. 1990−2014, Copyright MPI für Festkörperforschung Stuttgart 1997−2001, http://www.cpmd.org. (22) Branduardi, D.; Gervasio, F. L.; Parrinello, M. From A to B in Free Energy Space. J. Chem. Phys. 2007, 126, 054103−054112. (23) Wykoff, R. W. G. Crystal Structures; Wiley: New York, 1971; Vol. 1. (24) Csonka, G. I.; Perdew, J. P.; Ruzsinszky, A.; Philipsen, P. H. T.; Lebègue, S.; Paier, J.; Vydrov, O. A.; Á ngyán, J. G. Assessing the performance of recent density functionals for bulk solids. Phys. Rev. B 2009, 79, 155107−155120.

S Supporting Information *

■

Article

ACKNOWLEDGMENTS

We thank Alessandro Curioni, Jaap Kroes, and Fabio Trani for several useful discussions during the development of this project. Computer resources were awarded by CADMOS and by the Swiss National Supercomputing Centre - CSCS (under project ID s382). The financial support for CADMOS and for the Blue Gene systems is provided by the Canton of Geneva, the Canton Vaud, the Hans Wilsdorf foundation, the LouisJeantet foundation, the University of Geneva, the University of Lausanne, and the Ecole Polytechnique Fédérale de Lausanne. G



Article

(25) Haas, P.; Tran, F.; Blaha, P.; Schwarz, K.; Laskowski, R. Insight into the performance of GGA functionals for solid-state calculations. Phys. Rev. B 2009, 80, 195109−195121. (26) Rowe, J.; Nicklow, R.; Price, D.; Zanio, K. Lattice Dynamics of Cadmium Telluride. Phys. Rev. B 1974, 10, 671−675. (27) Tatarenko, S.; Bassani, F.; Klein, J. C.; Saminadayar, K.; Cibert, J.; Etgens, V. H. Surface Reconstructions of (001) CdTe and Their Role in the Dynamics of Evaporation and Molecular Beam Epitaxy Growth. J. Vac. Sci. Technol., A 1994, 12, 140−147. (28) Tatarenko, S.; Daudin, B.; Brun, D.; Etgens, V. H.; Veron, M. B. Cd and Te Desorption from (001), (111)B, and (110) CdTe Surfaces. Phys. Rev. B 1994, 50, 18479−18488. (29) Gaspard, J.; Raty, J.; Ceolin, R.; Bellissent, R. J. Non-Cryst. Solids 1996, 205−207, 75−78. (30) Prigent, G.; Bellissent, R.; Ceolin, R.; Fischer, H.; Gaspard, J. J. Non-Cryst. Solids 1999, 250−252, 297−300. (31) Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158−6170.

H


DFT-Derived Reactive Potentials for the Simulation of Activated

Recommend Documents