ABINIT Discussion Forums

Hi all,

I've been using ABINIT v5.8.4 to calculate the properties of a high pressure phase of CO2, space group P42/mnm.

I first optimized the strucutre (cell paramaters + atomic positions) for a given stress target
then I ran a response function calculation to get the phonon dispersion in the IBZ, following the indications given in the tutorials (file CO2ph2_13.in).



Everything looked to be working out right so far. I then merged the DDB and ran anaddb to get the phonon bands (file CO2ph2_15.in).



The results showed imaginary modes on the lower acoustic branch along Gamma-M(0.5,0.5,0). To make sure this was not due to convergence issues, I increased ecut from 35 Ha to 40 Ha but got the same result. I upload the figure, this is an eps file but added the ".in" extension to be able to upload it

Since I found the presence of these imaginary modes rather weird, I ran more rf calculations at 2 wavevectors along Gamma-M, q=[0.1,0.1,0] and q=[0.15,0.15,0]. The reuslts do not agree at all with those given by anaddb, in particular all frequencies are positive in these direct rf calculations. They are shown with "+" symobls on the figure.

I also made a rf calculation along Gamma-X(0,0.5,0),q=[0,0.15,0], and in this case the results are in perfect agreement with the output of anaddb.

Can someone help me understand what the problem is ?

Thanks in advance,
Frederic.

Hi,

Did you try removing the variable symdynmat from the input file for ANADDB?
The default (symdynmat 1) is better has it restores the symmetry of the system, even with Fourier interpolation.
I am not sure this will correct the behavior, I would need to see the outputs.

Hi,
Thanks for your suggestion. I m out of office at present, I will try it out as soon as I come back.
Frederic

Your approach is quite correct: check that the negative modes are "real" by doing explicit calculations at the given q-points. As you showed the modes should be positive. this means that the negative values are a result of the interpolation scheme. The straightforward solution is to increase the density of q-points. You will get better converged real-space interatomic force constants, and the negative dips will disappear.

You can also use an approximation, with the rifcsph anaddb input variable. It cuts off the IFC at the given radius, and will eventually eliminate the large, long-range IFC which are causing the dips, but also other oscillations and dips which may be physical. Start with a large value like 60 bohr, then reduce rifcsph progressively to see the effect. Again, to be used with caution, as it is a severe approximation (see the appendix of prb 78 045119 ).

matthieu

Thanls Matthieu.

I tried with symdinmat 1, it did not change anything.

I then tried to use decreasing values of rifcsph, it does affect the results below 20 Bohr but does not help to recover the correct behavior along Gamma-M and introduces more negative values.

Could you indicate how to generate a grid of denser q-points ?

Thanks,
Frederic

That's the safe way to go.

Run a ground state calculation (with nstep 1 and ecut as low as you want) with a k-point grid equal to the q-point grid you need, e.g.

prtvol 10
ngkpt 8 8 8
shiftk 0 0 0
kptopt 1

you will get from the output the irreducible k-points for your cell, which you can use as q-points for the phonons. Note that if you already have the DDB for a 4x4x4 grid you only have to do the _new_ q-points which were not in the previous grid. Add the new datasets on and merge the final DDBs for everybody. Obviously if you choose 6x6x6, say, you can't reuse the 0.25 0 0 q-point.

matthieu

Thanks matthieu.

Using a 8x8x8 grid adds a lot more q-points and thus as many more rf calculations.
Since the problem seems to be only for the Gamma-M branch, is it possible to only increase the number of q-points along this branch ?

Frederic

Another question: should I use the same number/grid of k-points as of q-points or these are completely unrelated ?

Frederic

Your q-grid and k-grid are better off if they are compatible: q is Gamma centered and a sub grid of the k.

However this is not obligatory, if it's not the case you need to calculate the wavefunctions at k+q explicitly (see tutorials).

matthieu