Skip to content

Commit

Permalink
[feat] improved standard behavior of block struct
Browse files Browse the repository at this point in the history 
* previously the default gf_struct_solver had keys up / down,
inconsistent with the default behavior after analyse_block_structure was
run: up_0 / down_0. Now the default solver structure always has the _0
in the key.
* old behavior resulted in error when analyse_block_structure was called
twice
* fixed analyse block structure tests with new changes
* to correctly use analyse_block_structure use now
extract_G_loc(transform_to_solver_blocks=False)
* changed density_matrix function to use directly extract_G_loc() if
using_gf is selected.
* print deprecation warning in density_matrix, same can be achieved via
extract_G_loc and [G.density() for G in Gloc]
* new function density_matrix_using_point_integration()
* fixed other tests accordingly
* fixed small bug in initial block structure regarding length of lists
  • Loading branch information
@the-hampel
the-hampel committed Feb 26, 2024
1 parent d4d2317 commit a485adc
Show file tree
Hide file tree
Showing 5 changed files with 166 additions and 140 deletions.
138 changes: 79 additions & 59 deletions python/triqs_dft_tools/sumk_dft.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,8 +168,8 @@ def __init__(self, hdf_file, h_field=0.0, mesh=None, beta=40, n_iw=1025, use_dft
# blocks possible
self.gf_struct_sumk = [[(sp, self.corr_shells[icrsh]['dim']) for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]]
for icrsh in range(self.n_corr_shells)]
# First set a standard gf_struct solver:
self.gf_struct_solver = [dict([(sp, self.corr_shells[self.inequiv_to_corr[ish]]['dim'])
# First set a standard gf_struct solver (add _0 here for consistency with analyse_block_structure):
self.gf_struct_solver = [dict([(sp+'_0', self.corr_shells[self.inequiv_to_corr[ish]]['dim'])
for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]])
for ish in range(self.n_inequiv_shells)]
# Set standard (identity) maps from gf_struct_sumk <->
Expand All @@ -180,12 +180,12 @@ def __init__(self, hdf_file, h_field=0.0, mesh=None, beta=40, n_iw=1025, use_dft
for ish in range(self.n_inequiv_shells)]
for ish in range(self.n_inequiv_shells):
for block, inner_dim in self.gf_struct_sumk[self.inequiv_to_corr[ish]]:
self.solver_to_sumk_block[ish][block] = block
self.solver_to_sumk_block[ish][block+'_0'] = block
for inner in range(inner_dim):
self.sumk_to_solver[ish][
(block, inner)] = (block, inner)
(block, inner)] = (block+'_0', inner)
self.solver_to_sumk[ish][
(block, inner)] = (block, inner)
(block+'_0', inner)] = (block, inner)
# assume no shells are degenerate
self.deg_shells = [[] for ish in range(self.n_inequiv_shells)]

Expand Down Expand Up @@ -1065,13 +1065,6 @@ def analyse_block_structure_from_gf(self, G, threshold=1.e-5, include_shells=Non

gf = self._get_hermitian_quantity_from_gf(G)

# initialize the variables
self.gf_struct_solver = [{} for ish in range(self.n_inequiv_shells)]
self.sumk_to_solver = [{} for ish in range(self.n_inequiv_shells)]
self.solver_to_sumk = [{} for ish in range(self.n_inequiv_shells)]
self.solver_to_sumk_block = [{}
for ish in range(self.n_inequiv_shells)]

# the maximum value of each matrix element of each block and shell
max_gf = [{name:np.max(np.abs(g.data),0) for name, g in gf[ish]} for ish in range(self.n_inequiv_shells)]

Expand All @@ -1080,6 +1073,10 @@ def analyse_block_structure_from_gf(self, G, threshold=1.e-5, include_shells=Non
include_shells = list(range(self.n_inequiv_shells))

for ish in include_shells:
self.gf_struct_solver[ish] = {}
self.sumk_to_solver[ish] = {}
self.solver_to_sumk[ish] = {}
self.solver_to_sumk_block[ish] = {}
for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim']

Expand Down Expand Up @@ -1455,21 +1452,12 @@ def calculate_diagonalization_matrix(self, prop_to_be_diagonal='eal', calc_in_so

return trans


def density_matrix(self, method='using_gf'):
"""Calculate density matrices in one of two ways.
Parameters
----------
method : string, optional
- if 'using_gf': First get lattice gf (g_loc is not set up), then density matrix.
It is useful for Hubbard I, and very quick.
No assumption on the hopping structure is made (ie diagonal or not).
- if 'using_point_integration': Only works for diagonal hopping matrix (true in wien2k).
beta : float, optional
Inverse temperature.
def density_matrix_using_point_integration(self):
"""
Calculate density matrices using point integration: Only works for
diagonal hopping matrix (true in wien2k). Consider using extract_G_loc
together with [G.density() for G in Gloc] instead. Returned density
matrix is always given in SumK block structure.
Returns
-------
Expand All @@ -1483,34 +1471,21 @@ def density_matrix(self, method='using_gf'):
[self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], complex)

ikarray = np.array(list(range(self.n_k)))
ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO]
for ik in mpi.slice_array(ikarray):

if method == "using_gf":

G_latt_iw = self.lattice_gf(ik=ik, mu=self.chemical_potential)
G_latt_iw *= self.bz_weights[ik]
dm = G_latt_iw.density()
MMat = [dm[sp] for sp in self.spin_block_names[self.SO]]

elif method == "using_point_integration":

ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO]
dims = {sp:self.n_orbitals[ik, ntoi[sp]] for sp in spn}
MMat = [np.zeros([dims[sp], dims[sp]], complex) for sp in spn]

