Ab initio calculation of transport and optical properties of aluminum: influence of simulation parameters

Ab initio calculation of transport and optical properties of aluminum: influence of simulation parameters

D.V. Knyazev d.v.knyazev@yandex.ru P.R. Levashov Joint Institute for High Temperatures RAS, Izhorskaya 13 bldg. 2, Moscow 125412, Russia Moscow Institute of Physics and Technology (State University), Institutskii per. 9, Dolgoprudny, Moscow Region, 141700, Russia
Abstract

This work is devoted to the ab initio calculation of transport and optical properties of aluminum. The calculation is based on the quantum molecular dynamics simulation, density functional theory and the Kubo-Greenwood formula. Mainly the calculations are performed for liquid aluminum at near-normal densities for the temperatures from melting up to 20000 K. The results on dynamic electrical conductivity, static electrical conductivity and thermal conductivity are obtained and compared with available reference and experimental data and the calculations of other authors. The influence of the technical parameters on the results is investigated in detail. The error of static electrical conductivity calculation is estimated to be about 20%; more accurate results require bigger number of atoms.

keywords:
quantum molecular dynamics, density functional theory, Kubo-Greenwood formula, transport and optical properties, aluminum
journal: Computational Materials Science

1 Introduction

The knowledge of the transport and optical properties of matter in the warm dense matter regime is crucial for many modern fundamental and applied researches: femtosecond laser heating Veysman et al. (2008), experiments with exploding wires DeSilva and Katsouros (1998), investigation of the metal-nonmetal transition Korobenko and Rakhel (2012).

Ab initio calculations of the transport and optical properties based on the quantum molecular dynamics (QMD) simulation, finite temperature density functional theory (FT-DFT) and the Kubo-Greenwood formula are widely used nowadays. The computations are performed for the substances of different types: liquid metals Desjarlais et al. (2002), hydrogen Collins et al. (2001), molten salts Silvestrelli et al. (1996), silica Laudernet et al. (2004) and plastics Lambert and Recoules (2012).

The paper Silvestrelli (1999) is one of the first works that applies the method specified to the calculation of transport and optical properties of aluminum. The works Desjarlais et al. (2002); Mazevet et al. (2005); Recoules et al. (2002) report optical properties and the static electrical conductivity of aluminum in the region on the phase diagram available with exploding wire experiments. The obtained results on the static electrical conductivity are compared with the exploding wire experiments DeSilva and Katsouros (1998); Recoules et al. (2002). Other papers Recoules and Crocombette (2005); Alemany et al. (2004); Knider et al. (2007) provide the comparison with slow thermophysical experiments on the optical properties, static electrical conductivity and thermal conductivity.

The number of atoms in the supercell is an important technical parameter during the calculation according to the method discussed. The number of atoms used in the papers on aluminum mentioned above is usually not larger than 108. Our previous paper Povarnitsyn et al. (2012) reports results on the optical properties of aluminum calculated with 108 atoms in the supercell. Only the work Alemany et al. (2004) demonstrates the computation with 205 atoms.

On the other hand, the most recent works Lambert et al. (2011) and Pozzo et al. (2011) show significant influence of the size effects on the transport and optical properties. The supercell should contain 1000 — 2000 atoms for convergence with the number of atoms to be achieved. The calculations in the papers Lambert et al. (2011) and Pozzo et al. (2011) are performed for hydrogen and sodium, respectively. The recent paper Alfè et al. (2012) shows that the increasing of the number of atoms up to 1024 is necessary to achieve well-converged results for solid iron.

In this paper we perform calculations for aluminum, taking into account three electrons per atom. Mainly we use 256 atoms in the supercell. The influence of the main technical parameters on the results is discussed. We vary the technical parameters and estimate the error of the static electrical conductivity calculation.

The transport and optical properties of aluminum were mainly calculated for the liquid phase at near-normal densities for the temperatures from melting up to 20000 K. The results are compared with available experimental and reference data and calculations by other authors.

Our paper is organized as follows. The method of calculation is described in section 2. Here the main technical parameters are introduced. The section 3 gives the numerical values of the technical parameters. The choice of these parameters is explained, the convergence according to these parameters is investigated, the error of the static electrical conductivity is estimated. The section 4 contains some technical information on the expression for the calculation of the dynamic Onsager coefficients. The main results on the optical properties, static electrical and thermal conductivity are described in section 5.

2 Method of Calculation

The calculation consists of three main parts: quantum molecular dynamics simulation, precise band structure resolution and calculation of transport and optical properties via the Kubo-Greenwood formula.

Initially atoms are placed to the supercell with periodic boundary conditions. Then during the QMD simulation at given density and temperature the ionic trajectories are calculated. The independent ionic configurations are selected from the equilibrium section of the QMD simulation for the further calculation of optical properties and thermal conductivity.

Precise band structure calculation is performed for the selected ionic configurations. As well as during the QMD simulation electronic structure is calculated via the density functional theory. Technical parameters, used during the precise band structure resolution lead to more accurate results than during the QMD simulation. The obtained wave functions are used to calculate the matrix elements of the nabla operator.

The real part of the dynamic electrical conductivity is calculated via the Kubo-Greenwood formula (9) using the matrix elements, energy eigenvalues and Fermi-weights obtained during the precise resolution of the band structure. The values of the dynamic electrical conductivity obtained for the different ionic configurations are averaged. The imaginary part of the dynamic electrical conductivity is reconstructed via the Kramers-Kronig transformation (11). If the real and imaginary parts of the dynamic electrical conductivity are known the optical properties may be calculated.

The dynamic Onsager coefficients are calculated via the Kubo-Greenwood formula (23). The values of the Onsager coefficients obtained for the different ionic configurations are averaged. The static Onsager coefficients are obtained by extrapolation of the dynamic to the zero frequency. Thermal conductivity is expressed via the static Onsager coefficients (21).

2.1 Quantum molecular dynamics simulation

During the QMD simulation the atoms at given density are placed to the cubic supercell with periodic boundary conditions.

At each step of the QMD simulation the electronic structure is calculated within the framework of density functional theory (DFT). The calculation is performed in the Born-Oppenheimer approximation: due to the great difference in masses between electrons and ions, the electrons immediately adjust to the current spatial positions of ions.

