ABINIT Discussion Forums

Dear users/developers,
I am trying to calculate the phonon dispersion of graphane.
Prior to this calculation I've made an optimization of the Lattice parameter, and Ecut.
But when it comes to phonon dispersion it is just terrible! the dispersion looks Ugly and I have a negative frequencies around -1000cm-1.
The input for the phonon calculations:
-----------------------------------------------------------------------
#tolvrs 1.0d-9
prtvol 3
ndtset 28
#jdtset 1 25 26 27 28
# ground state
#getwfk1 1
kptopt1 1
tolvrs1 1d-16
prtden1 1
getwfk1 0 # Cancel default
rfphon1 0 # Cancel default
nqpt1 0
# Response Function calculation
nqpt 1
rfphon 1
rfatpol 1 4
rfdir 1 1 1
getwfk 1
kptopt 1
tolvrs 1.0d-9

qpt2 0.00000000E+00 0.00000000E+00 0.00000000E+00
qpt3 6.66666667E-02 0.00000000E+00 0.00000000E+00
qpt4 1.33333333E-01 0.00000000E+00 0.00000000E+00
qpt5 2.00000000E-01 0.00000000E+00 0.00000000E+00
qpt6 2.66666667E-01 0.00000000E+00 0.00000000E+00
qpt7 3.33333333E-01 0.00000000E+00 0.00000000E+00
qpt8 4.00000000E-01 0.00000000E+00 0.00000000E+00
qpt9 4.66666667E-01 0.00000000E+00 0.00000000E+00
qpt10 1.33333333E-01 6.66666667E-02 0.00000000E+00
qpt11 2.00000000E-01 6.66666667E-02 0.00000000E+00
qpt12 2.66666667E-01 6.66666667E-02 0.00000000E+00
qpt13 3.33333333E-01 6.66666667E-02 0.00000000E+00
qpt14 4.00000000E-01 6.66666667E-02 0.00000000E+00
qpt15 4.66666667E-01 6.66666667E-02 0.00000000E+00
qpt16 -4.66666667E-01 6.66666667E-02 0.00000000E+00
qpt17 2.66666667E-01 1.33333333E-01 0.00000000E+00
qpt18 3.33333333E-01 1.33333333E-01 0.00000000E+00
qpt19 4.00000000E-01 1.33333333E-01 0.00000000E+00
qpt20 4.66666667E-01 1.33333333E-01 0.00000000E+00
qpt21 -4.66666667E-01 1.33333333E-01 0.00000000E+00
qpt22 4.00000000E-01 2.00000000E-01 0.00000000E+00
qpt23 4.66666667E-01 2.00000000E-01 0.00000000E+00
qpt24 -4.66666667E-01 2.00000000E-01 0.00000000E+00
qpt25 -4.00000000E-01 2.00000000E-01 0.00000000E+00
qpt26 -4.66666667E-01 2.66666667E-01 0.00000000E+00
qpt27 -4.00000000E-01 2.66666667E-01 0.00000000E+00
qpt27 -4.00000000E-01 2.66666667E-01 0.00000000E+00
qpt28 -3.33333333E-01 3.33333333E-01 0.00000000E+00
# SCF

acell 4.7471698113 4.7471698113 15

rprim 1.0000000000 0.0000000 0.000000
0.5000000000 0.8660254037844386 0.000000
0.0000000000 0.00000000 1.000000

#DAT : Definition of the types of atoms

ntypat 2
znucl 6 1

#DA : Definition of the atoms

natom 4
typat 1 1 2 2
xred 1/3 1/3 0.400000/15
-1/3 -1/3 -0.400000/15
1/3 1/3 2.500000/15
-1/3 -1/3 -2.500000/15

#DPWset : Definition of the plane wave basis set
occopt 1
nband 14
ecut 38
diemac 6.00

#Dkptgrid : Definition of the k-points grid

#kptopt 1
ngkpt 15 15 1
nshiftk 1
shiftk 0 0 0
#toldfe 1e-8
#prtdos 2
#prtden 1
#band
#DSCFproc : Definition of the SCF procedure
#prtdos 2
nstep 50
----------------------------------------------------
And the input for anaddb:
------------------------------------------------------
!Input file for the ifc code. Analysis of a restricted SiO2 DDB
!Will use only the Gamma point information : compute
!Compute LO-TO splitting for different directions,
!as well as the dielectric tensor at zero frequency

