shifted grids method for dipoles

[Xiaoming Wang]

Hi Davide,

Thanks for your suggestions. I figured out that it is not the q shift magnitude that matters as you mentioned. The problem is that by using SVD, under which the overlap matrix elements are in parallel transport gauge as in David Vanderbilt's book, the diagonal matrix elements are supposed to be near 1 and the off-diagonal ones are supposed to be close to 0! The ypp renormalization issue can be solved by commenting out the if clause.

Indeed, I'm more interested in the diagonal dipole matrix elements difference (A_cc - A_vv in the following fig). So with the shifted grid method implemented in Yambo, if we expand the u_nk+q to second order, we have

Screen Shot 2021-10-18 at 2.17.23 AM.png

The alpha and gamma can be obtained from a shifted grid run, while beta is calculated by another -q shift. Note that I comment out the symmetrization procedure of the matrix elements and set the loops over valence and conduction bands to all the bands while calculating the overlap matrix elements. The problem is that while the A_vc is quite robust, the A_cc - A_vv is however very sensitive to the q shift. Are these due to the numerical accuracy of the wave functions? Btw, I'm using the double precision version of yambo and the conv_thr for QE calculations is set to 1.d-13. Any comment is appreciated.

Best,
Xiaoming

[Xiaoming Wang]

Maybe the shift q is too small. So I tried to increase q to 0.01. But the dipole calculation stopped with the following error.

Code: Select all

[ERROR] STOP signal received while in[04] Dipoles
[ERROR] Reading File shiftx/SAVE//ns.wf_fragments_1_1; Variable WF_COMPONENTS_@_SP_POL1_K1_BAND_GRP_1; NetCDF: Start+count exceeds dimension bound

Any idea?

[Davide Sangalli]

Dear Xiaoming Wang

Thanks for your suggestions. I figured out that it is not the q shift magnitude that matters as you mentioned. The problem is that by using SVD, under which the overlap matrix elements are in parallel transport gauge as in David Vanderbilt's book, the diagonal matrix elements are supposed to be near 1 and the off-diagonal ones are supposed to be close to 0! The ypp renormalization issue can be solved by commenting out the if clause.

All makes sense to me.
For the renormalization I guess you mean these lines in ypp/k-points/k_grids.F

Code: Select all

       if (v_norm(U_q_shifts(i1,:))>1.E-4.or.v_norm(U_q_shifts(i1,:))<1.E-6) then
         q0 = U_q_shifts(i1,:)*1.E-4/v_norm(U_q_shifts(i1,:))
         call msg("s","Renormalizing shift to :",v_norm(q0))
       else
         q0 = U_q_shifts(i1,:)
       endif

and you directly set q0 = U_q_shifts(i1,:)
I think it was there since
- for q<1.E-6 single precision error starts to give issues in the off-diagonal dipoles
- for q>1.E-4 the distance between k and k+q starts to be too big to get a good approximation for the Taylor expansion
I would expect this range to be optimal also for what you need

Indeed, I'm more interested in the diagonal dipole matrix elements difference (A_cc - A_vv in the following fig). So with the shifted grid method implemented in Yambo, if we expand the u_nk+q to second order, we have
..
The alpha and gamma can be obtained from a shifted grid run, while beta is calculated by another -q shift. Note that I comment out the symmetrization procedure of the matrix elements and set the loops over valence and conduction bands to all the bands while calculating the overlap matrix elements. The problem is that while the A_vc is quite robust, the A_cc - A_vv is however very sensitive to the q shift. Are these due to the numerical accuracy of the wave functions? Btw, I'm using the double precision version of yambo and the conv_thr for QE calculations is set to 1.d-13. Any comment is appreciated.

I do not have specific comments. Double precision and conv_thr look good to me.

Maybe the shift q is too small. So I tried to increase q to 0.01. But the dipole calculation stopped with the following error.

Code: Select all

[ERROR] STOP signal received while in[04] Dipoles
[ERROR] Reading File shiftx/SAVE//ns.wf_fragments_1_1; Variable WF_COMPONENTS_@_SP_POL1_K1_BAND_GRP_1; NetCDF: Start+count exceeds dimension bound

Any idea?

Before loading the database of the shifted grids, the code checks that it has the same number of components.
Can it be that the number of bands differs between the two databases (shifted and non-shifted)?
Try to do

Code: Select all

ncdump -h SAVE/ns.wf_fragments_1_1

to check if the size are the same.

Best,
D.

[Xiaoming Wang]

Dear Davide,

