ABINIT, anaddb input variables:

Govern the imposition of the Acoustic Sum Rule (ASR).

• 0 => no ASR for interatomic force constants is imposed.
• 1 or 2 => the ASR for interatomic force constants is imposed by modifying the on-site interatomic force constants, in a symmetric way (asr=2), or in the more general case, unconstrained way (asr=1).

More detailed explanations: the total energy should be invariant under translation of the crystal as a whole. This would guarantee that the three lowest phonon modes at Gamma have zero frequency (Acoustic Sum Rule - ASR). Unfortunately, the way the DDB is generated (presence of a discrete grid of points for the evaluation of the exchange-correlation potential and energy) slightly breaks the translational invariance. Well, in some pathological cases, the breaking can be rather important.

Two quantities are affected: the interatomic forces (or dynamical matrices), and the effective charges. The ASR for the effective charges is called the charge neutrality sum rule, and will be dealt with by the variable chneut. The ASR for the interatomic forces can be restored, by modifying the interatomic force of the atom on itself, (called self-IFC), as soon as the dynamical matrix at Gamma is known. This quantity should be equal to minus the sum of all interatomic forces generated by all others atoms (action-reaction law!), which is determined by the dynamical matrix at Gamma.

So, if asr is non-zero, the correction to the self-force will be determined, and the self-force will be imposed to be consistent with the ASR. This correction will work if IFCs are computed (ifcflag/=0), as well as if the IFCs are not computed (ifcflag==0). In both cases, the phonon frequencies will not be the same as the ones determined by the output of abinit, RF case. If you want to check that the DDB is correct, by comparing phonon frequencies from abinit and anaddb, you should turn off both asr and chneut.

Until now, we have not explained the difference between asr=1 and asr=2. This is rather subtle. In some local low-symmetry cases (basically the effective charges should be anisotropic), when the dipole-dipole contribution is evaluated and subtracted, the ASR cannot be imposed without breaking the symmetry of the on-site interatomic forces. That explains why two options are given: the second case (asr=2, sym) does not entirely impose the ASR, but simply the part that keeps the on-site interatomic forces symmetric (which means that the acoustic frequencies do not go to zero exactly), the first case (asr=1, asym) imposes the ASR, but breaks the symmetry. asr=2 is to be preferred for the analysis of the interatomic force constant in real space, while asr=1 should be used to get the phonon band structure.

The actual numbers of the atoms for which the interatomic force constant have to be written and eventually analysed.

Allow setting the target band gap, in eV. (elphflag@anaddb=1).

Allows to specify the Bravais lattice of the crystal, in order to help to generate a grid of special q points.

• 1 => all the lattices (including FCC, BCC and hexagonal)
• 2 => specific for Face Centered lattices
• 3 => specific for Body Centered lattices
• 4 => specific for the Hexagonal lattice

Note that in the latter case, the rprim of the unit cell have to be 1.0 0.0 0.0 -.5 sqrt(3)/2 0.0 0.0 0.0 1.0 in order for the code to work properly.

Set the treatment of the Charge Neutrality requirement for the effective charges.

• chneut=0 => no ASR for effective charges is imposed
• chneut=1 => the ASR for effective charges is imposed by giving to each atom an equal portion of the missing charge. See Eq.(48) in Phys. Rev. B55, 10355 (1997).
• chneut=2 => the ASR for effective charges is imposed by giving to each atom a portion of the missing charge proportional to the screening charge already present. See Eq.(49) in Phys. Rev. B55, 10355 (1997).

Integer. Frequency-dependent dielectric tensor flag.