!Input file for the ifc code. Analysis of the BN-sheet DDB
!Flags
dipdip 0
eivec 1
ifcflag 1 ! Interatomic force constant flag
asr 1 ! Acoustic Sum Rule. 1 => imposed asymetrically
symdynmat 1

!Wavevector grid number 1 (coarse grid, from DDB)
brav 1 ! Bravais Lattice : 1-S.C., 2-F.C., 3-B.C., 4-Hex.)
ngqpt 15 15 1 ! Monkhorst-Pack indices
nqshft 1 ! number of q-points in repeated basic q-cell
q1shft 0 0 0.0

!Effective charges
chneut 1 ! Charge neutrality requirement for effective charges.

!Wavevector list number 1 (Reduced coordinates and normalization factor)
nph1l 239 ! number of phonons in list 1

qph1l

-0.33333 0.33333 0.00000 1.00000
-0.33003 0.33003 0.00000 1.00000
-0.32673 0.32673 0.00000 1.00000
-0.32343 0.32343 0.00000 1.00000
-0.32013 0.32013 0.00000 1.00000
-0.31683 0.31683 0.00000 1.00000
-0.31353 0.31353 0.00000 1.00000
-0.31023 0.31023 0.00000 1.00000
-0.30693 0.30693 0.00000 1.00000
-0.30363 0.30363 0.00000 1.00000
-0.30033 0.30033 0.00000 1.00000
-0.29703 0.29703 0.00000 1.00000
-0.29373 0.29373 0.00000 1.00000
-0.29043 0.29043 0.00000 1.00000
-0.28713 0.28713 0.00000 1.00000
-0.28383 0.28383 0.00000 1.00000
-0.28053 0.28053 0.00000 1.00000
-0.27723 0.27723 0.00000 1.00000
-0.27393 0.27393 0.00000 1.00000
-0.27063 0.27063 0.00000 1.00000
----------- ----------- --------- ----------
---------------------------------------------------------------------------
This example is for 15x15 grid for kpts and qpts.
But I also tried the 30x30 for k-point but still have the same dispersion.
Many thanks for your help,
Samir

Before doing a finite q calculation, it might be a good idea to non-scf calculation with commensurable q mesh with the k mesh in G.S.

In other words, add the following 7 lines before phonone calculation. You can see the the two parts have same qpt, but one is non-scf, the other is scf.
nqpt2 1
qpt2 0.0 0.0 0.1
getwfk2 1 #from SCF high-accuracy calculation
getden2 1 #from SCF high-accuracy calculation
kptopt2 3
tolwfr2 1.0d-18
iscf2 -2 # Prepare for finite q phonon

After that, do the phonon calculation

rfphon3 1
rfatpol3 1 2
rfdir3 1 1 1
nqpt3 1
qpt3 0.0 0.0 0.1
getwfk3 1
getwfq3 2
kptopt3 3
tolvrs3 1.0d-8
iscf3 5

Then, when doing the finite q calculation, you can call "getwfq" for the non-scf DDK.

I'm not totally sure about your concerns, but it is the general procedure to be helpful.

What's more, I believe at the beginning you don't want to calculate the entire dispersion, but few points for tests; extending to the full dispersion can be done later.

Dear Ludwig,
I am really grateful for your help.
I'll check what you're suggesting but I am not sure that this problem comes directly from the fact that the qpts is not commensurable with respect to kpts.
I am experiencing the same problems for graphene. Indeed, as you might know graphene is much simpler than graphane. So I decided to reproduce graphene's phonons first. Especially since the structure is very similar (they have the same hex space group)
The phonon dispersion for graphene is even worse! as you can see it in the attached plot (the path is G-K-M-G)!. The input:
# phonons
prtvol 3
ndtset 28
udtset 14 2
# ground state
#Dataset 1 : SCF
# Chksymbreak 0
# rfasr 2
getden 11
getwfk 11
kptopt 3
# Response Function calculation
# qpt generation
nqpt 1
qptopt 1
ngqpt 10 10 1
nshiftq 1
shiftq 0.0 0.0 0.0
iqpt:? 1
iqpt+? 1
#Dataset 11 : SCF
getden11 0
getwfk11 0
iscf11 7
nqpt11 0
prtden11 1
###############
#Dataset 2 : phonon at gamma
getwfq12 0
##################

