Symmetric DFT implementation

Author: Dario Fiore Mosca

.idea/.gitignore

@@ -72,22 +72,26 @@ def is_vasp_running(vasp_pid):
     pid_exists = mpi.bcast(pid_exists)
     return pid_exists
 
+
 def get_dft_energy():
     """
-    Reads energy from the last line of OSZICAR.
+    Reads DFT energy from the last line of Vasp's OSZICAR or from vasptriqs.h5
     """
-    with open('OSZICAR', 'r') as f:
-        nextline = f.readline()
-        while nextline.strip():
-            line = nextline
-            nextline = f.readline()
-#            print "OSZICAR: ", line[:-1]
+    h5_energy = False
+    if os.path.isfile('vaspout.h5'):
+        with HDFArchive('vaspout.h5', 'r') as h5:
+            if 'oszicar' in h5['intermediate/ion_dynamics']:
+                dft_energy = h5['intermediate/ion_dynamics/oszicar'][-1,1]
+                h5_energy = True
 
-    try:
+    # as backup use OSZICAR file
+    if not h5_energy:
+        with open('OSZICAR', 'r') as file:
+            nextline = file.readline()
+            while nextline.strip():
+                line = nextline
+                nextline = file.readline()
         dft_energy = float(line.split()[2])
-    except ValueError:
-        print("Cannot read energy from OSZICAR, setting it to zero")
-        dft_energy = 0.0
 
     return dft_energy
 
@@ -98,9 +102,8 @@ class bcolors:
     YELLOW = '\033[93m'
     RED = '\033[91m'
     ENDC = '\033[0m'
-#
+
 # Main self-consistent cycle
-#
 def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
     """
     """
@@ -117,15 +120,14 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
         mpi.barrier()
         while is_vasp_lock_present():
             time.sleep(1)
-#            if debug: print bcolors.YELLOW + " waiting: rank %s"%(mpi.rank) + bcolors.ENDC
+            if debug: print(bcolors.YELLOW + " waiting: rank %s"%(mpi.rank) + bcolors.ENDC)
             if not is_vasp_running(vasp_pid):
                 mpi.report("  VASP stopped")
                 vasp_running = False
                 break
 
-# Tell VASP to stop if the maximum number of iterations is reached
-        
-        
+        # Tell VASP to stop if the maximum number of iterations is reached
+
         if debug: print(bcolors.MAGENTA + "rank %s"%(mpi.rank) + bcolors.ENDC)
         err = 0
         exc = None
@@ -161,7 +163,7 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
         # electron.F around line 644
         iter_dft = 0
 
-        if vasp_version == 'standard':
+        if vasp_version == 'standard' or vasp_version == 'ncl':
             copyfile(src='GAMMA',dst='GAMMA_recent')
         while iter_dft < n_iter_dft:
             # insert recalculation of GAMMA here
@@ -190,7 +192,7 @@ def run_all(vasp_pid, dmft_cycle, cfg_file, n_iter, n_iter_dft, vasp_version):
                     vasp_running = False
                     break
             iter_dft += 1
-            if vasp_version == 'standard':
+            if vasp_version == 'standard' or vasp_version == 'ncl':
                 copyfile(src='GAMMA_recent',dst='GAMMA')
         iter += 1
         if iter == n_iter:
@@ -253,8 +255,8 @@ def main():
     except KeyError:
         vasp_version = 'standard'
 
-    if vasp_version != 'standard' and vasp_version != 'no_gamma_write':
-        raise Exception('vasp_version has to be standard or no_gamma_write')
+    #if vasp_version != 'standard' and vasp_version != 'no_gamma_write':
+    #    raise Exception('vasp_version has to be standard or no_gamma_write')
 
 #    if len(sys.argv) > 1:
 #        vasp_path = sys.argv[1]

.idea/dft_tools.src.iml

@@ -82,7 +82,7 @@ def __init__(self, vasp_dir, read_all=True, efermi_required=True):
             self.plocar = h5Plocar(h5path)
             self.poscar = h5Poscar(h5path)
             self.kpoints = h5Kpoints(h5path)