• 0 => No dielectric tensor is calculated.
• 1 => The frequency-dependent dielectric tensor is calculated. Requirements for preceding response-function DDB generation run: electric-field and full atomic-displacement responses. Set rfstrs = 1, 2, or 3 (preferably 3). Set rfatpol and rfdir to do a full calculation of phonons at Q=0. The frequencies are defined by the nfreq,frmin,frmax variables. Also, the generalized Lyddane-Sachs-Teller relation will be used as an independent check of the dielectric tensor at zero frequency (this for the directions defined in the phonon list 2. See nph2l.
• 2 => Only the electronic dielectric tensor is calculated. It corresponds to a zero-frequency homogeneous field, with quenched atomic positions. For large band gap materials, this quantity is measurable because the highest phonon frequency is on the order of a few tenths of eV, and the band gap is larger than 5eV. Requirements for preceding response-function DDB generation: electric-field response. Set rfelfd = 1 or 3 (preferably 3). Note that the same information on the electronic dielectric tensor will be printed in the .out file of the rfelfd run.
• 3 => Compute and print the relaxed-ion dielectric tensor. Requirements for preceding response-function DDB generation run: electric-field and full atomic-displacement responses. Set rfstrs = 1, 2, or 3 (preferably 3). Set rfatpol and rfdir to do a full calculation of phonons at Q=0 (needed because the inverse of force-constant tensor is required). Furthermore, in the anaddb input file the variable nph2l must be nonzero in order to initiate computation of atomic displacements. If only the dielectric response is needed it is sufficient to set nph2l to 1 and leave qph2l at its default value (the Gamma point). Note that the relaxed-ion dielectric tensor computed here can also be obtained as the zero-frequency limit of the frequency-dependent dielectric tensor using input variables dieflag=1 and frmin=0.0. (The results obtained using these two approaches should agree to good numerical precision.) The ability to compute and print the static dielectric tensor here is provided for completeness and consistency with the other tensor quantities that are computed in this section of the code.
• 4 => Calculate dielectric tensor of both relaxed ion and free stress. We need information of internal strain and elastic tensor (relaxed ion) in this computation. So please set: elaflag=2,3,4 or 5 and instrflag=1

• 0 => the dipole-dipole interaction is not handled separately in the treatment of the interatomic forces. This option is available for testing purposes or if effective charge and/or dielectric tensor is not available in the derivative database. It gives results much less accurate than dipdip=1.
• 1 => the dipole-dipole interaction is subtracted from the dynamical matrices before Fourier transform, so that only the short-range part is handled in real space. Of course, it is reintroduced analytically when the phonon spectrum is interpolated, or if the interatomic force constants have to be analysed in real space.

• 0 => do not write the phonon eigenvectors;
• 1 or 2 => write the phonon eigenvectors;
• 4 => generate output files for band2eps (drawing tool for the phonon band structure);

Flag for calculation of elastic and compliance tensors

• 0 => No elastic or compliance tensor will be calculated.
• 1 => Only clamped-ion elastic and compliance tensors will be calculated. Requirements for preceding response-function DDB generation run: Strain perturbation. Set rfstrs to 1, 2, or 3. Note that rfstrs=3 is recommended so that responses to both uniaxial and shear strains will be computed.
• 2 => Both relaxed- and clamped-ion elastic and compliance tensor will be calculated, but only the relaxed-ion quantities will be printed. The input variable instrflag should also be set to 1, because the internal-strain tensor is needed to compute the relaxed-ion corrections. Requirements for preceding response-function DDB generation run: Strain and atomic-displacement responses at Q=0. Set rfstrs = 1, 2, or 3 (preferably 3). Set rfatpol and rfdir to do a full calculation of phonons at Q=0 (needed because the inverse of force-constant tensor is required).
• 3 => Both relaxed and clamped-ion elastic and compliance tensors will be printed out. The input variable instrflag should also be set to 1. Requirements for preceding response-function DDB generation run: Same as for elaflag=2.
• 4 => Calculate the elastic and compliance tensors (relaxed ion) at fixed displacement field, the relaxed-ion tensors at fixed electric field will be printed out too, for comparison. When elaflag=4, we need the information of internal strain and relaxed-ion dielectric tensor to build the whole tensor, so we need set instrflag=1 and dieflag=3 or 4 .
• 5 => Calculate the relaxed ion elastic and compliance tensors, considering the stress left inside cell. At the same time, bare relaxed ion tensors will still be printed out for comparison. In this calculation, stress tensor is needed to compute the correction term, so one supposed to merge the first order derivative data base (DDB file) with the second order derivative data base (DDB file) into a new DDB file, which can contain both information. And the program will also check for the users.

Give the energy for the phonon frequency output (in the output file, not in the console log file, for which Hartree units are used).

• 0 => Hartree and cm-1;
• 1 => meV and Thz;
• 2 => Hartree, cm-1, meV, Thz, and Kelvin.

• 0 => no analysis of interatomic force constants;
• 1 => analysis of interatomic force constants.

• 0 => do all calculations directly from the DDB, without the use of the interatomic force constant.
• 1 => calculate and use the interatomic force constants for interpolating the phonon spectrum and dynamical matrices at every q wavevector, and eventually analyse the interatomic force constants, according to the informations given by atifc,dipdip,ifcana,ifcout,natifc,nsphere,rifcsph

Internal strain tensor flag.

• 0 => No internal-strain calculation.
• 1 => Print out both force-response and displacement-response internal-strain tensor. Requirements for preceding response-function DDB generation run: Strain and full atomic-displacement responses. Set rfstrs = 1, 2, or 3 (preferably 3). Set rfatpol and rfdir to do a full calculation of phonons at Q=0.

Indices of the elements of the strain tensor that are fixed during a structural relaxation at constrained polarisation:

• 0 => No elastic or compliance tensor will be calculated.
• 1 => Only clamped-ion elastic and compliance tensors will be calculated. Requirements for preceding response-function DDB generation run: Strain perturbation. Set rfstrs to 1, 2, or 3. Note that rfstrs>=3 is recommended so that responses to both uniaxial and shear strains will be computed.
• 2 => Both relaxed- and clamped-ion elastic and compliance tensor will be calculated, but only the relaxed-ion quantities will be printed. The input variable instrflag should also be set to 1, because the internal-strain tensor is needed to compute the relaxed-ion corrections. Requirements for preceding response-function DDB generation run: Strain and atomic-displacement responses at Q=0. Set rfstrs = 1, 2, or 3 (preferably 3). Set rfatpol and rfdir to do a full calculation of phonons at Q=0 (needed because the inverse of force-constant tensor is required).
• 3 => Both relaxed and clamped-ion elastic and compliance tensors will be printed out. The input variable instrflag should also be set to 1. Requirements for preceding response-function DDB generation run: Same as for elaflag=2'.
• 4 => Calculate the elastic and compliance tensors (relaxed ion) at fixed displacement field, the relaxed-ion tensors at fixed electric field will be printed out too, for comparison. When elaflag=4, we need the information of internal strain and relaxed-ion dielectric tensor to build the whole tensor, so we need set instrflag=1 and dieflag=3 or 4 .
• 5 => Calculate the relaxed ion elastic and compliance tensors, considering the stress left inside cell. At the same time, bare relaxed ion tensors will still be printed out for comparison. In this calculation, stress tensor is needed to compute the correction term, so one supposed to merge the first order derivative data base (DDB file) with the second order derivative data base (DDB file) into a new DDB file, which can contain both information. And the program will also check for the users.

Non-linear properties flag.

• 0 => do not compute non-linear properties ;
• 1 => the electrooptic tensor, Raman susceptibilities and non-linear optical susceptibilities are calculated;
• 2 => only the non-linear optical susceptibilities and first-order changes of the dielectric tensor induced by an atomic displacement are calculated;
• 3 => only the non-linear optical susceptibility is calculated.

The number of wavevectors in phonon list 2, defining the directions along which the non-analytical splitting of phonon frequencies at Gamma will be calculated. The actual values of the wavevector directions will be specified by qph2l. These are actually all wavectors at Gamma, but obtained by a limit along a different direction in the Brillouin-zone. It is important to note that non-analyticities in the dynamical matrices are present at Gamma, due to the long-range Coulomb forces. So, going to Gamma along different directions can give different results.

Number of atoms included in the cut-off sphere for interatomic force constant, see also the alternative rifcsph. If nsphere= 0: maximum extent allowed by the grid. If nsphere= -1: the code analyzes different values of nsphere and find the value that does not lead to unstable frequencies in a small sphere around Gamma. The truncated IFCs are then used for further post-processing. The results of the test are reported in the main output file. This option is useful to obtain a initial guess of nsphere: the value that leads to stable frequencies and gives linear dispersion for the acoustic modes around Gamma is usually smaller that the one reported by nsphere -1.

Flag for calculation of piezoelectric tensors

• 0 => No piezoelectric tensor will be calculated.
• 1 => Only the clamped-ion piezoelectric tensor is computed and printed. Requirements for preceding response-function DDB generation run: Strain and electric-field responses. For the electric-field part, one needs results from a prior 'ddk perturbation' run. Note that even if only a limited number of piezoelectric tensor terms are wanted (as determined by rfstrs and rfdir in this calculation) it is necessary to set rfdir = 1 1 1 in the d/dk calculation for most structures. The only obvious exception to this requirement is cases in which the primitive lattice vectors are all aligned with the cartesian axes. The code will omit terms in the output piezoelectric tensor for which the available d/dk set is incomplete. Thus: Set rfstrs to 1, 2, or 3i (preferably 3)
• 2 => Both relaxed- and clamped-ion elastic and compliance tensor will be calculated, but only the relaxed-ion quantities will be printed. The input variable instrflag should also be set to 1, because the internal-strain tensor is needed to compute the relaxed-ion corrections. Requirements for preceding response-function DDB generation run: Strain, electric-field and full atomic-displacement responses at Q=0. Set rfstrs = 1, 2, or 3 (preferably 3). Set rfelfd = 3. Set rfatpol and <rfdir to do a full calculation of phonons at Q=0 (needed because the inverse of force-constant tensor is required).
• 3 => Both relaxed and clamped-ion piezoelectric tensors will be printed out. The input variable instrflag should also be set to 1. Requirements for preceding response-function DDB generation run: Same as for piezoflag=2.
• 4 => Calculate the piezoelectric d tensor (relaxed ion). In order to calculate the piezoelectric d tensor, we need information of internal strain and elastic tensor (relaxed ion). So we should set elaflag= 2,3,4, or 5 and instrflag=1. The subroutine will also do a check for you, and print warning message without stopping even if flags were not correctly set.
• 5 => Calculate the piezoelectric g tensor (relaxed ion). In this computation, we need information of internal strain, elastic tensor (relaxed ion) and dielectric tensor (relaxed ion). So we should set: instrflag=1, elaflag=2,3,4 or 5, dieflag=3 or 4. The subroutine will also do a check for you, and print warning message without stopping even if flags were not correctly set.
• 6 => Calculate the piezoelectric h tensor (relaxed ion). In this calculation, we need information of internal strain and dielectric tensor (relaxed ion). So we need set: instrflag=1 and dieflag=3 or 4. The subroutine will also do a check for you, and print warning message without stopping even if flags were not correctly set.
• 7 => calculate all the possible piezoelectric tensors, including e (clamped and relaxed ion), d, g and h tensors. The flags should be set to satisfy the above rules from 1 to 6.

If activated, compute polarization in cartesian coordinates, and update lattice constants and atomic positions in order to perform a structural optimization at constrained polarization.

More detailed explanation: ANADDB can use the formalism described in Na Sai et al, PRB 66, 104108 (2002), to perform structural relaxations under the constraint that the polarization is equal to a value specified by the input variable targetpol. The user starts from a given configuration of a crystal and performs a ground-state calculation of the Hellman-Feynman forces and stresses and the Berry phase polarization as well as a linear response calculation of the whole matrix of second-order energy derivatives with respect to atomic displacement, strains and electric field.
In case polflag=1, ANADDB solves the linear system of equations (13) of the Na Sai paper, and computes new atomic positions (if relaxat=1) and lattice constant (if relaxstr=1). Then, the user uses these parameters to perform a new ground-state and linear-response calculation. This must be repeated until convergence is reached. The user can also fix some atomic positions, or strains, thanks to the input variables natfix, nstrfix,iatfix,istrfix.
In case both relaxat and relaxstr are 0, while polflag=1, ANADDB only computes the polarization in cartesian coordinates.

Flag to print out the Interatomic Force Constants in real space to a file. The available options are:

• 0 => do nothing (IFC are printed to the log file);
• 1 => write out the IFC in file ifcinfo.out (the name is fixed) to be used by AI2PS from John Rehr's group

• 0 => do not write the BoltzTraP input files;
• 1 => write out the input files for BoLTZTRaP code.

Flag to print out the DDB file interpolated with the Interatomic Force Constants.

The available options are:

• 0 => no output of DDB (default);
• 1 => Interpolate the DDB and write out the DDB and DDB.nc files.

The prtdos variable is used to calculate the phonon density of states, PHDOS, by Fourier interpolating the interatomic force constants on the (dense) q-mesh defined by ng2qpt. Note that the variable ifcflag must be set to 1 since the interatomic force constants are supposed to be known.

The available options are:

• 0 => no output of PHDOS (default);
• 1 => calculate PHDOS using the gaussian method and the broadening defined by dossmear.

Only for electron-phonon calculations. The available options are:

• 0 => do not write the Fermi Surface;
• 1 => write out the Fermi Surface in the BXSF format used by XCrySDen.

Further comments:

• 0 => do not write the mode-by-mode decomposition of the electrooptic tensor;
• 1 => write out the contribution of the individual zone-center phonon modes to the electrooptic tensor.

Only for electron-phonon calculations. This input variable is used to calculate the nesting function defined as: \chi_{nm}(q) = \sum_k \delta(\epsilon_{k,n}-epsilon_F) \delta(\epsilon_{k+q,m}-\epsilon_F). The nesting factor is calculated for every point of the k-grid employed during the previous GS calculation. The values are subsequently interpolated along the trajectory in q space defined by qpath, and written in the _NEST file using the X-Y format (prtnest=1). It is also possible to analyze the behavior of the function in the reciprocal unit cell saving the values in the NEST_XSF file that can be read using XCrySDen (prtnest=2). Note that in the present implementation what is really printed to file is the "total nesting" defined as \sum_{nm} \chi_{nm}(q). Limitations: the k-grid defined by kptrlatt must be orthogonal in reciprocal space, moreover off-diagonal elements are not allowed, i.e kptrlatt 4 0 0 0 4 0 0 0 4 is fine while kprtlatt = 1 0 0 0 1 1 0 -1 1 will not work.

• 0 => do not write the nesting function;
• 1 => write only the nesting function along the q-path in the X-Y format;
• 2 => write out the nesting function both in the X-Y and in the XSF format.

Only if ifcflag=1. This option specifies the file format for the phonon band structure. Possible values:

• 1 Write frequencies in xmgrace format. A file with extension `PHBANDS.agr` is produced. Use `xmgrace file_PHBANDS.agr` to visualize the data
• 2 Write frequencies in gnuplot format. The code produces a `PHBANDS.dat` file with the eigenvalues and a `PHBANDS.gnuplot` script. Use `gnuplot file_PHBANDS.gnuplot` to visualize the phonon band structure.

Only if ifcflag=1. The available options are:

• 0 => do not write the SR/LR decomposition of phonon frequencies;
• 1 => write out the SR/LR decomposition of the square of phonon frequencies at each q-point specified in qph1l.

This vector gives the shifts needed to define the coarse q-point grid.

a) Case nqshft=1 In general, 0.5 0.5 0.5 with the ngqpt's even will give very economical grids. On the other hand, is it sometimes better for phonons to have the Gamma point in the grid. In that case, 0.0 0.0 0.0 should be OK. For the hexagonal lattice, the above mentioned quantities become 0.0 0.0 0.5 and 0.0 0.0 0.0 .

