ABINIT Discussion Forums

Hi all,

I have a question about calculating the phonon DOS, including the perturbation from the E field. I'm following RF2 (https://docs.abinit.org/tutorial/rf2/). As they suggest, to generate my phonon wavevectors, qpt, I first run a GS calculation. I'm using a low-symmetry monoclinic cell, so I end up with 518 symmetry-inequivalent wavevectors. I don't want to run this calculation 518 times (so expensive!), and I'm surprised there isn't a faster way to do this. I know the code is doing a DFPT SCF calculation for each qpt, but I'm thinking that one could also just calculate the dynamical matrix at gamma, the Born Effective Charge, the high-frequency and static dielectric tensors, and could then calculate the long-range corrections to the dynamical matrix. Then you could diagonalize this dynamical matrix for each qpt to get the frequencies and integrate get the DOS. Is there already an option to do this in Abinit or a tutorial I have not seen that describes this process? All I know is that I can't do this calculation 520 times for several materials! So if anyone spots an error in my logic or has a solution, I would be very glad to know!
Cheers

Hi
I wonder What's the k grid you use to generate phonon wavevectors?
I have been calculated some low symmetry system but haven't seen such many wavevectors yet.

Andy

Hi Andy,

I run a GS calculation, and I set :
ngkpt 12 12 12
kptopt 1
nshiftk 1
shiftk 0 0 0
prtkpt 1
nstep 1
nline 1
prtvol 2 #more than first 50 of kpts


which makes the outfile contain all the kpoints used. This way, I generate all the symmetry inequivalent qpts needed to do a realistic integration over the Brillouin Zone.

Hi

In my experiment, I think the 12x12x12 is too big for q-grids since the 12x12x12 q-points equal to generate 12x12x12 supercell's force constant. It's not surprising that the calculation will be expensive.
I suggest you use the 2x2x2 at begin.
If you want to check the convergence of the result, than add the k-grid to check the convergence of the properties you interest in.

Best wishes,
Andy

Hi Andy,

Thanks for the response. I tried lowering the number of qpts, but I don't think this run is using a supercell technique, right? Abinit uses DFPT to get around the supercell, frozen-phonon method for phonons. So while I agree, 12x12x12 is too much for this calculation, isn't 2x2x2 very few points in the BZ, considering it's just the primitive cell in this calculation?

Cheers

Hi

Yes the DFPT method is not using a supercell technique directly, but the DFPT phonon calculation is to construct a set of Finite-wave-vector phonon calculations on the q-grid and use the discrete Fourier transform to gain the IFC in the real space to get the Phonon dispersion. So it doesn't matter if you are using the supercell or not. That's a advantage for the DFPT method i think.
You can find more information in Phys. Rev. B 55, 10355-10368 (1997)

The convergence for the phonon dispersion need to test depend on the system. Generally, the 2x2x2 q-grid is usually good enough for the discrete Fourier transform. if you want to get the phonon DOS. you can try to do the converge test by yourself.(I think for your system, the 2x2x2 q-grid still cost some time)

Cheers
Andy

Hi Andy,

Okay, I read that paper. So Abinit uses a Fourier transform to take values of the dynamical matrix on a grid around the BZ (2x2x2 for example) in order to calculate the force constants and other derivatives needed to calculated phonon frequencies. After this main Abinit run when you build up the ddb, you can sample many points in the BZ with your anaddb input file because all of the force constants and such are known, so you can calculate the modes anywhere in the BZ. So your qpts in the main run don't have to match the qph1l list of qpts in the anaddb run?

Is that correct?

Many thanks!
J

Hi

Yes the qph1l list of qpts doesn't need to match the the main run, but the anaddb input ngqpt need to be the q-grid you use.
Once you have the information at the q-grid, the rest of the phonon dispersion can be form by the Fourier transform. so you can choose any path in the Brillouin zone to form the dispersion you interested in.

Cheers
Andy

