Controls the range of atoms for which displacements will be considered in non-
linear computations (using the 2n+1 theorem), for the 1st perturbation.
May take values from 1 to natom, with d3e_pert1_atpol (1)<=
d3e_pert1_atpol (2). See rfatpol for additional details.

Gives the directions to be considered in non-linear computations (using the
2n+1 theorem), for the 1st perturbation.
The three elements corresponds to the three primitive vectors, either in real
space (atomic displacement), or in reciprocal space (electric field perturbation).
See rfdir for additional details.

Turns on electric field perturbation in non-linear computation, as 1st
perturbation. Actually, such calculations requires first the non-self-
consistent calculation of derivatives with respect to k, independently of the
electric field perturbation itself. See rfelfd for additional details.

Controls the range of atoms for which displacements will be considered in non-
linear computations (using the 2n+1 theorem), for the 2nd perturbation.
May take values from 1 to natom, with d3e_pert2_atpol (1) <= d3e_pert2_atpol (2).
See rfatpol for additional details.

Gives the directions to be considered in non-linear computations (using the
2n+1 theorem), for the 2nd perturbation.
The three elements corresponds to the three primitive vectors, either in real
space (atomic displacement), or in reciprocal space (electric field perturbation).
See rfdir for additional details.

Turns on electric field perturbation in non-linear computation, as 2nd
perturbation. Actually, such calculations requires first the non-self-
consistent calculation of derivatives with respect to k, independently of the
electric field perturbation itself. See rfelfd for additional details.

Controls the range of atoms for which displacements will be considered in non-
linear computations (using the 2n+1 theorem), for the 3rd perturbation.
May take values from 1 to natom, with d3e_pert3_atpol (1)<= d3e_pert3_atpol (2).
See rfatpol for additional details.

Gives the directions to be considered in non-linear computations (using the
2n+1 theorem), for the 3rd perturbation.
The three elements corresponds to the three primitive vectors, either in real
space (atomic displacement), or in reciprocal space (electric field perturbation).
See rfdir for additional details.

Turns on electric field perturbation in non-linear computation, as 3rd
perturbation. Actually, such calculations requires first the non-self-
consistent calculation of derivatives with respect to k, independently of the
electric field perturbation itself.
See rfelfd for additional details.

It is the value of the “scissors operator”, the shift of conduction band
eigenvalues, used in response function calculations.
Can be specified in Ha (the default), Ry, eV or Kelvin, since ecut has the
‘ENERGY‘ characteristics (1 Ha=27.2113845 eV).
Typical use is for response to electric field (rfelfd = 3), but NOT for d/dk
(rfelfd = 2) and phonon responses.

Turns on effective mass tensor calculations. Such calculations requires the
non-self-consistent calculation of derivatives with respect to k, in the same
dataset. It must therefore be used with rfelfd = 2 (or 1).

Note

At the present time, both norm-conserving (NC) and PAW calculations are
supported. Also, for PAW calculations only, nspinor == 2 and
pawspnorb == 1 (i.e. spin-orbit (SO) calculations) is supported. NC SO
calculations are NOT currently supported. Also, for both NC and PAW,
nspden/=1 and nsppol/=1 are NOT supported.

Mnemonics: EFfective MASs, BANDS to be treated.Mentioned in topic(s):topic_EffectiveMassVariable type: integerDimensions: (2,nkpt)Default value: The full range of band available in the calculation for each k-point.Only relevant if:efmas == 1

This variable controls the range of bands for which the effective mass is to
be calculated. If a band is degenerate, all other bands of the degenerate
group will automatically be treated, even if they were not part of the user specified range.

Allows the user to calculate the scalar effective mass of all bands specified
by efmas_bands along specific directions in reciprocal space. This is
particularly useful when considering degenerate bands, which are usually
warped, and thus cannot have their dispersion (hessian) and effective mass
expressed as a tensor. This allows the user to see the more complex angular
behavior of effective masses in these cases, for instance.

When efmas_calc_dirs == 0, no directions are read from the input file (using
efmas_dirs) and the effective masses along the 3 cartesian directions are
output by default.

When efmas_calc_dirs == 1, 2 or 3, efmas_n_dirs directions are read from
efmas_dirs, assuming cartesian, reduced or angular (\theta,\phi)
coordinates, respectively. In the case efmas_calc_dirs == 3, 2 real values
per directions are read, whereas 3 real values are read in the two other cases.

