QE and Hybrid functional for BerkeleyGW

[sergio_pi...@berkeley.edu]

I am a graduate student at uc berkeley in the chemistry department, and I'm trying to use BerkeleyGW to do a GW calculation of MnO using hybrid functionals (PBE0). I'll try to explain my issue in terms of the documentation diagram online 

I see from the "BerkeleyGW Help" page that the it looks like using HSE seems entirely out of the questions (https://groups.google.com/a/berkeleygw.org/forum/#!searchin/help/hybrid|sort:date/help/8n3jWu_mL84/HKsHV4UICwAJ) but the post seems to imply that PBE0 should be possible.  From the same post (and from trying) it is clear that I can't run a calculation='bands'  in quantum espresso (QE) for PBE0. So in order to produce the kohn sham states I have to save them directly from an SCF calculation to produce WFN_inner. In the same step I can save the local part of the exchange correlation potential into vxc.dat and the density into RHO. 

The manual states "For Hartree-Fock/hybrid functional, no epsmat/eps0mat files are needed" which is consistent with the fact that I will not be able to run a "bands" calculation to produce WFNq with PBE0.

At this point I have tried various inputs for the sigma step (because I assume my only issue is a syntax error on my part) but always ends in one of the following errors, 
" ERROR: qgrid must be specified for Hartree-Fock."
"forrtl: severe (24): end-of-file during read, unit -5, file Internal List-Directed Read"

I should state that I'm using version berkeleygw/2.1.knl available on nersc and running  sigma.cplx.x because I have spin polarized KS-DFT on the QE side of things.  Do you know why I am having trouble running this calculation? Below is a version of my input that I think is most correct. but I've done variations... 1) make k points and q points the same, 2) change qgrid array to have no spacing between number 3) added a rightmost column in qpoints equivalent to rightmost column of kpoints etc. 

number_bands 27

band_index_min 1

band_index_max 27

spin_index_min 1

spin_index_max 2

frequency_dependence -1

screening_semiconductor

begin kpoints

  0.000000000  0.000000000  0.000000000   1.0 1

  0.000000000  0.000000000  0.500000000   1.0 0

  0.000000000  0.500000000  0.500000000   1.0 0

  0.500000000  0.500000000  0.500000000   1.0 0

end

bare_exchange_fraction 0.75

begin qpoints

  0.000000000  0.000000000  0.001000000   1.0 

  0.000000000  0.000000000  0.501000000   1.0 

  0.000000000  0.500000000  0.001000000   1.0 

  0.000000000  0.500000000  0.501000000   1.0 

  0.500000000  0.500000000  0.001000000   1.0 

  0.500000000  0.500000000  0.501000000   1.0 

end

qgrid 1 0 0 0 0 0

Thank you in advance for your help

[Mauro Del Ben]

Hi Sergio,

I'm not an expert on HF calculations with sigma, but I think I can help here. 

If I understood right, the qpoints list should be the list of q-points you would

have used if you had to perform the epsilon calculation, but without the small

q0 shift. I think that you can just try to use the same list as that you have for the

kpoints.

qgrid should be the size of the grid you are using, for example qgrid 3 3 3 for a 3x3x3 grid.

Hope this helps

Best

-M

--
You received this message because you are subscribed to the Google Groups "BerkeleyGW Help" group.
To unsubscribe from this group and stop receiving emails from it, send an email to help+uns...@berkeleygw.org.
To view this discussion on the web visit https://groups.google.com/a/berkeleygw.org/d/msgid/help/8d55970d-3a85-4cf1-af87-60d2d1954a15%40berkeleygw.org.

[sergio_pi...@berkeley.edu]

Thank you for the prompt reply.

taking what you said into consideration I get further along, but run into a new issue. using the following input ....

number_bands 27

band_index_min 1

band_index_max 27

spin_index_min 1

spin_index_max 2

no_symmetries_q_grid

frequency_dependence -1

screening_semiconductor

bare_exchange_fraction 0.75