Thanks for your detailed reply. For the last point, the two wave functions did have different components but same bands. However, I have already commented out the checking procedure in the subroutine. Indeed, For smaller q shift of 0.001, the two components are also different, but it works. I don't know why for the q shift of 0.01 it doesn't.

Best,
Xiaoming

[Xiaoming Wang]

Btw, does the difference of the components matter? I mean which grid has larger components?

Best,
Xiaoming

[Davide Sangalli]

The number of components is defined by the DFT code to respect |k+G|^2<=E_wf_cutoff.
Now shifting the k-points by delta, it may happen that for some k_i, G_j for which |k_i+G_j|^2>E_wf_cutoff, G_j becomes allowed, i.e. |k_i+delta+G_j|^2<=E_wf_cutoff (or viceversa)
I would say it does not really matter, but you can eventually check that by playing with the cutoff.

D.

[Xiaoming Wang]

Thanks Davide. I would expect to have different components for larger shifts anyway. The problem is that sometimes it works, for example,

Code: Select all

 <---> loading info from shiftx
 <---> [WARNING] Shifted grid in folder "shiftx":change in number of WFs components
 <---> Grid 1 has shift vector        :  0.100000E-2   0.00000      0.00000    [c.c.]
 <---> loading info from shifty
 <---> [WARNING] Shifted grid in folder "shifty":change in number of WFs components
 <---> Grid 2 has shift vector        :   0.0000      0.10000E-02 -0.44452E-17 [c.c.]
 <---> loading info from shiftz
 <---> [WARNING] Shifted grid in folder "shiftz":change in number of WFs components
 <---> Grid 3 has shift vector        :   0.00000      0.00000     0.100000E-2 [c.c.]
 <09s> Dipoles shifted grids |########################################| [100%] 08s(E) 08s(X)
 <45s> Dipoles: P (T) |########################################| [100%] 35s(E) 35s(X)
 <45s> [DIP] Writing dipoles header
 <45s> [05] Optics
 <45s> [LA] SERIAL linear algebra
 <45s> [DIP] Checking dipoles header
 <45s> [X-CG] R(p) Tot o/o(of R):    712   7104    100
 <45s> [X] Upper matrix triangle filled
 <46s> Xo@q[1] |########################################| [100%] --(E) --(X)
 <46s> [06] Timing Overview
 <46s> [07] Memory Overview
 <46s> [08] Game Over & Game summary

but sometimes it does not.

Code: Select all

 <---> [04] Dipoles
 <---> [DIP] Checking dipoles header
 <---> [WARNING] DIPOLES database not correct or missing
 <---> loading info from shiftx
 <---> [WARNING] Shifted grid in folder "shiftx":change in number of WFs components
 <---> Grid 1 has shift vector        :  0.10000E-01  0.25237E-16  0.18354E-16 [c.c.]
 <---> loading info from shifty
 <---> [WARNING] Shifted grid in folder "shifty":change in number of WFs components
 <---> Grid 2 has shift vector        :   0.0000      0.10000E-01 -0.26767E-16 [c.c.]
 <---> loading info from shiftz
 <---> [WARNING] Shifted grid in folder "shiftz":change in number of WFs components
 <---> Grid 3 has shift vector        :  0.000000  0.000000  0.010000 [c.c.]
 <---> Dipoles shifted grids |                                        | [000%] --(E) --(X)

[ERROR] STOP signal received while in[04] Dipoles
[ERROR] Reading File shiftx/SAVE//ns.wf_fragments_1_1; Variable WF_COMPONENTS_@_SP_POL1_K1_BAND_GRP_1; NetCDF: Start+count exceeds dimension bound

So the error seems to point out that the components for shifted grid are larger than the unshifted one? If it is smaller, the problem would disappear?

Best,
Xiaoming

[Davide Sangalli]

Sorry, I did not reply anymore.
If you accept databases where the wave-functions have a different number of components, then you need handle the two different values.
In the subroutine

Code: Select all

  src/dipoles/DIPOLE_shifted_grids.F

1) define wf_ncx_shifted(3), change the set_parms subroutine to accept in input the iq value

Code: Select all

  call set_parms(iq,grid='shifted',command=' define')

and inside it

Code: Select all

  wf_ncx_shifted(iq) = wf_ncx

2) before the call to WF_load(WF,...) you need to set

Code: Select all

  wf_ncx=wf_ncx_save
  call to WF_load(WF,...)

and later, before the call to WF_load(WF_shifted,...)

Code: Select all

  wf_ncx=wf_ncx_shifted(iq)
  call to WF_load(WF_shifted,...)

After that, it should work.

Best,
D.