Activate (==1) or not (==0) the treatment of degenerate bands
(criterion efmas_deg_tol is used to determine whether bands are degenerate).
Also compute the transport equivalent effective mass (see [Mecholsky2014]).

For 2D or 1D systems, the band dispersion goes to 0 perpendicular to the
system, which causes the inverse effective mass to be singular, i.e. the
effective mass to be NaN. This keyword circumvents the problem by eliminating
the troublesome dimensions from the inverse effective mass.

In 2D, the Z axis is ignored and, in 1D, the Z and Y axis are ignored.

Also, note that in the 2D degenerate case, a subtlety arises: the ‘transport
equivalent’ effective mass does not determine the scale of the transport
tensors (conductivity and others). Therefore, for this specific case, the
factor by which these transport tensors should be scaled once determined from
the ‘transport equivalent’ effective mass tensor is output separately on the
line immediately after the effective mass.

When a band is degenerate, the usual definition of effective mass becomes
invalid. However, it is still possible to define a ‘transport equivalent mass
tensor’ that reproduces the contribution of the band to the conductivity
tensor. To obtain this tensor, an integration over the solid sphere is
required. The angular variables are sampled using efmas_ntheta points for the theta coordinate,
and twice efmas_ntheta points for the phi coordinate.
The default value gives a tensor accurate to the 4th decimal in Ge.

The variable elph2_imagden determines the imaginary shift of the
denominator of the sum-over-states in the perturbation,
(e_{nk}-e_{n'k'}+ielph2_imagden). One should use a width comparable with
the Debye frequency or the maximum phonon frequency.
Can be specified in Ha (the default), Ry, eV or Kelvin, since ecut has the
‘ENERGY‘ characteristics (1 Ha=27.2113845 eV).

The variable esmear determines the width of the functions approximating
the delta function, \delta(e_{nk}-e_{n'k'}), present in the expression of the
lifetimes. One should use a width comparable with the Debye frequency or the
maximum phonon frequency.
Can be specified in Ha (the default), Ry, eV or Kelvin, since ecut has the
‘ENERGY‘ characteristics (1 Ha=27.2113845 eV).

Can be used to suppress artificially the first-order change of Fermi energy,
in case of Response Function calculation for metals at Q=0. If the input variable
frzfermi is set to 1, this contribution is suppressed, even though this is incorrect.

If ieig2rf is greater then 0, the code will produce a file, named with the
trailing suffix _EIGR2D, containing the second-order electronic eigenvalues
for the perturbation. These files are used in the calculation of the thermal
correction to the electronic eigenvalues.

If ieig2rf is set to 1, the second-order electronic eigenvalues will be
calculated from the DFPT method (Sternheimer).

If ieig2rf is set to 2, the second-order electronic eigenvalues will be
calculated from the Allen-Cardona method (sum over states).

If ieig2rf is set to 3, the second-order electronic eigenvalues will be
calculated from the DFPT method (sum over states) but using a different part
of the code. This is equivalent to ieig2rf = 1 [debuging].

If ieig2rf is set to 4, the second-order electronic eigenvalues will be
calculated from the dynamical DFPT method (Sternheimer). The code will
generate _EIGR2D.nc files that contain the electron-phonon matrix element
squared on the space orthogonal to the active space. The code will also
produce _FAN.nc files that contain the electron-phonon matrix elements
squared.

If ieig2rf is set to 5, the second-order electronic eigenvalues will be
calculated from the dynamical DFPT method (Sternheimer). The code will
generate _EIGR2D.nc files that contain the electron-phonon matrix element
square on the space orthogonal to the active space. The code will also produce
_GKK.nc files that contain electron-phonon matrix elements. This option is
preferable for large system to ieig2rf = 4 as the GKK files take much
less disk space and memory (but run a little bit slower).

Note

ieig2rf = 4 and 5 can only be used if Abinit is compiled with NETCDF support.

Method of calculation of the 1st order XC potential in non-collinear DFPT
calculations. The possible values 1,2 and 3 correspond to the following
methods

If ixcrot=1, the spinor rotation matrix U at each fft point is not calculated explicitly. Instead the needed expressions involving U are derived based on the general properties of the U matrix.

If ixcrot=2, U is computed explicitly