#Dataset ?1 : non-SCF
iscf?1 -2
tolwfr?1 1.0d-20
#Dataset ?2 : RF
getwfq?2 -1
rfatpol?2 1 1
rfdir?2 1 1 1
rfphon?2 1
tolvrs?2 1.0d-16
#######################################################################
#Common input variables
#######################################################################
acell 4.6511685 4.6511685 14
xred 2/3 1/3 0.0
1/3 2/3 0.0
rprim 1.0 0.0 0.0 # Hex. primitive vectors(to be scaled by acell)
-0.5 0.8660254037844386 0.0
0.0 0.0 1.0
#Definition of the atom types
ntypat 1 # There are two types of atoms
znucl 6 # The keyword "zatnum" refers to the atomic number of the
# possible type(s) of atom. The pseudopotential(s)
# mentioned in the "files" file must correspond
# to the type(s) of atom.
#Definition of the atoms
natom 2
typat 1 1
ngkpt 10 10 1
nshiftk 1 # Use one copy of grid only (default)
shiftk 0.0 0.0 0.0 # This choice of origin for the k point grid
# preserves the hexagonal symmetry of the grid,
# which would be broken by the default choice.
nband 16
nbdbuf 4
occopt 3
tsmear 0.005
#Definition of the planewave basis set
ecut 35.0 # Minimal kinetic energy cut-off, in Hartree
#Definition of the SCF procedure
nstep 500 # Maximal number of SCF cycles
diemac 12.0 # Although this is not mandatory, it is worth to
######################################
you can check in a separate test (you can just generate kpts for the GS using 10x10x1 grid for kptopt 1 and 3 and compare), that the qpts can be expressed as in terms of the kpts differences.
I spent a lot a lot of time on this problem (I am really suffering!) and I have the feeling that the problem might come from the 6.8.2 version of abinit.
So, I don't really know if installing the latest version could me helpful!. I already tried different sampling of kpts/qpts.......without success. I don't know what are the other parameter that could be relevant for phonons.
I really appreciate your help/suggestions.
Many thanks,
Sam.

Hello, I went through the input file and seems the essential elements are already there. So maybe let's cut the problem into several:

1) is the convergence test good? (from observation it is already enough)

2) is the acell(3)=14 large enough? (seems enough, too)

3) is the structure fully relaxed? (check fcart and stress in the out file or log file)

4) based on this structure, can the electronic structure be correctly reproduced?

If you don't mind, could you please take a look back to confirm the above, so we can focus that the problem is only in the phonon part?

If all of the above are correct, since I saw some negative frequencies, it might be from the acoustic sum rule; without which, there are negative frequencies.

Also, the ZA mode should be proportional to q^2.

If it still does not work, my suggestion is to follow two ways:

1) get several q points manually (with commensurate k mesh of non-scf Ground State), and compare it.

2) Calculate the dynamical matrix of commensurate in-equivalent q points(which q points needed to be calculated? just do a simple scf calculation, and grep "kpt"), and interpolate them with the phonon path you need.

I'm sure that with the acoustic sum rule, enough large k mesh (all convergence), you are quite close to the final solution.

Best,
Ludwig