b) Case nqshft=2 The two q1shft vectors must form a BCC lattice. For example, use 0.0 0.0 0.0 and 0.5 0.5 0.5

c) Case nqshft=4 The four q1shft vectors must form a FCC lattice. For example, use 0.0 0.0 0.0 , 0.0 0.5 0.5 , 0.5 0.0 0.5 , 0.5 0.5 0.0 or 0.5 0.5 0.5 , 0.0 0.0 0.5 , 0.0 0.5 0.0 , 0.5 0.0 0.0 (the latter is referred to as shifted)

Further comments: by using this technique, it is possible to increase smoothly the number of q-points, at least less abruptly than relying on series of grids like (for the full cubic symmetry):
1x1x1 => (0 0 0)
2x2x2 (shifted) => (.25 .25 .25)
2x2x2 => 1x1x1 + (.5 0 0) (.5 .5 0) (.5 .5 0)
4x4x4 => 2x2x2 + (.25 0 0) (.25 .25 0) (.25 .5 0) (.25 .25 .25) (.25 .25 .5) (.25 .5 .5)
...

with respectively 1, 1, 4 and 10 q-points, corresponding to a number of points in the full BZ of 1, 8, 8 and 64. Indeed, the following grids are made available:
1x1x1 with nqshft=2 => (0 0 0) (.5 .5 .5)
1x1x1 with nqshft=4 => (0 0 0) (.5 .5 0)
1x1x1 with nqshft=4 (shifted) => (.5 0 0) (.5 .5 .5)
2x2x2 with nqshft=2 => 2x2x2 + (.25 .25 .25)
2x2x2 with nqshft=4 => 2x2x2 + (.25 .25 0) (.25 .25 .5)
2x2x2 with nqshft=4 (shifted) => (.25 0 0) (.25 .25 .25) (.5 .5 .25) (.25 .5 0)
...

