[fix] calc_density_correction: see issue #250

* fix not initialized variable Glatt_iw
* add simple run tests for all 4 DFT codes to test at least if
  calc_density_correction runs without errors
* small formatting fixes

python/triqs_dft_tools/sumk_dft.py

@@ -30,7 +30,7 @@ import triqs.utility.dichotomy as dichotomy
 from triqs.gf import *
 import triqs.utility.mpi as mpi
 from triqs.utility.comparison_tests import assert_arrays_are_close
-from h5 import *
+from h5 import HDFArchive
 from .symmetry import *
 from .block_structure import BlockStructure
 from .util import compute_DC_from_density
@@ -2128,7 +2128,7 @@ class SumkDFT(object):
 
         """
         #automatically set dm_type if required
-        if dm_type==None:
+        if dm_type is None:
             dm_type = self.dft_code
 
         assert dm_type in ('vasp', 'wien2k','elk', 'qe'), "'dm_type' must be either 'vasp', 'wienk', 'elk' or 'qe'"
@@ -2242,9 +2242,9 @@ class SumkDFT(object):
                     f1.write("%.14f\n" %
                              (self.chemical_potential / self.energy_unit))
                 # write beta in rydberg-1
-                f.write("%.14f\n" % (G_latt_iw.mesh.beta * self.energy_unit))
+                f.write("%.14f\n" % (self.mesh.beta * self.energy_unit))
                 if self.SP != 0:
-                    f1.write("%.14f\n" % (G_latt_iw.mesh.beta * self.energy_unit))
+                    f1.write("%.14f\n" % (self.mesh.beta * self.energy_unit))
 
                 if self.SP == 0:  # no spin-polarization
 
@@ -2314,7 +2314,7 @@ class SumkDFT(object):
                     n_spin_blocks = self.SP + 1 - self.SO
                     nbmax = np.max(self.n_orbitals)
         # output beta and mu in Hartrees
-                    beta = G_latt_iw.mesh.beta * self.energy_unit
+                    beta = self.mesh.beta * self.energy_unit
                     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))
@@ -2361,7 +2361,7 @@ class SumkDFT(object):
                         delta_N[ik, inu, imu] = valre + 1j*valim
             if mpi.is_master_node():
                 with HDFArchive(self.hdf_file, 'a') as ar:
-                    if not subgrp in ar:
+                    if subgrp not in ar:
                         ar.create_group(subgrp)
                     things_to_save = ['delta_N']
                     for it in things_to_save:

test/python/calc_density_correction.py

@@ -0,0 +1,92 @@
+# Copyright (c) 2013 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
+# Copyright (c) 2013 Centre national de la recherche scientifique (CNRS)
+# Copyright (c) 2019-2024 Simons Foundation
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You may obtain a copy of the License at
+#     https:#www.gnu.org/licenses/gpl-3.0.txt
+#
+# Authors: A. Hampel
+
+import shutil
+from triqs.gf import Gf, MeshImFreq
+from triqs_dft_tools.sumk_dft import SumkDFT
+
+# Simple run w/o error test for all DFT codes to write the density correction to a file
+# Comparison with real self energy against refrence data should be done in the future
+
+# define mesh for all calculations
+beta = 40
+mesh = MeshImFreq(beta, statistic='Fermion', n_iw=1024)
+
+# Wien2k test
+sumk = SumkDFT(hdf_file='SrVO3.ref.h5', mesh=mesh)
+
+Sigma_iw = [
+    sumk.block_structure.create_gf(ish=iineq, gf_function=Gf, space='solver', mesh=sumk.mesh) for iineq in range(sumk.n_inequiv_shells)
+]
+
+for iineq in range(sumk.n_inequiv_shells):
+    Sigma_iw[iineq] << 0.1 + 0.0j
+
+sumk.set_Sigma(Sigma_iw)
+
+deltaN, dens = sumk.calc_density_correction(dm_type='wien2k')
+
+####################################
+# Elk
+sumk = SumkDFT(hdf_file='elk/elk_convert/elk_convert.ref.h5', mesh=mesh)
+
+Sigma_iw = [
+    sumk.block_structure.create_gf(ish=iineq, gf_function=Gf, space='solver', mesh=sumk.mesh) for iineq in range(sumk.n_inequiv_shells)
+]
+
+for iineq in range(sumk.n_inequiv_shells):
+    Sigma_iw[iineq] << 0.1 + 0.0j
+
+sumk.set_Sigma(Sigma_iw)
+
+deltaN, dens = sumk.calc_density_correction(dm_type='elk')
+
+####################################
+# Vasp
+sumk = SumkDFT(hdf_file='plovasp/converter/lunio3.ref.h5', mesh=mesh)
+
+Sigma_iw = [
+    sumk.block_structure.create_gf(ish=iineq, gf_function=Gf, space='solver', mesh=sumk.mesh) for iineq in range(sumk.n_inequiv_shells)
+]
+
+for iineq in range(sumk.n_inequiv_shells):
+    Sigma_iw[iineq] << 0.1 + 0.0j
+
+sumk.set_Sigma(Sigma_iw)
+
+deltaN, dens, en_corr = sumk.calc_density_correction(dm_type='vasp')
+
+####################################
+# QE
+shutil.copy('w90_convert/SrVO3_col_blochbasis.ref.h5', 'w90_convert/SrVO3_col_blochbasis.test.h5')
+sumk = SumkDFT(hdf_file='w90_convert/SrVO3_col_blochbasis.test.h5', mesh=mesh)
+
+Sigma_iw = [
+    sumk.block_structure.create_gf(ish=iineq, gf_function=Gf, space='solver', mesh=sumk.mesh) for iineq in range(sumk.n_inequiv_shells)
+]
+
+for iineq in range(sumk.n_inequiv_shells):
+    Sigma_iw[iineq] << 0.1 + 0.0j
+
+sumk.set_Sigma(Sigma_iw)
+# little hack to speed up the calculation
+sumk.n_k = 1
+
+deltaN, dens, en_corr = sumk.calc_density_correction(dm_type='qe', filename='dump.h5')
+