- Email: [email protected]
Total energy pseudopotential calculation in local density approximation for BPO4
THEO CHEM Journal of Molecular Structure (Theochem) 426 (I 998) 24 I-247
ELSEVIER
Total energy pseudopotential calculation in local density approximation for BP04 V. Louis-Achille*, Commissariat
L. De Windt, M. Defranceschi
ci I’ Energie Atomique. IPSN/DPRE/SERGD.
BP&92265
Fontenav-aux-Roses
Cedex. France
Received 15October 1996; accepted 18 June 1997
Abstract We use the density functional pseudopotential method to describe successfully the ground state structural and electronic structure of solids as well as their elastic properties. Employing a primitive cell with periodic boundary conditions, and full relaxation of all internal coordinates, and using BP04 as a test case, we show that the above-mentioned can be extrapolated to
more sophisticated condensed matter systems. 0 1998 Elsevier Science B.V. Kqvwords:
Local density functional; Pseudopotential; Mineral phosphate; Monazite
1. Introduction Attempts to model the electronic structure of condensed matter systems have appeared in the literature for tens of years, but it is only recently that they have reached such a level that they can give predictive information accurate enough to be a useful tool for materials engineering. One can easily imagine that dealing with solids rapidly becomes an insuperable task as the size of the system increases if a bruteforce all-electron method is used. Consequently, all attempts to develop condensed matter modeling are based on approximations in order to reduce the computational effort. Two main categories of methods have been developed so far: those related to the so-called tight binding method (which can be related to semi-empirical methods) and those considered as ab initio methods because there is no need to calibrate the results with experimental data. Among * Corresponding
author.
the ab initio methods the density functional theory experiences a considerable amount of success in predicting ground state properties and gives the opportunity to extend the applicability of calculations of electronic structure of solids. Presently, the calculation of ground state properties for complex systems is limited to the availability of computer time and memory. Therefore, in order to reduce the computational time, a pseudopotential approximation is used to describe the core region. The selection of the pseudopotential is as crucial as the selection of the basis set in molecular quantum chemistry calculations. Nowadays, numerous codes are available for dealing with such problems; however, the use of the density functional pseudopotential approach is not straightforward: it necessitates a good skill. It is the purpose of the present paper to present a strategy to calculate structural energetical, and mechanical properties (bulk modulus, stress tensor components) using this type of method on BPOI as a test case. This is done without loss of generality, since BP04 is a simple phosphate
0166-1280/98/$19.00 0 1998 Elsevier Science B.V. All rights reserved PII SO166-
1280(97)00325-4
242
V. Louis-Achille
et al./Journal of Molecular Structure (Theochem) 426 (1998) 241-247
structure composed of light atoms and with a small unit cell, which is structurally related to monazites and apatites which are the structures we are interested in. In another paper of the same issue, the density functional pseudopotential approach is applied to two orthophosphates &PO4 and YP04 (structurally related to the monazites) and to the fluoroapatite Ca
IO(P04)6F2.
In a first step we present the theoretical techniques we used for the determination of the total energy and we detail the corrections used for calculating this energy and the stress tensor. Then, we compare two different procedures of geometry optimization which, starting with a trial structure, yield to the most stable structure at 0 K. Then, in Section 4, from the dependence of the total energy on lattice parameters just obtained, we determine the cohesive energy and the bulk modulus of the solids.
2. Methodology Our calculations are based on density functional theory and pseudopotentials performed in periodic boundary conditions with the occupied valence orbitals represented by a plane wave expansion. The calculations performed throughout the paper have been obtained using the Plane-Wave 3.0.0 code [l]. Density functional theory and the local density approximation (LDA), using the Ceperley and Alder potential as parametrized by Teter [2], are employed to model the electronic interactions. The minimization of the Kohn-Sham energy functional [3,4] with respect to the volume is realized with a conjugated gradient technique [5]. The advantage of such a method is that the full Hamiltonian matrix is not calculated, leading to a reduction of storage and an increase in speed of computation. One of the difficulties of applying LDA calculations to minerals comes from the computational resources needed. Typically, the unit cells are very large, since they contain numerous independent atoms, and the calculations require an enormous set of basis functions to solve the Kohr-Sham equations. In the present approach the problem is overcome by the use of pseudopotentials representing the nuclei and the core electrons: the original solid is replaced by pseudovalence electrons and a pseudoion core.
Consequently, a first step in a calculation of the pseudopotential. 2. I. Description
is the choice
of pseudoatoms
As direct expansions of the all-electron orbitals are very tightly bound to the nucleus, and the valence states require an enormous set of basis functions, and because the core orbitals are oscillatory functions of the atom, the use of a pseudopotential allows one to replace a deep all-electron potential with a much softer pseudopotential. This removes the chemically inert core electrons from the calculation and replaces the all-electron valence functions with smooth pseudowave functions that can be expanded with a moderate basis of plane-waves. The pseudopotential used in this work is a smooth norm-conserving pseudopotential generated using the procedure of Troullier and Martins [6]: it includes a non-local contribution expressed via the Kleinman-Bylander foundation [7]. Following Bloch’s theorem for a periodic solid, plane waves are used as basis functions to model the electron density outside the core region. Each electronic wave function is written as
where {G} are the reciprocal lattice vectors of the system, each wave function is labeled with a band index n and a Bloch wave vector k. The power of Bloch’s theorem is that it reduces the number of basis functions that must be included in the calculation of wave functions from a continuous to a discrete set. The electronic wave functions at a given k point that are very close together are almost identical, then the electronic wave function over a region of k space is represented by the wave function at a single k point. We used the method of Monkhorst and Pack to obtain a very accurate approximation to the electronic potential and to the contribution to the total energy from the filled band by calculating the electronic states at special k points in the Brillouin zone [8]. We tested that if we change the number of k points from three to six the energy changes about 10m3%. The coefficient for the plane waves with a small kinetic energy are typically more important, than those with a large kinetic energy. Hence, the plane
lC Louis-Achille
E -2056
et allJournal
‘SIP)
4
Fig. 1, Total energy evolution with respect to cut-off energy variations for three k points.
plane waves of highest values of k change the electronic wavefunctions only within the ion cores and not in the bonding regions. However, running energy calculations with a reduced cut-off leads to a ragged result in the energy curve, and here again a good compromise is to be looked for. Finally, it is worth stressing that one of the drawbacks of the use of plane wave bases is that, in performing a minimization with respect to the cell parameters, the number of plane waves involved in the calculation changes discontinuously with the size of the unit cell, thereby introducing discontinuities in the calculated energy and leading to spurious oscillations in the curve. In order to reduce these errors two types of correction have been considered: we have introduced energy and stress corrections which smooth the discontinuities due to the NP” change with the cell size and by the use of a reduced cut-off energy. They are detailed now. 2.2. Correction
wave basis set is truncated to include only plane waves that have kinetic energy of less than some particular cut-off energy E,,, [9] according to the equation (in a.u.) ilk + Cl2 < Ecu,
243
of Molecular Structure (Theochem) 426 (1998) 241-247
(2)
The number of plane waves denoted NP” depends on the chosen cut-off energy and on the size of the primitive cell. We studied the evolution of the energy with respect to the cut-off energy keeping the number of special k points fixed. Fig. 1 presents the result for three special k points. One can see on the curve, that the energy is well converged for a cut-off of 870 eV. The dilemma for choosing the computational parameters is that the larger the E,,, is, then the more accurate are the desired results and more memory and CPU time are needed. Thus, to find good parameters one has to look for a good compromise between accuracy and computational time and memory. In particular, the value of E,,, depends strongly on the choice of the pseudopotential. For instance, in another paper of this volume, the use of another pseudopotential necessitates an E,,, value of 1600 eV. However, it is worth pointing out that the convergence of energy differences occurs at a much smaller cut-off energy than that required for absolute convergence. This is easily understood if one remembers that the
to energy
According to the cut-off definition, the number of plane waves deduced from Eq. (2) increases with the cell volume. For each increase (decrease) of the cell volume, the G values in the reciprocal space decrease (increase) accordingly for a given cut-off, the number of plane waves defined by k, changes discontinuously. A new set (for each G vector) of (k + G) vectors might be added (removed) to the basis set, the total energy discontinuously decreases (increases) due to the added (removed) variational freedom in the wave function. This effect is more important when the cut-off energy is low. The energy correction factor accounts for the difference between the number of states in a basis set with infinitely large number of k points and the number of basis set used in the calculation [lo]. Then (in a.u.) E corrected
-E -
calculated
-
%caicuiatedNP” a WpW)
v,,llg,
(3)
and (4) where Np” is the geometric average of the number of plane waves for each k point calculated by the code, and V,,,i is the unit cell’s volume. The product V,,,,g,
244
V Louis-Achille
et al./Journal ofMolecular
corresponds to an ideal number of plane waves which takes no account of the discrete nature of the number of plane waves. To calculate the energy derivative, we performed three calculations at a fixed geometry including the lattice constants close to the expected value of final geometry; here it was the experimental geometry. We performed a first calculation with the chosen cut-off energy (870 eV), then two others with the energy cut-off increased and decreased by 3%. For each point the geometric average number of plane waves is different, and the total energy decreases with the increase of i@‘“. From this point, we calculate the total energy derivative vs. log Np”. The corrected energy is more physical than the original energy of the material unless a very high E,,, is used in the calculation. In addition, the energy difference in the corrected energy will always converge at a much smaller cut-off energy than that required for energy differences in the uncorrected energy. A possible way out is to work with a fixed number of plane waves; however, it is generally observed that choosing a constant cut-off energy leads to smaller errors in static equilibrium [ Ill.
Structure (Theochem) 426 (1998) 241-247
Fig. 2. The BPO, crystal structure
properties much closer to the experimental values than is otherwise possible at this cut-off energy. We can, therefore, perform accurate total energy calculations in order to determine total energy differences and physical properties using significantly smaller basis sets than is possible without applying the correction, allowing substantial savings in memory and computer time.
3. Optimization 2.3. Correction
to the stress tensor
We calculated the macroscopic stress tensor components. All the non-diagonal components are close to zero. The diagonal components are written as
cell(i)dJ%lculated uii = y V
procedures
d cell(i)
(5)
where i = x,y,z, cell(i) is an array containing the cell parameters and V is the volume of the cell. When the stress tensor is determined for different volumes, a problem similar to the one encountered for the energy calculation is met. We can then distinguish the calculated stress corresponding to a locally constant NP” and the corrected stress corresponding to a strictly constant E,,, [ Ill. This difference is defined as the Pulay stress. Thus
(6) The last term is obtained as above from three calculations of total energy with different Np” at a fixed volume. The inclusion of this correction yields physical
The BP04 structure crystallizes in the tetragonaltype structure according to the tetragonal space group 4,. The experimental cell parameters for a and c parameters are respectively 4.332 A and 6.640 A. In the tetragonal system, b is equal to c and the three angles are equal to 90”. The primitive cell holds six atoms, that is one BP04 formula [12]. BP04 is an ionocovalent-type solid. Fig. 2 represents the tetragonal multiple cell. 3.1. Optimization
with energy minimization
In order to obtain the equilibrium cell parameters, two parameters are sufficient to vary: the cell sides a and c, or equivalently the cell volume, and the ratio c/u. The main characteristic of our approach is that instead of optimizing the geometry only with energy minimization criterion, we mix two methods: the first one minimizes the total energy and the second one looks at the stress tensor for the determination of the c/a ratio. This allows us to save computational time without losing accuracy.
245
li. Louis-Achille et al/Journal of Molecular Structure (Theochem) 426 (1998) 241-247 E
ieV)
Table I Calculated
and experimental
cell parameters
-2ssa.ab t
2OSS.bO
7
“.,
-2093.aa
4
‘\;
.
Calculated Experimental
a (I\)
(’ (A,
cla
4.329 4.332
6.659 6.640
1.538 1.533
In the second step, the volume is kept equal to the value corresponding to the above-determined minimum and c/a is varied about 1.5% around the experimental value and the anisotropic part of the stress tensor is studied. As BP04 is a tetragonal structure, the anisotropic part is reduced to (u,~,+,) where z stands for the c direction and x for the a direction. The minimum of the total energy corresponds to the volume where the stress tensor is isotropic, that is the anisotropic part is zero [ 131. By a linear extrapolation we calculated this volume, then we deduced the c/a ratio. In Table 1 are reported the cell parameters and the corresponding c/a ratio from calculation after the correction and from experiment [ 121. The difference between the experimental and the calculated values is small, 0.07% for a, 0.3% for c and 0.3% for the c/a ratio. These results follow the general trends that LDA calculations give good equilibrium structures due to cancellation of errors. 3.2. Optimization
(b)
Fig. 3. Total energy evolution with respect to isotropic variations of the cell volume: (a) without the correction factor; (b) with the correction factor.
In the first step, the da ratio is kept constant and the parameters a and c are varied isotropically. For each volume the total energy is determined and the atomic coordinates are relaxed using the Hellman-Feynmann forces as a guide; the correction defined by Eq. (3) is applied on top of the minimization. In Fig. 3(a) and Fig. 3(b) are reported the curves of the evolution of the energy vs. the cell volume with and without energy correction respectively. The effect of the correction is twofold: it smooths out the curve and it shifts the minimum of the curve toward the origin.
with pressure
calculation
Another method allows us to calculate the cell parameters of the equilibrium. As mentioned previously, a system is considered as being in an equilibrium state if the anisotropic part of the stress tensor is zero. Thus, we can calculate the external pressure P needed for maintaining the crystal at the considered non-equilibrium volume which, for the tetragonal structure of BP04, is written as follows: P = - $(c& + ax,, + a,)
(7)
where uXX) u I)”u :Z are the diagonal components of the stress tensor [l 11. At the equilibrium state, since the anisotropic part of the stress is zero, uXXis equal to uZZ.In addition, the summation of uXX,uVV and cr, is equal to zero, that is the external pressure is zero. Using this method we realized only the first step of the last method. For each volume (the same as the previous one used for the energy minimization) we calculated the stress tensor
246
V. Louis-Achille
et d/Journal
of Molecular Structure (Theochem) 426 (1998) 241-247
IJ(GPaj
Table 2 Calculated
k
and experimental
Calculated Experimental
11‘I
117
I20
123
126
129
132
13s
we 64
cell parameters
a (‘Q
c (A)
cia
4.332 4.332
6.623 6.640
1.538 1.533
the c/a ratio. The results are given in Table 2. They are in good agreement with experiments, meaning that the method is also well adapted for the calculation of the structural parameters of equilibrium. Furthermore, results obtained with total energy and stress optimization procedures are consistent, meaning that the optimization and correction schemes have been performed properly. For all these reasons, we think that our approach is suitable for the geometry calculation of solids.
P( GPn)
4. Subsequent
f\
physical properties
Obtaining the equilibrium geometry with accuracy is not a purpose in itself. We are more particularly concerned with the calculation of physical properties (cohesive energy, bulk modulus) which can be derived from the dependence of the total energy on structural parameters. 4. I. Cohesive energy I 114
I
/
117
120
-7----7---j-4 123
126
129
132
135
VA')
@I
Fig. 4. Pressure evolution with respect to isotropic variations of the cell volume: (a) without the correction factor; (b) with the correction factor.
and then the pressure. The point where the pressure is zero gives the equilibrium volume. For each calculation, we apply the correction on the stress. Fig. 4(a) and Fig. 4(b) show the resulting curves without and with correction. The shape of the curve does not vary with the correction but, contrary to what happens for the energy, the zero is shifted towards greater values of the volume. Then, we use the above result according to the cancellation of the anisotropic part of the stress on
The cohesive energy of a solid is defined as the difference between the crystal total energy and the energy of the isolated atoms at zero temperature. Both the pseudopotential energies of the equilibrium crystal structures and of the isolated atoms were calculated using Plane-Wave. Plane-Wave atomic calculations are performed putting each atom in the center of a cubic box with a Pzzl periodicity in order to keep a spherical symmetry. We choose the ridge large enough so that the interactions between adjacent atoms are negligible. The calculated energy includes a spin-polarization correction. Thermal and vibrational contributions were not considered in the calculations of cohesive energy, since they are expected to be much smaller than the inherent errors from LDA. The experimental value of the cohesive energy is calculated with the reticular energy of BPOc the ionization potential of B and P and the electronic affinity of 0. The
V. Louis-Achille et al./Journal of Molecular Structure (TheochemJ 426 (I 998) 241-247
reticular energy is the difference between the crystal energy and the energy of the isolated gaseous ions. The experimental value of the cohesive energy is 47.6 eV/ BP04 [ 141. The calculated value is 5 1.5 eV/BPO+ We find a good agreement between the calculation and the experiment. An overestimation is usual within the LDA approximation. 4.2. Bulk modulus The bulk modulus is a measure of the resistance to an external pressure. For a crystal of volume V in an isotropic environment where the pressure has a given value, if the pressure is changed by an amount AI’, the volume of the crystal changes by an amount of AV. The isothermal bulk modulus B of the crystal is defined in a manner analogous to the definition of Young’s modulus and B = - V(AP/A V) for constant temperature. In the limiting case, where the pressure change has an infinitesimal value dP, B becomes -V(dP/dV). An estimate of the zero pressure isothermal bulk modulus and its pressure derivative was obtained using the total energy calculation technique. The values for the total energy as a function of the cell volume were fitted, using a least squares technique to Mumaghan’s equation of state [ 151:
1
E(V)= VB (v’vo)B” + 1 +const &I [ B;-1
(8)
where BO, and Bb are the bulk modulus and its pressure derivative at the equilibrium volume VO. The calculated bulk modulus is 28 GPa and its pressure derivative is 5. No experimental data of bulk modulus are available in the literature at present. The prediction of the bulk modulus is generally more difficult than that of lattice constants, since elastic properties need more accuracy. It is also very difficult to evaluate the error of these estimates. However, the uncertainty of B is related to those of the cell parameters and to the derivative Bb [ 161 and an uncertainty of 0.3% in lattice parameters and a Bb value of 5 leads to an uncertainty of 6% on BO.
247
simple crystal structure contains six atoms in the primitive cell including only light atoms. The method proposed in this paper has been shown to provide valuable results on the ground state structure of solids on the one hand, as well as on their macroscopic properties (e.g. cohesive energy, elastic properties) in this test case on the other hand. The introduction of energy and stress corrections allowed us to run with a reduced cut-off, thereby saving computational time. The computed parameters are in good agreement with the experimental data for both geometry optimization methods. Lattice parameters are within 0.3% of the experimental structure and cohesive energy is about 8%. Consequently, the bulk modulus is estimated within 6%. Moreover, the results obtained with total energy and stress optimization procedures are consistent, meaning that the optimization and correction schemes have been performed properly. The underlying strategy to calculate all these quantities has been developed, and the numerical tests performed on BP04 are promising and open up the prospects for future extension to more complex systems. Results along these lines are presented in another paper for other iono-covalent-type phosphates, but the strategy can also be extended to other types of mineral.
References [l] Plane-Wave 3.0.0, M.S.I., San Diego, 1994. [2] D.M. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (1980) 566. M.P. Teter, unpublished results, 1990. [3] P. Hohenberg, W. Kohn, Phys. Rev. B. 136 (1964) 864. [4] W. Kohn, L.J. Sham, Phys. Rev. A. 140 (1965) 1133. [5] M.P. Teter, M.C. Payne, D.C. Allan, Phys. Rev. B. 40 (1989) 12255. [6] [7] [8] [9] [IO]
[I l] [ 121
5. Conclusion
[13] [ 141
A local density functional pseudopotential approach was applied to the iono-covalent solid BPO?. Its
[15] [I61
N. Troullier, J.L. Martins, Phys. Rev. B. 43 ( 1991) 1993. L. Kleinman, P.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. H.J. Monkhorst, J.D. Pack, Phys. Rev. B. 13 (1976) 5188. M.C. Payne. M.P. Teter, D.C. Allan, T.A. Arias, J.D. Joannopopoulos, Rev. Mod. Phys. 64 (1992) 1045. G.P. Francis, M.C. Payne, J. Phys.: Condens. Matter 2 ( 1990) 4395. P. Dacosta, O.H. Nielsen, K. Kunc, J. Phys. C: 19 (1986) 3 163. R.W.G. Wychkoff, Crystal Structure, vol. 4, Interscience. New York, 1968. R.J. Needs, R.M. Martin, Phys. Rev. B 43 (1991) 1993. Handbook of Physical Constants, The Geological Society of America, 1996. F.D. Mumaghan, Proc. Natl. Acad. Sci. USA 30 (1994) 422. O.H. Neilsen, R.M. Martin, Phys. Rev. B 32 (1985) 3792.
ELSEVIER
Total energy pseudopotential calculation in local density approximation for BP04 V. Louis-Achille*, Commissariat
L. De Windt, M. Defranceschi
ci I’ Energie Atomique. IPSN/DPRE/SERGD.
BP&92265
Fontenav-aux-Roses
Cedex. France
Received 15October 1996; accepted 18 June 1997
Abstract We use the density functional pseudopotential method to describe successfully the ground state structural and electronic structure of solids as well as their elastic properties. Employing a primitive cell with periodic boundary conditions, and full relaxation of all internal coordinates, and using BP04 as a test case, we show that the above-mentioned can be extrapolated to
more sophisticated condensed matter systems. 0 1998 Elsevier Science B.V. Kqvwords:
Local density functional; Pseudopotential; Mineral phosphate; Monazite
1. Introduction Attempts to model the electronic structure of condensed matter systems have appeared in the literature for tens of years, but it is only recently that they have reached such a level that they can give predictive information accurate enough to be a useful tool for materials engineering. One can easily imagine that dealing with solids rapidly becomes an insuperable task as the size of the system increases if a bruteforce all-electron method is used. Consequently, all attempts to develop condensed matter modeling are based on approximations in order to reduce the computational effort. Two main categories of methods have been developed so far: those related to the so-called tight binding method (which can be related to semi-empirical methods) and those considered as ab initio methods because there is no need to calibrate the results with experimental data. Among * Corresponding
author.
the ab initio methods the density functional theory experiences a considerable amount of success in predicting ground state properties and gives the opportunity to extend the applicability of calculations of electronic structure of solids. Presently, the calculation of ground state properties for complex systems is limited to the availability of computer time and memory. Therefore, in order to reduce the computational time, a pseudopotential approximation is used to describe the core region. The selection of the pseudopotential is as crucial as the selection of the basis set in molecular quantum chemistry calculations. Nowadays, numerous codes are available for dealing with such problems; however, the use of the density functional pseudopotential approach is not straightforward: it necessitates a good skill. It is the purpose of the present paper to present a strategy to calculate structural energetical, and mechanical properties (bulk modulus, stress tensor components) using this type of method on BPOI as a test case. This is done without loss of generality, since BP04 is a simple phosphate
0166-1280/98/$19.00 0 1998 Elsevier Science B.V. All rights reserved PII SO166-
1280(97)00325-4
242
V. Louis-Achille
et al./Journal of Molecular Structure (Theochem) 426 (1998) 241-247
structure composed of light atoms and with a small unit cell, which is structurally related to monazites and apatites which are the structures we are interested in. In another paper of the same issue, the density functional pseudopotential approach is applied to two orthophosphates &PO4 and YP04 (structurally related to the monazites) and to the fluoroapatite Ca
IO(P04)6F2.
In a first step we present the theoretical techniques we used for the determination of the total energy and we detail the corrections used for calculating this energy and the stress tensor. Then, we compare two different procedures of geometry optimization which, starting with a trial structure, yield to the most stable structure at 0 K. Then, in Section 4, from the dependence of the total energy on lattice parameters just obtained, we determine the cohesive energy and the bulk modulus of the solids.
2. Methodology Our calculations are based on density functional theory and pseudopotentials performed in periodic boundary conditions with the occupied valence orbitals represented by a plane wave expansion. The calculations performed throughout the paper have been obtained using the Plane-Wave 3.0.0 code [l]. Density functional theory and the local density approximation (LDA), using the Ceperley and Alder potential as parametrized by Teter [2], are employed to model the electronic interactions. The minimization of the Kohn-Sham energy functional [3,4] with respect to the volume is realized with a conjugated gradient technique [5]. The advantage of such a method is that the full Hamiltonian matrix is not calculated, leading to a reduction of storage and an increase in speed of computation. One of the difficulties of applying LDA calculations to minerals comes from the computational resources needed. Typically, the unit cells are very large, since they contain numerous independent atoms, and the calculations require an enormous set of basis functions to solve the Kohr-Sham equations. In the present approach the problem is overcome by the use of pseudopotentials representing the nuclei and the core electrons: the original solid is replaced by pseudovalence electrons and a pseudoion core.
Consequently, a first step in a calculation of the pseudopotential. 2. I. Description
is the choice
of pseudoatoms
As direct expansions of the all-electron orbitals are very tightly bound to the nucleus, and the valence states require an enormous set of basis functions, and because the core orbitals are oscillatory functions of the atom, the use of a pseudopotential allows one to replace a deep all-electron potential with a much softer pseudopotential. This removes the chemically inert core electrons from the calculation and replaces the all-electron valence functions with smooth pseudowave functions that can be expanded with a moderate basis of plane-waves. The pseudopotential used in this work is a smooth norm-conserving pseudopotential generated using the procedure of Troullier and Martins [6]: it includes a non-local contribution expressed via the Kleinman-Bylander foundation [7]. Following Bloch’s theorem for a periodic solid, plane waves are used as basis functions to model the electron density outside the core region. Each electronic wave function is written as
where {G} are the reciprocal lattice vectors of the system, each wave function is labeled with a band index n and a Bloch wave vector k. The power of Bloch’s theorem is that it reduces the number of basis functions that must be included in the calculation of wave functions from a continuous to a discrete set. The electronic wave functions at a given k point that are very close together are almost identical, then the electronic wave function over a region of k space is represented by the wave function at a single k point. We used the method of Monkhorst and Pack to obtain a very accurate approximation to the electronic potential and to the contribution to the total energy from the filled band by calculating the electronic states at special k points in the Brillouin zone [8]. We tested that if we change the number of k points from three to six the energy changes about 10m3%. The coefficient for the plane waves with a small kinetic energy are typically more important, than those with a large kinetic energy. Hence, the plane
lC Louis-Achille
E -2056
et allJournal
‘SIP)
4
Fig. 1, Total energy evolution with respect to cut-off energy variations for three k points.
plane waves of highest values of k change the electronic wavefunctions only within the ion cores and not in the bonding regions. However, running energy calculations with a reduced cut-off leads to a ragged result in the energy curve, and here again a good compromise is to be looked for. Finally, it is worth stressing that one of the drawbacks of the use of plane wave bases is that, in performing a minimization with respect to the cell parameters, the number of plane waves involved in the calculation changes discontinuously with the size of the unit cell, thereby introducing discontinuities in the calculated energy and leading to spurious oscillations in the curve. In order to reduce these errors two types of correction have been considered: we have introduced energy and stress corrections which smooth the discontinuities due to the NP” change with the cell size and by the use of a reduced cut-off energy. They are detailed now. 2.2. Correction
wave basis set is truncated to include only plane waves that have kinetic energy of less than some particular cut-off energy E,,, [9] according to the equation (in a.u.) ilk + Cl2 < Ecu,
243
of Molecular Structure (Theochem) 426 (1998) 241-247
(2)
The number of plane waves denoted NP” depends on the chosen cut-off energy and on the size of the primitive cell. We studied the evolution of the energy with respect to the cut-off energy keeping the number of special k points fixed. Fig. 1 presents the result for three special k points. One can see on the curve, that the energy is well converged for a cut-off of 870 eV. The dilemma for choosing the computational parameters is that the larger the E,,, is, then the more accurate are the desired results and more memory and CPU time are needed. Thus, to find good parameters one has to look for a good compromise between accuracy and computational time and memory. In particular, the value of E,,, depends strongly on the choice of the pseudopotential. For instance, in another paper of this volume, the use of another pseudopotential necessitates an E,,, value of 1600 eV. However, it is worth pointing out that the convergence of energy differences occurs at a much smaller cut-off energy than that required for absolute convergence. This is easily understood if one remembers that the
to energy
According to the cut-off definition, the number of plane waves deduced from Eq. (2) increases with the cell volume. For each increase (decrease) of the cell volume, the G values in the reciprocal space decrease (increase) accordingly for a given cut-off, the number of plane waves defined by k, changes discontinuously. A new set (for each G vector) of (k + G) vectors might be added (removed) to the basis set, the total energy discontinuously decreases (increases) due to the added (removed) variational freedom in the wave function. This effect is more important when the cut-off energy is low. The energy correction factor accounts for the difference between the number of states in a basis set with infinitely large number of k points and the number of basis set used in the calculation [lo]. Then (in a.u.) E corrected
-E -
calculated
-
%caicuiatedNP” a WpW)
v,,llg,
(3)
and (4) where Np” is the geometric average of the number of plane waves for each k point calculated by the code, and V,,,i is the unit cell’s volume. The product V,,,,g,
244
V Louis-Achille
et al./Journal ofMolecular
corresponds to an ideal number of plane waves which takes no account of the discrete nature of the number of plane waves. To calculate the energy derivative, we performed three calculations at a fixed geometry including the lattice constants close to the expected value of final geometry; here it was the experimental geometry. We performed a first calculation with the chosen cut-off energy (870 eV), then two others with the energy cut-off increased and decreased by 3%. For each point the geometric average number of plane waves is different, and the total energy decreases with the increase of i@‘“. From this point, we calculate the total energy derivative vs. log Np”. The corrected energy is more physical than the original energy of the material unless a very high E,,, is used in the calculation. In addition, the energy difference in the corrected energy will always converge at a much smaller cut-off energy than that required for energy differences in the uncorrected energy. A possible way out is to work with a fixed number of plane waves; however, it is generally observed that choosing a constant cut-off energy leads to smaller errors in static equilibrium [ Ill.
Structure (Theochem) 426 (1998) 241-247
Fig. 2. The BPO, crystal structure
properties much closer to the experimental values than is otherwise possible at this cut-off energy. We can, therefore, perform accurate total energy calculations in order to determine total energy differences and physical properties using significantly smaller basis sets than is possible without applying the correction, allowing substantial savings in memory and computer time.
3. Optimization 2.3. Correction
to the stress tensor
We calculated the macroscopic stress tensor components. All the non-diagonal components are close to zero. The diagonal components are written as
cell(i)dJ%lculated uii = y V
procedures
d cell(i)
(5)
where i = x,y,z, cell(i) is an array containing the cell parameters and V is the volume of the cell. When the stress tensor is determined for different volumes, a problem similar to the one encountered for the energy calculation is met. We can then distinguish the calculated stress corresponding to a locally constant NP” and the corrected stress corresponding to a strictly constant E,,, [ Ill. This difference is defined as the Pulay stress. Thus
(6) The last term is obtained as above from three calculations of total energy with different Np” at a fixed volume. The inclusion of this correction yields physical
The BP04 structure crystallizes in the tetragonaltype structure according to the tetragonal space group 4,. The experimental cell parameters for a and c parameters are respectively 4.332 A and 6.640 A. In the tetragonal system, b is equal to c and the three angles are equal to 90”. The primitive cell holds six atoms, that is one BP04 formula [12]. BP04 is an ionocovalent-type solid. Fig. 2 represents the tetragonal multiple cell. 3.1. Optimization
with energy minimization
In order to obtain the equilibrium cell parameters, two parameters are sufficient to vary: the cell sides a and c, or equivalently the cell volume, and the ratio c/u. The main characteristic of our approach is that instead of optimizing the geometry only with energy minimization criterion, we mix two methods: the first one minimizes the total energy and the second one looks at the stress tensor for the determination of the c/a ratio. This allows us to save computational time without losing accuracy.
245
li. Louis-Achille et al/Journal of Molecular Structure (Theochem) 426 (1998) 241-247 E
ieV)
Table I Calculated
and experimental
cell parameters
-2ssa.ab t
2OSS.bO
7
“.,
-2093.aa
4
‘\;
.
Calculated Experimental
a (I\)
(’ (A,
cla
4.329 4.332
6.659 6.640
1.538 1.533
In the second step, the volume is kept equal to the value corresponding to the above-determined minimum and c/a is varied about 1.5% around the experimental value and the anisotropic part of the stress tensor is studied. As BP04 is a tetragonal structure, the anisotropic part is reduced to (u,~,+,) where z stands for the c direction and x for the a direction. The minimum of the total energy corresponds to the volume where the stress tensor is isotropic, that is the anisotropic part is zero [ 131. By a linear extrapolation we calculated this volume, then we deduced the c/a ratio. In Table 1 are reported the cell parameters and the corresponding c/a ratio from calculation after the correction and from experiment [ 121. The difference between the experimental and the calculated values is small, 0.07% for a, 0.3% for c and 0.3% for the c/a ratio. These results follow the general trends that LDA calculations give good equilibrium structures due to cancellation of errors. 3.2. Optimization
(b)
Fig. 3. Total energy evolution with respect to isotropic variations of the cell volume: (a) without the correction factor; (b) with the correction factor.
In the first step, the da ratio is kept constant and the parameters a and c are varied isotropically. For each volume the total energy is determined and the atomic coordinates are relaxed using the Hellman-Feynmann forces as a guide; the correction defined by Eq. (3) is applied on top of the minimization. In Fig. 3(a) and Fig. 3(b) are reported the curves of the evolution of the energy vs. the cell volume with and without energy correction respectively. The effect of the correction is twofold: it smooths out the curve and it shifts the minimum of the curve toward the origin.
with pressure
calculation
Another method allows us to calculate the cell parameters of the equilibrium. As mentioned previously, a system is considered as being in an equilibrium state if the anisotropic part of the stress tensor is zero. Thus, we can calculate the external pressure P needed for maintaining the crystal at the considered non-equilibrium volume which, for the tetragonal structure of BP04, is written as follows: P = - $(c& + ax,, + a,)
(7)
where uXX) u I)”u :Z are the diagonal components of the stress tensor [l 11. At the equilibrium state, since the anisotropic part of the stress is zero, uXXis equal to uZZ.In addition, the summation of uXX,uVV and cr, is equal to zero, that is the external pressure is zero. Using this method we realized only the first step of the last method. For each volume (the same as the previous one used for the energy minimization) we calculated the stress tensor
246
V. Louis-Achille
et d/Journal
of Molecular Structure (Theochem) 426 (1998) 241-247
IJ(GPaj
Table 2 Calculated
k
and experimental
Calculated Experimental
11‘I
117
I20
123
126
129
132
13s
we 64
cell parameters
a (‘Q
c (A)
cia
4.332 4.332
6.623 6.640
1.538 1.533
the c/a ratio. The results are given in Table 2. They are in good agreement with experiments, meaning that the method is also well adapted for the calculation of the structural parameters of equilibrium. Furthermore, results obtained with total energy and stress optimization procedures are consistent, meaning that the optimization and correction schemes have been performed properly. For all these reasons, we think that our approach is suitable for the geometry calculation of solids.
P( GPn)
4. Subsequent
f\
physical properties
Obtaining the equilibrium geometry with accuracy is not a purpose in itself. We are more particularly concerned with the calculation of physical properties (cohesive energy, bulk modulus) which can be derived from the dependence of the total energy on structural parameters. 4. I. Cohesive energy I 114
I
/
117
120
-7----7---j-4 123
126
129
132
135
VA')
@I
Fig. 4. Pressure evolution with respect to isotropic variations of the cell volume: (a) without the correction factor; (b) with the correction factor.
and then the pressure. The point where the pressure is zero gives the equilibrium volume. For each calculation, we apply the correction on the stress. Fig. 4(a) and Fig. 4(b) show the resulting curves without and with correction. The shape of the curve does not vary with the correction but, contrary to what happens for the energy, the zero is shifted towards greater values of the volume. Then, we use the above result according to the cancellation of the anisotropic part of the stress on
The cohesive energy of a solid is defined as the difference between the crystal total energy and the energy of the isolated atoms at zero temperature. Both the pseudopotential energies of the equilibrium crystal structures and of the isolated atoms were calculated using Plane-Wave. Plane-Wave atomic calculations are performed putting each atom in the center of a cubic box with a Pzzl periodicity in order to keep a spherical symmetry. We choose the ridge large enough so that the interactions between adjacent atoms are negligible. The calculated energy includes a spin-polarization correction. Thermal and vibrational contributions were not considered in the calculations of cohesive energy, since they are expected to be much smaller than the inherent errors from LDA. The experimental value of the cohesive energy is calculated with the reticular energy of BPOc the ionization potential of B and P and the electronic affinity of 0. The
V. Louis-Achille et al./Journal of Molecular Structure (TheochemJ 426 (I 998) 241-247
reticular energy is the difference between the crystal energy and the energy of the isolated gaseous ions. The experimental value of the cohesive energy is 47.6 eV/ BP04 [ 141. The calculated value is 5 1.5 eV/BPO+ We find a good agreement between the calculation and the experiment. An overestimation is usual within the LDA approximation. 4.2. Bulk modulus The bulk modulus is a measure of the resistance to an external pressure. For a crystal of volume V in an isotropic environment where the pressure has a given value, if the pressure is changed by an amount AI’, the volume of the crystal changes by an amount of AV. The isothermal bulk modulus B of the crystal is defined in a manner analogous to the definition of Young’s modulus and B = - V(AP/A V) for constant temperature. In the limiting case, where the pressure change has an infinitesimal value dP, B becomes -V(dP/dV). An estimate of the zero pressure isothermal bulk modulus and its pressure derivative was obtained using the total energy calculation technique. The values for the total energy as a function of the cell volume were fitted, using a least squares technique to Mumaghan’s equation of state [ 151:
1
E(V)= VB (v’vo)B” + 1 +const &I [ B;-1
(8)
where BO, and Bb are the bulk modulus and its pressure derivative at the equilibrium volume VO. The calculated bulk modulus is 28 GPa and its pressure derivative is 5. No experimental data of bulk modulus are available in the literature at present. The prediction of the bulk modulus is generally more difficult than that of lattice constants, since elastic properties need more accuracy. It is also very difficult to evaluate the error of these estimates. However, the uncertainty of B is related to those of the cell parameters and to the derivative Bb [ 161 and an uncertainty of 0.3% in lattice parameters and a Bb value of 5 leads to an uncertainty of 6% on BO.
247
simple crystal structure contains six atoms in the primitive cell including only light atoms. The method proposed in this paper has been shown to provide valuable results on the ground state structure of solids on the one hand, as well as on their macroscopic properties (e.g. cohesive energy, elastic properties) in this test case on the other hand. The introduction of energy and stress corrections allowed us to run with a reduced cut-off, thereby saving computational time. The computed parameters are in good agreement with the experimental data for both geometry optimization methods. Lattice parameters are within 0.3% of the experimental structure and cohesive energy is about 8%. Consequently, the bulk modulus is estimated within 6%. Moreover, the results obtained with total energy and stress optimization procedures are consistent, meaning that the optimization and correction schemes have been performed properly. The underlying strategy to calculate all these quantities has been developed, and the numerical tests performed on BP04 are promising and open up the prospects for future extension to more complex systems. Results along these lines are presented in another paper for other iono-covalent-type phosphates, but the strategy can also be extended to other types of mineral.
References [l] Plane-Wave 3.0.0, M.S.I., San Diego, 1994. [2] D.M. Ceperley, B.J. Alder, Phys. Rev. Lett. 45 (1980) 566. M.P. Teter, unpublished results, 1990. [3] P. Hohenberg, W. Kohn, Phys. Rev. B. 136 (1964) 864. [4] W. Kohn, L.J. Sham, Phys. Rev. A. 140 (1965) 1133. [5] M.P. Teter, M.C. Payne, D.C. Allan, Phys. Rev. B. 40 (1989) 12255. [6] [7] [8] [9] [IO]
[I l] [ 121
5. Conclusion
[13] [ 141
A local density functional pseudopotential approach was applied to the iono-covalent solid BPO?. Its
[15] [I61
N. Troullier, J.L. Martins, Phys. Rev. B. 43 ( 1991) 1993. L. Kleinman, P.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. H.J. Monkhorst, J.D. Pack, Phys. Rev. B. 13 (1976) 5188. M.C. Payne. M.P. Teter, D.C. Allan, T.A. Arias, J.D. Joannopopoulos, Rev. Mod. Phys. 64 (1992) 1045. G.P. Francis, M.C. Payne, J. Phys.: Condens. Matter 2 ( 1990) 4395. P. Dacosta, O.H. Nielsen, K. Kunc, J. Phys. C: 19 (1986) 3 163. R.W.G. Wychkoff, Crystal Structure, vol. 4, Interscience. New York, 1968. R.J. Needs, R.M. Martin, Phys. Rev. B 43 (1991) 1993. Handbook of Physical Constants, The Geological Society of America, 1996. F.D. Mumaghan, Proc. Natl. Acad. Sci. USA 30 (1994) 422. O.H. Neilsen, R.M. Martin, Phys. Rev. B 32 (1985) 3792.