Problem extracting the head of eps^-1 - possible bug?

Biller's picture
Submitted by Biller on Thu, 04/17/2014 - 06:32

Forums 

Bug reports

Hi,
I am trying to look at the behavior of eps^-1(q)_00 for a 2D system with a truncated Coulomb interaction
(as seen in http://dx.doi.org/10.1103/PhysRevB.73.233103 and in http://www.nersc.gov/assets/Uploads/6-boronnitride.pdf)

However, I believe that the my output for the head of epsilon after inversion, eps(1,1), is not always the value for G=(0,0,0), G'=(0,0,0)...

Aside from the current problem of extracting the value, would this be a problem in sigma.x calculations? or does sigma.x treat epsmat by the actual G,G' value order that is given in epsmat itself?

I have verified this for various q points, all while inverting with scaLAPACK, for different ecuts and nbands, in the following way:

* For calculations where epsilon.log contained the line 0 0 0 0 0 0 under the " g gp inverse epsilon" report,
I made a comparison with the value reported in the output (e.g. epsilon.out) and it was not the same number.

* For some of my epsilon.log files, the reported Head after inversion matched the first line under the "g gp inverse epsilon" line, but that was not necessarily 0 0 0 0 0 0 (the elements in the report are just the ones processed by the head node, right?)

*For one calculation where epsilon.log did not contain the 0 0 0 0 0 0 line, I converted epsmat to ascii and looked at the matrix elements.
The value reported for the head was indeed the first value, and the actual head was at index 119341.

(my thoughts about this: I can only think of this happening if the scaLAPACK grid was somehow mangled. I'm using ifort 13.1.3 and mkl version 11.0.5.192)

jdeslip's picture

Submitted by jdeslip on Sun, 05/04/2014 - 16:29

HI Ariel,

Sorry for delay on this. You are right that epsilon.log prints only values owned by mpi-task0. For small q ~ 0, usually the G=0 vector is the first in the list.

I believe the code assumes when printing the head of epsilon/inverse-epsilon.

If you convert epsmat file to ascii (as you have done) you can look at the index of of gvectors for each q-point. Is G=0 not the first element for your qpoint?

The only thing I can think of that would cause this is if you have a large q, but it shouldn't affect the epsmat file other than making the "head" as reported by output meaningless.

-Jack

Biller's picture

Submitted by Biller on Tue, 05/06/2014 - 03:52

Nope, when I converted eps0mat the first element was not the head of epsilon.
I should probably look into it more, but it would have to wait a few weeks, since right now I have a pressing deadline.

However, I would like to comment that currently I am doing real epsilon calculations, separating chi and eps runs and using cray's libsci scalapack for the inversion, since intel's MKL pdgesv does not play nice (at the moment?) with cray's mpi. My original report was based on intel scalapack on our local cluster. It looks like that for this procedure the first element in the epsmat.ascii for each qpoint is in fact the value reported in the "head after inversion" line.
I am still pretty sure that this is not necessarily the 0,0 term, but it is hard to tell for this system. For example, for q=0.5, 0.5, 0, the first term in epsilon.log is for (-1,-1,0),(-1,-1,0) and it is 0.85525673. The value reported in "head after" is 0.855256730391331. However, the value for (0,0,0)(0,0,0) in epsilon.log is 0.85525672.

Another (seemingly) related issue is the following:

I have the following output in eps0mat for a cell slab truncation:

Head After Inversion = 0.994992007118876

This is also the value reported in epsilon.log for the head, and indeed the first value in the ascii file where actual matrix elements are reported, BUT it is clearly different from what is reported on the sigma stdout:
grepping eps\(1,1\) from a sigma stdout gives me the following values for the head of eps:

q= 0.000 0.000 0.000 n= 2561 eps(1,1)= 0.90
q= 0.000 0.125 0.000 n= 2582 eps(1,1)= 0.80
q= 0.000 0.250 0.000 n= 2606 eps(1,1)= 0.81
q= 0.000 0.375 0.000 n= 2608 eps(1,1)= 0.83
q= 0.000 0.500 0.000 n= 2624 eps(1,1)= 0.85
q= 0.250 0.000 0.000 n= 2564 eps(1,1)= 0.83
q= 0.250 0.125 0.000 n= 2590 eps(1,1)= 0.81
q= 0.250 0.250 0.000 n= 2615 eps(1,1)= 0.81
q= 0.250 0.375 0.000 n= 2636 eps(1,1)= 0.83
q= 0.250 0.500 0.000 n= 2638 eps(1,1)= 0.85
q= 0.500 0.000 0.000 n= 2564 eps(1,1)= 0.83
q= 0.500 0.125 0.000 n= 2618 eps(1,1)= 0.82
q= 0.500 0.250 0.000 n= 2624 eps(1,1)= 0.82
q= 0.500 0.375 0.000 n= 2618 eps(1,1)= 0.84
q= 0.500 0.500 0.000 n= 2636 eps(1,1)= 0.86

I have verified this behavior for multiple eps0mats and sigma calculations. There are several eps values in the ascii file with ~0.90 but they are certainly not the head.
I wonder if this somehow interferes with the truncation subroutines?

Btw grepping "Head After" from an epsmat inversion:

Head After Inversion = 0.802935519963003
Head After Inversion = 0.811607388977761
Head After Inversion = 0.834019346139109
Head After Inversion = 0.853921404731958
Head After Inversion = 0.834748010044180
Head After Inversion = 0.808666157570722
Head After Inversion = 0.814029645536077
Head After Inversion = 0.834247220721728
Head After Inversion = 0.853656431680774
Head After Inversion = 0.829404324344164
Head After Inversion = 0.819099494346391
Head After Inversion = 0.821188042623688
Head After Inversion = 0.837326686837760
Head After Inversion = 0.855256730391331

It is notable there are some points other than the first which are in slight disagreement, but that is probably due to round-off.

jdeslip's picture

Submitted by jdeslip on Sun, 05/11/2014 - 11:37

Hi Ariel,

OK, so about your first issue, when q > 0.5, then it is true that what is report as the head may not be the head and may not match eps(1,1). This shouldn't really cause other issues.

About your second issue, I think what is going on for the q = 0 head as reported by sigma. I believe that what is reported is written after we go through the fix_wings (which actually "fixes" the head, too) which does the following for head: epstemp(i) = wcoul0/vcoul

Basically, we do cell-averaging (for 2D using a model epsilon(q) in the mini BZ for q=0.). We adjust the value of epsilon so that when we later multiply vcoul we get the correct averaged screened Wcoul we want.