The electronic structure is calculated by solving the Kohn-Sham equations at the finite temperature Martin (2004):

(1)
(2)
(3)
(4)

Here are the Kohn-Sham wave functions,  — the corresponding energy eigenvalues,  — the corresponding Fermi-weights,  — the electron density,  — the chemical potential,  — the number of electrons included in calculation,  — the Planck constant, and  — the electron mass and charge respectively.

is the exchange-correlation potential defined as a variation of exchange-correlation functional

(5)

is the ion potential. The ions are treated within the pseudopotential approach, which significantly reduces computational effort. The choice of the pseudopotential and exchange-correlation functional is discussed below (see 3.5).

The electron temperature is introduced by the parameter in the Fermi-Dirac distribution (3).

The solution to the Kohn-Sham equations (1)–(4) for a supercell with periodic boundary conditions is the energy eigenvalue and the corresponding wave function which is represented as the set of plane wave orbitals:

(6)

Here is the volume of the supercell, are the coefficients of the expansion over plane waves, G is so-called G-vector, its coordinates are multiples of reciprocal lattice basis vectors. The number of the G-vectors is controlled by the energy cut-off :

(7)

The energy cut-off is the important technical parameter; its choice is discussed below (see 3.4).

k-vectors in the expansion (6) are called k-points in the Brillouin zone. The number of k-points is also an important technical parameter, we briefly discuss its influence on the results in A.

Each solution to the Kohn-Sham equations (1)–(4) in the supercell with periodic boundary conditions is called a band in this paper. The number of bands included to the solving of the Kohn-Sham equations is a technical parameter, discussed below (see 3.6).

The forces acting on the ions are calculated during the QMD simulations via the Hellmann-Feynman theorem. Ions are treated classically. The Newton equations acting on the ions are solved using the Verlet algorithm. The temperature of ions is maintained via the Nose-Hoover thermostat.

The equilibrium section of the QMD simulation is used to choose the ionic configurations for the further precise resolution of the band structure. The choice of ionic configurations is discussed below (see 3.1).

2.2 Precise resolution of the band structure

The ionic configurations selected during the QMD simulation are used to calculate the band structure. For these configurations the Kohn-Sham equations (1)–(4) are solved one more time, but with the larger energy cut-off, number of bands, number of k-points in the Brillouin zone. The ions are not moved. This stage is referred to as the precise resolution of the band structure in this work.

The technical parameters used in the precise resolution of the band structure are discussed in section 3.

The Kohn-Sham wave functions are necessary to calculate the matrix elements of the nabla operator . Some technical details about the calculation of the matrix elements are discussed below (see 3.5).

2.3 Calculation of optical properties

The complex dynamic electrical conductivity is the complex coefficient between the current density and applied electric field of the frequency :

(8)

The real part of the dynamic electrical conductivity (also referred to as simply dynamic electrical conductivity in this work) is connected with the energy absorbed by the electrons. It is calculated via the Kubo-Greenwood formula for each ionic configuration:

(9)

Here sum is over all -points within the Brillouin zone (by ), all bands participating in the calculation (by and ), and three spatial dimensions (by ). is the weight of the particular k-point in the Brillouin zone. A simple, short and intuitively clear derivation of the Kubo-Greenwood formula for the dynamic electrical conductivity may be found in Moseley and Lukes (1978).

For the practical calculation the -function in the Kubo-Greenwood formula (9) should be broadened by the Gaussian function:

(10)

where is the broadening of the -function. The broadening is an important technical parameter, the choice of its value is described below (see 3.3).

The values of the dynamic electrical conductivity obtained for the different ionic configurations are averaged.

If the real part of the dynamic electrical conductivity is known, the imaginary part may be reconstructed via the Kramers-Kronig transformation Collins et al. (2001):

(11)

where stands for the principle value of the integral.

If the complex dynamic electrical conductivity is known, the following optical properties may be calculated:

complex dielectric function :

(12)

where is the dielectric permittivity of vacuum (SI units);

complex index of refraction :

(13)
(14)

reflectivity (at the normal incidence of laser radiation):

(15)

absorption coefficient:

(16)

2.4 Calculation of thermal conductivity

The Onsager coefficients (referred to as the static Onsager coefficients in this work) are introduced as the coefficients that connect the electric current density and heat current density with the applied constant in time electric field and temperature gradient :

(17)
(18)

The Onsager coefficient is the static electrical conductivity . The Onsager reciprocal relation is valid for the static coefficients and Ziman (1960):

(19)

The typical experiment on the thermal conductivity measurement is performed at the zero electric current density. Under this condition the thermal conductivity is defined as the factor between the heat current density and the applied temperature gradient:

(20)

Setting to zero electric current density in (17) thermal conductivity may be expressed from (18) and (20):

(21)

It may be shown Ashcroft and Mermin (1976), that at relatively low temperatures the second term in expression (21) (referred to as the thermoelectric term in this work) gives negligible contribution to the thermal conductivity and approximate expression is valid:

(22)

In this work we calculate thermal conductivity via the exact expression (21) and verify, whether approximate expression (22) is valid (see 5.2.3).

In this work we use the Kubo-Greenwood formula to calculate the dynamic Onsager coefficients :

(23)

The Onsager coefficient is the real part of the dynamic electrical conductivity and equation (23) reduces to equation (9). The other dynamic Onsager coefficients do not have direct physical meaning in this work and are used only to get the static Onsager coefficients by the extrapolation to the zero frequency. The static electrical conductivity is obtained by the extrapolation of to the zero frequency. The procedure of the extrapolation is described below (see 3.7).

The equation (23) gives the equal values of and at non-zero frequencies. Some technical details connected with this fact are discussed below (section 4).

The values of the dynamic Onsager coefficients obtained for the different ionic configurations are averaged.

The derivation of the equation (23) may be found in Holst et al. (2011).

Finally, the experimentally discovered Wiedemann-Franz law should be mentioned. According to it, the ratio of the thermal conductivity to the product of static electrical conductivity and temperature is constant:

(24)

