Phonon DOS [SOLVED]
Phonon DOS
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 lowsymmetry monoclinic cell, so I end up with 518 symmetryinequivalent 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 highfrequency and static dielectric tensors, and could then calculate the longrange 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
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 lowsymmetry monoclinic cell, so I end up with 518 symmetryinequivalent 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 highfrequency and static dielectric tensors, and could then calculate the longrange 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

 Posts: 45
 Joined: Tue Nov 26, 2019 12:50 am
Re: Phonon DOS
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
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
Re: Phonon DOS
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.
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.

 Posts: 45
 Joined: Tue Nov 26, 2019 12:50 am
Re: Phonon DOS
Hi
In my experiment, I think the 12x12x12 is too big for qgrids since the 12x12x12 qpoints 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 kgrid to check the convergence of the properties you interest in.
Best wishes,
Andy
In my experiment, I think the 12x12x12 is too big for qgrids since the 12x12x12 qpoints 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 kgrid to check the convergence of the properties you interest in.
Best wishes,
Andy
Re: Phonon DOS
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, frozenphonon 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
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, frozenphonon 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

 Posts: 45
 Joined: Tue Nov 26, 2019 12:50 am
Re: Phonon DOS
Hi
Yes the DFPT method is not using a supercell technique directly, but the DFPT phonon calculation is to construct a set of Finitewavevector phonon calculations on the qgrid 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, 1035510368 (1997)
The convergence for the phonon dispersion need to test depend on the system. Generally, the 2x2x2 qgrid 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 qgrid still cost some time)
Cheers
Andy
Yes the DFPT method is not using a supercell technique directly, but the DFPT phonon calculation is to construct a set of Finitewavevector phonon calculations on the qgrid 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, 1035510368 (1997)
The convergence for the phonon dispersion need to test depend on the system. Generally, the 2x2x2 qgrid 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 qgrid still cost some time)
Cheers
Andy
Re: Phonon DOS
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
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

 Posts: 45
 Joined: Tue Nov 26, 2019 12:50 am
Re: Phonon DOS
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 qgrid you use.
Once you have the information at the qgrid, 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
Yes the qph1l list of qpts doesn't need to match the the main run, but the anaddb input ngqpt need to be the qgrid you use.
Once you have the information at the qgrid, 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
Re: Phonon DOS
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 symmetryinequivalent 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.
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 symmetryinequivalent 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.

 Posts: 45
 Joined: Tue Nov 26, 2019 12:50 am
Re: Phonon DOS
Hi
From my experience, the ngkpt in the main run is okay if it is commensurate to the qgrids.
if you want to use the fine kgrid with only 2x2x2 qgrid, I suggest to use integral multiple ngkpt of the qgrids or it might cause the numerically error.(Image that the qpoint is not on the kgrid, 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 kgrid.
Cheers
Andy
From my experience, the ngkpt in the main run is okay if it is commensurate to the qgrids.
if you want to use the fine kgrid with only 2x2x2 qgrid, I suggest to use integral multiple ngkpt of the qgrids or it might cause the numerically error.(Image that the qpoint is not on the kgrid, 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 kgrid.
Cheers
Andy
Re: Phonon DOS
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 nonanalytic wavevectors generated by the wavevector grid(ng2qpt) near the zone center, which would be affected by the longrange 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 cm1, which means that the qvectors generated from ng2qpt would need to have components less than 163 cm1 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 cm1. If I use a 30x30x30 grid, the smallest wavevectors (besides gamma) will be so much larger than 163 cm1 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
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 nonanalytic wavevectors generated by the wavevector grid(ng2qpt) near the zone center, which would be affected by the longrange 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:
Code: Select all
!Flags
dieflag 1 ! Frequencydependent Dielectric tensor flag
ifcflag 1 ! Interatomic force constant flag
thmflag 1 ! Thermal flag. Gives Internal energy, entropy,
! heat capacity, phonon DOS, DebyeWaller factor)
prtdos 1
dosdeltae 5E05
!Effective charges
chneut 2 ! Charge neutrality requirement for effective charges.
! 2=> imposed with weights proportional to screening)
!Interatomic force constant info
dipdip 1 ! Dipoledipole 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
!Thermal information
nchan 1250 ! # of channels for the DOS with channel width 1 cm1
nwchan 5 ! # of different channel widths from this integer down to 1 cm1
thmtol 0.060 ! Tolerance on thermodynamical function fluctuations
!dosdeltae set to default of 1 cm1
!Wavevector grid number 1 (coarse grid, from DDB), should match grid used to generate qpts
brav 1
ngqpt 2 2 2 ! MonkhorstPack indices
nqshft 1 ! number of qpoints in repeated basic qcell
q1shft 0 0 0
!Wavevector grid number 2 (series of fine grids, extrapolated from intrat forces)
ng2qpt 30 30 30 ! 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
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 cm1, which means that the qvectors generated from ng2qpt would need to have components less than 163 cm1 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 cm1. If I use a 30x30x30 grid, the smallest wavevectors (besides gamma) will be so much larger than 163 cm1 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

 Posts: 45
 Joined: Tue Nov 26, 2019 12:50 am