If ixcrot=3, the brute force evaluation of the 1st order XC potential as a functional derivative is used. Rotation matrices are not computed.

In theory, all methods give identical results. However, due to different
implementation approaches, the round-off errors can lead to slight differences
intermediate and final results obtained using methods 1,2 and 3. The choice of
the method can also affect the convergence.

Note

For non-zero perturbation wavevector, only ixcrot=3 implementation is currently available.

Control the output of the non-linear implementation (only when usepead == 0).
The default value, nonlinear_info == 0 does nothing. If nonlinear_info == 1,
different contributions of 3rd derivatives of the energy are written in the
output file (non time consuming).

Higher values activate some internal tests for
checking the implementation correctness (time consuming, not useable in parallel).
If nonlinear_info == 2, same effect than 1 and tests are done in non-linear
(optdriver==5 and usepead == 0).
If nonlinear_info == 3, same effect than 1 and tests are done in rf2_init
(rf2_dkdk/=0 or rf2_dkde/=0).
If nonlinear_info == 4, same effect than 1 and tests are done in both non-linear and rf2_init.
A line containining “NOT PASSED” (and other information) is added to the output file
for each test that does not pass, otherwise nothing is printed. However, more information concerning
the tests is always printed in the standard output file.

The computation of third-order derivatives from the 2n+1 theorem requires the
first-order wavefunctions and densities obtained from a linear response
calculation. The standard approach in a linear response calculation is:

compute only the irreducible perturbations;

use symmetries to reduce the number of k-points for the k-point integration.

This approach cannot be applied, presently (v4.1), if the first-order
wavefunctions are to be used to compute third-order derivatives. First, for
electric fields, the code needs the derivatives along the three directions.
Still, in case of phonons, only the irreducible perturbations are required.
Second, for both electric fields and phonons, the wavefunctions must be
available in half the BZ (kptopt=2), or the full BZ (kptopt=3).
During the linear response calculation, in order to prepare a non-linear
calculation, one should put prepanl to 1 in order to force ABINIT to
compute the electric field perturbation along the three directions explicitly,
and to keep the full number of k-points.

In the case of a 2nd derivative of wavefunction (rf2_dkdk or rf2_dkde),
prepanl == 1 can be used in order to skip directions of perturbations that
will not be used by the non-linear routine (see rf2_dkdk for more details).

The calculation of electron-phonon coupling quantities requires the presence
of all the perturbations (all atoms in all directions) for the chosen set of
(irreducible) q-points. To impose this and prevent ABINIT from using symmetry
to reduce the number of perturbations, set prepgkk to 1. Use in
conjunction with prtgkk.

If 1, at the end of an effective mass calculation (efmas = 1), create a file *_EFMAS, that contains the generalized second-order k-derivatives, see Eq.(66) in [Laflamme2016], in view of further processing.

If set to 1, the output _1WF files will contain the full 1st-order wavefunctions, for both valence and conduction bands.
Otherwise, the _1WF files are not really 1st-order perturbed wavefunctions, but merely a set of perturbed wavefunctions that yield the correct perturbed density.
This is used when one expect to perform post-processing of the 1st-order wavefunctions.

If is equal to 1, activates computation of mixed second derivatives of wavefunctions with respect to
wavevector and electric field (ipert = natom+11 is activated). This is not strictly a response function
but is a needed auxiliary quantity in the calculations of 3rd-order derivatives of the energy
(non-linear response) if usepead == 0. The directions for the derivatives are determined by
rf2_pert1_dir, rf2_pert2_dir and prepanl in the same way than rf2_dkdk.

If is equal to 1, activates computation of second derivatives of wavefunctions with respect to
wavevectors (ipert = natom+10 is activated). This is not strictly a response function but is a needed
auxiliary quantity in the calculations of 3rd-order derivatives of the energy
(non-linear response) if usepead == 0. The directions for the derivatives are determined by
rf2_pert1_dir and rf2_pert2_dir and prepanl as the following:

The computation of the 2nd derivative of wavefunction with respect to “lambda_1” and “lambda_2” is computed if
if rf2_pert1_dir[idir1] AND rf2_pert2_dir[idir2] are equal to 1, where “idir1” (“idir2”) is direction of
the perturbation “lambda_1” (“lambda_2”).
If ALL directions are activated (default behavior) AND prepanl == 1, then the code automatically selects
only the directions that will be used by the non-linear routine (optdriver == 5) using crystal symmetries.

