problems regarding absorption calculation

Submitted by zwang9505 on Fri, 04/28/2017 - 07:02

Forums 

User questions

Hi all,

I’m now doing a GW calculation of system similar to MoS2 with a 2D hexagonal crystal structure. I’ve finished the epsilon, sigma and kernel calculation smoothly. However, problems occur during absorption calculation. I used an unreduced k point grid 6 * 6 * 1 for WFN_fi plus 0.001 shift in kx for WFNq_fi.

The error message is :
ERROR: q-shift vector may not be zero.

There is another warning message before that reads:
WARNING: checkbz: unfolded BZ from WFN_fi has missing k-points.

The run stops at here in absorption.out:
Highest occupied band (shifted fine grid) = 19
Shift vector : 0.00000 0.00000 0.00000 Length = 0.00000

My k points grid in Quantum Espresso for WFN_fi is
K_POINTS crystal
36
0.0000000 0.0000000 0.0000000 1.0
0.1666667 0.0000000 0.0000000 1.0
0.3333333 0.0000000 0.0000000 1.0
0.5000000 0.0000000 0.0000000 1.0
0.6666667 0.0000000 0.0000000 1.0
0.8333333 0.0000000 0.0000000 1.0
0.0000000 0.1666667 0.0000000 1.0
0.1666667 0.1666667 0.0000000 1.0
0.3333333 0.1666667 0.0000000 1.0
0.5000000 0.1666667 0.0000000 1.0
0.6666667 0.1666667 0.0000000 1.0
0.8333333 0.1666667 0.0000000 1.0
0.0000000 0.3333333 0.0000000 1.0
0.1666667 0.3333333 0.0000000 1.0
0.3333333 0.3333333 0.0000000 1.0
0.5000000 0.3333333 0.0000000 1.0
0.6666667 0.3333333 0.0000000 1.0
0.8333333 0.3333333 0.0000000 1.0
0.0000000 0.5000000 0.0000000 1.0
0.1666667 0.5000000 0.0000000 1.0
0.3333333 0.5000000 0.0000000 1.0
0.5000000 0.5000000 0.0000000 1.0
0.6666667 0.5000000 0.0000000 1.0
0.8333333 0.5000000 0.0000000 1.0
0.0000000 0.6666667 0.0000000 1.0
0.1666667 0.6666667 0.0000000 1.0
0.3333333 0.6666667 0.0000000 1.0
0.5000000 0.6666667 0.0000000 1.0
0.6666667 0.6666667 0.0000000 1.0
0.8333333 0.6666667 0.0000000 1.0
0.0000000 0.8333333 0.0000000 1.0
0.1666667 0.8333333 0.0000000 1.0
0.3333333 0.8333333 0.0000000 1.0
0.5000000 0.8333333 0.0000000 1.0
0.6666667 0.8333333 0.0000000 1.0
0.8333333 0.8333333 0.0000000 1.0

The corresponding pw2bgw input file reads:
&input_pw2bgw
prefix = 'InS'
outdir = './tmp'
real_or_complex = 2
wfng_flag = .true.
symm_type = 'hexagonal'
wfng_file = 'wfn.cplx'
wfng_kgrid = .true.
wfng_nk1 = 6
wfng_nk2 = 6
wfng_nk3 = 1
wfng_dk1 = 0
wfng_dk2 = 0
wfng_dk3 = 0
/

My k points grid in Quantum Espresso for WFN_fi is
K_POINTS crystal
36
0.0010000 0.0000000 0.0000000 1.0
0.1676667 0.0000000 0.0000000 1.0
0.3343333 0.0000000 0.0000000 1.0
0.5010000 0.0000000 0.0000000 1.0
0.6676667 0.0000000 0.0000000 1.0
0.8343333 0.0000000 0.0000000 1.0
0.0010000 0.1666667 0.0000000 1.0
0.1676667 0.1666667 0.0000000 1.0
0.3343333 0.1666667 0.0000000 1.0
0.5010000 0.1666667 0.0000000 1.0
0.6676667 0.1666667 0.0000000 1.0
0.8343333 0.1666667 0.0000000 1.0
0.0010000 0.3333333 0.0000000 1.0
0.1676667 0.3333333 0.0000000 1.0
0.3343333 0.3333333 0.0000000 1.0
0.5010000 0.3333333 0.0000000 1.0
0.6676667 0.3333333 0.0000000 1.0
0.8343333 0.3333333 0.0000000 1.0
0.0010000 0.5000000 0.0000000 1.0
0.1676667 0.5000000 0.0000000 1.0
0.3343333 0.5000000 0.0000000 1.0
0.5010000 0.5000000 0.0000000 1.0
0.6676667 0.5000000 0.0000000 1.0
0.8343333 0.5000000 0.0000000 1.0
0.0010000 0.6666667 0.0000000 1.0
0.1676667 0.6666667 0.0000000 1.0
0.3343333 0.6666667 0.0000000 1.0
0.5010000 0.6666667 0.0000000 1.0
0.6676667 0.6666667 0.0000000 1.0
0.8343333 0.6666667 0.0000000 1.0
0.0010000 0.8333333 0.0000000 1.0
0.1676667 0.8333333 0.0000000 1.0
0.3343333 0.8333333 0.0000000 1.0
0.5010000 0.8333333 0.0000000 1.0
0.6676667 0.8333333 0.0000000 1.0
0.8343333 0.8333333 0.0000000 1.0

