mirror of
https://github.com/triqs/dft_tools
synced 2024-09-12 05:38:31 +02:00
Corrections for new calc_dc and calc_mu
* typos and subbed check in spin polarized calculations for quantum espresso with a warning when computing the deltaN * fixed typos in comments * removed legacy mode maintaining only compatibility layer and switched to old (<3.10) python syntax * added target density output in mu finder for brent and newton, refactored tunit test for DC, changed some comments
This commit is contained in:
parent
d4f3c48784
commit
0d25aefc73
1
.gitignore
vendored
1
.gitignore
vendored
@ -1,3 +1,2 @@
|
|||||||
compile_commands.json
|
compile_commands.json
|
||||||
doc/cpp2rst_generated
|
doc/cpp2rst_generated
|
||||||
build/
|
|
||||||
|
@ -38,10 +38,11 @@ from warnings import warn
|
|||||||
from scipy import compress
|
from scipy import compress
|
||||||
from scipy.optimize import minimize, newton, brenth
|
from scipy.optimize import minimize, newton, brenth
|
||||||
|
|
||||||
def compute_DC_from_density(N_tot, U, J, N_spin = None, n_orbitals=5, method='cFLL'):
|
def compute_DC_from_density(N_tot, U, J, N_spin = None, n_orbitals=5, method='sFLL'):
|
||||||
"""
|
"""
|
||||||
Computes the double counting correction in terms using various methods
|
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
|
Parameters
|
||||||
----------
|
----------
|
||||||
@ -69,45 +70,48 @@ def compute_DC_from_density(N_tot, U, J, N_spin = None, n_orbitals=5, method='
|
|||||||
List of floats:
|
List of floats:
|
||||||
- DC_val: double counting potential
|
- DC_val: double counting potential
|
||||||
- E_val: double counting energy
|
- E_val: double counting energy
|
||||||
|
|
||||||
|
|
||||||
|
todo: See whether to move this to TRIQS directly instead of dft_tools
|
||||||
"""
|
"""
|
||||||
if N_spin is not None:
|
if N_spin is not None:
|
||||||
N_spin2 = N_tot-N_spin
|
N_spin2 = N_tot-N_spin
|
||||||
Mag = N_spin - N_spin2
|
Mag = N_spin - N_spin2
|
||||||
L_orbit = (n_orbitals-1)/2
|
L_orbit = (n_orbitals-1)/2
|
||||||
|
|
||||||
match method:
|
if method == 'cFLL':
|
||||||
case 'cFLL':
|
E_val = 0.5 * U * N_tot * (N_tot-1) - 0.5 * J * N_tot * (N_tot*0.5-1)
|
||||||
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)
|
||||||
DC_val = U * (N_tot-0.5) - J *(N_tot*0.5-0.5)
|
|
||||||
|
|
||||||
case 'sFLL':
|
elif method == 'sFLL':
|
||||||
assert N_spin is not None, "Spin density not given"
|
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
|
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)
|
DC_val = U * (N_tot-0.5) - J *(N_spin-0.5)
|
||||||
|
|
||||||
case 'cAMF':
|
elif method == 'cAMF':
|
||||||
E_val = +0.5 * U * N_tot **2
|
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)*N_tot**2
|
||||||
DC_val = U * N_tot - 0.5*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot
|
DC_val = U * N_tot - 0.5*(U+2*L_orbit*J)/(2*L_orbit+1)*N_tot
|
||||||
|
|
||||||
|
|
||||||
case 'sAMF':
|
elif method == 'sAMF':
|
||||||
assert N_spin is not None, "Spin density not given"
|
assert N_spin is not None, "Spin density not given"
|
||||||
E_val = 0.5 * U * N_tot **2
|
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)*N_tot**2
|
||||||
E_val -= 0.25*(U+2*L_orbit*J)/(2*L_orbit+1)*Mag**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
|
DC_val = U * N_tot - (U+2*L_orbit*J)/(2*L_orbit+1)*N_spin
|
||||||
|
|
||||||
case 'cHeld':
|
elif method == 'cHeld':
|
||||||
U_mean = (U + (n_orbitals-1)*(U-2*J)+(n_orbitals-1)*(U-3*J))/(2*n_orbitals-1)
|
# Valid for a Kanamori-type Hamiltonian where U'=U-2J
|
||||||
E_val = 0.5 * U_mean * N_tot * (N_tot - 1)
|
U_mean = (U + (n_orbitals-1)*(U-2*J)+(n_orbitals-1)*(U-3*J))/(2*n_orbitals-1)
|
||||||
DC_val = U_mean * (N_tot-0.5)
|
E_val = 0.5 * U_mean * N_tot * (N_tot - 1)
|
||||||
|
DC_val = U_mean * (N_tot-0.5)
|
||||||
|
|
||||||
case 'sHeld':
|
elif method == 'sHeld':
|
||||||
raise ValueError(f"Method sHeld not yet implemented")
|
raise ValueError(f"Method sHeld not yet implemented")
|
||||||
|
|
||||||
case _:
|
else:
|
||||||
raise ValueError(f"DC type {method} not supported")
|
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"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")
|
mpi.report(f"E_DC using the {method} method, E_DC = {E_val:.6f} eV")
|
||||||
@ -1735,7 +1739,7 @@ class SumkDFT(object):
|
|||||||
self.dc_energ = dc_energ
|
self.dc_energ = dc_energ
|
||||||
|
|
||||||
def calc_dc(self, dens_mat, orb=0, U_interact=None, J_hund=None,
|
def calc_dc(self, dens_mat, orb=0, U_interact=None, J_hund=None,
|
||||||
use_dc_formula=0, use_dc_value=None, transform=True, legacy=True):
|
use_dc_formula=0, use_dc_value=None, transform=True):
|
||||||
r"""
|
r"""
|
||||||
Calculate and set the double counting corrections.
|
Calculate and set the double counting corrections.
|
||||||
|
|
||||||
@ -1767,8 +1771,13 @@ class SumkDFT(object):
|
|||||||
Value of interaction parameter `U`.
|
Value of interaction parameter `U`.
|
||||||
J_hund : float, optional
|
J_hund : float, optional
|
||||||
Value of interaction parameter `J`.
|
Value of interaction parameter `J`.
|
||||||
use_dc_formula : int, optional
|
use_dc_formula : int or string, optional
|
||||||
Type of double-counting correction (see description).
|
Type of double-counting correction (see description of `compute_DC_from_density` above).
|
||||||
|
There is an interface with the legacy implementation which allows for the old convention:
|
||||||
|
* 0 -> 'sFLL' spin dependent fully localized limit
|
||||||
|
* 1 -> 'cHeld' spin independent Held formula
|
||||||
|
* 2 -> 'sAMF' spin dependent around-mean field approximation
|
||||||
|
|
||||||
use_dc_value : float, optional
|
use_dc_value : float, optional
|
||||||
Value of the double-counting correction. If specified
|
Value of the double-counting correction. If specified
|
||||||
`U_interact`, `J_hund` and `use_dc_formula` are ignored.
|
`U_interact`, `J_hund` and `use_dc_formula` are ignored.
|
||||||
@ -1807,68 +1816,21 @@ class SumkDFT(object):
|
|||||||
dim //= 2
|
dim //= 2
|
||||||
|
|
||||||
if use_dc_value is None:
|
if use_dc_value is None:
|
||||||
if legacy:
|
#For legacy compatibility
|
||||||
|
if use_dc_formula == 0:
|
||||||
if U_interact is None and J_hund is None:
|
mpi.report(f"Detected {use_dc_formula=}, changing to sFLL")
|
||||||
raise ValueError("set_dc: either provide U_interact and J_hund or set use_dc_value to dc value.")
|
use_dc_formula = "sFLL"
|
||||||
|
if use_dc_formula == 1:
|
||||||
if use_dc_formula == 0: # FLL
|
mpi.report(f"Detected {use_dc_formula=}, changing to cHeld")
|
||||||
|
use_dc_formula = "cHeld"
|
||||||
self.dc_energ[icrsh] = U_interact / \
|
if use_dc_formula == 2:
|
||||||
2.0 * Ncrtot * (Ncrtot - 1.0)
|
mpi.report(f"Detected {use_dc_formula=}, changing to sAMF")
|
||||||
for sp in spn:
|
use_dc_formula = "sAMF"
|
||||||
Uav = U_interact * (Ncrtot - 0.5) - \
|
|
||||||
J_hund * (Ncr[sp] - 0.5)
|
|
||||||
self.dc_imp[icrsh][sp] *= Uav
|
|
||||||
self.dc_energ[icrsh] -= J_hund / \
|
|
||||||
len(spn) * (Ncr[sp]) * (Ncr[sp] - 1.0)
|
|
||||||
mpi.report(
|
|
||||||
"DC for shell %(icrsh)i and block %(sp)s = %(Uav)f" % locals())
|
|
||||||
|
|
||||||
elif use_dc_formula == 1: # Held's formula, with U_interact the interorbital onsite interaction
|
|
||||||
|
|
||||||
self.dc_energ[icrsh] = (U_interact + (dim - 1) * (U_interact - 2.0 * J_hund) + (
|
|
||||||
dim - 1) * (U_interact - 3.0 * J_hund)) / (2 * dim - 1) / 2.0 * Ncrtot * (Ncrtot - 1.0)
|
|
||||||
for sp in spn:
|
|
||||||
Uav = (U_interact + (dim - 1) * (U_interact - 2.0 * J_hund) + (dim - 1)
|
|
||||||
* (U_interact - 3.0 * J_hund)) / (2 * dim - 1) * (Ncrtot - 0.5)
|
|
||||||
self.dc_imp[icrsh][sp] *= Uav
|
|
||||||
mpi.report(
|
|
||||||
"DC for shell %(icrsh)i and block %(sp)s = %(Uav)f" % locals())
|
|
||||||
|
|
||||||
elif use_dc_formula == 2: # AMF
|
|
||||||
|
|
||||||
self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot
|
|
||||||
for sp in spn:
|
|
||||||
Uav = U_interact * \
|
|
||||||
(Ncrtot - Ncr[sp] / dim) - \
|
|
||||||
J_hund * (Ncr[sp] - Ncr[sp] / dim)
|
|
||||||
self.dc_imp[icrsh][sp] *= Uav
|
|
||||||
self.dc_energ[
|
|
||||||
icrsh] -= (U_interact + (dim - 1) * J_hund) / dim / len(spn) * Ncr[sp] * Ncr[sp]
|
|
||||||
mpi.report(
|
|
||||||
"DC for shell %(icrsh)i and block %(sp)s = %(Uav)f" % locals())
|
|
||||||
|
|
||||||
mpi.report("DC energy for shell %s = %s" %
|
|
||||||
(icrsh, self.dc_energ[icrsh]))
|
|
||||||
|
|
||||||
else:
|
for sp in spn:
|
||||||
#For legacy compatibility
|
DC_val, E_val = compute_DC_from_density(N_tot=Ncrtot,U=U_interact, J=J_hund, n_orbitals=dim, N_spin=Ncr[sp], method=use_dc_formula)
|
||||||
match use_dc_formula:
|
self.dc_imp[icrsh][sp] *= DC_val
|
||||||
case 0:
|
self.dc_energ[icrsh] = E_val
|
||||||
mpi.report(f"Detected {use_dc_formula=}, changing to sFLL")
|
|
||||||
use_dc_formula = "sFLL"
|
|
||||||
case 1:
|
|
||||||
mpi.report(f"Detected {use_dc_formula=}, changing to cHeld")
|
|
||||||
use_dc_formula = "cHeld"
|
|
||||||
case 2:
|
|
||||||
mpi.report(f"Detected {use_dc_formula=}, changing to sAMF")
|
|
||||||
use_dc_formula = "sAMF"
|
|
||||||
|
|
||||||
for sp in spn:
|
|
||||||
DC_val, E_val = compute_DC_from_density(N_tot=Ncrtot,U=U_interact, J=J_hund, n_orbitals=dim, N_spin=Ncr[sp], method=use_dc_formula)
|
|
||||||
self.dc_imp[icrsh][sp] *= DC_val
|
|
||||||
self.dc_energ[icrsh] -= E_val
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -2048,7 +2010,6 @@ class SumkDFT(object):
|
|||||||
def calc_mu(self, precision=0.01, broadening=None, delta=0.5, max_loops=100, method="dichotomy"):
|
def calc_mu(self, precision=0.01, broadening=None, delta=0.5, max_loops=100, method="dichotomy"):
|
||||||
r"""
|
r"""
|
||||||
Searches for the chemical potential that gives the DFT total charge.
|
Searches for the chemical potential that gives the DFT total charge.
|
||||||
A simple bisection method is used.
|
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
@ -2061,7 +2022,7 @@ class SumkDFT(object):
|
|||||||
max_loops : int, optional
|
max_loops : int, optional
|
||||||
Number of dichotomy loops maximally performed.
|
Number of dichotomy loops maximally performed.
|
||||||
|
|
||||||
max_loops : string, optional
|
method : string, optional
|
||||||
Type of optimization used:
|
Type of optimization used:
|
||||||
* dichotomy: usual bisection algorithm from the TRIQS library
|
* dichotomy: usual bisection algorithm from the TRIQS library
|
||||||
* newton: newton method, faster convergence but more unstable
|
* newton: newton method, faster convergence but more unstable
|
||||||
@ -2116,52 +2077,50 @@ class SumkDFT(object):
|
|||||||
|
|
||||||
mpi.report("Trying out mu = {}".format(str(mu)))
|
mpi.report("Trying out mu = {}".format(str(mu)))
|
||||||
calc_dens = self.total_density(mu=mu, broadening=broadening).real - density
|
calc_dens = self.total_density(mu=mu, broadening=broadening).real - density
|
||||||
mpi.report("Delta to target density = {}".format(str(calc_dens)))
|
mpi.report(f"Target density = {density}; Delta to target = {calc_dens}")
|
||||||
return calc_dens
|
return calc_dens
|
||||||
|
|
||||||
#check for lowercase matching for the method variable
|
#check for lowercase matching for the method variable
|
||||||
match method.lower():
|
if method.lower()=="dichotomy":
|
||||||
|
mpi.report("\nSUMK calc_mu: Using dichtomy adjustment to find chemical potential\n")
|
||||||
|
self.chemical_potential = dichotomy.dichotomy(function=F_bisection,
|
||||||
|
x_init=self.chemical_potential, y_value=density,
|
||||||
|
precision_on_y=precision, delta_x=delta, max_loops=max_loops,
|
||||||
|
x_name="Chemical Potential", y_name="Total Density",
|
||||||
|
verbosity=3)[0]
|
||||||
|
elif method.lower()=="newton":
|
||||||
|
mpi.report("\nSUMK calc_mu: Using Newton method to find chemical potential\n")
|
||||||
|
self.chemical_potential = newton(func=F_optimize,
|
||||||
|
x0=self.chemical_potential,
|
||||||
|
tol=precision, maxiter=max_loops,
|
||||||
|
)
|
||||||
|
|
||||||
case "dichotomy":
|
elif method.lower()=="brent":
|
||||||
mpi.report("\nSUMK calc_mu: Using dichtomy adjustment to find chemical potential\n")
|
mpi.report("\nSUMK calc_mu: Using Brent method to find chemical potential")
|
||||||
self.chemical_potential = dichotomy.dichotomy(function=F_bisection,
|
mpi.report("SUMK calc_mu: Finding bounds \n")
|
||||||
x_init=self.chemical_potential, y_value=density,
|
|
||||||
precision_on_y=precision, delta_x=delta, max_loops=max_loops,
|
|
||||||
x_name="Chemical Potential", y_name="Total Density",
|
|
||||||
verbosity=3)[0]
|
|
||||||
case 'newton':
|
|
||||||
mpi.report("\nSUMK calc_mu: Using Newton method to find chemical potential\n")
|
|
||||||
self.chemical_potential = newton(func=F_optimize,
|
|
||||||
x0=self.chemical_potential,
|
|
||||||
tol=precision, maxiter=max_loops,
|
|
||||||
)
|
|
||||||
|
|
||||||
case 'brent':
|
|
||||||
mpi.report("\nSUMK calc_mu: Using Brent method to find chemical potential")
|
|
||||||
mpi.report("SUMK calc_mu: Finding bounds \n")
|
|
||||||
|
|
||||||
mu_guess_0, mu_guess_1 = find_bounds(function=F_optimize,
|
|
||||||
x_init=self.chemical_potential,
|
|
||||||
delta_x=delta, max_loops=max_loops,
|
|
||||||
)
|
|
||||||
mu_guess_1 += 0.01 #scrambles higher lying interval to avoid getting stuck
|
|
||||||
mpi.report("\nSUMK calc_mu: Searching root with Brent method\n")
|
|
||||||
self.chemical_potential = brenth(f=F_optimize,
|
|
||||||
a=mu_guess_0,
|
|
||||||
b=mu_guess_1,
|
|
||||||
xtol=precision, maxiter=max_loops,
|
|
||||||
)
|
|
||||||
|
|
||||||
case _:
|
mu_guess_0, mu_guess_1 = find_bounds(function=F_optimize,
|
||||||
raise ValueError(
|
x_init=self.chemical_potential,
|
||||||
f"SUMK calc_mu: The method selected: {method}, is not implemented\n",
|
delta_x=delta, max_loops=max_loops,
|
||||||
"""
|
)
|
||||||
Please check for typos or select one of the following:
|
mu_guess_1 += 0.01 #scrambles higher lying interval to avoid getting stuck
|
||||||
* dichotomy: usual bisection algorithm from the TRIQS library
|
mpi.report("\nSUMK calc_mu: Searching root with Brent method\n")
|
||||||
* newton: newton method, fastest convergence but more unstable
|
self.chemical_potential = brenth(f=F_optimize,
|
||||||
* brent: finds bounds and proceeds with hyperbolic brent method, a compromise between speed and ensuring convergence
|
a=mu_guess_0,
|
||||||
"""
|
b=mu_guess_1,
|
||||||
)
|
xtol=precision, maxiter=max_loops,
|
||||||
|
)
|
||||||
|
|
||||||
|
else:
|
||||||
|
raise ValueError(
|
||||||
|
f"SUMK calc_mu: The selected method: {method}, is not implemented\n",
|
||||||
|
"""
|
||||||
|
Please check for typos or select one of the following:
|
||||||
|
* dichotomy: usual bisection algorithm from the TRIQS library
|
||||||
|
* newton: newton method, fastest convergence but more unstable
|
||||||
|
* brent: finds bounds and proceeds with hyperbolic brent method, a compromise between speed and ensuring convergence
|
||||||
|
"""
|
||||||
|
)
|
||||||
|
|
||||||
return self.chemical_potential
|
return self.chemical_potential
|
||||||
|
|
||||||
@ -2414,7 +2373,8 @@ class SumkDFT(object):
|
|||||||
f.write("\n")
|
f.write("\n")
|
||||||
|
|
||||||
elif dm_type == 'qe':
|
elif dm_type == 'qe':
|
||||||
assert self.SP == 0, "Spin-polarized density matrix is not implemented"
|
if self.SP == 0:
|
||||||
|
mpi.report("SUMK calc_density_correction: WARNING! Averaging out spin-polarized correction in the density channel")
|
||||||
|
|
||||||
subgrp = 'dft_update'
|
subgrp = 'dft_update'
|
||||||
delta_N = np.zeros([self.n_k, max(self.n_orbitals[:,0]), max(self.n_orbitals[:,0])], dtype=complex)
|
delta_N = np.zeros([self.n_k, max(self.n_orbitals[:,0]), max(self.n_orbitals[:,0])], dtype=complex)
|
||||||
|
BIN
test/python/NiO.ref.h5
Executable file
BIN
test/python/NiO.ref.h5
Executable file
Binary file not shown.
BIN
test/python/NiO.ref.h5_copy
Executable file
BIN
test/python/NiO.ref.h5_copy
Executable file
Binary file not shown.
Binary file not shown.
BIN
test/python/SrVO3.ref.h5_copy
Normal file
BIN
test/python/SrVO3.ref.h5_copy
Normal file
Binary file not shown.
168
test/python/dc_test.py
Executable file
168
test/python/dc_test.py
Executable file
@ -0,0 +1,168 @@
|
|||||||
|
##########################################################################
|
||||||
|
#
|
||||||
|
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
||||||
|
#
|
||||||
|
# Copyright (C) 2023 by A. Carta
|
||||||
|
#
|
||||||
|
# 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
|
||||||
|
from triqs_dft_tools.sumk_dft import *
|
||||||
|
from triqs.gf import *
|
||||||
|
from h5 import HDFArchive
|
||||||
|
from triqs.operators.util import *
|
||||||
|
import triqs.utility.mpi as mpi
|
||||||
|
|
||||||
|
|
||||||
|
Uval = 5
|
||||||
|
Jval = 0.3
|
||||||
|
|
||||||
|
method_dict = {
|
||||||
|
"FLL":{
|
||||||
|
"numbering_convention":0,
|
||||||
|
"new_convention":"cFLL"
|
||||||
|
},
|
||||||
|
"AMF":{
|
||||||
|
"numbering_convention":2,
|
||||||
|
"new_convention":"cAMF"
|
||||||
|
},
|
||||||
|
"Held":{
|
||||||
|
"numbering_convention":1,
|
||||||
|
"new_convention":"cHeld"
|
||||||
|
},
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
def test_dc(method, method_dict, dens, Uval, Jval, filename):
|
||||||
|
|
||||||
|
dc_no = method_dict[method]["numbering_convention"]
|
||||||
|
dc_string = method_dict[method]["new_convention"]
|
||||||
|
|
||||||
|
mpi.report("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX")
|
||||||
|
mpi.report(f"\n Testing interface {method} \n")
|
||||||
|
mpi.report("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX")
|
||||||
|
|
||||||
|
mpi.report("\nTesting legacy compatibility layer:\n")
|
||||||
|
SK_compat.calc_dc(dens_mat=dens[0], U_interact=Uval, J_hund=Jval, use_dc_formula=dc_no )
|
||||||
|
mpi.report("Up DC matrix:")
|
||||||
|
mpi.report(SK_compat.dc_imp[0]['up'])
|
||||||
|
mpi.report(f"Double counting energy = {SK_compat.dc_energ} ")
|
||||||
|
|
||||||
|
mpi.report("\nTesting new dc interface:\n")
|
||||||
|
SK_new.calc_dc(dens_mat=dens[0], U_interact=Uval, J_hund=Jval, use_dc_formula=dc_string)
|
||||||
|
mpi.report("Up DC matrix:")
|
||||||
|
mpi.report(SK_new.dc_imp[0]['up'])
|
||||||
|
mpi.report(f"Double counting energy = {SK_new.dc_energ} ")
|
||||||
|
|
||||||
|
# Load previously computed DC values from h5 archive
|
||||||
|
R = HDFArchive(f'./{filename}', 'r')
|
||||||
|
dc_comp = R[f'DC_{method}_benchmark']['dc_imp']
|
||||||
|
en_comp = R[f'DC_{method}_benchmark']['dc_energ']
|
||||||
|
del R
|
||||||
|
|
||||||
|
mpi.report(f"\nAsserting comparison for method: {method}")
|
||||||
|
assert np.allclose(SK_compat.dc_imp[0]['up'], dc_comp, atol=1e-12), f"Assertion failed comparing legacy Vdc to reference, method: {method} "
|
||||||
|
assert np.allclose(SK_compat.dc_energ, en_comp, atol=1e-12), f"Assertion failed comparing legacy energy to reference, method {method} "
|
||||||
|
assert np.allclose(SK_new.dc_imp[0]['up'], dc_comp, atol=1e-12), f"Assertion failed comparing Vdc to reference, method: {method} "
|
||||||
|
assert np.allclose(SK_new.dc_energ, en_comp, atol=1e-12), f"Assertion failed comparing energy to reference, method: {method} "
|
||||||
|
mpi.report("Comparison with stored DC values successfull!\n")
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
# %% 5 orbitals testing
|
||||||
|
|
||||||
|
mpi.report("\n############################################")
|
||||||
|
mpi.report("############################################")
|
||||||
|
mpi.report(f"\n \n Starting tests for 5 orbitals \n \n")
|
||||||
|
mpi.report("############################################")
|
||||||
|
mpi.report("############################################")
|
||||||
|
dft_filename = "./NiO.ref"
|
||||||
|
use_blocks = False
|
||||||
|
SK_compat = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks)
|
||||||
|
SK_compat.set_mu(13.9)
|
||||||
|
SK_new = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks)
|
||||||
|
SK_new.set_mu(13.9)
|
||||||
|
|
||||||
|
|
||||||
|
icrsh = 0
|
||||||
|
dens = SK_compat.density_matrix()
|
||||||
|
|
||||||
|
with np.printoptions(precision=5):
|
||||||
|
for key in dens[0].keys():
|
||||||
|
mpi.report(f"{key} channel")
|
||||||
|
mpi.report(dens[0][key].real)
|
||||||
|
|
||||||
|
N_up = np.trace(dens[0]['up'].real)
|
||||||
|
N_down = np.trace(dens[0]['down'].real)
|
||||||
|
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")
|
||||||
|
|
||||||
|
#in case implementation changes, to write new testing data into archive
|
||||||
|
#R = HDFArchive('./NiO.ref.h5', 'a')
|
||||||
|
#R.create_group(f'DC_{method}_benchmark')
|
||||||
|
#R[f'DC_{method}_benchmark']['dc_imp']= SK_new.dc_imp[0]['up']
|
||||||
|
#R[f'DC_{method}_benchmark']['dc_energ']= SK_new.dc_energ
|
||||||
|
#del R
|
||||||
|
|
||||||
|
|
||||||
|
# 3 orbital testing
|
||||||
|
|
||||||
|
mpi.report("############################################")
|
||||||
|
mpi.report("############################################")
|
||||||
|
mpi.report(f"\n \n Starting tests for 3 orbitals \n \n")
|
||||||
|
mpi.report("############################################")
|
||||||
|
mpi.report("############################################")
|
||||||
|
dft_filename = "./SrVO3.ref"
|
||||||
|
use_blocks = False
|
||||||
|
SK_compat = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks)
|
||||||
|
SK_new = SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks)
|
||||||
|
|
||||||
|
icrsh = 0
|
||||||
|
dens = SK_compat.density_matrix()
|
||||||
|
|
||||||
|
|
||||||
|
with np.printoptions(precision=5):
|
||||||
|
for key in dens[0].keys():
|
||||||
|
mpi.report(f"{key} channel")
|
||||||
|
mpi.report(dens[0][key].real)
|
||||||
|
|
||||||
|
N_up = np.trace(dens[0]['up'].real)
|
||||||
|
N_down = np.trace(dens[0]['down'].real)
|
||||||
|
N_tot = N_up + N_down
|
||||||
|
|
||||||
|
mpi.report(f"{N_up=} ,{N_down=}, {N_tot=}\n")
|
||||||
|
|
||||||
|
Uval = 5
|
||||||
|
Jval = 0.3
|
||||||
|
|
||||||
|
for method in ["FLL", "AMF", "Held"]:
|
||||||
|
test_dc(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')
|
||||||
|
#R.create_group(f'DC_{method}_benchmark')
|
||||||
|
#R[f'DC_{method}_benchmark']['dc_imp']= SK_new.dc_imp[0]['up']
|
||||||
|
#R[f'DC_{method}_benchmark']['dc_energ']= SK_new.dc_energ
|
||||||
|
#del R
|
Loading…
Reference in New Issue
Block a user