Help file for the Anaddb utility of the ABINIT package.

• 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 garantee 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.

(NOTE : in order to confuse even more the situation, it seems that the acoustic phonon frequencies generated by the code for both the sym and asym options are exactly the same ... likely due to an extra symmetrisation in the diagonalisation routine. Of course, when the matrix at Gamma has been generated from IFCs coming from dynamical matrices none of which are Gamma, the breaking of the ASR is rather severe. In order to clear the situation, one should use a diagonalisation routine for non-hermitian matrices. So, at the present status of understanding, one should always use the asr=2 option ).



Go to the top | List of ANADDB input variables

• 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.

Warning : the generation of q-points in anaddb is rather old-fashioned, and should be replaced by routines used by the main abinit code.



Go to the top | List of ANADDB input variables

More detailed explanation : the sum of the effective charges in the unit cell should be equal to zero. It is not the case in the DDB, and this sum rule is sometimes strongly violated. In particular, this will make the lowest frequencies at Gamma non-zero. There is no "best" way of imposing the ASR on effective charges.



Go to the top | List of ANADDB input variables

If the analysis is activated, one get the trace of the matrices between pairs of atoms, if dipdip is 1, get also the trace of the short-range and electrostatic part, and calculate the ratio with the full matrice; then define a local coordinate reference (using the next-neighbour coordinates), and express the interatomic force constant matrix between pairs of atoms in that local coordinate reference (the first vector is along the bond; the second vector is along the perpendicular force exerted on the generic atom by a longitudinal displacement of the neighbouring atom - in case it does not vanish; the third vector is perpendicular to the two other) also calculate ratios with respect to the longitudinal force constant ( the (1,1) element of the matrix in local coordinates).



Go to the top | List of ANADDB input variables

The wavevectors in list 2 will be used to :
- generate and diagonalize a dynamical matrix, and print the corresponding eigenvalues.
- calculate the generalized Lyddane-Sachs-Teller relation. Note that if the three first numbers are zero, then the code will do a calculation at Gamma without non-analyticities.



Go to the top | List of ANADDB input variables

This number defines the atoms for which the short range part of the interatomic force constants, after imposition of the acoustic sum rule, will not be put to zero. This option is available for testing purposes (evaluate the range of the interatomic force constants), because the acoustic sum rule will be violated if some atoms are no more included in the inverse Fourier Transform.



Go to the top | List of ANADDB input variables

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 configurationof 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.

As described in the Na Sai's paper, it is important to use the finite difference expression of the ddk (berryopt=2 or -2) in the linear response calculation of the effective charges and the piezoelectric tensor.



Go to the top | List of ANADDB input variables

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.

• 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.

• 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.

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

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).

• 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.

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

• 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.

• 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

• 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 (if iavfrq=1)
• 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).

Input variables that may be needed if this flag is activated : dostol, nchan, ntemper, temperinc, tempermin, as well as the wavevector grid number 2 definition, ng2qpt, ngrids, q2shft.



Go to the top | List of ANADDB input variables