for isp, sp in enumerate(spn):
ind = ntoi[sp]
for inu in range(self.n_orbitals[ik, ind]):
# only works for diagonal hopping matrix (true in
# wien2k)
if (self.hopping[ik, ind, inu, inu] - self.h_field * (1 - 2 * isp)) < 0.0:
MMat[isp][inu, inu] = 1.0
else:
MMat[isp][inu, inu] = 0.0

else:
raise ValueError("density_matrix: the method '%s' is not supported." % method)
dims = {sp:self.n_orbitals[ik, ntoi[sp]] for sp in spn}
MMat = [np.zeros([dims[sp], dims[sp]], complex) for sp in spn]

for isp, sp in enumerate(spn):
ind = ntoi[sp]
for inu in range(self.n_orbitals[ik, ind]):
# only works for diagonal hopping matrix (true in
# wien2k)
if (self.hopping[ik, ind, inu, inu] - self.h_field * (1 - 2 * isp)) < 0.0:
MMat[isp][inu, inu] = 1.0
else:
MMat[isp][inu, inu] = 0.0

for icrsh in range(self.n_corr_shells):
for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
Expand All @@ -1519,11 +1494,7 @@ def density_matrix(self, method='using_gf'):
dim = self.corr_shells[icrsh]['dim']
n_orb = self.n_orbitals[ik, ind]
projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb]
if method == "using_gf":
dens_mat[icrsh][sp] += np.dot(np.dot(projmat, MMat[isp]),
projmat.transpose().conjugate())
elif method == "using_point_integration":
dens_mat[icrsh][sp] += self.bz_weights[ik] * np.dot(np.dot(projmat, MMat[isp]),
dens_mat[icrsh][sp] += self.bz_weights[ik] * np.dot(np.dot(projmat, MMat[isp]),
projmat.transpose().conjugate())

# get data from nodes:
Expand All @@ -1546,6 +1517,55 @@ def density_matrix(self, method='using_gf'):

return dens_mat


def density_matrix(self, method='using_gf', mu=None, with_Sigma=True, with_dc=True, broadening=None,
transform_to_solver_blocks=True, show_warnings=True):
"""Calculate density matrices in one of two ways.
Parameters
----------
method : string, optional
- if 'using_gf': First get lattice gf (g_loc is not set up), then density matrix.
It is useful for Hubbard I, and very quick.
No assumption on the hopping structure is made (ie diagonal or not).
- if 'using_point_integration': Only works for diagonal hopping matrix (true in wien2k).
mu : real, optional
Input chemical potential. If not provided the value of self.chemical_potential is used as mu.
with_Sigma : boolean, optional
If True then the local GF is calculated with the self-energy self.Sigma_imp.
with_dc : boolean, optional
If True then the double-counting correction is subtracted from the self-energy in calculating the GF.
broadening : float, optional
Imaginary shift for the axis along which the real-axis GF is calculated.
If not provided, broadening will be set to double of the distance between mesh points in 'mesh'.
Only relevant for real-frequency GF.
transform_to_solver_blocks : bool, optional
If True (default), the returned G_loc will be transformed to the block structure ``gf_struct_solver``,
else it will be in ``gf_struct_sumk``.
show_warnings : bool, optional
Displays warning messages during transformation
(Only effective if transform_to_solver_blocks = True
Returns
-------
dens_mat : list of dicts
Density matrix for each spin in each correlated shell.
"""

if method == "using_gf":
warn("WARNING: density_matrix: method 'using_gf' is deprecated. Use 'extract_G_loc' instead.")
Gloc = self.extract_G_loc(mu, with_Sigma, with_dc, broadening,
transform_to_solver_blocks, show_warnings)
dens_mat = [G.density() for G in Gloc]
elif method == "using_point_integration":
warn("WARNING: density_matrix: method 'using_point_integration' is deprecated. Use 'density_matrix_using_point_integration' instead. All additionally provided arguments are ignored.")
dens_mat = self.density_matrix_using_point_integration()
else:
raise ValueError("density_matrix: the method '%s' is not supported." % method)

return dens_mat

# For simple dft input, get crystal field splittings.
def eff_atomic_levels(self):
r"""
Expand Down
Loading

0 comments on commit a485adc

Please sign in to comment.