  Based on ab initio band structure calculations we formulate a general
theoretical method for description of the temperature dependence of electric
field gradient in solids. The method employs a procedure of averaging multipole
electron density component (l \neq 0) inside a sphere vibrating with the
nucleus at its center. As a result of averaging each Fourier component (K \neq
0) on the sphere is effectively reduced by the square root of the Debye-Waller
factor [exp(-W)]. The electric field gradient related to a sum of K-components
most frequently decreases with temperature (T), but under certain conditions
because of the interplay between terms of opposite signs it can also increase
with T. The method is applied to calculations of the temperature evolution of
the electric field gradients of pristine zinc and cadmium crystallized in the
hexagonal lattice. For calculations within our model of crucial importance is
the temperature dependence of mean-square displacements which can be taken from
experiment or obtained from the phonon modes in the harmonic approximation. For
the case of Zn we have used data obtained from single crystal x-ray
diffraction. In addition, for Zn and Cd we have calculated mean-square
displacements with the density functional perturbation treatment of the Quantum
Espresso package. With the experimental data for displacements in Zn our
calculations reproduce the temperature dependence of the electric field
gradient very accurately. Within the harmonic approximation of the Quantum
Espresso package the decrease of electric field gradients in Zn and Cd with
temperature is overestimated. Our calculations indicate that the anharmonic
effects are of considerable importance in the temperature dependence of
electric field gradients.
