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

 Posts: 43
 Joined: Wed Jul 20, 2022 9:36 am
Dealing with Coulomb divergencies
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 kpoint 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 Voronoicells 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): 13921403.
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 kpoint 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 Voronoicells 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): 13921403.
Franz Fischer
PhD student / IMPRSUFAST fellow
Institute of Physical Chemistry
University of Hamburg
PhD student / IMPRSUFAST fellow
Institute of Physical Chemistry
University of Hamburg
 Daniele Varsano
 Posts: 3979
 Joined: Tue Mar 17, 2009 2:23 pm
 Contact:
Re: Dealing with Coulomb divergencies
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_GS_G)
where S_Gamma is a sphere contained in R_G. The integral over the sphere is analytic, and the integral over R_GS_G is done stochastically.
This operation is done in: src/coulomb/rim_integrate_v.F
Best,
Daniele
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_GS_G)
where S_Gamma is a sphere contained in R_G. The integral over the sphere is analytic, and the integral over R_GS_G is done stochastically.
This operation is done in: src/coulomb/rim_integrate_v.F
Best,
Daniele
Dr. Daniele Varsano
S3CNR Institute of Nanoscience and MaX Center, Italy
MaX  Materials design at the Exascale
http://www.nano.cnr.it
http://www.maxcentre.eu/
S3CNR Institute of Nanoscience and MaX Center, Italy
MaX  Materials design at the Exascale
http://www.nano.cnr.it
http://www.maxcentre.eu/

 Posts: 43
 Joined: Wed Jul 20, 2022 9:36 am
Re: Dealing with Coulomb divergencies
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 qpoints 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 kgrid 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) semianalytically? 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
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 qpoints 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 kgrid 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) semianalytically? 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 / IMPRSUFAST fellow
Institute of Physical Chemistry
University of Hamburg
PhD student / IMPRSUFAST fellow
Institute of Physical Chemistry
University of Hamburg

 Posts: 8
 Joined: Thu Jun 24, 2021 9:10 am
Re: Dealing with Coulomb divergencies
Dear Franz,
I wonder how exactly these regions R_G are constructed. Do you use Voronoicells 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 qpoints 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 kgrid 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) semianalytically? 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 Wav 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
I wonder how exactly these regions R_G are constructed. Do you use Voronoicells 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 qpoints 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 kgrid 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) semianalytically? 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 Wav 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

 Posts: 43
 Joined: Wed Jul 20, 2022 9:36 am
Re: Dealing with Coulomb divergencies
Dear Alberto,
first of all thanks a lot for your fast and thorough reply.
I have some followup questions though.
Is k_grid_uc_vol the Brillouin zone volume or the kspace volume of the small cells constructed around each qpoint?
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
first of all thanks a lot for your fast and thorough reply.
I have some followup questions though.
Is k_grid_uc_vol the Brillouin zone volume or the kspace volume of the small cells constructed around each qpoint?
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 / IMPRSUFAST fellow
Institute of Physical Chemistry
University of Hamburg
PhD student / IMPRSUFAST fellow
Institute of Physical Chemistry
University of Hamburg

 Posts: 8
 Joined: Thu Jun 24, 2021 9:10 am
Re: Dealing with Coulomb divergencies
Dear Franz,
I am happy to help.
In the following the answers:
Is k_grid_uc_vol the Brillouin zone volume or the kspace volume of the small cells constructed around each qpoint?
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 nonperiodic 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 kpoints = 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
I am happy to help.
In the following the answers:
Is k_grid_uc_vol the Brillouin zone volume or the kspace volume of the small cells constructed around each qpoint?
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 nonperiodic 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 kpoints = 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

 Posts: 43
 Joined: Wed Jul 20, 2022 9:36 am
Re: Dealing with Coulomb divergencies
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 qvectors.
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 qpoint where it has a hyperlinear slope.
Thus, the average value in each qinterval 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): 13921403.
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 qvectors.
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 qpoint where it has a hyperlinear slope.
Thus, the average value in each qinterval 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): 13921403.
You do not have the required permissions to view the files attached to this post.
Franz Fischer
PhD student / IMPRSUFAST fellow
Institute of Physical Chemistry
University of Hamburg
PhD student / IMPRSUFAST fellow
Institute of Physical Chemistry
University of Hamburg
 Daniele Varsano
 Posts: 3979
 Joined: Tue Mar 17, 2009 2:23 pm
 Contact:
Re: Dealing with Coulomb divergencies
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 nonintegrated potential (noRIM) 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 miniBZ 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 nonRIM 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
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 nonintegrated potential (noRIM) 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 miniBZ 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 nonRIM 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
S3CNR Institute of Nanoscience and MaX Center, Italy
MaX  Materials design at the Exascale
http://www.nano.cnr.it
http://www.maxcentre.eu/
S3CNR Institute of Nanoscience and MaX Center, Italy
MaX  Materials design at the Exascale
http://www.nano.cnr.it
http://www.maxcentre.eu/

 Posts: 100
 Joined: Fri Mar 26, 2021 11:27 am
Re: Dealing with Coulomb divergencies
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 postprocessing 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
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 postprocessing 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
 Daniele Varsano
 Posts: 3979
 Joined: Tue Mar 17, 2009 2:23 pm
 Contact:
Re: Dealing with Coulomb divergencies
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 yambopy.
Best,
Daniele
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 yambopy.
Best,
Daniele
Dr. Daniele Varsano
S3CNR Institute of Nanoscience and MaX Center, Italy
MaX  Materials design at the Exascale
http://www.nano.cnr.it
http://www.maxcentre.eu/
S3CNR Institute of Nanoscience and MaX Center, Italy
MaX  Materials design at the Exascale
http://www.nano.cnr.it
http://www.maxcentre.eu/