where is the Boltzmann constant. The mentioned ratio is called the Lorenz number. Its value  WK mentioned in equation (24) is called the ideal value. In the work Chester and Thellung (1961) the Wiedemann-Franz law and the ideal value of the Lorenz number are explained theoretically under rather general assumptions. In this work we calculate the Lorenz number and thus verify the Wiedemann-Franz law.

In this work the QMD and band structure calculations are performed using the Vienna ab initio simulation package (VASP) Kresse and Hafner (1993, 1994); Kresse and Furthmüller (1996). We created a special parallel program module that computes the dynamic Onsager coefficients, extrapolates them to zero frequency and calculates the electrical and thermal conductivity. This module uses information on the wave functions or matrix elements from the VASP output (see 3.5).

3 Technical Details

This section is devoted to the dependence of the dynamic and static electrical conductivity on the technical parameters. The choice of these parameters will be demonstrated for liquid aluminum at temperature  K and density g/cm.

In this work we do not vary the number of k-points and use only -point in all calculations (see, however, A). In earlier work Desjarlais et al. (2002) it was shown that for aluminum plasma the increase of the number of k-points gave insignificant effect on the corresponding electrical conductivity values. Nevertheless, we plan to clarify this question in our future work in detail.

Currently the calculation for this point is performed with the following technical parameters: 256 atoms, 1500 steps of QMD simulation, the step of the QMD simulation is 2 fs. During the QMD simulation ultrasoft pseudopotential of Vanderbilt Vanderbilt (1990) (US-PP) and local density approximation (LDA) for exchange-correlation functional are used. QMD simulation is performed with only 1 k-point in the Brillouin zone (-point), energy cut-off is  eV, the number of bands is 500. Ionic configurations for the calculation of transport and optical properties are taken every 100 steps. The precise resolution of band structure is performed for each of these ionic configurations using the same pair US-PP, LDA as during the QMD simulation. Band structure is calculated with only 1 k-point in the Brillouin zone (-point), energy cut-off is  eV, the number of bands is 1300. The dynamic electrical conductivity and the dynamic Onsager coefficients are calculated for the frequencies from 0.005 eV up to 10 eV with the step 0.005 eV. The broadening of the -function in the Kubo-Greenwood formula is  eV. The static Onsager coefficients are obtained by the linear extrapolation of the dynamic values to zero frequency.

The dependence of the results on the technical parameters is examined as follows: one parameter is varied while others are kept fixed at the values mentioned above. Then the obtained information is used to estimate the error of static electrical conductivity calculation.

3.1 The ionic configurations

The typical dependence of the total energy of electrons and ions on the step number is shown in Fig. 1. At the start of simulation the atoms form the fcc lattice. The simulation is performed for 1500 steps, each step of simulation corresponds to 2 fs. Thus we observe the behavior of the system for 3 ps. During about 100 initial steps the lattice breaks down and energy comes to its equilibrium value. Then at the equilibrium section energy fluctuates around its average value. The configurations for the calculation of dynamic electrical conductivity are taken from this equilibrium section. The distance between neighbouring configurations is 100 steps.

Figure 1: (Color online) The dependence of the total energy of electrons and ions on the step number. Solid line is the dependence of energy on the step number; square points are configurations selected for the further precise resolution of the band structure.

To examine the convergence with the number of configurations we calculate the standard deviation of the mean conductivity value:

(25)

where is the value of the dynamic electrical conductivity for the ionic configuration, is the dynamic conductivity value averaged over configurations, is the number of ionic configurations. The ratio is less than 2% for all frequencies under consideration. Thus we can consider the statistical error of the conductivity calculation to be about 2%.

We can use expression (25) to estimate the statistical error only if the ionic configurations are independent. The significant decrease of the velocity-velocity autocorrelation function during the time between two subsequent configurations could be a good criterion of the independence of configurations. In this work, however, this criterion is not used.

3.2 The dependence on the number of atoms

The change of the number of atoms strongly influences the calculated transport and optical properties. Fig. 2 shows the dependence of the calculated dynamic electrical conductivity on the number of atoms in the supercell.

Figure 2: (Color online) The dependence of the calculated dynamic electrical conductivity on the number of atoms. Solid line — 108 atoms; dashed line — 256 atoms.

The change of the number of atoms from 108 to 256 leads to the significant changes in the dynamic electrical conductivity. For example for the frequency  eV the conductivity changes by the factor of 1.3. At the frequencies near to zero the change of the conductivity is about 10%. Evidently, even larger numbers of atoms should be taken to achieve the convergence. But at the current stage it looks more reasonable to take the results obtained with the largest (256) number of atoms.

It should be mentioned, that the convergence with the number of atoms is strongly connected with the dependence of the results on the broadening of the -function in the Kubo-Greenwood formula. This is to be discussed in the next section.

a) b)

Figure 3: (Color online) The dependence of the calculated dynamic electrical conductivity on the broadening of the -function in the Kubo-Greenwood formula. a) 256 atoms; b) 108 atoms.

3.3 The dependence on the broadening of the -function

The dependence of the calculated transport and optical properties on the broadening of the -function in the Kubo-Greenwood formula (10) is rather important (especially at low frequencies). The dependence of the dynamic electrical conductivity at low frequencies on the broadening of the -function is shown in Fig. 3a.

The Kubo-Greenwood formula takes into account the transitions between all bands. In the real substance the number of atoms is very large, the bands are very close to each other. Therefore a smooth curve of dynamic electrical conductivity is formed. The number of atoms in numerical calculation is rather small, and the transition between each two bands gives its own -peak to the curve of dynamic electrical conductivity. To form a smooth curve, each -peak is to be broadened by the Gaussian function (10). If the broadening is too small (0.02 eV in Fig. 3a), the transitions between bands are still noticeable and the curve is oscillating. If the broadening is too large (0.2 eV in Fig. 3a), too many transitions contribute to each point on the curve and the physical dependence on the frequency is smoothed.

We consider, that the convergence according to the broadening must occur as follows: when the broadening is diminished, at first, excessive smoothing of physical dependence disappears. Then for some value the further diminishing of the broadening does not affect the . For some range of the curves do not depend on the broadening and any curve from this range may be taken as final. And only at very small values of the oscillations appear.