-            self.eigenval = h5Eigenval(h5path)
+            self.eigenval = h5Eigenval(h5path, self.kpoints.ksymmap)
             self.doscar = h5Doscar(h5path)
         else:
             self.plocar = Plocar()
@@ -716,28 +716,34 @@ def __init__(self, h5path):
         # h5path = './vasptriqs.h5'
         with HDFArchive(h5path, 'a') as archive:
             kpoints = archive['results/electron_eigenvalues']
-            self.nktot = kpoints['kpoints']
-            self.kpts = kpoints['kpoint_coords']
-            self.kwghts = kpoints['kpoints_symmetry_weight']
+            self.nkred = kpoints['kpoints']
+            self.kpts = kpoints['kpoint_coords_full']
+            self.nktot = len(self.kpts)
+            self.kwghts = kpoints['kpoints_symmetry_weight_full']
+            self.ksymmap = kpoints['kpoints_symmetry_mapping']
+            self.ksymmap -= 1
             try:
                 self.ntet = kpoints['num_tetrahedra']
                 self.vtet = kpoints['volume_weight_tetrahedra']
                 self.itet = kpoints['coordinate_id_tetrahedra']
             except KeyError:
-                print("  No tetrahedron data found in vasptriqs.h5. Skipping...")
+                print("  No tetrahedron data found in vaspout.h5. Skipping...")
                 self.ntet = 0
 
         print()
+        print("   {0:>26} {1:d}".format("Reduced number of k-points:", self.nkred))
         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))
 
 
 class h5Eigenval:
 
-    def __init__(self, h5path):
+    def __init__(self, h5path, symmap):
         with HDFArchive(h5path, 'a') as archive:
             self.eigs = archive['results/electron_eigenvalues']['eigenvalues']
+            self.eigs = self.eigs[:, symmap, :]
             self.ferw = archive['results/electron_eigenvalues']['fermiweights']
+            self.ferw = self.ferw[:, symmap, :]
         # TODO Change the format in VASP to have [kpoints, bands, spin]
         self.eigs = np.transpose(self.eigs, (1, 2, 0))
         self.ferw = np.transpose(self.ferw, (1, 2, 0))

.idea/inspectionProfiles/Project_Default.xml

@@ -578,7 +578,8 @@ def lattice_gf(self, ik, mu=None, broadening=None, mesh=None, with_Sigma=True, w
                 mesh_values = self.mesh_values
 
         elif not mesh is None:
-            assert isinstance(mesh, (MeshReFreq, MeshDLRImFreq, MeshImFreq)),  "mesh must be a triqs MeshReFreq or MeshImFreq"
+            assert isinstance(mesh,
+                              (MeshReFreq, MeshDLRImFreq, MeshImFreq)), "mesh must be a triqs MeshReFreq or MeshImFreq"
             if isinstance(mesh, MeshImFreq):
                 mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh))
             elif isinstance(mesh, MeshDLRImFreq):
@@ -596,7 +597,7 @@ def lattice_gf(self, ik, mu=None, broadening=None, mesh=None, with_Sigma=True, w
             gf_struct = [(spn[isp], block_structure[isp])
                          for isp in range(self.n_spin_blocks[self.SO])]
             block_ind_list = [block for block, inner in gf_struct]
-            glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner),len(inner)])
+            glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner), len(inner)])
                              for block, inner in gf_struct]
             G_latt = BlockGf(name_list=block_ind_list,
                              block_list=glist(), make_copies=False)