begin kpoints

    0.000000    0.000000    0.000000   1.0 1 

   -0.000000   -0.000000    0.500000   1.0 0 

   -0.000000    0.500000   -0.000000   1.0 0 

   -0.000000    0.500000    0.500000   1.0 0 

    0.500000   -0.000000   -0.000000   1.0 0 

    0.500000   -0.000000    0.500000   1.0 0 

    0.500000    0.500000   -0.000000   1.0 0 

    0.500000    0.500000    0.500000   1.0 0 

end

begin qpoints

    0.000000    0.000000    0.000000   1.0 1 

   -0.000000   -0.000000    0.500000   1.0 0 

   -0.000000    0.500000   -0.000000   1.0 0 

   -0.000000    0.500000    0.500000   1.0 0 

    0.500000   -0.000000   -0.000000   1.0 0 

    0.500000   -0.000000    0.500000   1.0 0 

    0.500000    0.500000   -0.000000   1.0 0 

    0.500000    0.500000    0.500000   1.0 0 

end

qgrid 2 2 2

I reran the QE part to explicitly list the kpoints / qpoints using the K_points crystal card option. Now when I run the sigma part I get an error message saying ERROR: cannot find k-point in vxc.dat .
Using the wfn_rho_vxc_info.x tool I can see that these kpoints are in WFN file but the vxc.dat doesn't display kpoint info. WFN, RHO, VXC.dat all came from the same run of pw2bgw.x and I manually list out the kpoints for QE so I have no reason to believe that the k points in vxc.dat would be saved differently.

Maybe the error is due to how I'm running the QE calculation or my pp_in ?
for the QE side is there a way of saving the WFN for PBE0 other than at the end of calculation='scf' run? The examples for berkeleygw run calculation='bands' but this isn't an option for hybrid functionals.
for pp_in I am using    vxcg_flag = .true. and  vxcg_file = 'vxc.dat'   because I only want the local portion of the exchange correlation potential. The bare_exchange_fraction 0.75 option should take care of the nonlocal piece. Does this line of thought sound correct?

here is my pp_in

&input_pw2bgw

   prefix = 'MnO'

   real_or_complex = 2 

   wfng_flag = .true. 

   wfng_file = 'wfn'

   wfng_kgrid = .true.

   wfng_nk1 = 2 

   wfng_nk2 = 2 

   wfng_nk3 = 2 

   wfng_dk1 = 0.0 

   wfng_dk2 = 0.0 

   wfng_dk3 = 0.0 

   rhog_flag = .true.

   rhog_file = 'rho'

   vxcg_flag = .true.

   vxcg_file = 'vxc.dat'

   vxc_flag = .false.

   vxc_file = 'vxc.dat'

   vxc_diag_nmin = 1 

   vxc_diag_nmax = 27

   vxc_offdiag_nmin = 0 

   vxc_offdiag_nmax = 0 

/

and for completion my QE input

&control

   prefix = 'MnO'

   calculation = 'scf'

   restart_mode = 'from_scratch'

   wf_collect = .true.

   tstress = .false.

   tprnfor = .false.

   outdir = './'

   wfcdir = './'

   pseudo_dir = './'

/

&system

   ibrav = 0

   celldm(1) = 8.326099506257146

   nat    = 4

   nspin  = 2

   ntyp   = 3

   nbnd   = 28 

   ecutwfc = 300.0

   occupations               = 'fixed'

   starting_magnetization(1) = -0.5

   starting_magnetization(2) =  0.5

   tot_magnetization         = 0

   input_dft                 = 'pbe0'

nqx1=2

nqx2=2

nqx3=2

/

&electrons

   electron_maxstep = 500

   conv_thr = 1.0d-9

   mixing_mode = 'plain'

   mixing_beta = 0.3

   mixing_ndim = 8

   diagonalization = 'david'

   diago_david_ndim = 4

   diago_full_acc = .true.

/

CELL_PARAMETERS alat

    1.0        0.4945446  0.4945446

    0.4945446  1.0        0.4945446

    0.4945446  0.4945446  1.0

ATOMIC_SPECIES

   Mn1  54.900  Mn.RRKJ.upf

   Mn2  54.900  Mn.RRKJ.upf

     O  16.000   O.shin.ncpp