with respectively 2, 2, 2, 5, 6 and 4 q-points, corresponding to a number of points in the full BZ of 2, 4, 4, 16, 32 and 32.

For a FCC lattice, it is possible to sample only the Gamma point by using a 1x1x1 BCC sampling (nqshft=2).

Similar to q1shft, but for the series of fine grids.

List of phonon wavevector directions along which the non-analytical correction to the Gamma-point phonon frequencies will be calculated (for insulators). Four numbers, as for qph1l, but where the last one, that correspond to the normalisation factor, is 0.0 For the anaddb code, this has the meaning that the three previous values define a direction. The direction is in CARTESIAN COORDINATES, unlike the non-Gamma wavevectors defined in the first list of vectors...

Note that if the three first numbers are zero, then the code will do a calculation at Gamma without non-analyticities.

Govern the imposition of the sum-rule on the Raman tensors.
As in the case of the Born effective charges, the first-order derivatives of the linear dielectric susceptibility with respect to an atomic displacement must vanish when they are summed over all atoms. This sum rule is broken in most calculations. By putting ramansr equal to 1 or 2, this sum rule is imposed by giving each atom a part of the discrepancy.

• 0 => no sum rule is imposed;
• 1 => impose the sum rule on the Raman tensors, giving each atom an equal part of the discrepancy;
• 2 => impose the sum rule on the Raman tensors, giving each atom a part of the discrepancy proportional to the magnitude of its contribution to the Raman tensor.