@@ -610,7 +611,7 @@ def lattice_gf(self, ik, mu=None, broadening=None, mesh=None, with_Sigma=True, w
             ind = ntoi[spn[ibl]]
             n_orb = self.n_orbitals[ik, ind]
             if isinstance(mesh, (MeshImFreq, MeshDLRImFreq)):
-                gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl))
+                gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field * (1 - 2 * ibl))
                                     - self.hopping[ik, ind, 0:n_orb, 0:n_orb])
             else:
                 gf.data[:, :, :] = (idmat[ibl] *
@@ -655,8 +656,8 @@ def put_Sigma(self, Sigma_imp, transform_to_sumk_blocks=True):
 
         if (isinstance(self.mesh, (MeshImFreq, MeshDLRImFreq)) and
                 all(isinstance(gf.mesh, (MeshImFreq, MeshDLRImFreq)) and
-                isinstance(gf, Gf) and
-                gf.mesh == self.mesh for bname, gf in Sigma_imp[0])):
+                    isinstance(gf, Gf) and
+                    gf.mesh == self.mesh for bname, gf in Sigma_imp[0])):
             # Imaginary frequency Sigma:
             self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk')
                               for icrsh in range(self.n_corr_shells)]
@@ -692,7 +693,7 @@ def put_Sigma(self, Sigma_imp, transform_to_sumk_blocks=True):
                     self.max_band_energy - self.chemical_potential):
                 warn(
                     'The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f' % (
-                    mesh[0], mesh[-1], self.min_band_energy, self.max_band_energy))
+                        mesh[0], mesh[-1], self.min_band_energy, self.max_band_energy))
 
     def transform_to_sumk_blocks(self, Sigma_imp, Sigma_out=None):
         r""" transform Sigma from solver to sumk space
@@ -1843,7 +1844,8 @@ def add_dc(self):
             for bname, gf in sigma_minus_dc[icrsh]:
                 # Transform dc_imp to global coordinate system
                 if self.use_rotations:
-                    gf -= np.dot(self.rot_mat[icrsh], np.dot(self.dc_imp[icrsh][bname], self.rot_mat[icrsh].conjugate().transpose()))
+                    gf -= np.dot(self.rot_mat[icrsh],
+                                 np.dot(self.dc_imp[icrsh][bname], self.rot_mat[icrsh].conjugate().transpose()))
                 else:
                     gf -= self.dc_imp[icrsh][bname]
 
@@ -2205,7 +2207,8 @@ 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 = 'vaspgamma.h5'
+                filename = 'GAMMA'
             elif dm_type == 'elk':
                 filename = 'DMATDMFT.OUT'
             elif dm_type == 'qe':
@@ -2360,7 +2363,7 @@ def calc_density_correction(self, filename=None, dm_type=None, spinave=False, kp
                         vasp_h5['deltaN'] = deltaN
                 else:
                     with open(filename, 'w') as f:
-                        f.write(" %i  -1  ! Number of k-points, default number of bands\n" % len(kpts_to_write))
+                        f.write(" -1  -1  ! Number of k-points, default number of bands\n")  # % len(kpts_to_write))
                         for index, ik in enumerate(kpts_to_write):
                             ib1 = band_window[0][ik, 0]
                             ib2 = band_window[0][ik, 1]
@@ -2372,8 +2375,10 @@ def calc_density_correction(self, filename=None, dm_type=None, spinave=False, kp
                                         valim = (deltaN['ud'][ik][inu, imu].imag) / 1.0
                                         f.write(" %.14f  %.14f" % (valre, valim))
                                     else:
-                                        valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
-                                        valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
+                                        valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][
+                                            inu, imu].real) / 2.0
+                                        valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][
+                                            inu, imu].imag) / 2.0
                                         f.write(" %.14f  %.14f" % (valre, valim))
                                 f.write("\n")
 
@@ -2396,7 +2401,7 @@ def calc_density_correction(self, filename=None, dm_type=None, spinave=False, kp
                     mu = self.chemical_potential / self.energy_unit
                     # ouput n_k, nspin and max orbitals - a check
                     f.write(" %d  %d  %d  %.14f %.14f ! nkpt, nspin, nstmax, beta, mu\n" % (
-                    self.n_k, n_spin_blocks, nbmax, beta, mu))
+                        self.n_k, n_spin_blocks, nbmax, beta, mu))
                     for ik in range(self.n_k):
                         for ispn in range(n_spin_blocks):
                             # Determine the SO density matrix band indices from the spinor band indices