ATOMIC_POSITIONS crystal

   Mn1   0.0000   0.0000   0.0000

   Mn2   0.5000   0.5000   0.5000

     O   0.2500   0.2500   0.2500

     O   0.7500   0.7500   0.7500

K_POINTS crystal

  8

 0.0000000  0.0000000  0.0000000   0.1250000

 0.0000000  0.0000000  0.5000000   0.1250000

 0.0000000  0.5000000  0.0000000   0.1250000

 0.0000000  0.5000000  0.5000000   0.1250000

 0.5000000  0.0000000  0.0000000   0.1250000

 0.5000000  0.0000000  0.5000000   0.1250000

 0.5000000  0.5000000  0.0000000   0.1250000

 0.5000000  0.5000000  0.5000000   0.1250000

To unsubscribe from this group and stop receiving emails from it, send an email to he...@berkeleygw.org.

[Mauro Del Ben]

Hi Sergio,

You are using the same name for the VXC and vxc.dat

   vxcg_file = 'vxc.dat'

   vxc_flag = .false.

   vxc_file = 'vxc.dat'

change to:

   vxcg_file = 'VXC'

   vxc_flag = .true.

   vxc_file = 'vxc.dat'

Best

-M

To unsubscribe from this group and stop receiving emails from it, send an email to help+uns...@berkeleygw.org.
To view this discussion on the web visit https://groups.google.com/a/berkeleygw.org/d/msgid/help/12c8470b-e74d-40dc-b239-7e4a9199234c%40berkeleygw.org.

[sergio_pi...@berkeley.edu]

Thank you, that fixed it! I can run sigma now. 

One last question. I know I posted a link from 04/2019 that said range separated hybrid functionals such as HSE06 wouldn't be an option (https://groups.google.com/a/berkeleygw.org/forum/#!searchin/help/hybrid|sort:date/help/8n3jWu_mL84/HKsHV4UICwAJ ) but i found in the 2.1 manual which came out 06/2019 the following quote under the bare_exchange_fraction parameter 

Fraction of bare exchange. Set to 1.0 if you use the exchange-correlation matrix elements read from file vxc.dat. Set to 1.0 for local density functional, 0.0 for HF, 0.75 for PBE0, 0.80 for B3LYP if you use the local part of the exchange-correlation potential read from file VXC. For functionals such as HSE whose nonlocal part is not some fraction of bare exchange, use vxc.dat and not this option. This is set to 1.0 by default.

Does this mean I can use HSE from quantum espresso for GW if by simply using  the option    vxc_flag = .true.     vxc_file = 'vxc.dat'?

To view this discussion on the web visit https://groups.google.com/a/berkeleygw.org/d/msgid/help/12c8470b-e74d-40dc-b239-7e4a9199234c%40berkeleygw.org.

[Felipe H. da Jornada]

Hi Sergio,

You are right that, in principle, vxc.dat could contain all the required information for a calculation starting from an arbitrary hybrid functional, since it could in principle contain contributions from local and nonlocal potentials in the matrix elements. The problem is that, in practice, the wrapper pw2bgw currently does not compute the nonlocal contribution to the matrix elements in vxc.dat. My understanding is that other wrappers, such as JDFTx and Paratec, correctly compute these matrix elements including the non-local contribution due to the bare or modified exchange, but not pw2bgw.

The development team is also working hard to revamp pw2bgw and the way that such matrix elements are computed to easily support arbitrary exchange-correlation functionals. We don't have a date when this will be available, so you could look into these alternative codes as a workaround for now. I am asking the relevant developers to make sure these other wrappers correctly compute the nonlocal contribution into vxc.dat.

All the best,

Felipe

To unsubscribe from this group and stop receiving emails from it, send an email to help+uns...@berkeleygw.org.
To view this discussion on the web visit https://groups.google.com/a/berkeleygw.org/d/msgid/help/590329e0-9a43-44dc-b8fd-538c51fb01fc%40berkeleygw.org.

--

Felipe H. da Jornada

Assistant Professor

Department of Materials Science and Engineering

129 Durand Building

Stanford, CA 94305

Stanford University

