small adjustments to Vasp converter (closes #146)

* adding kpts_basis, kpts, and kpt_weights to h5 dft_input
* read these properties as optional arguments in Sumk_dft.py
* change accordingly the ref h5 files for vasp converter test
* soon all converters are demanted to store those properties
* bz_weights should then be replaced by kpt_weights
* closes PR #146

python/triqs_dft_tools/converters/plovasp/elstruct.py

@@ -118,25 +118,18 @@ class ElectronicStructure:
 #        assert natom_plo == self.natom, "PLOCAR is inconsistent with POSCAR (number of atoms)"
         self.structure = {'a_brav': vasp_data.poscar.a_brav}
         self.structure['nqtot'] = vasp_data.poscar.nq
+        self.structure['kpt_basis'] = vasp_data.poscar.kpt_basis
         self.structure['ntypes'] = vasp_data.poscar.ntypes
         self.structure['nq_types'] = vasp_data.poscar.nions
 # Concatenate coordinates grouped by type into one array
         self.structure['qcoords'] = np.vstack(vasp_data.poscar.q_types)
         self.structure['type_of_ion'] = vasp_data.poscar.type_of_ion
 
-        a = []
-        for ia in range(3):
-            a.append( vasp_data.poscar.a_brav[:,ia])
-        vol = np.dot(a[0],np.cross(a[1],a[2]))
-        b1 = 2.0*np.pi*np.cross(a[1],a[2])/vol
-        b2 = 2.0*np.pi*np.cross(a[2],a[0])/vol
-        b3 = 2.0*np.pi*np.cross(a[0],a[1])/vol
-        b = [b1,b2,b3]
         self.kmesh['kpoints_cart'] = 0.0 * self.kmesh['kpoints']
 
         for ik in range(self.nktot):
             for ii in range(3):
-                self.kmesh['kpoints_cart'][ik] += self.kmesh['kpoints'][ik,ii]*b[ii]
+                self.kmesh['kpoints_cart'][ik] += self.kmesh['kpoints'][ik,ii]*self.structure['kpt_basis'][:,ii]
 
 # FIXME: This can be removed if ion coordinates are stored in a continuous array
 ## Construct a map to access coordinates by index

python/triqs_dft_tools/converters/plovasp/plotools.py

@@ -310,6 +310,9 @@ def ctrl_output(conf_pars, el_struct, ng):
     head_dict['ngroups'] = ng
     head_dict['nk'] = el_struct.kmesh['nktot']
     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])
+    head_dict['kvec3'] = list(el_struct.structure['kpt_basis'][:,2])
     head_dict['nc_flag'] = 1 if el_struct.nc_flag else 0
 #    head_dict['efermi'] = conf_pars.general['efermi']  # We probably don't need Efermi
 

python/triqs_dft_tools/converters/plovasp/vaspio.py

@@ -337,12 +337,12 @@ class Poscar:
         - ntypes ([int]) : number of ion types
         - nions (int) : a list of number of ions of each type
         - a_brav (numpy.array((3, 3), dtype=float)) : lattice vectors
+        - kpt_basis (numpy.array((3, 3), dtype=float)) : reciprocal lattice vectors
         - q_types ([numpy.array((nions, 3), dtype=float)]) : a list of
           arrays each containing fractional coordinates of ions of a given type
     """
     def __init__(self):
         self.q_cart = None
-        self.b_rec = None
 
     def from_file(self, vasp_dir='./', poscar_filename='POSCAR'):
         """
@@ -409,8 +409,8 @@ class Poscar:
 
 # Check whether coordinates are cartesian or fractional
         cartesian = (sline[0].lower() in 'ck')
-        if cartesian:
-            brec = np.linalg.inv(self.a_brav.T)
+        # determine reciprocal basis in units of 2*pi
+        self.kpt_basis = np.linalg.inv(self.a_brav.T)
 
 # Read atomic positions
         self.q_types = []
@@ -424,7 +424,7 @@ class Poscar:
                 sline = readline_remove_comments()
                 qcoord = list(map(float, sline.split()[:3]))
                 if cartesian:
-                    qcoord = np.dot(brec, qcoord)
+                    qcoord = np.dot(self.kpt_basis, qcoord)
                 q_at_it[iq, :] = qcoord
 
             self.q_types.append(q_at_it)

python/triqs_dft_tools/converters/vasp.py

@@ -159,14 +159,23 @@ class VaspConverter(ConverterTools):
         SP = ctrl_head['ns'] - 1
         SO = ctrl_head['nc_flag']
 
+        # load reciprocal basis
+        kpt_basis = numpy.zeros((3,3))
+        kpt_basis[:,0] = ctrl_head['kvec1']
+        kpt_basis[:,1] = ctrl_head['kvec2']
+        kpt_basis[:,2] = ctrl_head['kvec3']
+
         kpts = numpy.zeros((n_k, 3))
         kpts_cart = numpy.zeros((n_k, 3))
         bz_weights = numpy.zeros(n_k)
+        kpt_weights = numpy.zeros(n_k)
         try:
             for ik in range(n_k):
                 kx, ky, kz = next(rf), next(rf), next(rf)
                 kpts[ik, :] = kx, ky, kz
                 bz_weights[ik] = next(rf)
+                # bz_weights soon to be removed, and replaced by kpt_weights
+                kpt_weights[ik] = bz_weights[ik]
             for ik in range(n_k):
                 kx, ky, kz = next(rf), next(rf), next(rf)
                 kpts_cart[ik, :] = kx, ky, kz
@@ -380,16 +389,18 @@ class VaspConverter(ConverterTools):
             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',
-                          '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','kpts','kpts_cart']
-            if self.proj_or_hk == 'hk' or True:
+                              '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',
+                              'kpts','kpt_weights', 'kpt_basis']
+            if self.proj_or_hk == 'hk' or self.proj_or_hk ==  True:
                 things_to_save.append('proj_mat_csc')
             for it in things_to_save: ar[self.dft_subgrp][it] = locals()[it]
 
             # Store Fermi weights to 'dft_misc_input'
             if not (self.misc_subgrp in ar): ar.create_group(self.misc_subgrp)
             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
 
         # Symmetries are used, so now convert symmetry information for *correlated* orbitals:

python/triqs_dft_tools/sumk_dft.py

@@ -108,7 +108,8 @@ class SumkDFT(object):
                 raise ValueError('ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:',self.values_not_read)
 
             # optional properties to load
-            optional_things_to_read = ['proj_mat_csc', 'proj_or_hk', 'kpts_cart']
+            # soon bz_weights is depraced and replaced by kpt_weights, kpts_basis and kpts will become required to read soon
+            optional_things_to_read = ['proj_mat_csc', 'proj_or_hk', 'kpts_basis','kpts','kpt_weights', 'kpt_basis']
             subgroup_present, self.optional_values_not_read = self.read_input_from_hdf(subgrp=self.dft_data, things_to_read=optional_things_to_read)
 
             if self.symm_op: