outsource calc_DC_from_density into util.py and cleanup

python/triqs_dft_tools/sumk_dft.py

@@ -33,93 +33,12 @@ from triqs.utility.comparison_tests import assert_arrays_are_close
 from h5 import *
 from .symmetry import *
 from .block_structure import BlockStructure
+from .util import compute_DC_from_density
 from itertools import product
 from warnings import warn
 from scipy import compress
 from scipy.optimize import minimize, newton, brenth
 
-def compute_DC_from_density(N_tot, U, J, N_spin = None,  n_orbitals=5,  method='sFLL'):
-    """
-    Computes the double counting correction using various methods.
-    For FLL and AMF DC the notations and equations from doi.org/10.1038/s41598-018-27731-4
-    are used, whereas for the Held DC the definitions from doi.org/10.1080/00018730701619647 are used.
-
-    Parameters
-    ----------
-    N_tot : float 
-        Total density of the impurity
-    N_spin : float , default = None
-        Spin density, defaults to N_tot*0.5 if not specified
-    U : float 
-        U value
-    J : float 
-        J value
-    n_orbitals : int, default = 5
-        Total number of orbitals
-    method : string, default = 'cFLL' 
-        possibilities:
-        -    cFLL: DC potential from Ryee for spin unpolarized DFT: (DOI: 10.1038/s41598-018-27731-4)
-        -    sFLL: same as above for spin polarized DFT
-        -    cAMF: around mean field
-        -    sAMF: spin polarized around mean field
-        -    cHeld: unpolarized Held's formula as reported in (DOI: 10.1103/PhysRevResearch.2.03308)
-        -    sHeld: NOT IMPLEMENTED
-    
-    Returns
-    -------
-    List of floats:
-        -   DC_val: double counting potential
-        -   E_val: double counting energy
-
-
-    todo: See whether to move this to TRIQS directly instead of dft_tools
-    """
-    if N_spin is not None:
-        N_spin2 = N_tot-N_spin
-        Mag = N_spin - N_spin2
-    L_orbit = (n_orbitals-1)/2 
-
-    if method == 'cFLL':
-        E_val = 0.5 * U * N_tot * (N_tot-1) - 0.5 * J * N_tot * (N_tot*0.5-1)
-        DC_val = U * (N_tot-0.5) - J *(N_tot*0.5-0.5)
-
-    elif method == 'sFLL':
-        assert N_spin is not None, "Spin density not given"
-        E_val = 0.5 * U * N_tot * (N_tot-1) - 0.5 * J * N_tot * (N_tot*0.5-1) - 0.25 * J * Mag**2
-        DC_val = U * (N_tot-0.5) - J *(N_spin-0.5)
-        
-    elif method == 'cAMF':
-        E_val = +0.5 * U * N_tot **2
-        E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot**2
-        DC_val = U * N_tot - 0.5*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot
-            
-
-    elif method == 'sAMF':
-        assert N_spin is not None, "Spin density not given"
-        E_val = 0.5 * U * N_tot **2
-        E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot**2
-        E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*Mag**2
-        DC_val = U * N_tot  - (U+2*L_orbit*J)/(2*L_orbit+1)*N_spin
-        
-    elif method == 'cHeld':
-        # Valid for a Kanamori-type Hamiltonian where U'=U-2J
-        U_mean = (U + (n_orbitals-1)*(U-2*J)+(n_orbitals-1)*(U-3*J))/(2*n_orbitals-1)
-        E_val = 0.5 * U_mean * N_tot * (N_tot - 1)
-        DC_val = U_mean * (N_tot-0.5)
-        
-    elif method == 'sHeld':
-        raise ValueError(f"Method sHeld not yet implemented")
-        
-    else:
-        raise ValueError(f"DC type {method} not supported")
-
-    mpi.report(f"DC potential computed using the {method} method, V_DC = {DC_val:.6f} eV")
-    mpi.report(f"E_DC using the {method} method, E_DC = {E_val:.6f} eV")
-    if 'Held' in method:
-        mpi.report(f"Held method for {n_orbitals} orbitals, computed U_mean={U_mean:.6f} eV")
-
-    return DC_val, E_val
-
 
 class SumkDFT(object):
     """This class provides a general SumK method for combining ab-initio code and triqs."""

python/triqs_dft_tools/util.py