https://jornada.stanford.edu/

[Andrew Wang]

Hi,

I'm trying to generate BGW input files from hybrid calculations as well and I had the same question; thanks for the conversation earlier - it has been really helpful in answering my question!

QE just released version 6.8 on Jul 19. From the update notes one of the changes is "Interface for BerkeleyGW extended to hybrid and meta-GGA functionals: Fangzhou Zhao (Berkeley)". Can I confirm if this means that the pw2bgw.x wrapper in this version has been updated to compute the non-local contributions to matrix elements in vxc.dat, such that it would be compatible with hybrid functionals (specifically HSE06)? 

Thanks!

Andrew

[Felipe Jornada]

Hi Andrew,

 

As far as I understand, the wrapper that comes with QE does handle spinors or hybrid functionals. You should copy the wrapper shipped with BerkeleyGW into the QE src/PP folder, and recompile QE. In the future, this wrapper will be shipped with QE out-of-the-box.

 

Best,

Felipe

 

–

Felipe H. da Jornada

Assistant Professor

Department of Materials Science and Engineering

 

129 Durand Building

Stanford, CA 94305

Stanford University

https://jornada.stanford.edu/

 

Pronouns: he/him/his

[Mauro Del Ben]

Hi Andrew, 

I think the wrapper shipped with QE 6.8 should work out-of-the-box for hybrid, spinors and meta-gga. 

I'm double checking with the other developers that everything is in place. For previous versions of QE, you 

can find the corresponding wrappers in the folders BerkeleyGW/MeanField/ESPRESSO/version-6.*

Best

-M

To view this discussion on the web visit https://groups.google.com/a/berkeleygw.org/d/msgid/help/6CCA2426-1DC6-4283-8002-E12C04CAC853%40stanford.edu.

[Fangzhou Zhao]

Hi Andrew,

Yes, the pw2bgw.x with QE 6.8 is compatible with hybrid functionals (including HSE06) and also metaGGA functionals. It does not compute the vxc.dat with non-local contributions of the hybrid functional directly, but it compute the rest part (energy eigenvalue - vxc) which does not have the non-local contributions. Please also refer to the README in /MeanField/ESPRESSO/version-6.4 or version 6.7. 

Best,

Fangzhou

[Andrew Wang]

Hi All,

Thanks so much for the responses. 

There seems to be a problem with the kih.dat file generated by pw2bgw.x in QE 6.8...

When I run the sigma code, with the generated kih file (file here), the program runs a bit but then crashes with 

"forrtl: severe (59): list-directed I/O syntax error, unit 876, file /gpfs/fs0/scratch/a/aspuru/awang95/berkeleygw/MnO/kih.dat", after 

 00:15:58   Dealing with k =  0.000000  0.000000  0.125000                1 / 28

================================================================================

 Reading kih.dat

After some searching I learned that it could be how the file was structured. The generated file has 4 columns in certain sections and 5 columns in others. For the parts with 4 columns, I added 0's after the 2nd column and made the file 5 columns ("kih_modified.dat") 

And then it worked; I was able to get after the "dealing with k" part

Number of q-points in the irreducible BZ(k) (nrq): 120

Started calculating Sigma with 960 block(s) at 00:16:31.

[ 00:16:32 |   0% ] block   1 / 960

So I guess maybe there be problems with the contents of the kih.dat file generated by pw2bgw that makes it currently not readable by sigma...?

Would you mind helping me with this?

Thanks!

Andrew

[Zhenglu Li]

Hi Andrew,

Thanks for sending the file and catching the problem!  Yes, this is due to Fortran using formatted output.  For example, the second line in your file reads:

      1       1-6642.841075585   -0.000000000

The second the third numbers need to have some spaces in between, but they are connected due to formatted output. That mean, "1-6642.841075585" should really be "1   -6642.841075585"

For now you may manually fix this issue.  We will work on a new update to fix this problem.  Thanks again for pointing this out!

Best,

Zhenglu

To view this discussion on the web visit https://groups.google.com/a/berkeleygw.org/d/msgid/help/0d248883-2287-450d-b582-43198b90977an%40berkeleygw.org.