Dear Ludwig,
Thank you so much for your help.
I think we can skip points 1 to 4 concerning the electronic structure of graphene. I already checked this and the DOS and the Band-structure are Ok.
Yes, I forced the ASR too in anaddb input (I checked both asr 1 and 2 but I got the same result I sent you).
Concerning phonons calculations I did not understand both points 1 and 2!
I don't understand what do you mean by:
1) get several q points manually (with commensurate k mesh of non-scf Ground State), and compare it.
for the input I already sent you I generated separately the following qpts:
0.00000000E+00 0.00000000E+00 0.00000000E+00
1.00000000E-01 0.00000000E+00 0.00000000E+00
2.00000000E-01 0.00000000E+00 0.00000000E+00
3.00000000E-01 0.00000000E+00 0.00000000E+00
4.00000000E-01 0.00000000E+00 0.00000000E+00
5.00000000E-01 0.00000000E+00 0.00000000E+00
1.00000000E-01 1.00000000E-01 0.00000000E+00
2.00000000E-01 1.00000000E-01 0.00000000E+00
3.00000000E-01 1.00000000E-01 0.00000000E+00
4.00000000E-01 1.00000000E-01 0.00000000E+00
2.00000000E-01 2.00000000E-01 0.00000000E+00
3.00000000E-01 2.00000000E-01 0.00000000E+00
4.00000000E-01 2.00000000E-01 0.00000000E+00
3.00000000E-01 3.00000000E-01 0.00000000E+00
By just using a usual input (GS) with ngkpt 10 10 1 and kptopt 1. And I checked for the same ngkpt but for kptopt 3 (as it should be for kpts when doing phonon calculations) and I found that the qpts (kptopt 1) can be expressed in terms of kpts (kptopt 3). I don't understand what do you mean by generate manually? and which should be compared with which?
For the second point concerning the dynamical matrix I don't see why I should calculate q point that are equivalent (by symmetry). I can understand that I have probably to check inequivalent q-point separately and then interpolate but there should not be a way to interpolate equivalent points since they have the same frequencies!. But again this is what I am doing so far. If I make in the input nqptopt 1 means that the 14 points I gave previously are inequivalent and then at the anaddb stage they will be interpolated according to the indicated path!. Probably I am missing something but that is my understanding!.
I would be happy if you could explain more these two points.
Many many thanks,
Sam

Hello Sam:

Sorry for the possible confusion. I mean, maybe we can first reduce the possible problematic range.

For example, if at certain q point, the phonon frequency is less than 0, we may first only check that single point (say, qx qy qz):

First, a scf calculation taking advantage of full symmetry:
kptopt1 1
tolvrs1 1.0d-18
iscf1 5
prtden1 1

Then, a non-scf calculation at finite q, and use entire B.Z

nqpt2 1
qpt2 qx qy qz
getwfk2 1
getden2 1
kptopt2 3
tolwfr2 1.0d-18
iscf2 -2

Finally, call the WFK of dataset1, and WFQ of dataset2

rfphon3 1
rfatpol3 1 2
rfdir3 1 1 1
nqpt3 1
qpt3 qx qy qz
getwfk3 1
getwfq3 2
kptopt3 3
tolvrs3 1.0d-8
iscf3 5

Then, finally at the output file, you can
grep -A3 "Phonon" filename.out
to see Phonon frequencies directly.

Can this get correct positive freq behavior?

This does not use futher mrgddb and anaddb, only abinit; also it does not interpolate, only use a finite q G.S., and call it to calculate phonon.
In other words, this is the most reduced version to calculate phonon at finite q.

If you still cannot get phonon correctly, then it might be a problem beyond our scope.

Dear Ludwig,
I really appreciate your help.
You're right!. I am surprised that I have already troubles during the first GS calculation Check.
scfcge : ERROR -
Potential-based CG line minimization not converged after 13 restarts.
Action : read the eventual warnings about lack of convergence.
Some might be relevant. Otherwise, raise nband
As you know nband is already enough (nband 16) so I will try other iscf (for example 7) and other k-point sampling other then 10x10x1.
What do you think?
Thanks,
Sam

Dear Ludwig,
I used iscf 7 instead of 5 and the three steps was done successfully as you suggested. The main output data are as follow:
loper3 : total number of electrons, from k and k+q
fully or partially occupied states are 8.000000E+00 and 8.000000E+00.
Phonon frequencies in cm-1 :
- 4.660460E+01 2.099668E+02 3.282793E+02 9.046315E+02 1.518869E+03
- 1.578322E+03
the 2DTE is also converged:
iter 2DEtotal(Ha) deltaE(Ha) residm vres2
-ETOT 1 23.630088896782 -2.389E+02 1.594E+00 2.949E+04
ETOT 2 11.757348810888 -1.187E+01 4.160E-02 4.168E+02
ETOT 3 11.559484232802 -1.979E-01 6.882E-04 1.498E+00
ETOT 4 11.559173911248 -3.103E-04 3.210E-05 2.370E-02
ETOT 5 11.559171583725 -2.328E-06 5.216E-07 5.066E-04
ETOT 6 11.559171555644 -2.808E-08 1.397E-08 5.825E-05
ETOT 7 11.559171535485 -2.016E-08 7.738E-10 2.206E-05
ETOT 8 11.559171521792 -1.369E-08 2.074E-11 2.351E-08
ETOT 9 11.559171521788 -4.178E-12 9.164E-13 6.588E-10
and the the dynamical matrix is calculated to be:
Dynamical matrix, in cartesian coordinates,
if specified in the inputs, asr has been imposed
j1 j2 matrix element
dir pert dir pert real part imaginary part