Gives the directions of the 1st perturbation to be considered when solving the 2nd order Sternheimer equation.
The three elements corresponds to the three primitive vectors, either in real
space (phonon calculations), or in reciprocal space (\,d/ \,d k, homogeneous
electric field, homogeneous magnetic field calculations).
If equal to 1, the 2nd order wavefunctions, as defined by rf2_dkdk or rf2_dkde, are computed for the
corresponding direction. If 0, this direction is not considered.
See rf2_dkdk for more details.

Gives the directions of the 2nd perturbation to be considered when solving the 2nd order Sternheimer equation.
The three elements corresponds to the three primitive vectors, either in real
space (phonon calculations), or in reciprocal space (\,d/ \,d k, homogeneous
electric field, homogeneous magnetic field calculations).
If equal to 1, the 2nd order wavefunctions, as defined by rf2_dkdk or rf2_dkde, are computed for the
corresponding direction. If 0, this direction is not considered.
See rf2_dkdk for more details.

Control the evaluation of the acoustic sum rule in effective charges and
dynamical matrix at Gamma within a response function calculation (not active
at the level of producing the DDB, but at the level of the phonon
eigenfrequencies output).

Control the range of atoms for which displacements will be considered in
phonon calculations (atomic polarizations).
These values are only relevant to phonon response function calculations.
May take values from 1 to natom, with rfatpol(1)<=rfatpol(2).
The atoms to be moved will be defined by the
do-loop variable iatpol:

For the calculation of a full dynamical matrix, use rfatpol(1)=1 and
rfatpol(2)=natom, together with rfdir 1 1 1. For selected
elements of the dynamical matrix, use different values of rfatpol and/or
rfdir. The name ‘iatpol’ is used for the part of the internal variable
ipert when it runs from 1 to natom. The internal variable ipert can also
assume values larger than natom, denoting perturbations of electric field
or stress type (see the DFPT help file).

Activates computation of derivatives of ground state wavefunctions with
respect to wavevectors. This is not strictly a response function but is a
needed auxiliary quantity in the electric field calculations (see rfelfd).
The directions for the derivatives are determined by rfdir.

0 → no derivative calculation

1 → calculation of first derivatives of wavefunctions with respect to k points (\,d/ \,d k calculation).
The exact same functionality is provided by rfelfd = 2.

Gives the directions to be considered for response function calculations (also
for the Berry phase computation of the polarization, see the berryopt
input variable).
The three elements corresponds to the three primitive vectors, either in real
space (phonon calculations), or in reciprocal space (\,d/ \,d k, homogeneous
electric field, homogeneous magnetic field calculations). So, they generate a
basis for the generation of the dynamical matrix or the macroscopic dielectric
tensor or magnetic susceptibility and magnetic shielding, or the effective charge tensors.
If equal to 1, response functions, as defined by rfddk, rfelfd,
rfphon, rfdir and rfatpol, are to be computed for the
corresponding direction. If 0, this direction should not be considered.

Turns on electric field response function calculations. Actually, such
calculations requires first the non-self-consistent calculation of derivatives
with respect to k, independently of the electric field perturbation itself.

0 → no electric field perturbation

1 → full calculation, with first the derivative of ground-state wavefunction with respect to k (\,d / \,d k calculation), by a non-self-consistent calculation, then the generation of the first-order response to an homogeneous electric field

2 → only the derivative of ground-state wavefunctions with respect to k

3 → only the generation of the first-order response to the electric field, assuming that the data on derivative of ground-state wavefunction with respect to k is available on disk.

Note

Because the tolerances to be used for derivatives or homogeneous
electric field are different, one often does the calculation of derivatives in
a separate dataset, followed by calculation of electric field response as well
as phonon.
The options 2 and 3 proves useful in that context; also, in case a scissor
shift is to be used, it is usually not applied for the \,d / \,d k response).

rfmagn allows one to run response function calculations with respect to
external magnetic field if set to 1. Currently, orbital magnetism is not taken into
account and the perturbing potential has Zeeman form.

Selects method used in response function calculations. Presently, only abs(rfmeth) = 1 is
allowed. This corresponds to storing matrix elements of the 2DTE computed using non-stationary expressions,
instead of stationary ones.

