{'fraction_non_alphanumeric': 0.038336663336663336, 'fraction_numerical': 0.013736263736263736, 'mean_word_length': 4.63305978898007, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The virtual crystal approximation revisited: Application to dielectric and piezoelectric properties of perovskites

We present an approach to the implementation of the virtual crystal approximation (VCA) for the study of properties of solid solutions in the context of density-functional methods. Our approach can easily be applied to any type of pseudopotential, and also has the advantage that it can be used to obtain estimates of the atomic forces that would arise if the real atoms were present, thus giving insight into the expected displacements in the real alloy. We have applied this VCA technique within the Vanderbilt ultrasoft-pseudopotential scheme to predict dielectric and piezoelectric properties of the Pb(Zr[0.5],Ti[0.5])O[3] solid solution in its paraelectric and ferroelectric phases, respectively. Comparison with calculations performed on ordered alloy supercells and with data on parents compounds demonstrates the adequacy of using the VCA for perovskite solid solutions. In particular, the VCA approach reproduces the anomalous Born effective charges and the large value of the piezoelectric coefficients.


I. INTRODUCTION
The application of first-principles electronic bandstructure methods to the study of disordered alloys and solid solutions requires some approximation for the treatment of the alloy disorder. A "direct" approach is to make use of the supercell approximation, i.e., to study one or more disordered configurations in a supercell with artificially imposed periodic boundary conditions. Such calculations generally require the use of very large supercells in order to mimic the distribution of local chemical environments, and tend to be computationally very demanding. A much simpler and computationally less expensive approach is to employ the virtual crystal approximation (VCA) [1], in which one studies a crystal with the primitive periodicity, but composed of fictitious "virtual" atoms that interpolate between the behavior of the atoms in the parent compounds. This technique has seen wide use in band-structure calculations [2][3][4][5][6][7][8][9][10][11]. Another possible approach would be to make use of the coherent potential approximation (CPA) [12], but unfortunately the CPA is generally not well suited for use in first-principles total-energy methods. A different way to go beyond the VCA is to carry out a systematic perturbation expansion in the difference between the true and VCA potentials, an approach that is sometimes referred to as "computational alchemy" [2][3][4]. However, this method is much more complicated than the usual VCA, requiring the use of density-functional linear-response techniques.
Clearly the VCA has the advantages of simplicity and computational efficiency, if two possible concerns can be addressed. First and foremost is the question of the accuracy of the VCA approximation. Previous work has demonstrated good accuracy for the VCA in some semiconductor and ferromagnetic materials [2][3][4][5][6][7], but it was found to be inadequate for an accurate treatment of the electronic structure of some unusual semiconductor systems [8][9][10]. Until the recent pioneering work of Ramer and Rappe [11], nothing was known about the ability of the VCA to describe the properties of an important class of materials, the ferroelectric perovskite solid solutions. Their work strongly suggests that these alloys are good candidates for modeling with the VCA, since it reproduces the strain-induced transitions of ordered supercells of Pb(Zr 0.5 Ti 0.5 )O 3 . However, it is not known whether the VCA is good enough to predict the anomalous dielectric and piezoelectric properties of perovskite solid solutions.
A second concern is more technical. By its nature, the VCA is closely tied to the pseudopotential approximation. Indeed, unless pseudopotentials are used, it is hopeless to apply the VCA to the usual case of isoelectronic substitution (i.e., atoms belonging to the same column but different rows of the Periodic Table). However, as pseudopotential methods have matured, it has become less obvious what is the correct or optimal way to implement the VCA. For the case of local pseudopotentials, the implementation is straightforward [8]: the potential of the virtual system made from the (A 1−x B x )C alloy is generated simply by compositionally averaging the potentials of the parent AC and BC compounds, In practice this is usually done in Fourier space by averaging V AC (G) and V BC (G). In the case of semilocal (e.g, Hamann-Schlüter-Chiang [13]) pseudopotentials, a similar averaging of the radial potentials V A,l (r) and V B,l (r) can be done separately in each angular momentum channel l. However, with the fully non-local Kleinman-Bylander type separable pseudopotentials [14] that are most commonly used in the current generation of electronic-structure calculations, the implementation of the VCA is neither straightforward nor unique. For example, Ramer and Rappe [11] discuss four different ways of implementing the VCA for such pseudopotentials, each of them providing different physical results. Similarly, the best way of applying the VCA to the case of ultrasoft pseudopotentials [15] is less obvious still. The purpose of the present work is to report progress in addressing both of the above concerns. Taking them in reverse order, we first present a first-principles VCA approach which is easily implemented for any type of pseudopotential. The method is demonstrated and tested in the context of calculations on Pb(Zr x Ti 1−x )O 3 (PZT), an important perovskite solid solution. Our approach also has the advantage that it can easily be used to obtain estimates of the atomic forces that would arise if the real atoms were present, thus giving insight into the expected displacements in the real alloy. Such information, which for example is highly relevant to many properties of ferroelectric systems, is not provided by the usual VCA techniques.
Second, we use our new approach to evaluate the quality of the VCA approximation for predicting the unusual dielectric and piezoelectric properties in ferroelectric perovskite alloys. Perovskite compounds are known to exhibit anomalous dielectric properties. For example, they display anomalously large values of the Born effective charges, resulting from hybridization between the transition-metal d and oxygen 2p orbitals [16,17]. Similarly, piezoelectric coefficients are large, compared to other classes of materials, both because of the large Born effective charges and because of the large microscopic reaction of the internal atomic coordinates to macroscopic strain [18][19][20]. Using Pb(Zr x Ti 1−x )O 3 with x=0.5 for our test system, Born effective charges and piezoelectric coefficients are calculated using our new VCA approach together with the modern theory of polarization [21,22]. We find that the VCA can be used with fair confidence to predict dielectric and piezoelectric properties of Pb(Zr 0.5 Ti 0.5 )O 3 alloys. As a matter of fact, our VCA technique yields large Born effective charges that are very nearly equal to the average between the effective charges of the parent compounds. Furthermore, comparison with calculations performed on ordered Pb(Zr 0.5 Ti 0.5 )O 3 supercells demonstrates the ability of the VCA to mimic piezoelectric coefficients of perovskite solid solutions.
The paper is organized as follows. In Sec. II, we implement our VCA approach in the context of densityfunctional theory, emphasizing the advantages of the new approach. Section III reports the predictions of this new VCA technique for the Born effective charges and piezoelectric coefficients of the Pb(Zr 0.5 Ti 0.5 )O 3 solid solution in its paraelectric and ferroelectric phases, respectively.
We conclude in Sec. IV with a discussion of perspectives and future directions. The Appendix contains details about the implementation of our VCA technique within the Vanderbilt ultrasoft-pseudopotential scheme [15].

II. THE VCA IMPLEMENTATION
Within a pseudopotential approach to density-functional theory, the total energy of N v valence electrons can be written in terms of the one-particle wavefunctions φ i as where R I is the location of the site I, and V I ps are the pseudopotentials. Here, n(r) is the electron density, E XC is the exchange and correlation energy, and U ({R I }) is the ion-ion interaction energy. A local pseudopotential takes the form V ext (r, r ′ ) = V ext (r) δ(r − r ′ ), while a non-local pseudopotential is written as a sum of projectors.
In either case, it is possible to derive a "VCA" operator equation by simply averaging the pseudopotentials of the alloyed elements on site I, where, e.g., A=Ti and B=Zr in Pb(Zr Then V ext can be written where V α ps is the pseudopotential for an atom of type α and w I α is a "weight" which specifies the statistical composition on site I. In cubic Pb(Zr x Ti 1−x )O 3 , for example, we set w = 1 for Pb and w = 0 otherwise at the cube-corner site; w = 1 for oxygen and w = 0 otherwise at the three face-center sites; and w = x and w = 1 − x respectively for Zr and Ti (and zero otherwise) at the cube-center site. In other words, we think of this crystal as composed of six atoms in the primitive cell: the usual one Pb and three oxygen atoms, and two "ghost" atoms (Zr and Ti) sharing the same lattice site (and having weights x and 1−x respectively). We then treat this unit cell containing six "atoms" in the usual first-principles pseudopotential approach, solving the Kohn-Sham equations in the presence of the potential given by Eq. (6).
This approach has the advantage of requiring only very slight modifications to the usual first-principles pseudopotential code. One simply inputs the weight w, along with the position and atom type, of each "atom" in the unit cell. These weights are then used in just a few places, e.g., in the construction of the total external potential (6). Perhaps the only subtlety is in the treatment of the Ewald energy [23] in the ion-ion interaction term U appearing in Eq. (2). Here, we clearly have to prevent any Coulomb interaction between two "ghost atoms" on the same site, or else the Ewald energy would be infinite. In fact, the Ewald energy that we calculate is that of a crystal having valence charge on site I (so that, in the case of isoelectronic substitution, this is just the usual Ewald energy). In practice, as long as the same gaussian width is used for all species when splitting the real-space and reciprocal-space Ewald contributions, this can be done very simply by replacing z α v → w α z α v inside the program and deleting the infinite onsite interaction terms that would occur in the real-space Ewald sum.
There are three definite advantages to this new VCA approach. First of all, it is extremely easy to implement, as already indicated; only the minor modifications of Eqs. (6) and (7) have to be implemented when starting from a conventional first-principles pseudopotential code. Second, in contrast to the approach of Ref. [11], there is no need to generate a pseudopotential for each virtual atom. Here, the pseudotentials are created once and for all for each true atomic species; only the weights w I change when dealing with a new composition of the solid solution. Third, the alloyed elements are still considered as separate atomic species (with corresponding weights), rather than creating a single virtual atom as a whole.
This last point is less trivial than it might seem. Because the two "ghost atoms" on a site are considered as separate "atoms," one can consider responses to the displacement of just one of these "atoms" alone. In fact, the program automatically reports the forces on all the "atoms" in the unit cell, including those on the two ghost atoms separately. For structures of low symmetry, these forces need not be the same [24]. These forces can provide some hints about the atomic distortions that would occur in the true disordered materials. The magnitudes of the force differences also act as a kind of internal diagnostic for the appropriateness of using the VCA for the system of interest. A large magnitude would suggest that the VCA approximation is not expected to be very accurate, while a smaller value implies that the VCA can be used with fair confidence to mimic structural properties of the disordered alloy under consideration.
It is also straightforward, using our approach, to use finite-difference methods to evaluate other kinds of response to the separate ghost-atom displacements. For example, the Born effective charges can be defined as the first-order polarization changes with respect to firstorder sublattice displacements. As we will see in Sec. III, the contribution of each ghost atom to the Born effective charge of the whole virtual atom can thus easily be calculated. To do the finite-difference calculation, we simply compute the change in polarization as the two ghost atoms are displaced to slightly different positions. Of course, to be meaningful, any real physical quantity (e.g., the derivative of some observable with respect to displacement) ultimately has to be evaluated at the configuration of identical ghost-atom positions.
The discussion above has been limited to normconserving pseudopotentials, but the extension to the case of Vanderbilt ultrasoft pseudopotentials [15] is fairly straightforward. This extension is discussed in the Appendix. In fact, all of our tests presented below have been carried out within the ultrasoft formulation.

III. APPLICATION TO PEROVSKITE SOLID SOLUTIONS
A. Born effective charges of the paraelectric Pb(Zr0.5Ti0.5)O3 alloy Our first goal is to determine the dynamical effective charges of the Pb(Zr 0.50 Ti 0.50 )O 3 solid solution in its paraelectric phase, as predicted by the VCA. This alloy is usually denoted as PZT. For this purpose, we first perform local-density approximation (LDA) [25] calculations within the Vanderbilt ultrasoft-pseudopotential scheme on the cubic perovskite structure, using our VCA technique. As detailed in Ref. [26], a conjugate-gradient technique is used to minimize the Kohn-Sham energy functional. The Pb 5d, Pb 6s, Pb 6p, Zr 4s, Zr 4p, Zr 4d, Zr 5s, Ti 3s, Ti 3p, Ti 3d, Ti 4s, O 2s and O 2p electrons are treated as valence electrons. A weight w of 1 is assigned to Pb and oxygen atoms on their corresponding sites, while w = 0.5 for both Ti and Zr at the cube-center site [see Eqs. (6) and (7)]. Consequently, the VCA calculation includes 44 electrons per cell. We use the Ceperley-Alder exchange and correlation [27] as parameterized by Perdew and Zunger [28]. A (6,6,6) Monkhorst-Pack mesh [29] is used in order to provide converged results [26]. The lattice parameter a 0 is fully optimized by minimizing the total energy, and is found to be very well described by Vegard's law. In other words, a 0 is very nearly equal to the compositional average between the lattice constants of pure PbTiO 3 and pure PbZrO 3 given in Ref. [26].
To mimic the paraelectric phase, the atoms are kept in the ideal cubic positions. The dynamical effective charge Z * 33,α of each real and ghost atom is then calculated by using the formula where dP z is the change in polarization along the zdirection induced by the displacements du z,α of the α atoms along the z-direction, and w α refers to the weight assigned to the atoms. We allow the two ghost atoms to be at different atomic positions in order to compute the contribution of each of them on the Born effective charge of the whole virtual atom. We follow the procedure introduced in Ref. [21] which consists in directly calculating the spontaneous polarization as a Berry phase of the Bloch states. Technically, we use roughly 660 Bloch states to assure a good convergence of the effective charges. Table I shows the Born effective charges of the different atoms, as well as the compositional average between the Z * 33 of pure paraelectric PbTiO 3 and PbZrO 3 as given in Ref. [16]. In this Table, the averaged transition-metal atom interpolating between Zr and Ti is referred to as B , and the oxygen atoms are grouped into two kinds: those denoted O 3 , located between two B atoms along the z direction; and those denoted O 1 , located between two B atoms in the perpendicular directions [19,30]. One can first see that the VCA reproduces very well the effective charges of Pb, B and O atoms found in paraelectric PZT solid solutions [31,32]: Z * 33 is around 4, −2.5, and −5.5 for Pb, O 1 and O 3 atoms, respectively, while the effective charge of the B atom is close to 6.5. In fact, Table I shows that the VCA approximation essentially averages the effective charges of the parent compounds for Pb, O 1 , O 3 and B . The VCA approach is thus able to mimic the weak and subtle interaction between the d orbitals of the transition-metal atom, treated as a single atom, and the O 2p orbitals. This interaction is responsible for the anomalous effective charges of both O 3 and B atoms [16]. On the other hand, Table I   We now apply the VCA procedure to determine the piezoelectric coefficients e ij of the ferroelectric tetragonal P4mm ground state of the Pb(Zr 0.50 Ti 0.50 )O 3 solid solution as predicted by the VCA. We use here the lattice constants minimizing the total energy of the ordered supercell of Ref. [20]; that is, the lattice parameter a 0 and the tetragonal axial ratio c/a are equal to 3.99 A and 1.0345, respectively. The atomic displacements are fully optimized by minimizing the total energy and the Hellmann-Feynman forces, the latter being smaller than 0.05 eV/Å at convergence. During these minimizations, the ghost atoms can move but always share the same position. The Hellmann-Feynman force on the virtual B atom is simply the sum of the forces on the Zr and Ti "ghost" atoms, i.e., having a weight of 0.5 in Eqs. (6) and (7). In the ferroelectric VCA ground state of PZT, the force on the Ti atom is along the polarization direction, i.e., along the z-axis. This force is exactly opposite to the force on the Zr atom. The magnitude of these forces is found to be 1.1 eV/Å, which is less than three times larger than the force used to get convergent results in Ref. [33]. This indicates that the VCA approach can be used with some confidence to describe the properties of PZT alloys. As a matter of fact, Table II shows that the VCA can reproduce remarkably well the atomic displacements leading to the appearance of ferroelectricity in the ordered Pb(Zr 0.50 Ti 0.50 )O 3 alloy [20].
Once the ferroelectric ground state is determined, the modern theory of polarization [21,22] is used to calculate the piezoelectric coefficient of PZT within our VCA procedure. More precisely, the piezoelectric coefficients e ij can be computed via [34] where Ω is the cell volume and α = 1, 2, 3 runs over the three real-space lattice vectors R α and reciprocal lattice vectors G α , and η j is the macroscopic strain. Eq. (9) has recently been derived in order to make the piezoelectric coefficients independent of the choice of branch of the Berry phase [34]. At the same time, Eq. (9) automatically eliminates of the so-called "improper" terms [19] as required to correctly predict the piezoelectric coefficients [34]. Technically, Eq. (9) is evaluated by finite differences between two strained configurations: first that of the ferroelectric ground state, and then for an additional 1% strain relative to this ground state. In the second run, the relative atomic coordinates naturally have to be reoptimized in response to the applied strain.
The "clamped-ion" or "homogeneous-strain" contribution e 33,c is given by Eq. (9) for the case i = j = 3 and evaluated at vanishing internal strain (that is, without allowing the additional relaxation of the relative atomic coordinates that would be induced by the strain). e 33,c reflects electronic effects, measuring the extent to which the Wannier centers fail to follow the homogeneous strain. The "internal-strain" part e 33,i measures just those contributions to the piezoelectric response coming from internal distortions, i.e., reflecting the extent to which the ions fail to follow the homogeneous strain. In practice, e 33 and e 33,c are computed, and e 33,i is then obtained from their difference.
The results for e 33 , e 33,c and e 33,i , as predicted by the VCA, are shown in Table III and compared with various calculations on ordered supercells. One can first notice that the VCA is able to reproduce not only the magnitudes but also the signs of the piezoelectric coefficients of ordered supercells. In particular, the clamped-ion contribution, which is negative and independent of the ordering, is very well described by the VCA approximation. Similarly, the VCA results for e 33,i and e 33 lie between those of the different ordered supercells, for which a larger spread exists. Overall, the results shown in Tables II and III confirm the adequacy of the VCA to mimic ferroelectric properties of PZT.

IV. CONCLUSIONS
In summary, we have developed a new first-principles virtual crystal approach. This method (1) is easy to implement, (2) does not require the generation of pseudopotentials for each alloy composition, and (3) its outputs, via the computation of the Hellmann-Feynman forces on the "ghost" alloyed elements, provide a hint about the ability of the VCA to mimic properties of the disordered alloys under consideration. This technique has been applied, within the Vanderbilt ultrasoft-pseudopotential scheme [15], to predict dielectric and piezoelectric properties of the Pb(Zr 0.5 Ti 0.5 )O 3 solid solution in its paraelectric and ferroelectric phase, respectively. Comparison with calculations performed on ordered supercells and with data on parent compounds demonstrates the adequacy of using the VCA for such properties, at least for isoelectronic perosvkite solid solutions. More work is needed to assess the ability of VCA to describe properties of heterovalent alloys [31,36], i.e., systems in which the alloyed elements belong to different columns of the periodic table (e.g., Pb(Sc 0.5 Nb 0.5 )O 3 or Pb(Mg The present study also strongly suggests that firstprinciples derived effective-Hamiltonian methods, already available for simple perovskite systems [37][38][39][40][41][42][43], can be used with confidence to predict finite-temperature properties of perovskite (isoelectronic) solid solutions, by modeling these alloys within the VCA approach. In this appendix, we indicate how Eq. (6) can be realized when using the Vanderbilt's ultrasoft pseudopotentials scheme [15]. In this approach, the total energy of N v valence electrons described by the monoelectronic wavefunctions φ i is given by where the external potential V ext contains a local part V ion loc and a fully nonlocal part V NL , The local part contains local ionic contributions V ion loc (r) = I V ion,I loc (|r − R I |) , while the fully nonlocal part is given by The functions β I n as well as the coefficients D where the augmentation functions Q I nm (r) are also provided by the pseudopotentials and are strictly localized in the core regions. The ultrasoft pseudopotential is thus fully determined by the functions V ion,I loc , Q I nm and β I n , and by the scalar D (0) nm,I . The algorithm used to generate these quantities is described in Refs [15,44]. The wavefunctions φ i are eigensolutions of: where S is an hermitian overlap operator given by S = 1 + nm,I q I nm |β I n β I m | with q I nm = drQ I nm (r), and where H = −∇ 2 + V eff + nm,I D I nm |β I n β I m | .
Here V eff is the screened effective local potential V eff (r) = δE tot δn(r) = V ion loc (r) + dr ′ n(r ′ ) |r − r ′ | + µ xc (r) (19) with µ xc = δE XC [n]/δn(r), and Eq. (6), which is the fundamental equation underlying our VCA approach, can thus be realized by simply replacing three ionic quantities provided by the ultrasoft pseudopotentials, namely V ion,α loc , D nm,α , and Q α nm , by their product with the corresponding atomic weight w α :