The corresponding pw2bgw input file reads:
&input_pw2bgw
prefix = 'InS'
outdir = './tmp'
real_or_complex = 2
wfng_flag = .true.
symm_type = 'hexagonal'
wfng_file = 'wfn.cplx'
wfng_kgrid = .true.
wfng_nk1 = 6
wfng_nk2 = 6
wfng_nk3 = 1
wfng_dk1 = 0.006
wfng_dk2 = 0
wfng_dk3 = 0
/

The absorption.inp file is:
number_val_bands_fine 4
number_val_bands_coarse 6
number_cond_bands_fine 4
number_cond_bands_coarse 6
use_velocity
energy_resolution 0.15
use_symmetries_coarse_grid
no_symmetries_fine_grid
no_symmetries_shifted_grid
screening_semiconductor
cell_slab_truncation
eqp_co_corrections

Any help would be appreciated!

Thanks,
Zhan

Submitted by babarker on Fri, 04/28/2017 - 13:03

Hello Zhan,

I would recommend using a random shift to generate the k-points for both WFN_fi and WFNq_fi. For example, in the "examples" directory, we see for Si, from the input file for pw2bgw,

WFN_fi:
wfng_nk1 = 10
wfng_nk2 = 10
wfng_nk3 = 10
wfng_dk1 = 0.47
wfng_dk2 = 0.37
wfng_dk3 = 0.31

WFNq_fi:
wfng_nk1 = 10
wfng_nk2 = 10
wfng_nk3 = 10
wfng_dk1 = 0.47
wfng_dk2 = 0.37
wfng_dk3 = 0.32

This way, both wavefunctions have grids that are not able to be reduced by symmetries, and they are slightly shifted (in k3) from each other. For a 2D system, the dk3 shift will be 0.0 for both files, as you did. If the system is symmetric in k1 and k2, pick either one for the shift. If the system is not symmetric, it makes the most sense to do two different absorption calculations for the different polarization choices.

Also, in the absorption.inp file, make sure to use

use_symmetries_coarse_grid
no_symmetries_fine_grid
no_symmetries_shifted_grid

And I take it that your choice of 6x6x6 fine grid sampling is strictly for testing purposes? 2D materials require very, very fine BZ sampling to converge, due to the structure of the dielectric matrix. (For instance, see Qiu, Jornada, and Louie, PRB 93, p. 235435).

Best,
Brad

Submitted by zwang9505 on Sat, 04/29/2017 - 06:26

Hi Brad,

Thanks! I used a random shift to generate the k grid like you said, and it worked out well.

And it is true that I used the 6*6*1 fine grid for testing purposes. However, I do have another question regarding this issue.

I understand that the convergence in absorption calculation requires a rather fine k grid, in fact, 72*72*1 seems to be good in MoS2 systems according to Qiu, Felipe and Louie PRL 111 (21), 216805. And for epsilon and sigma calculation, the convergence requires a k grid as dense as 12*12*1 plus ecutoff 500eV (about 40 Ry) and number of bands 1000 or more. However, according to my calculation, the epsilon and sigma calculation requires large available memory and calculation time. Would you please tell me if there is any way to reduce the memory cost or speed up the calculation? Thanks.

Best,
Zhan

Submitted by dyq on Sat, 04/29/2017 - 21:09

Hi Zhan,

A couple comments:

1. k-grid convergence for 2D systems is very tricky. There is an erratum to Qiu, Felipe and Louie, PRL (2013), which indicates that you need about 18x18x1 k-points in epsilon and sigma and for the coarse grid in the BSE steps to converge the gap to about 100 meV. Unfortunately, the standard interpolation scheme used in the BerkeleyGW for the BSE step does not work too well for 2D systems. It might work for you depending on the material you're studying, but you should carefully converge the BSE with respect to both the coarse and fine grids. See Qiu, Jornada, Louie, PRB (2016) for details.

2. Some easy options for saving memory and calculation time:
- Using a smaller vacuum (~15 Angstroms) will allow you to use fewer empty states. See Qiu, Jornada, Louie, PRB (2016) for details.
- The BSE typically converges much more quickly with respect to screened cutoff and bands included in the dielectric matrix than the sigma step. For MoS2, we find that ~5 Ry and ~300 bands works fine for BSE.
- Try using a scissor shift for the GW bandstructure, so you don't have to calculate every k point in the Sigma step.

3. Some options for saving memory and calculation which may require coding on your part:
- We have developed some new interpolation and subsampling schemes specifically for reducing computational cost in 2D and 1D systems. Unfortunately, they will not be available in the public BerkeleyGW code until after the next version is released. However, depending on your familiarity with coding you can try implementing something similar yourself. See Jornada, Qiu, and Louie, PRB (2017) for details.
- To reduce the cost of summing over many conduction bands, you can replace the higher energy conduction bands with a smaller number of effective conduction bands. See this paper for details: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5104980/ .

4. In your original question, the error message "ERROR: q-shift vector may not be zero" suggests that the q-shift of 0.001 may be reported incorrectly in the header of your WFNq_fi file. To check the header, use the command wfn_rho_wxc_info.x. You can also double-check that the line "q_shift 0.001 0.0 0.0" is reported correctly in your absorption.inp file.

Hope that helps!

Best,
Diana

Submitted by zwang9505 on Sat, 04/29/2017 - 22:21

Hi Diana,

Thanks very much for your detailed comments! I will try some of these approaches in the following calculations.

Best,
Zhan