use IBZ kpoint list for writing GAMMA vaspgamma.h5

* use IBZ kpoint list for writing GAMMA vaspgamma.h5
* make vaspgamma.h5 default
* write n_k_ibz into dft_misc_input so in can be used inside of sumk
* write n_k_ibz into seed.ctrl json file
* change ref h5 files for VASP tests by adding n_k_ibz

Author: the-hampel

python/triqs_dft_tools/converters/plovasp/elstruct.py

@@ -39,6 +39,7 @@ class ElectronicStructure:
 
     - *natom* (int) : total number of atoms
     - *nktot* (int) : total number of `k`-points
+    - *nkibz* (int) : number of `k`-points in IBZ
     - *nband* (int) : total number of bands
     - *nspin* (int) : spin-polarization
     - *nc_flag* (True/False) : non-collinearity flag
@@ -58,8 +59,9 @@ def __init__(self, vasp_data):
         self.natom = vasp_data.poscar.nq
         self.type_of_ion = vasp_data.poscar.type_of_ion
         self.nktot = vasp_data.kpoints.nktot
+        self.nkibz = vasp_data.kpoints.nkibz
 
-        self.kmesh = {'nktot': self.nktot}
+        self.kmesh = {'nktot': self.nktot, 'nkibz': self.nkibz}
         self.kmesh['kpoints'] = vasp_data.kpoints.kpts
         # VASP.6.
         self.nc_flag = vasp_data.plocar.nc_flag
@@ -86,7 +88,7 @@ def __init__(self, vasp_data):
 
 # Check that the number of k-points is the same in all files
         _, ns_plo, nk_plo, nb_plo = vasp_data.plocar.plo.shape
-        assert nk_plo == self.nktot, "PLOCAR is inconsistent with IBZKPT (number of k-points)"
+        assert nk_plo == self.nktot, "PLOCAR is inconsistent with IBZKPT (number of k-points). If you run VASP with symmetry make sure to use the h5 interface of the converter, i.e. have the locproj information written to vaspout.h5"
 
 # FIXME: Reading from EIGENVAL is obsolete and should be
 #        removed completely.

python/triqs_dft_tools/converters/plovasp/plotools.py

@@ -289,6 +289,8 @@ def ctrl_output(conf_pars, el_struct, ng):
 
      * *nk*: number of `k`-points
 
+     * *nkibz*: number of `k`-points in IBZ
+
      * *ns*: number of spin channels
 
      * *nc_flag*: collinear/noncollinear case (False/True)
@@ -309,6 +311,7 @@ def ctrl_output(conf_pars, el_struct, ng):
 # Construct the header dictionary
     head_dict['ngroups'] = ng
     head_dict['nk'] = el_struct.kmesh['nktot']
+    head_dict['nkibz'] = el_struct.kmesh['nkibz']
     head_dict['ns'] = el_struct.nspin
     head_dict['kvec1'] = list(el_struct.structure['kpt_basis'][:,0])
     head_dict['kvec2'] = list(el_struct.structure['kpt_basis'][:,1])

python/triqs_dft_tools/converters/plovasp/vaspio.py

@@ -40,6 +40,7 @@
 import numpy as np
 import re
 import os
+import time
 from h5 import HDFArchive
 
 log = logging.getLogger('plovasp.vaspio')
@@ -79,6 +80,10 @@ def __init__(self, vasp_dir, read_all=True, efermi_required=True):
         if vasph5:
             log.warning("Reading from vaspout.h5")
             h5path = os.path.join(vasp_dir, 'vaspout.h5')
+            # give VASP some time to write the file
+            with HDFArchive(h5path, 'r') as archive:
+                if 'locproj' not in archive['results']:
+                    time.sleep(2)
             self.plocar = h5Plocar(h5path)
             self.poscar = h5Poscar(h5path)
             self.kpoints = h5Kpoints(h5path)
@@ -401,7 +406,8 @@ class Kpoints:
     Class describing k-points and optionally tetrahedra.
 
     Properties:
