ABINIT Discussion Forums

I am working on how to calculate the second-order thermal corrections to eigenvalues. Following Doctor Boulanger's lecture, I obtain the *_TBS file. However, all the values are zero. I guess that I may ignore the use of foo_EIGR2D. But there is no such utility in abinit. So I turn to guys for help!
Thanks for your time and I quite appreciated it!
Best Wishes!

Hi,

I see this is an old post but I am having the same problem.

Best,
Jean Paul

Dear Jean Paul,

If your version of Abinit is recent enough, you can use the python script to post-process your _EIGR2D and _FAN files.

The scripts can be found in abinit/scripts/post_processing/temperature_para.py
This script requires abinit/scripts/post_processing/rf_mods.py

You can then run the script and follow the instructions.

Best,

Samuel

Dear Samuel,

Thanks for your reply and sorry for the delay. I guess I was expecting to receive an email when there was a new reply to the post (I wish I had seen your reply earlier because I have been stuck with this forever). I checked the control panel of the forum but I am not sure I can do that.

I am interested in calculating the thermal corrections to the electronic eigenvalues. Do I need to have ifcflag 1 as it says in the tutorial?

When I select ifcflag 1 in the anaddb.in file, I get the message:

The dynamical matrix number 2 cannot be built,
since no blok with wavevector 0.000000E+00 1.250000E-01 1.250000E-01

When I select ifcflag 0, then the calculation finishes. But I get all 0. I guess I don't need the forces anyway? Or I do, and that's why I'm getting all 0? The q points I have are the ones obtained by using kptopt 1, with a 8 8 8 grid, shifted to 0 0 0. This corresponds to 29 points. Being asked for extra q points suggests I would need the 512 points that are obtained without using any symmetries, as opposed to 29? Maybe abinit expects to have all 512 q points when calculating forces and nsym=1 in the ddb file?

I was now actually trying to use all 512 points to see if I could avoid these problems. But when I do so, I get the wrong Bravais lattice for each q point. I get messages like:

The Bravais lattice determined from modified primitive
vectors, bravais(1)= 2, has a lower symmetry than before,
but is still more symmetric than the real one, iholohedry= 1
obtained by taking into account the atomic positions.

and then for example