Okay, final question - so does that mean that nkpt in the main run has to have the same number of points as the qpts?

In the RF2 tutorial, ngkpt = 4 4 4, and it looks like the qpts are determined by getting the symmetry inequivalent points from a 4x4x4 mesh. I can't tell if this is a coincidence, but the tutorial says this in the first input file trf2_1.in:

#Complete set of symmetry-inequivalent qpt chosen to be commensurate
# with kpt mesh so that only one set of GS wave functions is needed.

To be commensurate, could you also do nkpt 12 1 2 12 (what I usually do for a run to converge) in the main run and use the qpts generated from a (2x2x2)run, because 12/2 = integer? Or must they be exactly same? As in ngkpt (l m n) needs to match the (n l m) used to generate the qpts, and also match ngqpt in the anaddb file . But ng2qpt can be whatever you want it to be.

Hi

From my experience, the ngkpt in the main run is okay if it is commensurate to the q-grids.
if you want to use the fine k-grid with only 2x2x2 q-grid, I suggest to use integral multiple ngkpt of the q-grids or it might cause the numerically error.(Image that the q-point is not on the k-grid, so the tiny error will be accumulation) The details should be tested case by case.

The only thing I can say that is your idea is correct from my perspective, but it's not a guarantee. Perhaps you could wait for the expert to answer about this one.
BTW, you should also check the converge with respect to the k-grid.

Cheers
Andy

Hi Andy,

Thanks so much for your help! I finally got a beautiful DOS. One question though, I noticed that there's no option to include the direction of the non-analytic wavevectors generated by the wavevector grid(ng2qpt) near the zone center, which would be affected by the long-range force constants. Usually, Abinit needs me to enter qph2l and nph21 to give the direction of the NA wavevector, and that's how I know it's being taken into account.

In my DOS, I have more lower frequency states than I expected. I wonder if it's because all the points in the near the zone center that should get corrections aren't getting a boost.

Here's my anaddb input file: I included nph2l and qph2l, but I don't think they're being used for the DOS calculation. Rather, they're for the phonon spectrum.

Does anyone know if Abinit can include NA corrections in the DOS calculation?

*******EDIT***********
I've considered the fact that the LST corrections really only matter when q < omega_transverse/speed of light. My w_T at gamma = 163 cm-1, which means that the q-vectors generated from ng2qpt would need to have components less than 163 cm-1 for the boost to occur. My BZ is almost cubic with lattice constant 4 Ang, so the zone edges are +/- pi/4 Ang ~ 7.85E7 cm-1. If I use a 30x30x30 grid, the smallest wavevectors (besides gamma) will be so much larger than 163 cm-1 that it doesn't matter. So only one qpt is missing the LST boost, gamma, and I don't think this is causing the larger than expected DOS for low frequencies. If you spot an error in my logic, please let me know.
***********************

Cheers
J

Hi...just some comment and different view

I cannot tell if your idea is correct since i am not familiar with the PDOS calculation.
But i think you can do the simple test by using the small grid to prove your point.

Actually I think that the much lower frequency states is not causing by the missing NA correction, but i am not so sure
I also wondering the answer, hoping that someone can answer the question.

Cheers,
Andy

Dear Andy and Jerkov,
The phonon DOS integrates the TO modes only at Gamma, but for large grid of q-point this should converge since the weight of Gamma point will be reduced by all the other q-points that contain the LO branch.
Best wishes,
Eric

Great! All good then. You can mark this as solved.

Hi Eric and Andy and whoever else is reading,

I'm working with a crystal that is stable at 0 K and with very small forces order ~1e-8 Hartree/Bohr. Yet I have a few imaginary frequencies in my BZ. 11 out of the 846 values when I choose to have nwchan = 5 and nchan = 1250 are imaginary.