@@ -0,0 +1,107 @@
+##########################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2023 by A. Carta, A. Hampel
+#
+# TRIQS 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.
+#
+# TRIQS 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 should have received a copy of the GNU General Public License along with
+# TRIQS. If not, see <http://www.gnu.org/licenses/>.
+#
+##########################################################################
+
+import numpy as np
+
+import triqs.utility.mpi as mpi
+
+def compute_DC_from_density(N_tot, U, J, N_spin=None,  n_orbitals=5,  method='sFLL'):
+    """
+    Computes the double counting correction using various methods.
+    For FLL and AMF DC the notations and equations from doi.org/10.1038/s41598-018-27731-4
+    are used, whereas for the Held DC the definitions from doi.org/10.1080/00018730701619647 are used.
+
+    Parameters
+    ----------
+    N_tot : float 
+        Total density of the impurity
+    N_spin : float , default = None
+        Spin density, defaults to N_tot*0.5 if not specified
+    U : float 
+        U value
+    J : float 
+        J value
+    n_orbitals : int, default = 5
+        Total number of orbitals
+    method : string, default = 'cFLL' 
+        possibilities:
+        -    cFLL: DC potential from Ryee for spin unpolarized DFT: (DOI: 10.1038/s41598-018-27731-4)
+        -    sFLL: same as above for spin polarized DFT
+        -    cAMF: around mean field
+        -    sAMF: spin polarized around mean field
+        -    cHeld: unpolarized Held's formula as reported in (DOI: 10.1103/PhysRevResearch.2.03308)
+        -    sHeld: NOT IMPLEMENTED
+    
+    Returns
+    -------
+    List of floats:
+        -   DC_val: double counting potential
+        -   E_val: double counting energy
+
+
+    todo: 
+        - See whether to move this to TRIQS directly instead of dft_tools
+        - allow as input full density matrix to allow orbital dependent DC
+    """
+    if N_spin is not None:
+        N_spin2 = N_tot-N_spin
+        Mag = N_spin - N_spin2
+    L_orbit = (n_orbitals-1)/2
+
+    if method == 'cFLL':
+        E_val = 0.5 * U * N_tot * (N_tot-1) - 0.5 * J * N_tot * (N_tot*0.5-1)
+        DC_val = U * (N_tot-0.5) - J * (N_tot*0.5-0.5)
+
+    elif method == 'sFLL':
+        assert N_spin is not None, "Spin density not given"
+        E_val = 0.5 * U * N_tot * (N_tot-1) - 0.5 * J * N_tot * (N_tot*0.5-1) - 0.25 * J * Mag**2
+        DC_val = U * (N_tot-0.5) - J * (N_spin-0.5)
+
+    elif method == 'cAMF':
+        E_val = +0.5 * U * N_tot ** 2
+        E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot**2
+        DC_val = U * N_tot - 0.5*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot
+
+    elif method == 'sAMF':
+        assert N_spin is not None, "Spin density not given"
+        E_val = 0.5 * U * N_tot ** 2
+        E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot**2
+        E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*Mag**2
+        DC_val = U * N_tot - (U+2*L_orbit*J)/(2*L_orbit+1)*N_spin
+
+    elif method == 'cHeld':
+        # Valid for a Kanamori-type Hamiltonian where U'=U-2J
+        U_mean = (U + (n_orbitals-1)*(U-2*J)+(n_orbitals-1)*(U-3*J))/(2*n_orbitals-1)
+        E_val = 0.5 * U_mean * N_tot * (N_tot - 1)
+        DC_val = U_mean * (N_tot-0.5)
+
+    elif method == 'sHeld':
+        raise ValueError(f"Method sHeld not yet implemented")
+
+    else:
+        raise ValueError(f"DC type {method} not supported")
+
+    mpi.report(f"DC potential computed using the {method} method, V_DC = {DC_val:.6f} eV")
+    mpi.report(f"E_DC using the {method} method, E_DC = {E_val:.6f} eV")
+    if 'Held' in method:
+        mpi.report(f"Held method for {n_orbitals} orbitals, computed U_mean={U_mean:.6f} eV")
+
+    return DC_val, E_val

test/python/dc_test.py

@@ -23,6 +23,7 @@
 
 import numpy as np
 from triqs_dft_tools.sumk_dft import *
+from triqs_dft_tools.util import compute_DC_from_density
 from triqs.gf import  *
 from h5 import HDFArchive
 from triqs.operators.util import *
@@ -50,7 +51,7 @@ method_dict = {
 
 
 
-def test_dc(method, method_dict, dens, Uval, Jval, filename):
+def test_dc(SK_compat, SK_new, method, method_dict, dens, Uval, Jval, filename):
 
     dc_no = method_dict[method]["numbering_convention"]
     dc_string = method_dict[method]["new_convention"]
@@ -117,7 +118,7 @@ N_tot = N_up + N_down
 mpi.report(f"{N_up=} ,{N_down=}, {N_tot=}\n")
 
 for method in ["FLL", "AMF", "Held"]: 
-    test_dc(method, method_dict, dens, Uval, Jval, filename = f"{dft_filename}.h5")
+    test_dc(SK_compat, SK_new, method, method_dict, dens, Uval, Jval, filename = f"{dft_filename}.h5")
 
     #in case implementation changes, to write new testing data into archive
     #R = HDFArchive('./NiO.ref.h5', 'a')
@@ -158,7 +159,7 @@ Uval = 5
 Jval = 0.3
 
 for method in ["FLL", "AMF", "Held"]: 
-    test_dc(method, method_dict, dens, Uval, Jval, filename = f"{dft_filename}.h5" )
+    test_dc(SK_compat, SK_new, method, method_dict, dens, Uval, Jval, filename = f"{dft_filename}.h5" )
     
     #in case implementation changes, to write new testing data into archive
     #R = HDFArchive(f'./{dft_filename}.h5', 'a')