[Zhenglu Li]

Hi Andrew,

Probably the best solution is, in specifying the bands of interest (those to be calculated by Sigma), you do not need to include bands from the beginning.  The first few bands seem to be deep core states that generally are not main research objects (which are usually near Fermi level).  If you set up the band indices skipping those bands, then the calculations should proceed normally.

Best,

Zhenglu

[Andrew Wang]

ok thanks so much!

[Andrew Wang]

sorry another question: when I use pw2bgw I haven't really been able to use the vxc_diag_nmax and vxc_offdiag_nmax flags...I seem to always get the message 

WARNING: resetting diag_nmax to max number of bands    15

WARNING: resetting offdiag_nmax to max number of bands    15

This message doesn't show up in the .out file, but rather shows up in the terminal as I run the program interactively. I've checked the kih.dat file and it seems that the bands do indeed only go up to 15 (see the link to kih files earlier)

I've tried to get around this limit on bands by generating the WFN file first and then using SAPO to generate another WFN with more bands (218 bands) and have so far been able to successfully complete the epsilon calculation, but I've been stuck at the sigma step. I can only set band_index_max = 15 even with the WFN from SAPO with 218 bands. When I set band_index_max = 218, even with number_diag 218 and the corresponding band indices under begin_diag (i.e. 1...218), I get the following error: 

read_matrix_elements: not enough diagonals present

I've tried but haven't been able to get around this...would you mind helping me with this? Sorry if any ignorance shows through on the related background physics; I'm a chemist by training but learning how to do GW/BSE calculations for my research including the theory but my progress in understanding has been quite slow...

I'm working with MnO with the HSE06 xc functional for the mean field input; here's a folder containing the relevant files leading up to the sigma calculation.

thanks!

 

[Fangzhou Zhao]

Hi Andrew,

This error message means that the number of bands in your DFT run with davidson diagonalization is still 15, so the pw2bgw code cannot calculate the kih.dat matrix elements from only 15 davidson diagonalized wavefunctions. The later error "read_matrix_elements: not enough diagonals present" is due to the previous error that you do not have kih.dat matrix elements up to 218 bands outputted. 

To my understanding, the SAPO code may generate wavefunctions with index higher than 15 with different format than wavefunctions generated by QE, so maybe only VXC is supported when you are doing SAPO to get more wavefunctions. I am not sure if SAPO will be compatible with the usage of vxc.dat or kih.dat which directly computes the matrix elements in pw2bgw. SAPO may be compatible with VXC, but you cannot use VXC for hybrid functional calculations. 

Best,

Fangzhou

--

Fangzhou Zhao

Ph.D. 

Department of Physics

University of California, Berkeley

[Andrew Wang]

sorry just an addition: the "read_matrix_elements: not enough diagonals present" goes away if I set band_index_nmin=1 and band_index_nmax=15 and the calculation seems to be able to proceed normally, but I'm just worried that this might restrict the number of bands involved to only 15, but we need a lot of bands for GW...?

thanks

[Andrew Wang]

Okay I see thanks a lot!

[Andrew Wang]

Hi All,

I've been trying to generate more bands for a WFN generated from a spin polarized HSE06 QE run of MnO (k-grid 8x8x8), which has a default of 15 bands (ecutwfc = 45, ecutrho = ecutfock = 180). With the default number of bands, I was able to successfully run epsilon and sigma. However, when I try to use Parabands to generate additional bands, I get the following error in epsilon:

ERROR: WFN ifmin/ifmax fields are inconsistent

I understand that this comes from the energies of conduction/valence bands not making sense. The new WFN after Parabands somehow has the conduction bands below the valence bands (VB minimum ~20 eV, CB minimum ~17.5 eV) while the original WFN directly from pw2bgw had reasonable values (I checked with wfn_rho_vxc_info.x and CB minimum ~ 12.9 eV,  VB minimum ~14.32 eV).

What could be the reason for this? Would it be something wrong with my scf calculation or my Parabands input? The relevant input files can be found here. I've done my best to figure things out on my own but am stuck...would really appreciate any help with this. 

Thanks a lot in advance,

Andrew