I was going through my output files, and in my log file I found the suggestion to increase the number of bands in the GS calculation. I noticed that the tutorial (https://docs.abinit.org/tutorial/rf2/) works with only the valence bands being possibly occupied (valence of Al:5, As: 3, so (3+5)/2 = 4) for all datasets, including the GS. I found this unusual. Is this recommended generally or just for expediency in this example?
Cheers,
J

Dear Jerkov,
For insulators, the number of bands necessary can be set to the occupied bands only (and with an insulator occopt value).
For metals, you have to be sure the number of bands is large enough, but this is a convergence test to do. Is your system metallic? If so, this can be one of the different sources of the problem.
Note that you can set a metallic occupation (and with nband>occupied-bands) for an insulator too if the tsmear is not set too large w.r.t. the band gap value but the number of unoccupied bands should not affect the result.
Did you impose the acoustic sum rule in the post-processing analysis of the phonons with anaddb?
Best wishes,
Eric

Hi Eric,

I believe I imposed ASR. The Gamma point itself is stable, but not wavevectors near it. Here's my anaddb input file:

!Input file for the anaddb code. Analysis of the SiO2 DDB

!Flags
dieflag 1 ! Frequency-dependent Dielectric tensor flag
ifcflag 1 ! Interatomic force constant flag
thmflag 1 ! Thermal flag. Gives Internal energy, entropy,
! heat capacity, phonon DOS, Debye-Waller factor)
prtdos 2
dosdeltae 4E-06

!Effective charges
chneut 2 ! Charge neutrality requirement for effective charges.
! 2=> imposed with weights proportional to screening)

!Interatomic force constant info
dipdip 1 ! Dipole-dipole interaction treatment
ifcana 1 ! Analysis of the IFCs
natifc 5 ! Number of atoms in the cell for which ifc's are analysed
atifc 1 2 3 4 5 ! List of atoms
asr 1

!Thermal information
nchan 1250 ! # of channels for the DOS with channel width 1 cm-1
nwchan 5 ! # of different channel widths from this integer down to 1 cm-1
thmtol 0.060 ! Tolerance on thermodynamical function fluctuations
!dosdeltae set to default of 1 cm-1

!Wavevector grid number 1 (coarse grid, from DDB), should match grid used to generate qpts
brav 1
ngqpt 2 2 2 ! Monkhorst-Pack indices
nqshft 1 ! number of q-points in repeated basic q-cell
q1shft 0 0 0

!Wavevector grid number 2 (series of fine grids, extrapolated from intrat forces)
ng2qpt 100 100 100 ! sample the BZ up to ngqpt2
ngrids 3 ! number of grids of increasing size
q2shft 0.0 0.0 0.0

nph2l 1 ! number of phonons in list 2
qph2l 1.0 1.0 1.0 0.0
symdynmat 0

Cheers,
J

Dear Jerkov,
OK, you have included the asr, it is not the problem then.
If only the points in between the main q-points at which you did the DFPT are giving strange negative results it might means that the q-grid is not large enough to have a good interpolation. You used 2x2x2 q-grid, which night not be large enough to get the right frequencies at some q-points, like the ones close to Gamma.
Best wishes,
Eric

I believe I imposed ASR. The Gamma point itself is stable, but not wavevectors near it. Here's my anaddb input file:

This may be an artefact of the Fourier interpolation if the DFPT q-mesh is too coarse and/or you are dealing with tricky systems.

A possibile solution (not the most elegant, I would say) is to cutoff the IFCs(R) in R-space.
See https://docs.abinit.org/variables/anaddb/#nsphere , in particular nsphere = -1.
For a more elegant treatment see https://arxiv.org/pdf/2004.08875.pdf (not yet available in the official version).

BTW: I saw that you're using the "old algorithm" to compute the phonon DOS with the bin method.
There's a more advanced and accurate method based on the tetrahedron method (see https://docs.abinit.org/variables/anaddb/#prtdos).

Hope it helps.

Hi everyone,

The issue is resolved. I used more qpts (effectively a 3x3x3 supercell), and the instability disappeared. Thanks for all the help!

J