Dealing with Coulomb divergencies

Various technical topics such as parallelism and efficiency, netCDF problems, the Yambo code structure itself, are posted here.

Moderators: Davide Sangalli, andrea.ferretti, myrta gruning, andrea marini, Daniele Varsano, Conor Hogan, Nicola Spallanzani

Franz Fischer
Posts: 38
Joined: Wed Jul 20, 2022 9:36 am

Dealing with Coulomb divergencies

Post by Franz Fischer » Tue Feb 07, 2023 4:39 pm

Dear Developers,

in both of your papers you present the random integration method as a way to remove the divergencies arising from the Coulomb potential (at G = 0) [1].
You subdivide the Brillouin zone in small regions R_G centered at each k-point where the unification of all R_G's is the Brillouin zone volume and then you kind of average the Coulomb interaction in these regions to get a suitable value at G = 0.
I wonder how exactly these regions R_G are constructed. Do you use Voronoi-cells or rectangles?

Maybe you could even guide me to the corresponding file in your source code, so I can have a look by myself.
Thanks a lot in advance!


[1] Marini, Andrea, et al. "Yambo: an ab initio tool for excited state calculations." Computer Physics Communications 180.8 (2009): 1392-1403.
Franz Fischer
PhD student / IMPRS-UFAST fellow
Institute of Physical Chemistry
University of Hamburg

User avatar
Daniele Varsano
Posts: 3773
Joined: Tue Mar 17, 2009 2:23 pm
Contact:

Re: Dealing with Coulomb divergencies

Post by Daniele Varsano » Wed Feb 08, 2023 9:29 am

Dear Franz,


The integral is done via a Monte Carlo stochastic method (eq. 30 of the paper you mention). and the cell shape around gamma R_Gamma depends on the symmetry of the BZ and it is defined by your sampling.
The integral is done for all Rq the BZ is divided (BZ=URq) and in particular for R_gamma where the potential diverges R_gamma is divided as:
R_G=S_gamma U (R_G-S_G)
where S_Gamma is a sphere contained in R_G. The integral over the sphere is analytic, and the integral over R_G-S_G is done stochastically.
This operation is done in: src/coulomb/rim_integrate_v.F
Best,

Daniele
Dr. Daniele Varsano
S3-CNR Institute of Nanoscience and MaX Center, Italy
MaX - Materials design at the Exascale
http://www.nano.cnr.it
http://www.max-centre.eu/

Franz Fischer
Posts: 38
Joined: Wed Jul 20, 2022 9:36 am

Re: Dealing with Coulomb divergencies

Post by Franz Fischer » Wed Feb 08, 2023 3:47 pm

Dear Daniele,

thanks for your fast reply.
I read the source code (src/coulomb/rim_integrate_v.F) you send me and I have some further questions.

1) Is N_out the number of random q-points used in the RIM (i.e. the same as the variable RIM_n_rand_pts)?
2) Why is there a factor of 8 in the rfac defined in line no. 43? Is this because you only integrate along the positive part of the axes of the k-grid volume?
3) What exactly is q0_def_norm (line no. 146)?
4) In line no. 161 are you computing the V(q=0, G .neq. 0) semi-analytically? It seems that this is a different expression than the previous one which was for V(q=0, G=0) and V(q .neq. 0, G .neq. 0).

I am only interested in the part that is related to 2D systems (i.e. if cut_is_slab is true).

Thanks in advance.

Best,
Franz
Franz Fischer
PhD student / IMPRS-UFAST fellow
Institute of Physical Chemistry
University of Hamburg

Alberto Guandalini
Posts: 8
Joined: Thu Jun 24, 2021 9:10 am

Re: Dealing with Coulomb divergencies

Post by Alberto Guandalini » Thu Feb 09, 2023 4:26 pm

Dear Franz,

|I wonder how exactly these regions R_G are constructed. Do you use Voronoi-cells or rectangles?

I do not remember exactly what Voroni cells are, but we do the following steps:
1) we select a set of random points in reciprocal lattice units around the point q = 0 (from -1 to +1 for all the three coordinates);
2) we measure the distance between the random point and the q=0 and its nearest neighbour points in iku (internal crystal coordinates);
3) if the minimum distance is the one with q = 0, we take the random point, otherwise we discard it;

In this way, we sample the Brillouin zone of the Monkhorst pack grid of the q grid, called mini BZ in the paper (red area in Fig. 2).