Select a particular set of Data Blocks in the DDB. (PRESENTLY, ONLY OPTION 1 IS AVAILABLE)

• 1 => Blocks obtained by a non-stationary formulation.
• 2 => Blocks obtained by a stationary formulation.

Cut-off radius for the sphere for interatomic force constant, see also the alternative nsphere. If rifcsph= 0: maximum extent allowed by the grid .

Select some parts of the effective charge tensor. (This is done after the application or non-application of the ASR for effective charges). The transformed effective charges are then used for all the subsequent calculations.

• 0 => The effective charge tensor is left as it is.
• 1 => For each atom, the effective charge tensor is made isotropic, by calculating the trace of the matrix, dividing it by 3, and using this number in a diagonal effective charge tensor.
• 2 => For each atom, the effective charge tensor is made symmetric, by simply averaging on symmetrical elements.

Note: this is for analysis the effect of anisotropy in the effective charge. The result with non-zero selectz are unphysical.

Flag controlling the Fermi surface integration technique used for electron-phonon quantities.

• 0 = tetrahedron method (no adjustable parameter)
• 1 = Gaussian smearing (see elphsmear)
• 2 = uniformly weighted band window between ep_b_min and ep_b_max, for all k-points

Lowest temperature (Kelvin) at which the thermodynamical quantities have to be evaluated. Cannot be zero when thmflag is 1.