At rather high frequencies (above 1 eV in Fig. 3a) the situation is similar to the described above. But at low frequencies the situation is different. When the broadening is diminished the oscillations appear earlier than the convergence of the dynamic electrical conductivity is reached.

The convergence with the broadening can be improved by the increasing of the number of atoms. When the number of atoms is relatively big, it becomes possible to take smaller value of to reach the convergence according to this parameter and avoid the appearance of oscillations. Fig. 3b confirms this statement: the curves for the same values are less smooth than in Fig. 3a.

At the currently available number of atoms (256, Fig. 3a) the error connected with the choice of is not less than 3% (the difference between the curves for  eV and  eV).

3.4 Convergence according to the energy cut-off

The value of energy cut-off  eV is used during the precise resolution of the band structure. The increasing of the energy cut-off to the value of 400 eV leads to the negligible (less than 1%) change of dynamic electrical conductivity within the whole range of frequencies under consideration. On the other hand, the diminishing of the energy cut-off to the value of 100 eV leads to the excessive (up to 13%) changes in the dynamic electrical conductivity. Therefore the value  eV is optimal during the precise resolution of the band structure.

During the QMD simulation the smaller value of the energy cut-off  eV is used. The increasing of the energy cut-off to the value 200 eV leads to the changes of dynamic electrical conductivity not larger than 3% in the whole range of frequencies under consideration. On the other hand, the QMD simulation with the energy cut-off 200 eV requires significant computation time.

Therefore the error connected with energy cut-off convergence is estimated as 3%.

3.5 The dependence on the pseudopotential and exchange-correlation functional

The ultrasoft pseudopotential of Vanderbilt Vanderbilt (1990) (US-PP) in pair with local density approximation of the exchange-correlation functional (LDA) is used during QMD simulation in this work. US-PP was chosen due to its high computational efficiency.

We also use the same pair US-PP, LDA during the precise resolution of the band structure. This leads to some difficulties. At this stage we get pseudo wave functions from VASP, use them to calculate pseudo matrix elements , and then use them to calculate the dynamic electrical conductivity via the Kubo-Greenwood formula (9). The recalculation of pseudo matrix elements to all electron matrix elements is a difficult procedure for US-PP (as mentioned in Desjarlais et al. (2002)), which is not currently implemented in the VASP package. We use parallel program module created for this work to calculate pseudo matrix elements and dynamic electrical conductivity. Currently the recalculation of pseudo matrix element to all electron is omitted.

An alternative approach is to use projected augmented-wave (PAW) pseudopotential Kresse and Joubert (1999) during the precise resolution of the band structure. The VASP package provides the module that calculates all electron matrix elements . Unfortunately in the current VASP version this module is not parallelized. Then we use matrix elements from VASP to calculate dynamic electrical conductivity via the Kubo-Greenwood formula. In this work we perform one calculation with PAW pseudopotential in pair with Perdew, Burke, Ernzerhof (PBE) exchange-correlation functional Perdew et al. (1996).

The replacement of US-PP, LDA by the PAW, PBE pair gives the change of the dynamic electrical conductivity of not more than 4% in the whole range of frequencies. The largest change (4%) is obtained for the frequencies near to zero frequency. The small difference in results makes it possible to use US-PP in spite of difficulties described above.

Figure 4: (Color online) The dependence of the band number (dashed line) and Fermi-weight (solid line with points) on the energy counted from the chemical potential for one ionic configuration. The points on the curve show the discrete bands.

3.6 The choice of the number of bands

During the QMD simulation all the bands with significant Fermi-weights should be taken into account. If some of the occupied bands are neglected, their contribution to electron density will be missed and ionic configurations will be calculated inaccurately. We consider the bands with the Fermi-weights larger than (standard VASP precision during output of the Fermi-weights) to be occupied. Fig. 4 shows the dependence of the band number and Fermi-weight on the energy counted from the chemical potential for one ionic configuration. The Fermi-weights become less than at  eV. This corresponds to 500 bands.

Precise resolution of the band structure for the further calculation of the dynamic electrical conductivity requires the larger number of bands. In this work the dynamic electrical conductivity is calculated for the frequencies up to  eV. The Kubo-Greenwood formula (9) shows that the bands higher than occupied bands by the energy  eV give contribution to the value of conductivity . Because of the broadening of the -function this quantity should be additionally increased by  eV. For the example shown in Fig. 4 the bands with the energies up to  eV should be taken into account. This corresponds to 1300 bands.

3.7 Extrapolation to the zero frequency

The static value of electrical conductivity is obtained by the extrapolation of dynamic values to the zero frequency. The procedure of interpolation is shown in Fig. 5. The dynamic electrical conductivity is calculated on the dense frequency mesh (in this work the frequency step is 0.005 eV). If the number of atoms is large enough and the broadening of the -function is chosen correctly the dynamic electrical conductivity at low frequencies is a smooth function of frequency. The linear extrapolation based on the two points nearest to the zero frequency gives the value of the static conductivity. The error of the extrapolation is negligibly small, and the error of the static value is completely determined by the error of the dynamic values used for extrapolation.

Figure 5: (Color online) The extrapolation of the dynamic electrical conductivity to the zero frequency. Square points — calculated values of the dynamic electrical conductivity; dashed line — linear extrapolation to the zero frequency.

3.8 Estimation of static electrical conductivity calculation error

If the influence of different technical parameters on the results is examined it is possible to estimate the calculation error.

The largest contribution to the error of static electrical conductivity calculation is given by the change of the number of atoms. When the number of atoms is increased from 108 to 256 the static electrical conductivity grows by 10%. The behavior of the static electrical conductivity under the further increasing of the number of atoms is not clear. But taking into account available data the error connected with the number of atoms may be estimated as +10%.

The replacement of the US-PP, LDA pair by the PAW, PBE pair leads to the decrease of the static electrical conductivity by 4%. The error connected with the pseudopotential and exchange-correlation functional choice can be estimated as .

The error connected with the energy cut-off choice is about .

The choice of the broadening of the -function yields the error +3%.

The statistical error is .