N.B. for 2D samplings, we use the same procedure as in 3D (for ease of implementation), but in the end we drop the value of
each coordinate of the random points to zero along the non parallel direction. In this way we manage to sample only the 2D BZ.

In the following, the answers of your other questions:
1) Is N_out the number of random q-points used in the RIM (i.e. the same as the variable RIM_n_rand_pts)?
N_out is the number of random points BEFORE the selection over the distance (point 3).
In the end, it is 8 times the number of random points selected.

2) Why is there a factor of 8 in the rfac defined in line no. 43? Is this because you only integrate along the positive part of the axes of the k-grid volume?
That is because N_out is 8 times the number of selected random points. See the answer above.

3) What exactly is q0_def_norm (line no. 146)?
q0_def_norm = 10^-5. That is the regularization of the q=0 term, which cannot be exactly zero due to numerical issues.

4) In line no. 161 are you computing the V(q=0, G .neq. 0) semi-analytically? It seems that this is a different expression than the previous one which was for V(q=0, G=0) and V(q .neq. 0, G .neq. 0).
Please, the next time copy and paste directly the piece of code. In this way it will be easier to understand.

If you refer to:

| do i2=2,i2max,1
| r1=iku_v_norm(qr(:,i1))*iku_v_norm(g_vec(i2,:)+qr(:,i1))
| RIM_acc(i2)=RIM_acc(i2)+rfac/r1
| enddo

that part of the code is deprecated since the W-av method circumvent the calculation of the wings of the symmetrized Coulomb potential v(q,G,G') = sqrt(v(q,G))*sqrt(v(q,G')).
You may see in a previous part of the code:

| if (RIM_is_diagonal) i2max=1

By default, RIM_is_diagonal == True. Than, the for cycle does nothing.

Best,
Alberto

Franz Fischer
Posts: 38
Joined: Wed Jul 20, 2022 9:36 am

Re: Dealing with Coulomb divergencies

Post by Franz Fischer » Mon Feb 13, 2023 4:27 pm

Dear Alberto,

first of all thanks a lot for your fast and thorough reply.
I have some follow-up questions though.

Is k_grid_uc_vol the Brillouin zone volume or the k-space volume of the small cells constructed around each q-point?
How does k_grid_uc_vol relate to the Recip. lattice volume you print in your r_setup files?

Where does this volume factor actually come from and why is there a factor of (2 pi)^-3 in rfac? It is defined as rfac=8._DP*real(k_grid_uc_vol,DP)/real(N_out,DP)/(2._DP*pi_DP)**3.

And finally, where does the factor of 2 come from in this line: RIM_acc(1)=RIM_acc(1)+2._DP*real(pre_factor,DP)*rfac/r1

Thanks a lot for your patience.

Best,
Franz
Franz Fischer
PhD student / IMPRS-UFAST fellow
Institute of Physical Chemistry
University of Hamburg

Alberto Guandalini
Posts: 8
Joined: Thu Jun 24, 2021 9:10 am

Re: Dealing with Coulomb divergencies

Post by Alberto Guandalini » Tue Feb 14, 2023 9:23 am

Dear Franz,
I am happy to help.

In the following the answers:

|Is k_grid_uc_vol the Brillouin zone volume or the k-space volume of the small cells constructed around each q-point?
k_grid_uc_vol is the volume associated to each k point of the BZ, thus the volume in which the average is done.
For 2D systems, the average is done in a correspondent area (volume divided by the length of the BZ in the non-periodic direction).

|How does k_grid_uc_vol relate to the Recip. lattice volume you print in your r_setup files?
You can check in your calculation that k_grid_uc_vol * number of k-points = volume of the BZ.

|Where does this volume factor actually come from and why is there a factor of (2 pi)^-3 in rfac? It is defined as rfac=8._DP*real(k_grid_uc_vol,DP)/real(N_out,DP)/(2._DP*pi_DP)**3.
When you do an average, you must divide by the volume of integration.
Regarding the (2pi)^{-3}, I do not remember exactly, but it think there should be a (2 pi)^-3 due to the integration in Fourier space.
In any case, the RIM is used for example in QP_ppa_coshex.F. I advise you to check also the factors in that soubroutine, e.g. in a plasmon pole GW calculations.
Maybe some coefficients are added so to be simplified in a subsequent step for the sake of compatibility with the version of the code without the RIM.