-        - nktot (int) : total number of k-points in the IBZ
+        - nktot (int) : total number of k-points in the BZ
+        - nkibz (int) : number of k-points in the IBZ
         - kpts (numpy.array((nktot, 3), dtype=float)) : k-point vectors (fractional coordinates)
         - ntet (int) : total number of k-point tetrahedra
         - itet (numpy.array((ntet, 5), dtype=float) : array of tetrahedra
@@ -412,6 +418,7 @@ class Kpoints:
     def __init__(self):
         self.kpts = None
         self.nktot = None
+        self.nkibz = None
         self.kwghts = None
 
     #
@@ -441,6 +448,9 @@ def from_file(self, vasp_dir='./', ibz_filename='IBZKPT'):
         #   Number of k-points
         line = next(ibz_file)
         self.nktot = int(line.strip().split()[0])
+        # when reading from IBZKPT file we do not know the full number of k-points, i.e.
+        # works only ISYM=-1
+        self.nkibz = self.nktot
 
         print()
         print("   {0:>26} {1:d}".format("Total number of k-points:", self.nktot))
@@ -716,7 +726,7 @@ def __init__(self, h5path):
         # h5path = './vasptriqs.h5'
         with HDFArchive(h5path, 'a') as archive:
             kpoints = archive['results/electron_eigenvalues']
-            self.nkred = kpoints['kpoints']
+            self.nkibz = kpoints['kpoints']
             self.kpts = kpoints['kpoint_coords_full']
             self.nktot = len(self.kpts)
             self.kwghts = kpoints['kpoints_symmetry_weight_full']
@@ -731,7 +741,7 @@ def __init__(self, h5path):
                 self.ntet = 0
 
         print()
-        print("   {0:>26} {1:d}".format("Reduced number of k-points:", self.nkred))
+        print("   {0:>26} {1:d}".format("Reduced number of k-points:", self.nkibz))
         print("   {0:>26} {1:d}".format("Total number of k-points:", self.nktot))
         print("   {0:>26} {1:d}".format("Total number of tetrahedra:", self.ntet))
 

python/triqs_dft_tools/converters/vasp.py

@@ -158,6 +158,8 @@ def convert_dft_input(self):
 
         ng = ctrl_head['ngroups']
         n_k = ctrl_head['nk']
+        n_k_ibz = ctrl_head['nkibz']
+
 # Note the difference in name conventions!
         SP = ctrl_head['ns'] - 1
         SO = ctrl_head['nc_flag']
@@ -387,7 +389,7 @@ def convert_dft_input(self):
         with HDFArchive(self.hdf_file,'a') as ar:
             if not (self.dft_subgrp in ar): ar.create_group(self.dft_subgrp)
             # The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten!
-            things_to_save = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required',
+            things_to_save = ['energy_unit','n_k', 'k_dep_projection','SP','SO','charge_below','density_required',
                               'symm_op','n_shells','shells','n_corr_shells','corr_shells','use_rotations','rot_mat',
                               'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights',
                               'hopping','n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr','proj_or_hk',
@@ -401,6 +403,8 @@ def convert_dft_input(self):
             ar[self.misc_subgrp]['dft_fermi_weights'] = f_weights
             ar[self.misc_subgrp]['kpts_cart'] = kpts_cart
             ar[self.misc_subgrp]['band_window'] = band_window
+            if n_k_ibz is not None:
+                ar[self.misc_subgrp]['n_k_ibz'] = n_k_ibz
 
         # Symmetries are used, so now convert symmetry information for *correlated* orbitals:
         self.convert_symmetry_input(ctrl_head, orbits=self.corr_shells, symm_subgrp=self.symmcorr_subgrp)

python/triqs_dft_tools/sumk_dft.py

@@ -2207,8 +2207,7 @@ def calc_density_correction(self, filename=None, dm_type=None, spinave=False, kp
                 filename = 'dens_mat.dat'
             elif dm_type == 'vasp':
                 # use new h5 interface to vasp by default, if not wanted specify dm_type='vasp' + filename='GAMMA'
-                # filename = 'vaspgamma.h5'
-                filename = 'GAMMA'
+                filename = 'vaspgamma.h5'
             elif dm_type == 'elk':
                 filename = 'DMATDMFT.OUT'
             elif dm_type == 'qe':
@@ -2229,12 +2228,17 @@ def calc_density_correction(self, filename=None, dm_type=None, spinave=False, kp
         if dm_type in ['vasp', 'qe']:
             fermi_weights = 0
             band_window = 0
+            n_k_ibz = self.n_k
             if mpi.is_master_node():
                 with HDFArchive(self.hdf_file, 'r') as ar:
                     fermi_weights = ar['dft_misc_input']['dft_fermi_weights']
                     band_window = ar['dft_misc_input']['band_window']
+
+                    if 'n_k_ibz' in ar['dft_misc_input']:
+                        n_k_ibz = ar['dft_misc_input']['n_k_ibz']
             fermi_weights = mpi.bcast(fermi_weights)
             band_window = mpi.bcast(band_window)
+            n_k_ibz = mpi.bcast(n_k_ibz)
 
             # Convert Fermi weights to a density matrix
             dens_mat_dft = {}
@@ -2349,8 +2353,13 @@ def calc_density_correction(self, filename=None, dm_type=None, spinave=False, kp
                             fout.write("\n")
                         fout.close()
         elif dm_type == 'vasp':
-            if kpts_to_write is None:
+            if kpts_to_write is None and self.n_k == n_k_ibz:
                 kpts_to_write = np.arange(self.n_k)
+            elif kpts_to_write is None and n_k_ibz < self.n_k:
+                # If the number of IBZ k-points is less than the total number of k-points,
+                # then we only write the IBZ k-points, which are the first n_k_ibz k-points
+                # in VASP
+                kpts_to_write = np.arange(n_k_ibz)
             else:
                 assert np.min(kpts_to_write) >= 0 and np.max(kpts_to_write) < self.n_k
 
@@ -2359,8 +2368,12 @@ def calc_density_correction(self, filename=None, dm_type=None, spinave=False, kp
             if mpi.is_master_node():
                 if filename == 'vaspgamma.h5':
                     with HDFArchive('vaspgamma.h5', 'w') as vasp_h5:
-                        vasp_h5['band_window'] = band_window
-                        vasp_h5['deltaN'] = deltaN
+                        # only store the ibz kpoints in the h5
+                        bnd_win_towrite = [band_window[0][:n_k_ibz,:]]
+                        vasp_h5['band_window'] = bnd_win_towrite
+                        vasp_h5.create_group('deltaN')
+                        vasp_h5['deltaN']['up'] = deltaN['up'][:n_k_ibz]
+                        vasp_h5['deltaN']['down'] = deltaN['down'][:n_k_ibz]
                 else:
                     with open(filename, 'w') as f:
                         f.write(" -1  -1  ! Number of k-points, default number of bands\n")  # % len(kpts_to_write))