The error connected with the extrapolation of the dynamic electrical conductivity to zero frequency is negligibly small.

Unfortunately, currently the dependence of the results on the number of k-points is not examined because of too large computational requirements (see, however, the analysis of computational errors in A).

We estimate total error by the simple algebraic sum of errors from different sources. Thus we estimate the range where the true value of static electrical conductivity may be. The factors increasing the value of static electrical conductivity give the total error of . The factors decreasing the value of static electrical conductivity give the total error of .

We consider that the precision of the calculation of static electrical conductivity can be improved by further increasing of the number of atoms and by the examination of the convergence with the number of k-points.

Figure 6: (Color online) The dynamic Onsager coefficients and calculated according to the asymmetric formula (27) for the different broadenings of the -function. Filled symbols — , empty symbols — .

4 Dynamic Onsager Coefficients

In this section we discuss different expressions for calculation of the dynamic Onsager coefficients. The comparison is performed for liquid aluminum at  K and  g/cm. The technical parameters are the same as described in section 3, except for the broadening of the -function  eV.

The Kubo-Greenwood formula for the calculation of the dynamic Onsager coefficients (23) is derived in the paper Holst et al. (2011). The Onsager coefficient is the dynamic electrical conductivity . Other dynamic Onsager coefficients in this work have no direct physical meaning and are used only to obtain static values by the extrapolation to the zero frequency. In the limit of the zero frequency the Kubo-Greenwood formula (23) reduces to:

(26)

Here the expression in the Kubo-Greenwood formula (23) tends to in the limit of zero frequency. Also due to the -function in the formula (26) the values and are equal and there is no difference between using the factor of or .

If the dynamic Onsager coefficients are treated only as some expressions which give correct values in the limit of zero frequency, different expressions for the dynamic Onsager coefficients may be chosen. In paper Recoules and Crocombette (2005) the asymmetric expression is used:

(27)

that gives correct zero frequency limit (26).

The Onsager coefficients and calculated according to the formula (27) are shown in Fig. 6. The coefficients are not equal at non-zero frequencies. At the reasonable (see 3.3) broadenings of the -function 0.1 eV and 0.05 eV and coefficients diverge at the zero frequency limit. And only at very small broadening 0.01 eV two curves come to the close values at the zero frequency limit. Such a small value of broadening can not be used because of too large oscillations on the (see 3.3).

Figure 7: (Color online) The frequency dependence of the dynamic Onsager coefficients calculated according to the symmetric expression (23).

It can be easily seen from Fig. 6 that the curves and at the zero frequency limit have approximately the same values and opposite signs. If added and divided by two they compensate each other and form the smooth curve at the zero frequency limit. The half-sum gives the same value as the symmetric expression (23). The values of the Onsager coefficients calculated according to the symmetric expression (23) for  eV are shown in Fig. 7. In this case the curve is smooth at the zero frequency limit.

As concerns the Onsager coefficient, its frequency dependences calculated according to the symmetric (23) and asymmetric (27) expression are shown in Fig. 8. In the case of the coefficient, the asymmetric curve is smooth at the limit of zero frequency. The symmetric and asymmetric expressions give the close static values, although the symmetric value is 20% higher than asymmetric one.

Thus, both symmetric (23) and asymmetric (27) expressions for the dynamic Onsager coefficients theoretically give the correct zero frequency limit (26). But the symmetric expression is more numerically stable and gives no divergence at zero frequency limit at the numerical calculation of and . Also it is more justified theoretically Holst et al. (2011). In this work we use only the symmetric expression (23) to calculate the dynamic Onsager coefficients.

Figure 8: (Color online) The frequency dependence of the Onsager coefficient calculated according to the symmetric (solid line, equation (23)) and asymmetric (dashed line, equation (27)) expressions. The broadening of the -function is 0.07 eV.

5 Results and Discussion

5.1 Dynamic electrical conductivity and optical properties

5.1.1 The comparison with other calculations

We perform the calculation at the same temperatures and densities as in the well-known paper of Desjarlais Desjarlais et al. (2002) to test our computational scheme. The technical parameters of calculation are taken close to those of Desjarlais. The results are presented in Fig. 9. The good agreement means that the calculation performed previously may be well reproduced using the same computational package and the close values of the parameters.

a) b)

Figure 9: (Color online) The dynamic electrical conductivity of aluminum in comparison with calculation Desjarlais et al. (2002). Solid line — this work, dashed line — calculation Desjarlais et al. (2002). The calculation parameters: a) 32 atoms; PAW pseudopotential; PBE exchange-correlation functional; 1 k-point in the Brillouin zone; energy cut-off 200 eV; the broadening of the -function 0.1 eV; 15 ionic configurations; 1500 steps; 1 step – 2 fs. b) 108 atoms; US-PP pseudopotential; LDA exchange-correlation functional; 1 k-point in the Brillouin zone; energy cut-off 100 eV; the broadening of the -function 0.1 eV; 15 ionic configurations; 1500 steps; 1 step – 2 fs.

5.1.2 Normal isochor

The results on the dynamic electrical conductivity at the normal-density isochor are shown in Fig. 10. The frequency dependence of the dynamic electrical conductivity has the Drude-like shape. The dynamic electrical conductivity decreases with the temperature growth at the low frequencies, and, on the contrary increases at high frequencies. Also, it is interesting to mention, that the dynamic electrical conductivity does not depend on the temperature at  eV.

Figure 10: (Color online) The frequency dependences of the dynamic electrical conductivity at the normal isochor for different temperatures. The calculation parameters are the same as described in section 3.

5.1.3 Dynamic Onsager coefficients

In this work only the frequency dependence of the Onsager coefficient has the direct physical meaning of the real part of dynamic electrical conductivity . The frequency dependences of the other Onsager coefficients are used only to obtain the static values by the extrapolation to the zero frequency. Nevertheless, we consider that it may be useful to show also the curves of the dynamic Onsager coefficients to understand the procedure of the extrapolation to the zero frequency.

Fig. 11a shows the frequency dependences of the Onsager coefficient at density  g/cm for different temperatures from 1000 K up to 10000 K. The static value is the approximate value of the thermal conductivity (neglecting the thermoelectric term, (22)). The shape of the curve changes significantly with the temperature increasing. Nevertheless the static values form the smooth dependence on the temperature (see 5.2.2).