|And finally, where does the factor of 2 come from in this line: RIM_acc(1)=RIM_acc(1)+2._DP*real(pre_factor,DP)*rfac/r1
In principle there should be a 4pi factor (vcol = 4pi/q^2), but the 2pi has not been included as it is simplified with another 2pi later in the calculation (mainly due to historical reasons).

Cheers,
Alberto

Franz Fischer
Posts: 38
Joined: Wed Jul 20, 2022 9:36 am

Re: Dealing with Coulomb divergencies

Post by Franz Fischer » Wed Feb 15, 2023 2:03 pm

Dear Alberto,

thanks again for your help.

I implemented the RIM in our code and plotted the Coulomb potential with and without the RIM over the norm of the q-vectors.
Comparing the results to Fig. 5 (right) in this paper [1] leaves me wondering.
For me it seems that the RIM Coulomb potential is always equal or a little bit larger than the bare Coulomb potential (see attached figure).
This seems very intuitive, as you are averaging the potential around each q-point where it has a hyper-linear slope.
Thus, the average value in each q-interval around q will be larger than the value at q itself.

Is this correct? Are the discrepancies coming from the different dimensionalities of the systems (1D molecule in [1] and 2D MoS2 slab for me)?

Thanks a lot in advance.

Best,
Franz

[1] Marini, Andrea, et al. "Yambo: an ab initio tool for excited state calculations." Computer Physics Communications 180.8 (2009): 1392-1403.
You do not have the required permissions to view the files attached to this post.
Franz Fischer
PhD student / IMPRS-UFAST fellow
Institute of Physical Chemistry
University of Hamburg

User avatar
Daniele Varsano
Posts: 3773
Joined: Tue Mar 17, 2009 2:23 pm
Contact:

Re: Dealing with Coulomb divergencies

Post by Daniele Varsano » Sun Feb 19, 2023 4:16 pm

Dear Franz,

the integral in the reference refers to a 3D coulomb potential 1/q^2 with a 1D sampling. The BZ shape if I remember well is orthorombic.
The non-integrated potential (no-RIM) is divergent in q=0 as it is calculated as \int d^3q 1/|q|^2 ~ \sum_q_i 1/|q_i|^2 Vol_qi while the one calculated by Montecarlo
do converge in 3 dimension, so it seems reasonable to me that it has a lower value for small q. As you say, the dimension, and the shape of the mini-BZ can be the reason of the discrepancy you get. The one shown in the figure is a 3D elongated BZ with 1D k point sampling, the non-RIM method is particularly bad in this case as you are sampling the 1/q^2 function only in one dimension so highly overestimating the integral for small q.

Best,
Daniele
Dr. Daniele Varsano
S3-CNR Institute of Nanoscience and MaX Center, Italy
MaX - Materials design at the Exascale
http://www.nano.cnr.it
http://www.max-centre.eu/

Quxiao
Posts: 96
Joined: Fri Mar 26, 2021 11:27 am

Re: Dealing with Coulomb divergencies

Post by Quxiao » Mon Dec 04, 2023 8:29 am

Dear Daniele:
When I finish calculating BSE, how can I obtain the Coulomb divergence integral at q=0, which is W (q=0)? What code do I need to modify in yambo? Can post-processing like ypp achieve similar things? If so, I also want to obtain the W matrix to further analyze and observe the physical processes involved.
I am very grateful for your help,Thanks!
Best wishes!
Quxiao
Quxiao in BIT,calculate the exciton

User avatar
Daniele Varsano
Posts: 3773
Joined: Tue Mar 17, 2009 2:23 pm
Contact:

Re: Dealing with Coulomb divergencies

Post by Daniele Varsano » Mon Dec 04, 2023 12:33 pm

Dear Quxiao,

it is not very clear what quantity you want to look at. Can you be more explicit?
If you want to look at the W(q=0), this is calculated for a small value of q (q=10^-5). In the ndb.em1s (or ndb.pp) is contained in the
vX elements and form that you can build the W. You can read them by using yambo-py.

Best,
Daniele
Dr. Daniele Varsano
S3-CNR Institute of Nanoscience and MaX Center, Italy
MaX - Materials design at the Exascale
http://www.nano.cnr.it
http://www.max-centre.eu/

Post Reply