Thermal_supercell defines the real space supercell in which a thermalized atomic configuration should be produced, following the prescription of Zacharias and Giustino (PRB 94 075125 (2016)). The displacements are chosen for each phonon mode according to a temperature, and the displacements are alternated in sign/direction to obtain maximal compensation of the linear electron phonon coupling. In this way in the PRB dielectric properties at finite T can be obtained from a single supercell calculation instead of lots of MD and configuration averaging.

The supercell vectors are not constrained to be collinear with the normal lattice vectors: this effect is obtained by using a diagonal Thermal_supercell. The lines of the matrix describe the linear combination of the primitive cell lattice vectors yielding the supercell vectors, as for kptrlatt.

For the moment this feature is under development and it looks like the relative phases of the displacements are not fixed properly yet... (Aug 2017)

Flag controlling the calculation of thermal quantities.

• When thmflag == 1, the code will compute, using the histogram method:
  • the normalized phonon DOS
  • the phonon internal energy, free energy, entropy, constant volume heat capacity as a function of the temperature
  • the Debye-Waller factors (tensors) for each atom, as a function of the temperature
  • the mean-square velocity tensor for each atom, as a function of temperature
  • the "average frequency" as a function of the temperature
• When thmflag == 2, all the phonon frequencies for the q points in the second grid are printed.
• When thmflag == 3, 5 or 7, the thermal corrections to the electronic eigenvalues are calculated. If thmflag==3, the list of phonon wavevector from the first list is used (with equal weight for all wavevectors in this list), while if thmflag==5 or 7, the first grid of wavevectors is used, possibly folded to the irreducible Brillouin Zone if symmetry operations are present, or if they are recomputed (this happens for thmflag==7).
• When thmflag == 4 or 6, the temperature broadening (electron lifetime) of the electronic eigenvalues is calculated. If thmflag==4, the list of phonon wavevector from the first list is used (with equal weight for all wavevectors in this list), while if thmflag==6, the first grid of wavevectors is used, possibly folded to the irreducible Brillouin Zone if symmetry operations are present or if they are recomputed (this happens for thmflag==8).

This variable activates the calculation of the speed of sound (requires ifcflag = 1). The first entry of the array defines the radius of the small sphere around the Gamma point (Bohr-1). The second entry gives the absolute tolerance in kilometer/second. The speed of sound is evaluated by performing a spherical average on the small sphere using Lebedev-Laikov grids (typical values for q-radius: 0.1 Bohr-1) The number of radial points is increased until the integration converges twice withing the tolerance specified by the user (typical values for tolkms: 0.05 km/s).