Re: Phonon DOS
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
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
Re: Phonon DOS
Dear Andy and Jerkov,
The phonon DOS integrates the TO modes only at Gamma, but for large grid of qpoint this should converge since the weight of Gamma point will be reduced by all the other qpoints that contain the LO branch.
Best wishes,
Eric
The phonon DOS integrates the TO modes only at Gamma, but for large grid of qpoint this should converge since the weight of Gamma point will be reduced by all the other qpoints that contain the LO branch.
Best wishes,
Eric
Re: Phonon DOS
Great! All good then. You can mark this as solved.
Re: Phonon DOS
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 ~1e8 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
I'm working with a crystal that is stable at 0 K and with very small forces order ~1e8 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
Re: Phonon DOS
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>occupiedbands) 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 postprocessing analysis of the phonons with anaddb?
Best wishes,
Eric
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>occupiedbands) 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 postprocessing analysis of the phonons with anaddb?
Best wishes,
Eric
Re: Phonon DOS
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 ! Frequencydependent Dielectric tensor flag
ifcflag 1 ! Interatomic force constant flag
thmflag 1 ! Thermal flag. Gives Internal energy, entropy,
! heat capacity, phonon DOS, DebyeWaller factor)
prtdos 2
dosdeltae 4E06
!Effective charges
chneut 2 ! Charge neutrality requirement for effective charges.
! 2=> imposed with weights proportional to screening)
!Interatomic force constant info
dipdip 1 ! Dipoledipole 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 cm1
nwchan 5 ! # of different channel widths from this integer down to 1 cm1
thmtol 0.060 ! Tolerance on thermodynamical function fluctuations
!dosdeltae set to default of 1 cm1
!Wavevector grid number 1 (coarse grid, from DDB), should match grid used to generate qpts
brav 1
ngqpt 2 2 2 ! MonkhorstPack indices
nqshft 1 ! number of qpoints in repeated basic qcell
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
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 ! Frequencydependent Dielectric tensor flag
ifcflag 1 ! Interatomic force constant flag
thmflag 1 ! Thermal flag. Gives Internal energy, entropy,
! heat capacity, phonon DOS, DebyeWaller factor)
prtdos 2
dosdeltae 4E06
!Effective charges
chneut 2 ! Charge neutrality requirement for effective charges.
! 2=> imposed with weights proportional to screening)
!Interatomic force constant info
dipdip 1 ! Dipoledipole 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 cm1
nwchan 5 ! # of different channel widths from this integer down to 1 cm1
thmtol 0.060 ! Tolerance on thermodynamical function fluctuations
!dosdeltae set to default of 1 cm1
!Wavevector grid number 1 (coarse grid, from DDB), should match grid used to generate qpts
brav 1
ngqpt 2 2 2 ! MonkhorstPack indices
nqshft 1 ! number of qpoints in repeated basic qcell
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
Re: Phonon DOS
Dear Jerkov,
OK, you have included the asr, it is not the problem then.
If only the points in between the main qpoints at which you did the DFPT are giving strange negative results it might means that the qgrid is not large enough to have a good interpolation. You used 2x2x2 qgrid, which night not be large enough to get the right frequencies at some qpoints, like the ones close to Gamma.
Best wishes,
Eric
OK, you have included the asr, it is not the problem then.
If only the points in between the main qpoints at which you did the DFPT are giving strange negative results it might means that the qgrid is not large enough to have a good interpolation. You used 2x2x2 qgrid, which night not be large enough to get the right frequencies at some qpoints, like the ones close to Gamma.
Best wishes,
Eric
Re: Phonon DOS [SOLVED]
This may be an artefact of the Fourier interpolation if the DFPT qmesh is too coarse and/or you are dealing with tricky systems.I believe I imposed ASR. The Gamma point itself is stable, but not wavevectors near it. Here's my anaddb input file:
A possibile solution (not the most elegant, I would say) is to cutoff the IFCs(R) in Rspace.
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.
Re: Phonon DOS
Hi everyone,
The issue is resolved. I used more qpts (effectively a 3x3x3 supercell), and the instability disappeared. Thanks for all the help!
J
The issue is resolved. I used more qpts (effectively a 3x3x3 supercell), and the instability disappeared. Thanks for all the help!
J