diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py index 14d16bf0..0a509415 100644 --- a/python/triqs_dft_tools/sumk_dft.py +++ b/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.""" diff --git a/python/triqs_dft_tools/util.py b/python/triqs_dft_tools/util.py new file mode 100644 index 00000000..59fc7cd7 --- /dev/null +++ b/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 . +# +########################################################################## + +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 diff --git a/test/python/dc_test.py b/test/python/dc_test.py index ba003355..adf527e7 100755 --- a/test/python/dc_test.py +++ b/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')