Fig. 11b presents the frequency dependences of the Onsager coefficient at the isochor  g/cm for the temperatures from 1000 K up to 10000 K. The static value is necessary to calculate thermoelectric term in the expression (21). The shape of the curve changes significantly with the temperature growth. The static values form the nonmonotonic dependence on the temperature. Also for all the temperatures less than 10000 K the thermoelectric term is negligibly small (see 5.2.3). These two facts mean that additional investigation should be performed to get precise values .

a) b)

Figure 11: (Color online) The frequency dependences of the Onsager coefficients (a) and (b) for aluminum,  g/cm. The technical parameters of the calculation are the same as in section 3 (except for the broadening of the -function at temperatures up to 2000 K,  eV).

5.2 Static electrical and thermal conductivity

5.2.1 Normal isobar

The calculation for liquid aluminum at the normal isobar for the temperatures from 973 K up to 1473 K was made to perform comparison with the reference data. The values of the density of liquid aluminum at the normal isobar, necessary during the QMD simulation are taken from Lide (2005). The results of the calculation in comparison with reference data are shown in Fig. 12.

The reference density of liquid aluminum at normal pressure and temperature  K is  g/cm. The error estimation in section 3 is performed for temperature  K and density  g/cm, that only slightly differs from the density at the normal isobar. Therefore, this error estimation is applicable for the results at the normal isobar, too.

The discrepancy between the calculated values and reference data is . The factors which increase the static electrical conductivity yield total error . The factors which decrease the static electrical conductivity yield total error . Therefore the reason of the discrepancy is still not clarified. Nevertheless, negative slopes of the static electrical conductivity curves are close for the calculated and reference data (see Fig.12).

The static electrical conductivity undergoes the most significant changes under the varying of the number of atoms (see 3.8). The result for less number of atoms — 108 is also shown in Fig. 12, the difference from the reference data is . Nevertheless we think it is more correct to consider the result with the larger number of atoms as preliminary in spite of its larger discrepancy with reference data. We consider that the further increasing of the number of atoms and the investigation of the convergence according to the number of k-points may improve the results.

Figure 12: (Color online) The comparison of the calculated static electrical conductivity on the normal isobar with the reference data Grigoriev and Meylikhov (1991). Circle with error bars — the calculated point  K,  g/cm, the error estimation for this point is performed in section 3, 256 atoms. Diamond — the calculated point at the normal isobar with the decreased number of atoms — 108. Squares — the calculated points at the normal isobar, 256 atoms. Triangles — reference data Grigoriev and Meylikhov (1991). The technical parameters are the same as described in section 3 (except for the point with 108 atoms).

a) b)

Figure 13: (Color online) The calculated static electrical conductivity (a) and thermal conductivity (b) in comparison with the results of other authors and experimental data along the isochor  g/cm. Filled triangles — the calculation of this work; filled squares — calculation Recoules and Crocombette (2005); empty square — experiment Rhim and Ishikawa (1998). The calculation parameters are the same as in section 3 (except for the broadening of the -function at temperatures up to 2000 K,  eV).

5.2.2 The isochor 2.35 g/cm

To compare our results with the results of other authors the static electrical conductivity and thermal conductivity of liquid aluminum at density  g/cm were calculated. The temperature was varied from 1000 K up to 10000 K. The results of the calculation in comparison with the results of other authors Recoules and Crocombette (2005) and experimental data Rhim and Ishikawa (1998) are shown in Fig. 13.

The static electrical conductivity decreases with the temperature growth, whereas thermal conductivity increases. For all the temperatures up to 10000 K the thermoelectric term in the equation (21) gives small (not larger than 2%) contribution to the thermal conductivity and the approximate formula (22) is valid. The calculated values of the static electrical conductivity are larger than in the calculation Recoules and Crocombette (2005) and experiment Rhim and Ishikawa (1998). We explain this in the same way as for the normal isobar: we use the larger number of atoms (256 atoms) than in the work Recoules and Crocombette (2005) (108 atoms). So the static electrical conductivity grows away from the calculation with the smaller number of atoms Recoules and Crocombette (2005) and experiment. At relatively low temperatures (up to 2000 K) the calculation conditions are close to that at the normal isobar, so the error estimation of the section 3 is applicable here.

As far as the thermal conductivity is concerned its values are also higher than in the calculation Recoules and Crocombette (2005) and experiment Rhim and Ishikawa (1998). The dependence of the thermal conductivity on the technical parameters is not investigated in this work and we can only suppose that the total error of the thermal conductivity calculation is close to that of the static electrical conductivity.

The calculated values of the Lorenz number in comparison with the calculation Recoules and Crocombette (2005) and ideal value  WK are shown in Fig. 14. The difference between the calculated Lorenz number, results of Recoules and Crocombette (2005) and ideal value is not larger than the estimated calculation error . The value of the Lorenz number for the experimental data Rhim and Ishikawa (1998) is ideal because in that work the thermal conductivity was recalculated from the electrical conductivity via the Wiedemann-Franz law.

Figure 14: (Color online) The Lorenz number at the isochor  g/cm. Triangles — calculation of this work; squares with error bars — calculation Recoules and Crocombette (2005); dashed line — the ideal value  WK.

5.2.3 Normal isochor

The results of the calculation at the normal isochor  g/cm and the isochor  g/cm are shown in Fig. 15.

Qualitatively the dependences of static electrical conductivity and thermal conductivity on the temperature at the density  g/cm have the same shape as at the density  g/cm. The static electrical conductivity and thermal conductivity decrease with the growth of density.

The dependences of the thermal conductivity on the temperature calculated according to the approximate expression (22) and the exact expression (21) taking into account the thermoelectric term are shown in Fig. 16. The relative contribution of the thermoelectric term increases with the temperature growth. At the temperatures less than 1000 K its contribution is less than 1% and approximate expression (22) is valid. At the maximum temperature under consideration  K the relative contribution of the thermoelectric term is 10%.

a) b)