1 1 1 1 0.5765623597 0.0000000000
1 1 2 1 0.0243877579 -0.0000000130
1 1 3 1 -0.0000003641 0.0000000919
1 1 1 2 -0.5177520210 -0.1304228436
1 1 2 2 -0.0050870905 -0.0326155268
1 1 3 2 0.0000000000 0.0000000000
which means that the frequencies are positive. As this frequencies are calculated for q: 0.1 0.0 0.0 I compared them with the calculation I've done for the complete q-points (for which I got a negative frequencies).
loper3 : total number of electrons, from k and k+q
fully or partially occupied states are 8.000000E+00 and 9.388526E+00.
Initialisation of the first-order wave-functions :
Phonon frequencies in cm-1 :
- -2.576720E+03 -2.230206E+03 -2.178868E+03 -2.152276E+03 -1.351607E+03
- -1.318281E+03
The 2DTE:
ETOT 25 -506.70773439633 -4.547E-12 2.353E-20 1.265E-13
ETOT 26 -506.70773439633 -3.183E-12 1.245E-20 3.292E-14
ETOT 27 -506.70773439632 1.273E-11 4.967E-21 3.145E-15
ETOT 28 -506.70773439631 1.000E-11 1.031E-21 2.390E-16
ETOT 29 -506.70773439631 -4.547E-12 7.155E-23 6.962E-18
Dynamical matrix, in cartesian coordinates,
if specified in the inputs, asr has been imposed
j1 j2 matrix element
dir pert dir pert real part imaginary part

1 1 1 1 -1.4766922071 0.0000000000
1 1 2 1 0.0426574927 -0.0001824885
1 1 3 1 -0.0022453447 -0.0121550290
1 1 1 2 -0.6628288078 -0.1260569249
1 1 2 2 0.0228799948 0.0269300226
1 1 3 2 0.0000000000 0.0000000000
As you see the output data is completely different!. I am really confused. If I look at the previous input (with the total q-point) The second phonon wave-vector (0.1 0.0 0.0) is calculated in a very similar way as you suggested with the only difference that the calculation is done in three steps instead of two with the third step being performed in a self-consistent way(iscf 7 which takes a lot of time when done for all q-points). Does it mean that my problem come from the Non-SCF calculation?. If so there is no many choices in abinit because we have two methods only iscf -2 and -3 (which I think I already checked).
Thanks a lot,
Sam

Dear Sam:

Thanks for the post. I see that the ETOT is well converged. How about "2DEtotal" and "non-var. 2DEtotal" (they are listed right after the ETOT)? Are they the same?

If yes, this means that the Abinit itself is making correct result (positive), however, the Anaddb, which interpolates various matrices, give bad answer in certain direction. The acoustic sum rule is a usual source for negative freq problem, so you may even drop it (and charge neutrality) to see whether there is a real difference.

If not, there is one likely problem, that you may have to increase the k-mesh, to obtain more dynamical matrices of more in-equivalent K-points. Then you have more accuracy to do the dynamical matrix interpolation. One more further step can be done, to incorporate the DDB of the q=0.1 0.0 0.0(your example) calculated directly, to the mrgddb all in-equivalent K-points.

As to the non-SCF procedure, as long as we control the "tolwfr" low enough, and from a highly- accurate ground state, it should be OK.

Best,
Ljludwig

Dear Ludwig,
The ETOT values I provided are those of 2DTE...and you see that both (for the case where frequencies are not negative and the case where they are not) converge but the converged values are not the same!. Not only the 2DTE are different but also the dynamical matrix. This is because we did iscf 7 (with one q-point) and in the wrong calculation iscf is -2.........to me this is the only difference. Probably there is a problem for iscf -2 in this version (6.8.2).
For the anaddb is still early to judge whether the problem comes from there or not. The negative frequencies I provided they are obtained at the level of abinit already (the first run for phonons).
Thanks any way,
Sam.

