It will be easier to discover the
present file with the help of the tutorial
(see lesson 5).
It is worthwhile to print this help file, for ease of reading.
ABINIT can give direct access to various response functions, that are second derivatives of total energy (2DTE) with respect to different perturbations. Presently, these perturbations can be of three types : (1) phonons, (2) static homogeneous electric field, and (3) strain. The physical properties connected to 2DTE with respect to perturbations (1) and (2) are the phonon dynamical matrices, the dielectric tensor, and the Born effective charges, while the additional strain perturbation (3), mixed with phonon and electric field leads to elastic constant, internal strain, and piezoelectricity.
More functionalities of the computation of responses should be implemented sooner or later. The alchemical perturbation is one of the candidate, as well as the homogeneous magnetic field perturbation. Some third derivatives of the total energy (3DTE) are also implemented. The 3DTE might give phonon-phonon coupling, non-linear electric response, anharmonic elastic constants, gruneisen parameters,...
The basic quantities that ABINIT will compute are the FIRST-order derivatives of the wavefunctions (1WF) with respect to these perturbations. The later calculation of the 2DTE and 3DTE from these 1WF is an easy computational task : the construction of the 2DTE with respect to perturbations j1 and j2 involves mainly evaluating matrix elements between the 1WF of j1 and/or the 1WF of j2. More detail on this technique can be found in X. Gonze, Phys. Rev. B55, 10337 (1997) and X. Gonze and C. Lee, Phys. Rev. B55, 10355 (1997).
The calculation of the 1WF for a particular perturbation is done using a variational principle and an algorithm rather similar to the one used to find the ground-state (unperturbed) wavefunctions. Thus, lot of technical details and parameters are the same for both ground-state and response-function calculations. This justifies the development of one unique code for these two classes of properties : many of the routines of abinit are common in these calculations, or had to be duplicated, but with relatively small modifications.
The ABINIT code performs a rather primitive analysis of the calculated 2DTEs. For example, it gives the phonon frequencies, electronic dielectric tensor and effective charges. But the main output of the code is the Derivative DataBase (DDB) : a file that contains the set of all 2DTEs and 3DTEs calculated by the code. This DDB can be manipulated by the MRGDDB code, and fully analyzed by the Anaddb code. See the corresponding help files (Mrgddb help file, and Anaddb help file).
The perturbation of the phonon type is the displacement of one atom along one of the axis of the unit cell, by a unit of length (in reduced coordinates). It is characterized by two integer numbers and one wavevector. The two integer numbers are the number of the moved atom, which will be noted "ipert", and the number of the axis of the unit cell, which will be noted "idir". "ipert" for phonon perturbation can have values between 1 and natom, "idir" can have values between 1 and 3. The set of all possible phonon perturbations for one wavevector has thus (3*natom) elements. From this basis set, all phonons can be constructed by linear combination. The set of atoms to be moved in one dataset of ABINIT is determined by the input variable rfatpol. The set of directions to be considered in one dataset of ABINIT is determined by the input variable rfdir. The wavevector to be considered in one dataset of ABINIT is determined by the input variables nqpt, qpt, and qptnrm.
If you follow the tutorial, you should go back to the tutorial window now.
The perturbations of the electric field type are
Note : the ddk perturbation is defined as a the derivative with respect to k in reduced coordinates; this is equivalent to applying a linear perturbation of strength 2*pi along the conjugate direction in real space. This statement comes from the derivation of the phase factor exp(i*2*pi*k*r) with respect to k in reduced coordinates.
These electric-type perturbations are also characterized by two numbers : "ipert" being natom+1 for the ddk perturbation and natom+2 for the electric field, and "idir" being 1, 2 or 3, as for phonon perturbations. Although the possibility of electric field characterized by a non-zero wavevector is envisioned for a future version of the code, at present only homogeneous fields are considered. So the wavevector of the electric field type perturbations is Gamma (q=0).
The perturbations of the strain type are either an uniaxial strain
or a shear strain. The strain perturbations are considered
in cartesian coordinates (x,y,z).
They are characterized by two numbers,
with "ipert" being natom+3 for the uniaxial strains, and natom+4
for the shear strains, and "idir" describe the particular component.
Explicitely, for uniaxial strains, idir=1 gives the xx strain perturbation,
idir=2 gives the yy strain perturbation, idir=3 gives the zz strain perturbation,
while for shear strains, idir=1 gives the yz strain perturbation,
idir=2 gives the xz perturbation, and idir=3 gives the xy perturbation.
Note that the "rigid-atom" elastic constants,
as output of ABINIT, are those obtained with frozen internal
coordinates. The internal coordinate relaxation, needed to give
"physical" elastic constants should be handled through the knowledge
of the internal strain and dynamical matrix at Gamma, in ANADDB.
WARNING : as of v4.2, the strain perturbation is a very
recent addition to ABINIT (thanks to D.R. Hamann). Some
caution is in order ...
In v4.2, no automatic computation
of the "physical" elastic constant is provided, although
all the ingredients for such a calculation exist. Of course,
if all the internal coordinate are fixed by symmetry, all the
internal strains vanish, and the "rigid-atom" and "physical" elastic
constants are equal.
Limitations of the present implementation :
We also define the index of the perturbation, called "pertcase", equal to idir + 3*ipert. Accordingly, pertcase runs from 1 to 3*(natom+4), and will be needed to identify output and input files, see section 6 .
To summarize, the perturbations are characterized by two numbers, "ipert" from 1 to natom+4, and "idir", from 1 to 3, as well as one wavevector (that is gamma when a non-phonon perturbation is considered); a number called "pertcase" combines "ipert" and "idir", and runs from 1 to 3*(natom+4).
The 2DTE, being derivative of the total energy with respect to two perturbations, will be characterized by two sets of (idir,ipert), or by two pertcase numbers, while 3DTE will need three such sets or pertcase numbers. In addition they will depend on one wavevector (for 2DTE) or two wavevectors (for 3DTE).
In the present (non-stationary) implementation of the 2DTE, the first pertcase corresponds to the perturbation that gives the derivative of the potential, and the second pertcase corresponds to the perturbation that gives the derivative of the wavefunctions.
The same 'files' file as for GS calculations is used for RF calculations. Actually, in the multi-dataset mode, one will be able to make in one ABINIT run, ground-state computations as well as response-function computations, so that the 'files' file must be the same.... The 'input' file will have many common input variables for these different cases, but also some separate ones.
Two ground-state wavefunction files might be needed :
Note that the second file (k+q) will be identical to the first one (k), in the case of zero-wavevector perturbation. Also, if the q-wavevector maps the k-point grid onto itself (q being the difference between points that belong to the grid), all the information on the k+q grid is contained in the k grid, a fact that ABINIT is able to exploit.
One should have a look at the input variables irdwfk, and irdwfq, for a description of the ground-state wavefunction file names generated from the root names provided in the 'files' file. In the multi-dataset mode, the following input variables will be relevant : getwfk, and getwfq. The file names of the ground-state wavefunction file follows the same convention as for the ground-state case. Thus, the corresponding section of the abinis_help file can be read, if needed.
In the case of an electric field perturbation, the output 1WF of the corresponding ddk perturbation is needed as input. If the option rfelfd=1 is used, then the code will take care of doing first the derivative dk perturbation calculation, then write the 1WF at the correct place, as an output file, then begin the homogeneous field perturbation calculation. Usually, the use of rfelfd=1 is not recommended, as the ddk computation is the most often done with different parameters as the electric field perturbation.
The nomenclature of first-order wavefunction files is also given in the abinis_help file, but it is worth to specify it in more detail here. The root name is formed from the string of character in the third line of the 'files' file (for an input file) or the fourth line of the 'files' file (for an output file), that is complemented, in the multi-dataset mode, by '_DS' and the dataset pertcase (the index of the perturbation). Then, the string '_1WF' is added, followed by pertcase. This gives, e.g., for a 'files' file whose fourth line is '/tmp/o', for the dataset number 3, and a perturbation corresponding to the displacement of the second atom in the x direction (pertcase=4), the following name of the corresponding 1st-order wavefunction output file:
/tmp/o_DS3_1WF4Such a file might be used as input of another computation, or of another dataset.
The relevant input variables are : ird1wf, and irdddk, as well as get1wf, and getddk, in the multi-dataset mode.
In order to understand correctly the behaviour of response-function runs, some information on the use of symmetries must be given.
Some perturbations (including their wavevector) may be invariant for some symmetries. ABINIT is able to use symmetries to skip perturbations of which a symmetric has already been calculated. ABINIT is also able to use the symmetries that keeps the perturbations invariant to reduce the number of k points needed for the sampling of electronic wavefunctions in the Brillouin zone (although this feature is not optimal yet). There is one exception to this, the ddk perturbation, for which the spatial symmetries cannot be used yet.
In any case, unlike for the ground-state, the input k-point set for response function should NOT have been decreased by using spatial symmetries, prior to the loop over perturbations (see section 4). Only the time-reversal symmetry, retained by calculations at q=0, ought to be used to decrease this input k-point set. Considering each perturbation in turn, ABINIT will be able to select from this non-reduced set of k-points, the proper k-point set, automatically, by using the symmetries that leave each perturbation invariant.
Accordingly, the preferred way to generate the k-point grid is of course to use the ngkpt or kptrlatt input variables, with different values of kptopt :
In agreement with the information provided in the previous sections, different cases can be distinguished.
When one considers the response to an atomic displacement with q=0, the following procedure is suggested :
When one considers the response to an electric field (with q=0), the following procedure is suggested :
When one considers the response to an atomic displacement in the special case where q connects k-points that both belong to the special k-point grid, the following procedure is suggested :
When one considers the response to an atomic displacement for a general q point, the following procedure is suggested :
Of course, these different steps can be combined when a set of responses is looked for. In particular, the computations of responses at gamma, in the case where the full dynamical matrix as well as the dielectric tensor and the Born effective charges are needed, can be combined as follows :
Note that the error in the 2DTE is linear in the ground-state wavefunction error (unlike the error due to the 1WFs). Moreover, a large prefactor is associated with this source of error (it can even cause cause the unstability of the SCF procedure). As a consequence, the convergence of the ground-state wavefunction should be very good. The same is true at the level of the ddk wavefunctions.
As mentioned in the introduction, inside the response-function part of the code, the calculation proceeds in two steps : first the calculation of the first order derivative of the wavefunctions (1WF), then the combinations of these 1WF to build the 2DTE and 3DTE.
In an initialisation part, the input file and all the ground-state files are read, and a few basic quantities are constructed.
In the first part, every perturbation is examined, one at a time, separately :
A subset of the ABINIT input variables have a modified meaning or a modified behaviour in case of RF calculations. Here is the list of these input variables, together with the variables that applies only to RF computations. Note that the code will do a RF calculation (optdriver=1) when one of rfphon or rfelfd is non-zero.
Output from the code goes to several places listed below.
This file is the same as the log file of the abinit code when computing ground state (GS) results. Eventually, the output of datasets related to response functions will be intertwinned with those concerned with ground-state case. The purpose of this file is the same as in the GS case, and the use of error messages is unchanged.
This file is the same as the main output file of the abinit code when computing ground state (GS) results. Eventually, the output of datasets related to response functions will be intertwinned with those concerned with ground-state case. We explain here the parts related to the RF computation.
The initialisation part is the same as for the GS. So, the reader is adviced to read the section 6.2 the abinis_help file, as well as the first paragraph of the section 6.3 of this file. Afterwards, the content of the main output file differs a bit...
The main output file reports on the initialisation of the ground-state wavefunctions at k+q, then the loop on the perturbations begins. For each perturbation, there is :
After the computation of each perturbation, the output file reports on the 2DTE matrix elements. This part is not executed if the only perturbation is the derivative dk perturbation. It will give the following information :
For this last analysis, the code assumes that the whole set of perturbations in a class has been calculated, either all the phonon perturbations or the homogeneous electric field perturbation, or both. If this is not true, the code will give results that may be wrong in the case that the reduced system of coordinates is not cartesian (for the dynamical matrix, the effective charge tensor of the dielectric matrix), or in all case wrong (the phonon frequencies); also the code will put zero in the matrix elements that have not been calculated. A Warning message is issued if the above information cannot be trusted.
Finally, the code provides the timing information.
6.3. The first-order wavefunction (1WF) files
These are unformatted data files containing the planewaves coefficients of all the wavefunctions, written in the following format. First, the header (see section 6.4 of the abinis_help), followed by
bantot=0 <-- counts over all bands
do isppol=1,nsppol
do ikpt=1,nkpt
write(unit) npw1,nband1 <-- for each k point
do iband=1,nband1
write(unit) (eigen1(iband+bantot),iband=1,2*nband) <-- the first-order eigenvalue changes
bantot=bantot+2*nband
write(unit) (cg(ii+...),ii=1,2*npw1) <-- single band and k
enddo
enddo
enddo
6.4. The first-order density files
They consist of the header lines, followed by
write(unit) (rhor1(ir),ir=1,cplex*ngfft(1)*ngfft(2)*ngfft(3))
Here, rhor1 is the change of electron density in electrons/bohr^3. The parameter cplex is 1 when q=0 and 2 when q/=0 . Indeed, for q=0, the density change is a real quantity, while it is complex in general when q/=0 .
6.5. The derivative database (DDB)
It is made of two parts. The first should allow one to unambiguously identify the run that has generated the DDB, while the second part contains the 2DTE, grouped by blocks of data.
Note that the DDB output of the ABINIT code can be merged with other DDBs as described in the Mrgddb help file.
The first part contains :
The second part contains :
The elements of a 2DTE block are described by 4 integers and 2 real numbers. The 2 first integers define a first perturbation in the form (idir1,ipert1), the two next define a second perturbation in the form (idir2,ipert2). The matrix element corresponds to the derivative of the total energy with respect to the parameters corresponding to these perturbations. The real numbers are the real and imaginary parts of the 2DTE. Sometimes, the code uses spatial symmetries, the time-reversal symmetry, or even the permutation of first and second perturbations to deduce the value of non-computed matrix elements. This behaviour might be improved, as it is somtimes confusing ...
If you follow the tutorial, you should go back to the tutorial window now.
It is possible to get from the RF calculations essentially EXACT derivatives of the total energy with respect to perturbations. There is a published account of this fact in Phys. Rev. A 52, 1096 (1995). An agreement of 8 digits or more was obtained in the comparison with finite-difference derivatives of GS data.
The accuracy of these calculations are thus entirely determined by the input parameters the user choose for the RF run, and the preparatory GS runs.
We will now review the convergence parameters, usually the same as for the GS calculations, indicating only the specific features related to RF calculations.
Input parameters that could influence the accuracy of the calculation are :
The input parameters ecut, ecutsm, ixc, intxc, the set of k-wavevectors, as well as the related variables have to be the SAME in the ground-state calculations that go before a RF run and this RF run.
Namely : do not try to use ground-state wavefunction files produced with ecut=10Ha for a RF run with ecut=20Ha. You cannot even use newsp for translating the WF from ecut=10Ha to ecut=20Ha and use this latter file as an input of the RF run. In some cases, the code will complain and stop, but in other cases, it will simply produce garbage !
If the value of ngfft is input by hand, its value must also be equal in the GS and RF cases. ALWAYS use ngfft large enough to have boxcut=2 or larger, in order to avoid any FFT filter error. In the GS case, boxcut as small as 1.5 could be allowed in some cases. It is not allowed with RF calculations, because these are more sensitive to that error.
The convergence tests with respect to ecut, and the k-point grid should be done carefully. This was already emphasized in the abinis_help file and is re-emphasized here. The user should test the convergence DIRECTLY on the property he or she is interested in. For example, if the user wants a phonon frequency accurate to 10 cm^-1, he/she could be lead to do a full calculation (GS+RF) of phonons at 30Ha, then another full calculation at 35Ha, then another at 40Ha... It is an error to rely on tolerance on the total energy (for example 1mHa/atom) or geometry (accuracy of one part per thousand on the bond lengths) to draw 'a priori' conclusions on the convergence of other quantities, and not monitor the convergence of these directly. To be clear : if phonon frequencies are needed, check the convergence of phonon frequencies !
The user should note that for bands with very small occupancy in the metallic case as well as unoccupied bands for insulators, the ground state run preceeding response function runs will not necessarily converge these wavefunctions using usual ground-state tests such as toldfe or (better) tolvrs. To be sure that inaccuracies are not introduced into the response function calculation by poorly converged unoccupied bands, a separate run starting from a saved charge density (prtden=1 in the self-consistent run) and using iscf=-2 and tolwfr may be needed.