Figure 15: (Color online) The calculated static electrical conductivity (a) and thermal conductivity (b) of aluminum. The technical parameters of the calculation are the same as described in section 3 (except for the broadening of the -function at the temperatures up to 2000 K,  eV).
Figure 16: (Color online) The influence of the thermoelectric term on the thermal conductivity. Filled squares — , thermal conductivity calculated according to the approximate expression (22). Empty squares — thermal conductivity calculated according to the exact expression (21).

6 Conclusions

1. The influence of the main technical parameters on the results was investigated for liquid aluminum at  K and  g/cm. It was shown, that the current number of atoms — 256 is still too small for the convergence according to the number of atoms to be achieved. The error of the calculation of the static electrical conductivity was estimated to be about 22%. The discrepancy with the reference data is about 25%. Further increasing of the number of atoms and the investigation of the convergence with the number of k-points is necessary.

2. The expressions for the dynamic Onsager coefficients from the articles Holst et al. (2011) and Recoules and Crocombette (2005) were compared. The symmetric expression from the article Holst et al. (2011) being theoretically justified was found out to be more numerically stable.

3. The calculations of the dynamic and static electrical conductivity as well as thermal conductivity were performed for liquid aluminum at near-normal densities for the temperatures form melting up to 20000 K. The results are in satisfactory agreement with reference and experimental data and the calculations of other authors.

Acknowledgements

This work was supported by the FAIR-Russia Research Center Grant for master students (2011-2012), FAIR-Russia Research Center Grant for PhD students (2012-2013) and the Russian Foundation for Basic Research, grants 13-02-91057-CNRS and 13-08-01179.

Appendix A Influence of the number of k-points

During the preparation of this article we obtained some additional information concerning the dependence of the results on the number of k-points in the Brillouin zone during precise resolution of band structure. The calculation was performed for liquid aluminum for the same -point as in section 3:  g/cm and  K. Due to huge computational time it is difficult to conduct this investigation for the supercell with 256 atoms. Nevertheless, computations were performed for 108 atoms in the supercell. All other technical parameters are the same as described in section 3.

In this Appendix firstly the meshes in the Brillouin zone are described. Then the procedure of reduction of the number of k-points due to symmetry reasons is introduced which allows to shorten significantly computational time with nearly negligible influence on the results. Further the obtained results are discussed. Though the convergence with the number of k-points was not achieved, the effects of the increment of the number of k-points may be observed. The changes in the static electrical conductivity caused by the introduction of denser meshes in the Brillouin zone are estimated.

a.1 k-point meshes

Two types of k-point meshes based on the Monkhorst-Pack scheme Monkhorst and Pack (1976) were used.

The first scheme is the -centered Monkhorst-Pack scheme. The coordinates of k-points in the Brillouin zone for the cubic supercell with the side are given by:

(28)

where , and run over all integer values in the range:

(29)

Here is called the size of the k-point mesh.

Both for even and odd mesh sizes the -point is included in the -centered Monkhorst-Pack scheme. The scheme is symmetric with respect to the -point for odd mesh sizes. However it is asymmetric in case of even mesh sizes because  ,  , and   coordinates are not included.

The second scheme is the original Monkhorst-Pack scheme, which is introduced for even meshes as:

(30)

where , and run over all integer values in the range:

Thus the scheme becomes symmetric with respect to the -point, but the -point itself is not included in the scheme. The original Monkhorst-Pack scheme for odd mesh sizes coincides with the -centered one.

Further the original Monkhorst-Pack scheme for even and odd meshes is designated as MP, and the -centered scheme for even meshes is contracted as MP. Both MP and MP schemes are described in the VASP manual vaspmanual (). In both schemes the weights of all k-points are set equal, with the unit value of the weights’ sum.

a.2 Reduction of the number of k-points

The number of k-points in the Brillouin zone may be reduced due to symmetry reasons. This enables the calculation with the sufficient numbers of atoms and ionic configurations on the one hand, and with dense enough meshes in the Brillouin zone on the other hand. The procedure of the reduction of the number of k-points is close to that mentioned in paper Pozzo et al. (2011).

The k-points with the coordinates differing only by the sign are reduced to one. For example, k-points and are merged into one k-point with the weight equal to the sum of weights of the initial k-points.

The k-points differing only by the order of coordinates are also reduced to one. For example, k-points and are merged into one k-point with the weight equal to the sum of weights of the initial k-points.

The absence of the particular direction during the simulation in the liquid phase may be mentioned as the support of the procedure described. However, in this paper the validity of this method is checked numerically, by the comparison of the results obtained for full and reduced meshes (see Fig. 17).

The reduction of the number of k-points is performed both for MP and MP meshes.

a.3 Obtained results

The obtained results on static electrical conductivity are presented in Fig. 17.

Figure 17: (Color online) The dependence of static electrical conductivity on the number of k-points. Here is the size of the mesh. Filled square — full MP scheme, empty squares — reduced MP scheme. Filled triangle — full MP scheme, empty triangles — reduced MP scheme. The dependence on the number of k-points is investigated for 108 atoms in the supercell. Other simulation parameters are the same as in section 3. The result for the calculation with 256 atoms and -point only is shown for comparison — filled circle.

Firstly the procedure of the reduction of k-points was validated. The results obtained both for full and reduced meshes are shown. The comparison is performed both for MP and MP schemes. The replacement of the full scheme by the reduced one yields less than 0.5% change in static electrical conductivity. The changes in dynamic electrical conductivity are less than 2% in the whole range of frequencies under consideration. This proves that the described procedure of reduction of the number of k-points is reliable enough and enables investigation of the influence of the mesh size on results.

The difference between the results obtained using MP and MP schemes is shown for reduced and meshes. The changes are noticeable for a small mesh size and become less as the mesh size grows.

The values of static electrical conductivity obtained for all k-point meshes are less than 14% larger than the value obtained for the -point only. This value is close to the effect of the replacement of 108-atom supercell by one with 256 atoms. The result for 256 atoms and -point only is also presented in Fig. 17.