Dear Sam:

I believe that you are correct. Before doing any interpolation, the in-equivalent q points themselves already have negative frequencies. Since they come from a scf ground state and non-scf DDK, we may have to set both ground state and ddk tolwfr to a same, very low value, and then calculate dynamical matrices.

iscf2 -3 #ddk calculation
kptopt2 2
rfphon2 0
rfelfd2 2
tolwfr2 1.0d-25

This might be the final idea I have; if it still doesn't work, maybe the version may have to be changed.
Best,

Dear Ludwig,
I am very happy to tell you that I finally calculated the phonon dispersion of graphene (have a look at the attached plot).
I have just to add kptopt1 1. which means that the problem at Gamma (negative frequency) comes from the fact that kptopt was set to 3 (including Gamma point).
I am targeting the graphane (hydrogenated graphene) phonon dispersion. I did the calculation for graphene because it is simpler and very similar to graphane (since we have the same hex structure with two carbon and hydrogen sub-lattices). So I proceeded in the same way as you suggested, just for one k-point near the Gamma point, for-example. And with the same graphene phonon input I am still having some problems!.
The frequencies near the Gamma point are usually negative and when I get read of them I don't have the correct number of non-zero modes at the Gamma point!.
For example, in the letter (Appl. Phys. Lett. 98, 051914 2011) there is 9 non-zero modes while I always find only 6 non-zero modes (around the Gamma point)!!! So I am wondering what is the parameter that control the number of modes (at Gamma)?
I've already made a full geometry optimization and relaxation before the phonon calculation.
Is there any special setting for phonon-calculations in a semi-conducting material!
Your help was really useful for me, and hope that we could solve the problem for graphane as well.
Thanks,
Sam

Dear Sam:

Glad to hear that your problem is solved. I still don't get it- at the very beginning of your post, isn't kptopt1 already added? Could you please provide more detail, such as input file for future reference of this forum? Many thanks.

Are you sure this new hydrogenated graphene is thermodynamically stable? If we look at the procedure how phonon is calculated, then it is simply displacement of atoms. The extent of displacement is denoted by wavevector. If the structure is not stable, then the distorted, displaced structure is possible to have lower energy the the "ground state", especially near Gamma. Since Gamma is long-wavelength, such distortion is only small shift to "equilibrium position". If the "equilibrium" ground state is not Real ground state with stability, then negative frequency occurs.

The modes are the same in different k points- some are degenerate. It depends on the atoms of the system. For instance the simplest one, 6 branches are 3 acoustic and 3 optic. 3 modes are Z direction, transverse , longitudinal. Assume you have 4 atoms in one u.c. (simple example), A, A, B,C, unlike graphene (A, A), changing A to another A gets back to same result, here two A are no longer the same (which is definitely related to symmetry). The AB, AC difference is then revealed as another 3 branches.

As to semiconductor, I don't think there is special setting. Mainly, the occupation of electrons are in k-space, not real space, so it should affect little to atom movement. But there does have something special for the k-mesh. Different crystal symmetry can perform better with a proper choice of k-grids.

So far as I see you might want to compare the Raman with and without the additional layer. That's a good idea.

Cheers!

Dear Ludwig,
Sorry if my statement was not clear!
Indeed, I have just to add a line (with the complete qpts-data input): kptop11 1 as follows:
# phonons
prtvol 3
ndtset 28
udtset 14 2
# ground state
#Dataset 1 : SCF
#Chksymbreak 0
#rfasr 2
getden 11
getwfk 11
kptopt 3
# Response Function calculation
# qpt generation
nqpt 1
qptopt 1
ngqpt 10 10 1
nshiftq 1
shiftq 0.0 0.0 0.0
iqpt:? 1
iqpt+? 1
#Dataset 11 : SCF
getden11 0
getwfk11 0
iscf11 7
nqpt11 0
prtden11 1
--> kptopt11 1
----
-------------------------
If you compare this input to the previous one (the one with negative G-frequency) you find that the only difference is that for the previous setting I had: kptopt 3, but for the correct one (the present one): kptopt 3 with kptopt11 1 (which correspond to the Gamma point).
Thanks,
Sam