DATASET 91 : space group P1 (# 1); Bravais aP (primitive triclinic)

but that is wrong. The calculation wasn't finishing. Now it finished for the first 100 points but I can't trust it since the Bravais lattice is not the right one.

--------------------------------------------------------
Other parameters in the input file are (as in the last part of the tutorial of response function 2):
brav 2
ngqpt 8 8 8
nqshft 1
q1shft 3*0.0

ng2qpt 20 20 20 # I don't know what are typical good values for this...
ngrids 5 # Nor this
q2shft 3*0.0

ntemper 10
temperinc 10.
tempermin 10.

symdynmat 0

For nchan, nwchan and thmtol I just left the default values.
----------------------

I will try the scripts now but I would really appreciate if you could give me some feedback on the above issues.

Regards,
Jean Paul

Hi!

Maybe I can give you an answer at least to a part of your question.

charlesdarnay wrote:The dynamical matrix number 2 cannot be built,
since no blok with wavevector 0.000000E+00 1.250000E-01 1.250000E-01

This message says that you merged DDBs for wrong q-points (at least one q-point). Please, check whether the list of q-points in your response function calculation is the same as the one you got from your trial calculation for 8x8x8 mesh. Pay special attention to minuses in the first column.

charlesdarnay wrote:When I select ifcflag 0, then the calculation finishes. But I get all 0. I guess I don't need the forces anyway? Or I do, and that's why I'm getting all 0? The q points I have are the ones obtained by using kptopt 1, with a 8 8 8 grid, shifted to 0 0 0. This corresponds to 29 points.

When you put ifcflag to 0, anaddb did not use merged DDB for analysis, but each DDB for itself, regardless the q-point it corresponds to. But you probably need the merged DDB for thermal corrections. I guess you have to solve the above problem first.

Best regards,

Igor

Hi! Thanks for the prompt reply.

The q points that I have are the ones generated with kpoint 1, a 8 8 8 grid and nshift 0. The points are the following:

qpt3 0 0 0
qpt4 1.25000000E-01 0.00000000E+00 0.00000000E+00
qpt5 2.50000000E-01 0.00000000E+00 0.00000000E+00
qpt6 3.75000000E-01 0.00000000E+00 0.00000000E+00
qpt7 5.00000000E-01 0.00000000E+00 0.00000000E+00
qpt8 1.25000000E-01 1.25000000E-01 0.00000000E+00
qpt9 2.50000000E-01 1.25000000E-01 0.00000000E+00
qpt10 3.75000000E-01 1.25000000E-01 0.00000000E+00
qpt11 5.00000000E-01 1.25000000E-01 0.00000000E+00
qpt12 -3.75000000E-01 1.25000000E-01 0.00000000E+00
qpt13 -2.50000000E-01 1.25000000E-01 0.00000000E+00
qpt14 -1.25000000E-01 1.25000000E-01 0.00000000E+00
qpt15 2.50000000E-01 2.50000000E-01 0.00000000E+00
qpt16 3.75000000E-01 2.50000000E-01 0.00000000E+00
qpt17 5.00000000E-01 2.50000000E-01 0.00000000E+00
qpt18 -3.75000000E-01 2.50000000E-01 0.00000000E+00
qpt19 -2.50000000E-01 2.50000000E-01 0.00000000E+00
qpt20 3.75000000E-01 3.75000000E-01 0.00000000E+00
qpt21 5.00000000E-01 3.75000000E-01 0.00000000E+00
qpt22 -3.75000000E-01 3.75000000E-01 0.00000000E+00
qpt23 5.00000000E-01 5.00000000E-01 0.00000000E+00
qpt24 3.75000000E-01 2.50000000E-01 1.25000000E-01
qpt25 5.00000000E-01 2.50000000E-01 1.25000000E-01
qpt26 -3.75000000E-01 2.50000000E-01 1.25000000E-01
qpt27 5.00000000E-01 3.75000000E-01 1.25000000E-01
qpt28 -3.75000000E-01 3.75000000E-01 1.25000000E-01
qpt29 -2.50000000E-01 3.75000000E-01 1.25000000E-01
qpt30 -3.75000000E-01 5.00000000E-01 1.25000000E-01
qpt31 -2.50000000E-01 5.00000000E-01 2.50000000E-01

That's 29 points. You can see 0 0.125 0.125 is not one of those. If I instead use kptopt 3, that gives 512 points (I am following the tutorial, generating the points with the trf2_2.in file (although actually that prints only the first 50 points, I got the 512 points from some other output file)). And 0 0.125 0.125 is one of those 512 points. That's why I was wondering, as the symmetries are not implemented in a thermal calculation (that's the error message I get when I don't use nsym 1 in the input file), so the individual ddb files have nsym=1 and thus the merged ddb file has nsym=1 as well, if maybe then I have to use all 512 q points.

I am now trying to use Samuel's script, but I am having trouble getting the netCDF4 module...I am pretty new to Linux.

Best,
Jean Paul

Hi!

If you need nsym=1, then I suppose that you have to use all 512 k-points (unfortunately). I hope some of the elders will confirm this.

In that case, try to change brav in your anaddb input file. Try first brav=1. And maybe you'll have to play with shifts, too.

Igor

I think I will try now just using the 512 points. And yes, it would be great if someone could confirm this. For phonons one has to set kptopt 3, but doesn't need to set nsym to 1. I should think more about this, but why is this? Anaddb seemed to be working just fine with the 29 points (I was assuming that doing the thermal calculation was just a matter of adding the ieig2rf line). However, as I said before, I get messages of the type:

--- !COMMENT
message: |
The Bravais lattice determined from modified primitive
vectors, bravais(1)= 2, has a lower symmetry than before,
but is still more symmetric than the real one, iholohedry= 1
obtained by taking into account the atomic positions.
...
symlatt : the Bravais lattice is hR (rhombohedral)

So "the bravais lattice is hR" is wrong (it's zincblende, spgroup=216 actually). However, then it says

space group P1 (# 1); Bravais aP (primitive triclinic)

which I guess it's actually correct (as opposed to what I said earlier), since in the input file I set nsym=1. So I'm not sure what to make out of this. So (i) Presumably Abinit tries to find the possible symmetries automatically, regardless of the symmetry one establishes in the input? Maybe trying to check if the input file information is consistent? However, the "real" symmetry is not rhombohedral...

(ii) Also, I only just found out that by default istwfk is 0 for the ground state, using time reversal symmetry, but 1 for response function calculations, not using time reversal symmetry. I don't recall seeing this in the tutorial. So there are 260 points in the ground state, but 512 in the other datasets. The calculation finishes anyway. This is actually a more general question, not directly related with thermal corrections, but if the calculation finishes, can one trust it in general? Should I change istwfk to 1 in the ground state? I think I remember reading in some tutorial that Abinit will try to finish the calculation in spite of the warnings. I had gotten the impression it is fine, but now I am wondering about this again.

(iii) I should probably ask this in a Linux forum (so sorry if this question is out of bounds; I'll keep it brief), but I have problems with installing the netCDF4 module. Basically I seem to get all links to the netCDF4 python, but in the readme they say to make sure I have installed the netcdf-4! I basically found a version with a readme that has a link to ftp://ftp.unidata.ucar.edu/pub/netcdf/netcdf-4, so I will try that next :/.

Best,
Jean Paul

Hi!

charlesdarnay wrote:I think I will try now just using the 512 points. And yes, it would be great if someone could confirm this. For phonons one has to set kptopt 3, but doesn't need to set nsym to 1. I should think more about this, but why is this?

Using symmetries is not appropriate for RF, thus kptopt=3 and nsym=1. This is coded somewhere in the routines and has reasons described (probably) in the papers by X. Gonze (and others) from mid 90's.

charlesdarnay wrote:Anaddb seemed to be working just fine with the 29 points (I was assuming that doing the thermal calculation was just a matter of adding the ieig2rf line). However, as I said before, I get messages of the type:

--- !COMMENT
message: |
The Bravais lattice determined from modified primitive
vectors, bravais(1)= 2, has a lower symmetry than before,
but is still more symmetric than the real one, iholohedry= 1
obtained by taking into account the atomic positions.
...
symlatt : the Bravais lattice is hR (rhombohedral)

So "the bravais lattice is hR" is wrong (it's zincblende, spgroup=216 actually). However, then it says

space group P1 (# 1); Bravais aP (primitive triclinic)

which I guess it's actually correct (as opposed to what I said earlier), since in the input file I set nsym=1. So I'm not sure what to make out of this. So (i) Presumably Abinit tries to find the possible symmetries automatically, regardless of the symmetry one establishes in the input? Maybe trying to check if the input file information is consistent? However, the "real" symmetry is not rhombohedral...

I don't know what to say about the COMMENT message, but the space group info you got (P1) is consistent with the choice of input variables from your input file (nsym=1). Abinit tries to find out symmetries based on the input file variables (rprim, acell, xred, spgroup, brvltt,...), but I think it takes none for nsym=1.

charlesdarnay wrote:(ii) Also, I only just found out that by default istwfk is 0 for the ground state, using time reversal symmetry, but 1 for response function calculations, not using time reversal symmetry. I don't recall seeing this in the tutorial. So there are 260 points in the ground state, but 512 in the other datasets. The calculation finishes anyway. This is actually a more general question, not directly related with thermal corrections, but if the calculation finishes, can one trust it in general? Should I change istwfk to 1 in the ground state? I think I remember reading in some tutorial that Abinit will try to finish the calculation in spite of the warnings. I had gotten the impression it is fine, but now I am wondering about this again.

As in the first answer, symmetries are not appropriate for RF, contrary for GS. A part of the answer can be found at

http://www.abinit.org/documentation/helpfiles/for-v7.10/input_variables/vardev.html#istwfk

Since you can use symmetries in GS part of your calculation, istwfk=0, while in RF part it is forced to 1. Yes, you can trust it, since using symmetries only reduce the duration of calculation (need confirmation here). You can put istwfk=1 in GS, but I think you will get the same result and your calculation will only last longer (need confirmation here). Yes, Abinit will finish with warnings. One has to careful with them, thou, meaning one should understand whether they influence the final result or not in some way.

charlesdarnay wrote:(iii) I should probably ask this in a Linux forum (so sorry if this question is out of bounds; I'll keep it brief), but I have problems with installing the netCDF4 module. Basically I seem to get all links to the netCDF4 python, but in the readme they say to make sure I have installed the netcdf-4! I basically found a version with a readme that has a link to ftp://ftp.unidata.ucar.edu/pub/netcdf/netcdf-4, so I will try that next :/.

Dunno anything about that.

Igor

Hi again.

Thanks Igor for keeping up with the follow ups.

So I could complete the calculation for the 512 points. For some reason when I divided the points into different sets (different .in files), to have each of them take less amount on the cluster, the calculation would only finish for the first set. I wonder if having the gamma point in the set for example, somehow changes what abinit does for the following q points? In any case, I had more time in the cluster than I thought of, and eventually I ran the calculation for all of the points together and it finished. I am just writing this in case someone stumbles upon a similar problem in the future.

I am using a cutoff of 1100 ev and a 8 8 8 k-grid shifted to 0. Each of the EIGR2D files has around 16 MB. So it's around 8GB in total. The merged file, using mrgddb, has 4.4 GB. Can this be a problem for abinit? I am just wondering because when I try to open it with emacs, I get "Buffer exceeds maximum size". Anyway, the error I get is:

BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= EXIT CODE: 24

A few lines above, it says

forrtl: severe (24): end-of-file during read, unit 2, file /address_to_file/2_3.eigr2d.out

2_3.eigr2f.out is the merged file of the _EIGR2D files.

-----------------------------------------
I have the following files:

*anaddb.files:
naddb.in
anaddb.out
2_3.ddb.out
band2eps1
2_3.eigr2d.out
temp_gap_therm1
temp_gap_dummy1

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

!Flags
ifcflag 1 ! Interatomic force constant flag
thmflag 3 ! Thermodynamical properties flag
eivec 4

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

!Effective charges
asr 1 ! Acoustic Sum Rule. 1 => imposed asymetrically
chneut 1 ! Charge neutrality requirement for effective charges.

!Interatomic force constant info
# dipdip 1 ! Dipole-dipole interaction treatment

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

!Thermal information
ntemper 10 ! Number of temperatures
temperinc 10. ! Increment of temperature in K for temperature dependency
tempermin 10. ! Minimal temperature in Kelvin

-------------------------
It seems I had installed the netCDF4 module, but there were lots of problems when doing make check. So I still can't use temperature_para.py.

So hopefully (well, or not) I just made a silly mistake in the input file. Or rather, hopefully this has an easy solution...

Jean Paul

Hi!

charlesdarnay wrote:So I could complete the calculation for the 512 points. For some reason when I divided the points into different sets (different .in files), to have each of them take less amount on the cluster, the calculation would only finish for the first set. I wonder if having the gamma point in the set for example, somehow changes what abinit does for the following q points?

You have to do the same ground state calc for each such division of q-points. Gamma point in other divisions is not essential concerning the pure phonon calculation. So, the other division calcs did not finish because of some other reason.

charlesdarnay wrote:I am using a cutoff of 1100 ev and a 8 8 8 k-grid shifted to 0. Each of the EIGR2D files has around 16 MB. So it's around 8GB in total. The merged file, using mrgddb, has 4.4 GB. Can this be a problem for abinit? I am just wondering because when I try to open it with emacs, I get "Buffer exceeds maximum size".

Not to my experience. But maybe it is a problem for your local machine which you use to open these files. These are big files.

Igor

Hi Igor,

Yeah, I did the same ground state calculation for each of the separate files. So I don't know why it wasn't finishing. Anyway, I have the calculation done now.

I tried using anaddb in my computer as opposed to the cluster and, as expected, the calculation also just ended. It finishes in "The interatomic forces have been obtained". No thermal calculation is carried out (there are no TBS, PDS nor G2F files). So as anaddb is still not working, I will add the following paragraph...

I still haven't figured out the netCDF4 issue to use temperature_para.py. I guess this is a silly python question. But first I got the message "No module named rf_mods". Abinit has a file rf_mods.py which I copied to the folder. And then the terminal wasn't complaining anymore about this module. Then I got the "ImportError: No module named netCDF4". Abinit has a file named netcdf.py (copying it to the working directory doesn't work though) in /scripts/post_processing/appa/gui. I was wondering if this file has anything to do with the issue...

Thanks,
Jean Paul

Hi!

I'm not sure, but this might be useful. From input variables I see

http://www.abinit.org/documentation/helpfiles/for-v7.10/users/anaddb_help.html#thmflag wrote:If thmflag==3, the list of phonon wavevector from the first list is used (with equal weight for all wavevectors in this list)

WARNING : The use of symmetries for the temperature dependence of the eigenenergies is tricky ! It can only be valid for the k points that respect the symmetries (i.e. the Gamma point), provided one also averages over the degenerate states.
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.

But in your anaddb input I cannot see anywhere the list of q-points

http://www.abinit.org/documentation/helpfiles/for-v7.10/users/anaddb_help.html#qph1l

charlesdarnay wrote:anaddb.in:
!Input file for the anaddb code. Analysis of the SiO2 DDB

!Flags
ifcflag 1 ! Interatomic force constant flag
thmflag 3 ! Thermodynamical properties flag
eivec 4

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

!Effective charges
asr 1 ! Acoustic Sum Rule. 1 => imposed asymetrically
chneut 1 ! Charge neutrality requirement for effective charges.

!Interatomic force constant info
# dipdip 1 ! Dipole-dipole interaction treatment

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

!Thermal information
ntemper 10 ! Number of temperatures
temperinc 10. ! Increment of temperature in K for temperature dependency
tempermin 10. ! Minimal temperature in Kelvin

For netcdf, maybe you should repost that issue to another forum section. Someone proficient in that is maybe not looking at this section.

Igor

Hi Igor,

Thanks again for offering help!

So it says that the variables that may be needed are: dostol, nchan, ntemper, temperinc, tempermin, ng2qpt, ngrids, q2shft.

I didn't specify dostol, because it has a default value anyway. And I did specify the other variables. Also, I don't get any error saying that I should specify certain variable. I think I don't have to specify qph1l? It's not one of the variables mentioned above. I think I would have to specify it if I want to plot the phonon band structure. The help file says "List of nph1l wavevectors, at which the phonon frequencies will be interpolated". When I was using 29 q points and getting all 0, the file had a list of values for each k point in the ddb file. And I think I wasn't specifying nph1l (as in the tutorial, I was getting some temperature independent and dependent "results" (all 0s). Putting the cart before the horse, I don't know what those values are. There are 8 values for each k point, and then 4 columns for the T indep part, and 2 for the T dep part...). Anyway, I tried it, but the calculation stops at the same place.

When I select ifcflag 1, I get extra lines of the form:

phonon wavelength (reduced coordinates) , norm, and energies in hartree
0.00 0.00 0.00 1.00
0.000000000E+00 0.000000000E+00 0.000000000E+00 2.325213299E-03 2.325226953E-03
2.325232007E-03
Zero Point Motion energy (sum of freqs/2)= 3.487836129E-03

but the log file again finishes in:

inprep8 : open file 2_3.eigr2d.out
inprep8 : nband(1)= 14
occ 0.20000000000000D+01 0.20000000000000D+01 0.20000000000000D+01
0.20000000000000D+01 0.20000000000000D+01 0.20000000000000D+01
0.20000000000000D+01 0.20000000000000D+01 0.20000000000000D+01
0.00000000000000D+00 0.00000000000000D+00 0.00000000000000D+00
0.00000000000000D+00 0.00000000000000D+00

So it seems there is a problem with the eigr2d.out file? Do I really need ifcflag?

Before, in the cluster, I was getting the message "forrtl: severe (24): end-of-file during read, unit 2, file /address_to_file/2_3.eigr2d.out". And I'm not even getting that (and there is no Abinit error or warning message) now in my computer..

Summary:
There is no error or warning message in the log file.
The calculation just finishes after the above lines that start with "inprep8 : open file 2_3.eigr2d.out ... "

Help.

Jean Paul

PS: By the way, http://www.abinit.org/documentation/hel ... elp.html#1 says the ddb file has derivatives of the total energy with respect to phonons (in this case). So the DDB has the interatomic force constant matrix, as (81) in Baroni, Gironcoli, Dal Corso, REVIEWS OF MODERN PHYSICS, VOLUME 73, APRIL 2001? Or does it give some portion of (86)? I was wondering, since I want to know what do the eigr2d files have. Is it $\epsilon^{(2)}$ as in Boulanger's tutorial? Rather, I guess it's the derivative of the energy with respect to the density (Debye + Fan).

Hopefully someone can help with this simple question.

In the tutorial https://wiki.abinit.org/doku.php?id=howto:temperature which unfortunately I wasn't aware of initially, it seems that using the temperatura_para.py script, symmetries are implemented. For example, in temperature3.files, there are files corresponding to 8 q points, instead of 24.

I was trying to follow that tutorial using anaddb because I don't have Abinit installed with netCDF4...But, anaddb says I have files missing in the DDB file (it does work for one q point as in the tutorial, but the tutorial doesn't explain how it would work for more q points using anaddb).

So, symmetries are not implemented in anaddb, but they are implemented in the temperature_para.py script? Or is there a way to work around it with anaddb? Is that the main advantage of the script?

Thank you,
Jean Paul

Anaddb has symmetries implemented. But usually you have different _DDB files for your q points.
Therefore you need to use mrgddb executable to merge all the _DDB files before using anaddb.
If your DDB file does not include a q point that is in your anaddb input file, then it will complain.

Cheers'

Jordan

Hi Jordan,

Thanks for the reply. I know I have to use mrgddb to merge the different q points. So I'm not sure I understand your answer.

Are you saying that I could only use 8 q points with anaddb for ngqpt 4 4 4, but that those 8 q points are not the same that abinit generates?

I mean, it seems I can only use 8 q points with the temperature_para.py script. But I need 24 points if I use anaddb.

Regards,
Jean Paul