However, the increment of the number of k-points hardly improves the convergence with the broadening of the -function, described in subsection 3.3. The dependence of the obtained dynamic electrical conductivity on the broadening of the -function at low frequencies is shown in Fig. 18. The results obtained for 108 atoms and reduced MP mesh (Fig. 18) are close to those obtained for 108 atoms and -point only (Fig. 3b). In both cases when the broadening is decreased the oscillations appear earlier than the convergence of dynamic electrical conductivity is achieved. Compared to this, increment of the number of atoms up to 256 (Fig. 3a) leads to more significant improvements of the convergence with the broadening . It is consistent with the statement of paper Lambert et al. (2011), that the convergence with the broadening of the -function may not be improved by the increment of the number of k-points.

Figure 18: (Color online) The dependence of the dynamic electrical conductivity on the broadening of the -function. The number of atoms in the supercell is 108, the reduced MP mesh of k-points is used. Other simulation parameters are the same as described in section 3.

References

  • Veysman et al. (2008) M. E. Veysman, M. B. Agranat, N. E. Andreev, S. I. Ashitkov, V. E. Fortov, K. V. Khishchenko, O. F. Kostenko, P. R. Levashov, A. V. Ovchinnikov, D. S. Sitnikov, J. Phys. B: At. Mol. Opt. Phys. 41 (2008) 125704.
  • DeSilva and Katsouros (1998) A. W. DeSilva, J. D. Katsouros, Phys. Rev. E 57 (1998) 5945–5951.
  • Korobenko and Rakhel (2012) V. N. Korobenko, A. D. Rakhel, Phys. Rev. B 85 (2012) 014208.
  • Desjarlais et al. (2002) M. P. Desjarlais, J. D. Kress, L. A. Collins, Phys. Rev. E 66 (2002) 025401.
  • Collins et al. (2001) L. A. Collins, S. R. Bickham, J. D. Kress, S. Mazevet, T. J. Lenosky, N. J. Troullier, W. Windl, Phys. Rev. B 63 (2001) 184110.
  • Silvestrelli et al. (1996) P. L. Silvestrelli, A. Alavi, M. Parrinello, D. Frenkel, Phys. Rev. B 53 (1996) 12750–12760.
  • Laudernet et al. (2004) Y. Laudernet, J. Clérouin, S. Mazevet, Phys. Rev. B 70 (2004) 165108.
  • Lambert and Recoules (2012) F. Lambert, V. Recoules, Phys. Rev. E 86 (2012) 026405.
  • Silvestrelli (1999) P. L. Silvestrelli, Phys. Rev. B 60 (1999) 16382–16388.
  • Mazevet et al. (2005) S. Mazevet, M. P. Desjarlais, L. A. Collins, J. D. Kress, N. H. Magee, Phys. Rev. E 71 (2005) 016409.
  • Recoules et al. (2002) V. Recoules, P. Renaudin, J. Clérouin, P. Noiret, G. Zérah, Phys. Rev. E 66 (2002) 056412.
  • Recoules and Crocombette (2005) V. Recoules, J.-P. Crocombette, Phys. Rev. B 72 (2005) 104202.
  • Alemany et al. (2004) M. M. G. Alemany, L. J. Gallego, D. J. González, Phys. Rev. B 70 (2004) 134206.
  • Knider et al. (2007) F. Knider, J. Hugel, A. V. Postnikov, J. Phys.: Condens. Matter 19 (2007) 196105.
  • Povarnitsyn et al. (2012) M. E. Povarnitsyn, D. V. Knyazev, P. R. Levashov, Contrib. Plasma Phys. 52 (2012) 145–148.
  • Lambert et al. (2011) F. Lambert, V. Recoules, A. Decoster, J. Clérouin, M. Desjarlais, Phys. Plasmas 18 (2011) 056306.
  • Pozzo et al. (2011) M. Pozzo, M. P. Desjarlais, D. Alfè, Phys. Rev. B 84 (2011) 054203.
  • Alfè et al. (2012) D. Alfè, M. Pozzo, M. P. Desjarlais, Phys. Rev. B 85 (2012) 024102.
  • Martin (2004) R. M. Martin, Electronic Structure. Basic Theory and Practical Methods, University Press, Cambridge, 2004.
  • Moseley and Lukes (1978) L. L. Moseley, T. Lukes, Am. J. Phys. 46 (1978) 676–677.
  • Ziman (1960) J. M. Ziman, Electrons and Phonons: The Theory of Transport Phenomena in Solids, Clarendon Press, Oxford, 1960.
  • Ashcroft and Mermin (1976) N. W. Ashcroft, N. D. Mermin, Solid State Physics, Holt, Rinehart, and Winston, New York, 1976.
  • Holst et al. (2011) B. Holst, M. French, R. Redmer, Phys. Rev. B 83 (2011) 235120.
  • Chester and Thellung (1961) G. V. Chester, A. Thellung, Proc. Phys. Soc. 77 (1961) 1005–1013.
  • Kresse and Hafner (1993) G. Kresse, J. Hafner, Phys. Rev. B 47 (1993) 558–561.
  • Kresse and Hafner (1994) G. Kresse, J. Hafner, Phys. Rev. B 49 (1994) 14251–14269.
  • Kresse and Furthmüller (1996) G. Kresse, J. Furthmüller, Phys. Rev. B 54 (1996) 11169–11186.
  • Vanderbilt (1990) D. Vanderbilt, Phys. Rev. B 41 (1990) 7892–7895.
  • Kresse and Joubert (1999) G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758–1775.
  • Perdew et al. (1996) J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868.
  • Lide (2005) D. R. Lide (Ed.), CRC Handbook of Chemistry and Physics, Internet Version 2005, CRC Press, Boca Raton, FL, 85th edition, 2005.
  • Grigoriev and Meylikhov (1991) I. S. Grigoriev, E. Z. Meylikhov (Eds.), Physical Quantities, Energoatomizdat, Moscow, 1991 [in Russian].
  • Rhim and Ishikawa (1998) W.-K. Rhim, T. Ishikawa, Rev. Sci. Instrum. 69 (1998) 3628–3633.
  • Monkhorst and Pack (1976) H. J. Monkhorst, J. D. Pack, Phys. Rev. B 13 (1976) 5188–5192.
  • (35) http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
229066
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description