The difference between positive and negative values is rather technical. Very often, the symmetries can be used in such a way that
some matrix elements can be proven to be zero even without doing any computation. Positive values of rfmeth activate
this use of symmetries, while it is denied when rfmeth is negative. There is an indirect additional outcome of this,
as a symmetrization of the whole 2DTE is sometimes rendered possible when the additional knowledge of the zero matrix elements
is available. Thus, the results obtained for positive and negative values of rfmeth might slightly differ for non-zero elements of the 2DTE,
if they are computed in both cases.

Available to the developpers, to activate the use of ipert=natom+6 and
ipert=natom+7, two sets of perturbations that the developpers can define.

0 → no computations for ipert=natom+6 or ipert=natom+7

1 → response with respect to perturbation natom+6 will be computed

2 → response with respect to perturbation natom+7 will be computed

3 → responses with respect to perturbations natom+6 and natom+7 will be computed

Important

In order to define and use correctly the new perturbations, the developper
might have to include code lines or additional routines at the level of the
following routines: dfpt_cgwf.F90, dfpt_dyout.F90, dfpt_symph.F90,
dfpt_dyout.F90, dfpt_etot.F90, littlegroup_pert.F90, dfpt_looppert.F90,
dfpt_mkcor.F90, dfpt_nstdy.F90, dfpt_nstwf.F90, respfn.F90, dfpt_scfcv.F90,
irreducible_set_pert.F90, dfpt_vloca.F90, dfpt_vtorho.F90, dfpt_vtowfk.F90. In
these routines, the developper should pay a particular attention to the rfpert
array, defined in the routine respfn (in m_respfn_driver.F90), as well as to the ipert local
variable.

When smdelta in non-zero, it will trigger the calculation of the imaginary
part of the second-order electronic eigenvalues, which can be related to the
electronic lifetimes. The delta function is evaluated using:

when smdelta == 3, Cold smearing by Marzari using the parameter a=-0.8165 (monotonic function in the tail): as 2 but different a

when smdelta == 4, Smearing of Methfessel and Paxton ([Methfessel1989]) with Hermite polynomial of degree 2, corresponding to “Cold smearing” of N. Marzari with a=0 (so, same smeared delta function as smdelta=2, with different a).

The Matrix to be diagonalized in the Casida framework (see [Casida1995])
is a NxN matrix, where, by default, N is the product of the number of occupied
states by the number of unoccupied states. The input variable td_maxene
allows one to diminish N: it selects only the pairs of occupied and unoccupied
states for which the Kohn-Sham energy difference is less than td_maxene.
The default value 0.0 means that all pairs are taken into account.
See td_mexcit for an alternative way to decrease N.

The Matrix to be diagonalized in the Casida framework (see [Casida1995])
is a NxN matrix, where, by default, N is the product of the number of occupied
states by the number of unoccupied states. The input variable td_mexcit
allows one to diminish N: it selects the first td_mexcit pairs of occupied and
unoccupied states, ordered with respect to increasing Kohn-Sham energy difference.
However, when td_mexcit is zero, all pairs are allowed.
See td_maxene for an alternative way to decrease N.

If tim1rev is equal to 1, the Sternheimer equation is solved simultaneously at
+q and -q perturbation wavevectors. The first order potential at -q is taken
to be equal to the Hermitian conjugate of the first order potential at +q.
The wavefunctions from both +q and -q are then combined to generate the first order density.
Relevant in the case of magnetic field perturbation (but will be relevant also in case of non-zero frequency DFPT, when implemented).

Determine which non-linear implementation is used. If usepead=1, the Perturbation
Expansion After Discretization formalism is used, as in [Veithen2005].
In that method, the electric field is treated numerically, i.e the k-space
gradient operator appearing in the expression of the electric field potential
is discretized (see Eq.7 and 10 of [Veithen2005]).
If usepead=0, the electric field is treated analytically, leading to a better
k-points convergence. Furthermore, the current implementation is compatible with PAW
pseudopentials, while usepead=1 is not. The drawback of the analytical method
is one has to solve a second order Sternheimer equation before actually computing
third derivatives of the energy, using rf2_dkdk and rf2_dkde.
This is not the most time-consumming part though.
Look at the inputs of related tests in the testsuite to see examples of the workflow.