Fixing TB dipoles, cpot and development of polarizability.

[rreho]

This issue is meant to be a thread for discussion about the calculation of the Coulomb potential within the analytical model approach. This problem concerns also the class TB_dipoles. @stevenbos123 recently worked on the issue via a series of commits. The main one is this one 574ca03 . I agree with the changes done in get_eps in wann_H2P.py. Regarding the class TB_dipoles, I find the vectorization of the arrays useful to improve speed-up and reduce the time spent on the double for loop over transition indices.

• I noticed that we use self.nb = nc+nv and self.bse_nb = bse_nv + bse_nc. I am ok with it, but I remind that the number of bands might be different the the number of wannier projections and it could be a potential bug.

• I notice that we added a new flag with_bse. I am ok with it as it makes the code clear. However, it's an overkill as the presence or not of h2peigvec is enough to determine wether we want to call get_dipoles or get_dipoles_bse. I understand that additional control is useful.

• I checked these lines:

          eigvec_v = self.eigvec[k_indices, :, v_indices]
  
          # Compute dipoles
          for dim in range(dim_hlm):
              # Compute the dot product
              dot_product = np.einsum('ij,ijk->ik', eigvec_v, self.hlm[k_indices, :, :, dim])
  
              # Compute the vdot and multiply with GR
              dipoles[:, dim] = GR * np.sum(np.conjugate(eigvec_c) * dot_product, axis=1)

  And realized that those work when the table comes from a wann_model.Table. However, in wann_H2P, the user can specify bse_nv and bse_nc which change the indexing of valence band. This should be a new issue and we should discuss a general solution about the transitions table. I went through the np.einsum logic and the way you build the auxiliary arrays. Seems ok to me. However, I am doubting the division by Ha2eV**3. The reason should be that there are two. The dipoles should have dimension of Bohr. The numerator should be in units of eV**2*Bohr while the denominator in eV**2 which cancel out.

• The logic of nq_double is not the best.

• Why are we not vectorizing get_dipoles_bse?

Lastly. I tested the absorption plot. The shape of real and imaginary part of eps are reasonable, while the magnited (1e-12) is totally wrong. I suspect units. TODO

• Delete the conversion HA2EV**3 in get_dipoles
• vectorize the other routines
• Fix units eps

[rreho]

In devel-tbwannier-modelcpot I re-implemented from scratch the class TB_dipoles. This is the relevant commit cd0b6f6 . I tested it and the code runs with a model Coulomb potentials. Note that I want to keep that k and q grid. It looks a bit useless, but from our point of view this code should compute the excitonic overlap matrix and therefore I need to have these two independent grids.

TODO @stevenbos123

• Do a general testing
• Improve the routine get_eps and class TB_dipoles with broadcasting and, perhaps, chunking
• Test on WS2 monolayer with a small grid. Benchmark your results with WantiBexos.

[rreho]

@stevenbos123 spotted a bug where all the F_kcv are the same. Action needed.

[stevenbos123]

related to the bug where all the 'F_kcv' are the same I pushed this change: #622d8b0 . 'F_kcv' should be oscillator strength of k,c,v for each transition, n. Instead it was oscillator strength of k,c,v for the transition that accompanies this k,c,v. Hence the values where all the same for each element of 'F_kcv'.

[rreho]

The bug spotted above is real. The fix works. The code runs up to get_eps. The values of eps are not "small" anymore and the shape looks reasonable.

TODO

• Benchmark against known results.
• Improve the speed as build_H2P_fromcpot is quite slow.

[stevenbos123]

Benchmark of WS2 6x6x1 kgrid. Comparing our get_eps() implementation to calculate the IPA absorption, with the values acquired from wantibexos. The coulomb potential used is the V2DT2 model potential, with values: v0 = 0.00 lc = 8.0 w = 0.0 r0=1.0 eta=0.08

The plot is normalized to the max of the respective values. From wantibexos the imaginary part of the sp_diel_xx.dat file is plotted which is the single particle dielectric function \epsilon.

As can be seen in the plot, the second peak position coincides with the result from wantibexos, but the first peak position is not.

To investigate this further I plotted the F_kcv both from our implementation and from wantibexos sp_optics.dat again renormalized:

[stevenbos123]

When comparing our results from the dielectric function, with those with wantibexos, I noticed one fundamental difference between the two implementations. I will explain the difference by using the example of setting nc=1 and nv=1 in our yambopy code and doing the same for Wantibexos by specifying NBANDSC= 1 NBANDSV= 1. The way wantibexos uses this input variables is that the first band up or down from the fermi energy is considered to be the valence or conduction band used for the calculation. Now the tricky part is, this is not the same band for each kpoint, so the index of the conduction band might chance for the kpoints. This of course physically makes sense as you are interested in the minimal band gap for each kpoint. Now this is where our implementation differs. In our implementation the band index does not change for the k-points. In other words the bands are not sorted for each k-point on energy. So we do not consider the minimal gap possible for each k-point. This difference of course becomes less of a problem when you consider more bands.

Now a question to ask ourselves, should we keep our implementation the way it is, or copy the wantibexos way? In this case the yambo transition table now does not make sense anymore.

[stevenbos123]

When comparing our results from the dielectric function, with those with wantibexos, I noticed one fundamental difference between the two implementations. I will explain the difference by using the example of setting nc=1 and nv=1 in our yambopy code and doing the same for Wantibexos by specifying NBANDSC= 1 NBANDSV= 1. The way wantibexos uses this input variables is that the first band up or down from the fermi energy is considered to be the valence or conduction band used for the calculation. Now the tricky part is, this is not the same band for each kpoint, so the index of the conduction band might chance for the kpoints. This of course physically makes sense as you are interested in the minimal band gap for each kpoint. Now this is where our implementation differs. In our implementation the band index does not change for the k-points. In other words the bands are not sorted for each k-point on energy. So we do not consider the minimal gap possible for each k-point. This difference of course becomes less of a problem when you consider more bands.

Now a question to ask ourselves, should we keep our implementation the way it is, or copy the wantibexos way? In this case the yambo transition table now does not make sense anymore.

I looked at the eigenvalues from h2p and model, and they are both sorted per k-point. This implies we compute the same thing as wantibexos does, just our band indexing does not change with it.