LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2025-05-01 12:44:59 +02:00
Code Issues Releases Wiki Activity

revert auto formatting from prev commits

Browse Source
This commit is contained in:
the-hampel 2025-01-09 14:11:50 +01:00
parent de701d17af
commit 2b2b86846c
1 changed files with 174 additions and 241 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

415
python/triqs_dft_tools/sumk_dft.py
View File

@ -1,3 +1,4 @@
########################################################################## ##########################################################################
# #
# TRIQS: a Toolbox for Research in Interacting Quantum Systems # TRIQS: a Toolbox for Research in Interacting Quantum Systems
@ -23,7 +24,6 @@
General SumK class and helper functions for combining ab-initio code and triqs General SumK class and helper functions for combining ab-initio code and triqs
""" """
import os
from types import * from types import *
import numpy as np import numpy as np
import triqs.utility.dichotomy as dichotomy import triqs.utility.dichotomy as dichotomy
@ -47,7 +47,7 @@ class SumkDFT(object):
def __init__(self, hdf_file, h_field=0.0, mesh=None, beta=40, n_iw=1025, use_dft_blocks=False, def __init__(self, hdf_file, h_field=0.0, mesh=None, beta=40, n_iw=1025, use_dft_blocks=False,
dft_data='dft_input', symmcorr_data='dft_symmcorr_input', parproj_data='dft_parproj_input', dft_data='dft_input', symmcorr_data='dft_symmcorr_input', parproj_data='dft_parproj_input',
symmpar_data='dft_symmpar_input', bands_data='dft_bands_input', transp_data='dft_transp_input', symmpar_data='dft_symmpar_input', bands_data='dft_bands_input', transp_data='dft_transp_input',
misc_data='dft_misc_input', bc_data='dft_bandchar_input', cont_data='dft_contours_input'): misc_data='dft_misc_input',bc_data='dft_bandchar_input',cont_data='dft_contours_input'):
r""" r"""
Initialises the class from data previously stored into an hdf5 archive. Initialises the class from data previously stored into an hdf5 archive.
@ -121,33 +121,26 @@ class SumkDFT(object):
self.block_structure = BlockStructure() self.block_structure = BlockStructure()
# Read input from HDF: # Read input from HDF:
req_things_to_read = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', req_things_to_read = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required',
'density_required', 'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat',
'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'bz_weights', 'hopping',
'rot_mat', 'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr']
'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'bz_weights',
'hopping',
'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr']
self.subgroup_present, self.values_not_read = self.read_input_from_hdf( self.subgroup_present, self.values_not_read = self.read_input_from_hdf(
subgrp=self.dft_data, things_to_read=req_things_to_read) subgrp=self.dft_data, things_to_read=req_things_to_read)
# test if all required properties have been found # test if all required properties have been found
if len(self.values_not_read) > 0 and mpi.is_master_node: if len(self.values_not_read) > 0 and mpi.is_master_node:
raise ValueError( raise ValueError('ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:', self.values_not_read)
'ERROR: One or more necessary SumK input properties have not been found in the given h5 archive:',
self.values_not_read)
# optional properties to load # optional properties to load
# soon bz_weights is depraced and replaced by kpt_weights, kpts_basis and kpts will become required to read soon # soon bz_weights is depraced and replaced by kpt_weights, kpts_basis and kpts will become required to read soon
optional_things_to_read = ['proj_mat_csc', 'proj_or_hk', 'kpt_basis', 'kpts', 'kpt_weights', 'dft_code'] optional_things_to_read = ['proj_mat_csc', 'proj_or_hk', 'kpt_basis', 'kpts', 'kpt_weights', 'dft_code']
subgroup_present, self.optional_values_not_read = self.read_input_from_hdf(subgrp=self.dft_data, subgroup_present, self.optional_values_not_read = self.read_input_from_hdf(subgrp=self.dft_data, things_to_read=optional_things_to_read)
things_to_read=optional_things_to_read)
# warning if dft_code was not read (old h5 structure) # warning if dft_code was not read (old h5 structure)
if 'dft_code' in self.optional_values_not_read: if 'dft_code' in self.optional_values_not_read:
self.dft_code = None self.dft_code = None
mpi.report( mpi.report('\nWarning: old h5 archive without dft_code input flag detected. Please specify sumk.dft_code manually!\n')
'\nWarning: old h5 archive without dft_code input flag detected. Please specify sumk.dft_code manually!\n')
if self.symm_op: if self.symm_op:
self.symmcorr = Symmetry(hdf_file, subgroup=self.symmcorr_data) self.symmcorr = Symmetry(hdf_file, subgroup=self.symmcorr_data)
@ -170,13 +163,11 @@ class SumkDFT(object):
# GF structure used for the local things in the k sums # GF structure used for the local things in the k sums
# Most general form allowing for all hybridisation, i.e. largest # Most general form allowing for all hybridisation, i.e. largest
# blocks possible # blocks possible
self.gf_struct_sumk = [ self.gf_struct_sumk = [[(sp, self.corr_shells[icrsh]['dim']) for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]]
[(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)]
for icrsh in range(self.n_corr_shells)]
# First set a standard gf_struct solver (add _0 here for consistency with analyse_block_structure): # 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']) self.gf_struct_solver = [dict([(sp+'_0', self.corr_shells[self.inequiv_to_corr[ish]]['dim'])
for sp in for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]])
self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]])
for ish in range(self.n_inequiv_shells)] for ish in range(self.n_inequiv_shells)]
# Set standard (identity) maps from gf_struct_sumk <-> # Set standard (identity) maps from gf_struct_sumk <->
# gf_struct_solver # gf_struct_solver
@ -186,12 +177,12 @@ class SumkDFT(object):
for ish in range(self.n_inequiv_shells)] for ish in range(self.n_inequiv_shells)]
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]]: for block, inner_dim in self.gf_struct_sumk[self.inequiv_to_corr[ish]]:
self.solver_to_sumk_block[ish][block + '_0'] = block self.solver_to_sumk_block[ish][block+'_0'] = block
for inner in range(inner_dim): for inner in range(inner_dim):
self.sumk_to_solver[ish][ self.sumk_to_solver[ish][
(block, inner)] = (block + '_0', inner) (block, inner)] = (block+'_0', inner)
self.solver_to_sumk[ish][ self.solver_to_sumk[ish][
(block + '_0', inner)] = (block, inner) (block+'_0', inner)] = (block, inner)
# assume no shells are degenerate # assume no shells are degenerate
self.deg_shells = [[] for ish in range(self.n_inequiv_shells)] self.deg_shells = [[] for ish in range(self.n_inequiv_shells)]
@ -213,9 +204,9 @@ class SumkDFT(object):
self.min_band_energy = None self.min_band_energy = None
self.max_band_energy = None self.max_band_energy = None
################ ################
# hdf5 FUNCTIONS # hdf5 FUNCTIONS
################ ################
def read_input_from_hdf(self, subgrp, things_to_read): def read_input_from_hdf(self, subgrp, things_to_read):
r""" r"""
@ -286,8 +277,8 @@ class SumkDFT(object):
with HDFArchive(self.hdf_file, 'a') as ar: with HDFArchive(self.hdf_file, 'a') as ar:
if not subgrp in ar: ar.create_group(subgrp) if not subgrp in ar: ar.create_group(subgrp)
for it in things_to_save: for it in things_to_save:
if it in ["gf_struct_sumk", "gf_struct_solver", if it in [ "gf_struct_sumk", "gf_struct_solver",
"solver_to_sumk", "sumk_to_solver", "solver_to_sumk_block"]: "solver_to_sumk", "sumk_to_solver", "solver_to_sumk_block"]:
warn("It is not recommended to save '{}' individually. Save 'block_structure' instead.".format(it)) warn("It is not recommended to save '{}' individually. Save 'block_structure' instead.".format(it))
try: try:
ar[subgrp][it] = getattr(self, it) ar[subgrp][it] = getattr(self, it)
@ -324,9 +315,9 @@ class SumkDFT(object):
raise ValueError("load: %s not found, and so not loaded." % it) raise ValueError("load: %s not found, and so not loaded." % it)
return list_to_return return list_to_return
################ ################
# CORE FUNCTIONS # CORE FUNCTIONS
################ ################
def downfold(self, ik, ish, bname, gf_to_downfold, gf_inp, shells='corr', ir=None): def downfold(self, ik, ish, bname, gf_to_downfold, gf_inp, shells='corr', ir=None):
r""" r"""
@ -486,8 +477,7 @@ class SumkDFT(object):
gf_rotated.from_L_G_R(rot_mat[ish].conjugate( gf_rotated.from_L_G_R(rot_mat[ish].conjugate(
), gf_rotated, rot_mat[ish].transpose()) ), gf_rotated, rot_mat[ish].transpose())
else: else:
gf_rotated.from_L_G_R(rot_mat[ish], gf_rotated, rot_mat[ gf_rotated.from_L_G_R(rot_mat[ish], gf_rotated, rot_mat[ish].conjugate().transpose())
ish].conjugate().transpose())
elif direction == 'toLocal': elif direction == 'toLocal':
@ -543,15 +533,14 @@ class SumkDFT(object):
broadening = 2.0 * ((mesh.w_max - mesh.w_min) / (len(mesh) - 1)) broadening = 2.0 * ((mesh.w_max - mesh.w_min) / (len(mesh) - 1))
# Check if G_latt is present # Check if G_latt is present
set_up_G_latt = False # Assume not set_up_G_latt = False # Assume not
if not hasattr(self, "G_latt"): if not hasattr(self, "G_latt" ):
# Need to create G_latt_(i)w # Need to create G_latt_(i)w
set_up_G_latt = True set_up_G_latt = True
else: # Check that existing GF is consistent else: # Check that existing GF is consistent
G_latt = self.G_latt G_latt = self.G_latt
GFsize = [gf.target_shape[0] for bname, gf in G_latt] GFsize = [gf.target_shape[0] for bname, gf in G_latt]
unchangedsize = all([self.n_orbitals[ik, ntoi[spn[isp]]] == GFsize[ unchangedsize = all([self.n_orbitals[ik, ntoi[spn[isp]]] == GFsize[isp] for isp in range(self.n_spin_blocks[self.SO])])
isp] for isp in range(self.n_spin_blocks[self.SO])])
if (not mesh is None) or (not unchangedsize): if (not mesh is None) or (not unchangedsize):
set_up_G_latt = True set_up_G_latt = True
@ -578,8 +567,7 @@ class SumkDFT(object):
mesh_values = self.mesh_values mesh_values = self.mesh_values
elif not mesh is None: elif not mesh is None:
assert isinstance(mesh, assert isinstance(mesh, (MeshReFreq, MeshDLRImFreq, MeshImFreq)), "mesh must be a triqs MeshReFreq or MeshImFreq"
(MeshReFreq, MeshDLRImFreq, MeshImFreq)), "mesh must be a triqs MeshReFreq or MeshImFreq"
if isinstance(mesh, MeshImFreq): if isinstance(mesh, MeshImFreq):
mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh)) mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh))
elif isinstance(mesh, MeshDLRImFreq): elif isinstance(mesh, MeshDLRImFreq):
@ -597,7 +585,7 @@ class SumkDFT(object):
gf_struct = [(spn[isp], block_structure[isp]) gf_struct = [(spn[isp], block_structure[isp])
for isp in range(self.n_spin_blocks[self.SO])] for isp in range(self.n_spin_blocks[self.SO])]
block_ind_list = [block for block, inner in gf_struct] block_ind_list = [block for block, inner in gf_struct]
glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner), len(inner)]) glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner),len(inner)])
for block, inner in gf_struct] for block, inner in gf_struct]
G_latt = BlockGf(name_list=block_ind_list, G_latt = BlockGf(name_list=block_ind_list,
block_list=glist(), make_copies=False) block_list=glist(), make_copies=False)
@ -611,11 +599,11 @@ class SumkDFT(object):
ind = ntoi[spn[ibl]] ind = ntoi[spn[ibl]]
n_orb = self.n_orbitals[ik, ind] n_orb = self.n_orbitals[ik, ind]
if isinstance(mesh, (MeshImFreq, MeshDLRImFreq)): if isinstance(mesh, (MeshImFreq, MeshDLRImFreq)):
gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field * (1 - 2 * ibl)) gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl))
- self.hopping[ik, ind, 0:n_orb, 0:n_orb]) - self.hopping[ik, ind, 0:n_orb, 0:n_orb])
else: else:
gf.data[:, :, :] = (idmat[ibl] * gf.data[:, :, :] = (idmat[ibl] *
(mesh_values[:, None, None] + mu + self.h_field * (1 - 2 * ibl) + 1j * broadening) (mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl) + 1j*broadening)
- self.hopping[ik, ind, 0:n_orb, 0:n_orb]) - self.hopping[ik, ind, 0:n_orb, 0:n_orb])
if with_Sigma: if with_Sigma:
@ -649,26 +637,23 @@ class SumkDFT(object):
if transform_to_sumk_blocks: if transform_to_sumk_blocks:
Sigma_imp = self.transform_to_sumk_blocks(Sigma_imp) Sigma_imp = self.transform_to_sumk_blocks(Sigma_imp)
assert isinstance(Sigma_imp, list), \ assert isinstance(Sigma_imp, list),\
"put_Sigma: Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!" "put_Sigma: Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!"
assert len(Sigma_imp) == self.n_corr_shells, \ assert len(Sigma_imp) == self.n_corr_shells,\
"put_Sigma: give exactly one Sigma for each corr. shell!" "put_Sigma: give exactly one Sigma for each corr. shell!"
if (isinstance(self.mesh, (MeshImFreq, MeshDLRImFreq)) and if (isinstance(self.mesh, (MeshImFreq, MeshDLRImFreq)) and
all(isinstance(gf.mesh, (MeshImFreq, MeshDLRImFreq)) and all(isinstance(gf.mesh, (MeshImFreq, MeshDLRImFreq)) and
isinstance(gf, Gf) and isinstance(gf, Gf) and
gf.mesh == self.mesh for bname, gf in Sigma_imp[0])): gf.mesh == self.mesh for bname, gf in Sigma_imp[0])):
# Imaginary frequency Sigma: # Imaginary frequency Sigma:
self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk') self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk')
for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells)]
SK_Sigma_imp = self.Sigma_imp SK_Sigma_imp = self.Sigma_imp
elif isinstance(self.mesh, MeshReFreq) and all( elif isinstance(self.mesh, MeshReFreq) and all(isinstance(gf, Gf) and isinstance(gf.mesh, MeshReFreq) and gf.mesh == self.mesh for bname, gf in Sigma_imp[0]):
isinstance(gf, Gf) and isinstance(gf.mesh, MeshReFreq) and gf.mesh == self.mesh for bname, gf in
Sigma_imp[0]):
# Real frequency Sigma: # Real frequency Sigma:
self.Sigma_imp = [ self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, gf_function=Gf, space='sumk')
self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, gf_function=Gf, space='sumk') for icrsh in range(self.n_corr_shells)]
for icrsh in range(self.n_corr_shells)]
SK_Sigma_imp = self.Sigma_imp SK_Sigma_imp = self.Sigma_imp
else: else:
@ -684,16 +669,13 @@ class SumkDFT(object):
else: else:
gf << Sigma_imp[icrsh][bname] gf << Sigma_imp[icrsh][bname]
# warning if real frequency self energy is within the bounds of the band energies #warning if real frequency self energy is within the bounds of the band energies
if isinstance(self.mesh, MeshReFreq): if isinstance(self.mesh, MeshReFreq):
if self.min_band_energy is None or self.max_band_energy is None: if self.min_band_energy is None or self.max_band_energy is None:
self.calculate_min_max_band_energies() self.calculate_min_max_band_energies()
mesh = np.linspace(self.mesh.w_min, self.mesh.w_max, len(self.mesh)) mesh = np.linspace(self.mesh.w_min, self.mesh.w_max, len(self.mesh))
if mesh[0] > (self.min_band_energy - self.chemical_potential) or mesh[-1] < ( if mesh[0] > (self.min_band_energy - self.chemical_potential) or mesh[-1] < (self.max_band_energy - self.chemical_potential):
self.max_band_energy - self.chemical_potential): warn('The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f'%(mesh[0], mesh[-1], self.min_band_energy, self.max_band_energy))
warn(
'The given Sigma is on a mesh which does not cover the band energy range. The Sigma MeshReFreq runs from %f to %f, while the band energy (minus the chemical potential) runs from %f to %f' % (
mesh[0], mesh[-1], self.min_band_energy, self.max_band_energy))
def transform_to_sumk_blocks(self, Sigma_imp, Sigma_out=None): def transform_to_sumk_blocks(self, Sigma_imp, Sigma_out=None):
r""" transform Sigma from solver to sumk space r""" transform Sigma from solver to sumk space
@ -708,14 +690,13 @@ class SumkDFT(object):
according to ``gf_struct_sumk``; if None, it will be created according to ``gf_struct_sumk``; if None, it will be created
""" """
assert isinstance(Sigma_imp, list), \ assert isinstance(Sigma_imp, list),\
"transform_to_sumk_blocks: Sigma_imp has to be a list of Sigmas for the inequivalent correlated shells, even if it is of length 1!" "transform_to_sumk_blocks: Sigma_imp has to be a list of Sigmas for the inequivalent correlated shells, even if it is of length 1!"
assert len(Sigma_imp) == self.n_inequiv_shells, \ assert len(Sigma_imp) == self.n_inequiv_shells,\
"transform_to_sumk_blocks: give exactly one Sigma for each inequivalent corr. shell!" "transform_to_sumk_blocks: give exactly one Sigma for each inequivalent corr. shell!"
if Sigma_out is None: if Sigma_out is None:
Sigma_out = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[self.corr_to_inequiv[icrsh]].mesh, Sigma_out = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[self.corr_to_inequiv[icrsh]].mesh, space='sumk')
space='sumk')
for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells)]
else: else:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
@ -824,7 +805,7 @@ class SumkDFT(object):
return G_loc return G_loc
def transform_to_solver_blocks(self, G_loc, G_out=None, show_warnings=True): def transform_to_solver_blocks(self, G_loc, G_out=None, show_warnings = True):
""" transform G_loc from sumk to solver space """ transform G_loc from sumk to solver space
Parameters Parameters
@ -862,7 +843,7 @@ class SumkDFT(object):
ish_to=ish, ish_to=ish,
space_from='sumk', space_from='sumk',
G_out=G_out[ish], G_out=G_out[ish],
show_warnings=show_warnings) show_warnings = show_warnings)
# return only the inequivalent shells: # return only the inequivalent shells:
return G_out return G_out
@ -892,12 +873,10 @@ class SumkDFT(object):
""" """
if dm is None: if dm is None:
warn( warn("WARNING: No density matrix given. Calculating density matrix with default parameters. This will be deprecated in future releases.")
"WARNING: No density matrix given. Calculating density matrix with default parameters. This will be deprecated in future releases.")
dm = self.density_matrix(method='using_gf', transform_to_solver_blocks=False) dm = self.density_matrix(method='using_gf', transform_to_solver_blocks=False)
assert len( assert len(dm) == self.n_corr_shells, "The number of density matrices must be equal to the number of correlated shells."
dm) == self.n_corr_shells, "The number of density matrices must be equal to the number of correlated shells."
dens_mat = [dm[self.inequiv_to_corr[ish]] dens_mat = [dm[self.inequiv_to_corr[ish]]
for ish in range(self.n_inequiv_shells)] for ish in range(self.n_inequiv_shells)]
if hloc is None: if hloc is None:
@ -912,8 +891,7 @@ class SumkDFT(object):
self.solver_to_sumk_block[ish] = {} self.solver_to_sumk_block[ish] = {}
for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]: for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
assert sp in dens_mat[ assert sp in dens_mat[ish], f"The density matrix does not contain the block {sp}. Is the input dm given in sumk block structure?"
ish], f"The density matrix does not contain the block {sp}. Is the input dm given in sumk block structure?"
n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim'] n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim']
# gives an index list of entries larger that threshold # gives an index list of entries larger that threshold
dmbool = (abs(dens_mat[ish][sp]) > threshold) dmbool = (abs(dens_mat[ish][sp]) > threshold)
@ -933,10 +911,10 @@ class SumkDFT(object):
pair = offdiag.pop() pair = offdiag.pop()
for b1, b2 in product(blocs, blocs): for b1, b2 in product(blocs, blocs):
if (pair[0] in b1) and (pair[1] in b2): if (pair[0] in b1) and (pair[1] in b2):
if blocs.index(b1) != blocs.index(b2): # In separate blocks? if blocs.index(b1) != blocs.index(b2): # In separate blocks?
# Merge two blocks # Merge two blocks
b1.extend(blocs.pop(blocs.index(b2))) b1.extend(blocs.pop(blocs.index(b2)))
break # Move on to next pair in offdiag break # Move on to next pair in offdiag
# Set the gf_struct for the solver accordingly # Set the gf_struct for the solver accordingly
num_blocs = len(blocs) num_blocs = len(blocs)
@ -1012,10 +990,9 @@ class SumkDFT(object):
""" """
# make a GfImTime from the supplied GfImFreq # make a GfImTime from the supplied GfImFreq
if all(isinstance(g_sh.mesh, MeshImFreq) for g_sh in G): if all(isinstance(g_sh.mesh, MeshImFreq) for g_sh in G):
gf = [BlockGf(name_block_generator=[(name, GfImTime(beta=block.mesh.beta, gf = [BlockGf(name_block_generator = [(name, GfImTime(beta=block.mesh.beta,
indices=block.indices, n_points=len(block.mesh) + 1)) indices=block.indices,n_points=len(block.mesh)+1)) for name, block in g_sh],
for name, block in g_sh], make_copies=False) for g_sh in G]
make_copies=False) for g_sh in G]
for ish in range(len(gf)): for ish in range(len(gf)):
for name, g in gf[ish]: for name, g in gf[ish]:
g.set_from_fourier(G[ish][name]) g.set_from_fourier(G[ish][name])
@ -1027,7 +1004,7 @@ class SumkDFT(object):
gf = [g_sh.copy() for g_sh in G] gf = [g_sh.copy() for g_sh in G]
for ish in range(len(gf)): for ish in range(len(gf)):
for name, g in gf[ish]: for name, g in gf[ish]:
g << 1.0j * (g - g.conjugate().transpose()) / 2.0 / np.pi g << 1.0j*(g-g.conjugate().transpose())/2.0/np.pi
elif all(isinstance(g_sh.mesh, MeshReTime) for g_sh in G): elif all(isinstance(g_sh.mesh, MeshReTime) for g_sh in G):
def get_delta_from_mesh(mesh): def get_delta_from_mesh(mesh):
w0 = None w0 = None
@ -1035,23 +1012,24 @@ class SumkDFT(object):
if w0 is None: if w0 is None:
w0 = w w0 = w
else: else:
return w - w0 return w-w0
gf = [BlockGf(name_block_generator = [(name, GfReFreq(
gf = [BlockGf(name_block_generator=[(name, GfReFreq( window=(-np.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh)),
window=(-np.pi * (len(block.mesh) - 1) / (len(block.mesh) * get_delta_from_mesh(block.mesh)), np.pi*(len(block.mesh)-1) / (len(block.mesh)*get_delta_from_mesh(block.mesh))),
np.pi * (len(block.mesh) - 1) / (len(block.mesh) * get_delta_from_mesh(block.mesh))),
n_points=len(block.mesh), indices=block.indices)) for name, block in g_sh], make_copies=False) n_points=len(block.mesh), indices=block.indices)) for name, block in g_sh], make_copies=False)
for g_sh in G] for g_sh in G]
for ish in range(len(gf)): for ish in range(len(gf)):
for name, g in gf[ish]: for name, g in gf[ish]:
g.set_from_fourier(G[ish][name]) g.set_from_fourier(G[ish][name])
g << 1.0j * (g - g.conjugate().transpose()) / 2.0 / np.pi g << 1.0j*(g-g.conjugate().transpose())/2.0/np.pi
else: else:
raise Exception("G must be a list of BlockGf of either GfImFreq, GfImTime, GfReFreq or GfReTime") raise Exception("G must be a list of BlockGf of either GfImFreq, GfImTime, GfReFreq or GfReTime")
return gf return gf
def analyse_block_structure_from_gf(self, G, threshold=1.e-5, include_shells=None, analyse_deg_shells=True):
def analyse_block_structure_from_gf(self, G, threshold=1.e-5, include_shells=None, analyse_deg_shells = True):
r""" r"""
Determines the block structure of local Green's functions by analysing Determines the block structure of local Green's functions by analysing
the structure of the corresponding non-interacting Green's function. the structure of the corresponding non-interacting Green's function.
@ -1084,8 +1062,7 @@ class SumkDFT(object):
""" """
assert isinstance(G, list), "G must be a list (with elements for each correlated shell)" assert isinstance(G, list), "G must be a list (with elements for each correlated shell)"
assert len( assert len(G) == self.n_corr_shells, "G must have one element for each correlated shell, run extract_G_loc with transform_to_solver_blocks=False to get the correct G"
G) == self.n_corr_shells, "G must have one element for each correlated shell, run extract_G_loc with transform_to_solver_blocks=False to get the correct G"
gf = self._get_hermitian_quantity_from_gf(G) gf = self._get_hermitian_quantity_from_gf(G)
@ -1099,12 +1076,11 @@ class SumkDFT(object):
self.solver_to_sumk[ish] = {} self.solver_to_sumk[ish] = {}
self.solver_to_sumk_block[ish] = {} self.solver_to_sumk_block[ish] = {}
for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]: for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
assert sp in gf[self.inequiv_to_corr[ assert sp in gf[self.inequiv_to_corr[ish]].indices, f"The Green's function does not contain the block {sp}. Is the input G given in sumk block structure?"
ish]].indices, f"The Green's function does not contain the block {sp}. Is the input G given in sumk block structure?"
n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim'] n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim']
# gives an index list of entries larger that threshold # gives an index list of entries larger that threshold
max_gf = np.max(np.abs(gf[self.inequiv_to_corr[ish]][sp].data), 0) max_gf = np.max(np.abs(gf[self.inequiv_to_corr[ish]][sp].data),0)
maxgf_bool = (max_gf > threshold) maxgf_bool = (max_gf > threshold)
# Determine off-diagonal entries in upper triangular part of the # Determine off-diagonal entries in upper triangular part of the
@ -1121,10 +1097,10 @@ class SumkDFT(object):
pair = offdiag.pop() pair = offdiag.pop()
for b1, b2 in product(blocs, blocs): for b1, b2 in product(blocs, blocs):
if (pair[0] in b1) and (pair[1] in b2): if (pair[0] in b1) and (pair[1] in b2):
if blocs.index(b1) != blocs.index(b2): # In separate blocks? if blocs.index(b1) != blocs.index(b2): # In separate blocks?
# Merge two blocks # Merge two blocks
b1.extend(blocs.pop(blocs.index(b2))) b1.extend(blocs.pop(blocs.index(b2)))
break # Move on to next pair in offdiag break # Move on to next pair in offdiag
# Set the gf_struct for the solver accordingly # Set the gf_struct for the solver accordingly
num_blocs = len(blocs) num_blocs = len(blocs)
@ -1150,13 +1126,13 @@ class SumkDFT(object):
# transform G to the new structure # transform G to the new structure
full_structure = BlockStructure.full_structure( full_structure = BlockStructure.full_structure(
[{sp: self.corr_shells[self.inequiv_to_corr[ish]]['dim'] [{sp: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 sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]}
for ish in range(self.n_inequiv_shells)], self.corr_to_inequiv) for ish in range(self.n_inequiv_shells)],self.corr_to_inequiv)
G_transformed = [ G_transformed = [
self.block_structure.convert_gf(G[ish], self.block_structure.convert_gf(G[ish],
full_structure, ish, mesh=G[ish].mesh.copy(), show_warnings=threshold, full_structure, ish, mesh=G[ish].mesh.copy(), show_warnings=threshold,
gf_function=type(G[ish]._first()), space_from='sumk', space_to='solver') gf_function=type(G[ish]._first()), space_from='sumk', space_to='solver')
for ish in range(self.n_inequiv_shells)] for ish in range(self.n_inequiv_shells)]
if analyse_deg_shells: if analyse_deg_shells:
@ -1217,7 +1193,7 @@ class SumkDFT(object):
for ish in include_shells: for ish in include_shells:
for block1 in self.gf_struct_solver[ish].keys(): for block1 in self.gf_struct_solver[ish].keys():
for block2 in self.gf_struct_solver[ish].keys(): for block2 in self.gf_struct_solver[ish].keys():
if block1 == block2: continue if block1==block2: continue
# check if the blocks are already present in the deg_shells # check if the blocks are already present in the deg_shells
ind1 = -1 ind1 = -1
@ -1250,9 +1226,9 @@ class SumkDFT(object):
e1 = np.linalg.eigvalsh(gf1.data[0]) e1 = np.linalg.eigvalsh(gf1.data[0])
e2 = np.linalg.eigvalsh(gf2.data[0]) e2 = np.linalg.eigvalsh(gf2.data[0])
if np.any(abs(e1 - e2) > threshold): continue if np.any(abs(e1-e2) > threshold): continue
for conjugate in [False, True]: for conjugate in [False,True]:
if conjugate: if conjugate:
gf2 = gf2.conjugate() gf2 = gf2.conjugate()
@ -1269,15 +1245,13 @@ class SumkDFT(object):
# product to get a linear problem, which consists # product to get a linear problem, which consists
# of finding the null space of M vec T = 0. # of finding the null space of M vec T = 0.
M = np.kron(np.eye(*gf1.target_shape), gf2.data[0]) - np.kron(gf1.data[0].transpose(), M = np.kron(np.eye(*gf1.target_shape),gf2.data[0])-np.kron(gf1.data[0].transpose(),np.eye(*gf1.target_shape))
np.eye(*gf1.target_shape))
N = null(M, threshold) N = null(M, threshold)
# now we get the intersection of the null spaces # now we get the intersection of the null spaces
# of all values of tau # of all values of tau
for i in range(1, len(gf1.data)): for i in range(1,len(gf1.data)):
M = np.kron(np.eye(*gf1.target_shape), gf2.data[i]) - np.kron(gf1.data[i].transpose(), M = np.kron(np.eye(*gf1.target_shape),gf2.data[i])-np.kron(gf1.data[i].transpose(),np.eye(*gf1.target_shape))
np.eye(*gf1.target_shape))
# transform M into current null space # transform M into current null space
M = np.dot(M, N) M = np.dot(M, N)
N = np.dot(N, null(M, threshold)) N = np.dot(N, null(M, threshold))
@ -1311,7 +1285,6 @@ class SumkDFT(object):
Then, it calculates the chi2 of Then, it calculates the chi2 of
sum y[i] N[:,:,i] y[j].conjugate() N[:,:,j].conjugate().transpose() - eye. sum y[i] N[:,:,i] y[j].conjugate() N[:,:,j].conjugate().transpose() - eye.
""" """
def chi2(y): def chi2(y):
# reinterpret y as complex number # reinterpret y as complex number
y = y.view(complex) y = y.view(complex)
@ -1319,7 +1292,7 @@ class SumkDFT(object):
for a in range(Z.shape[0]): for a in range(Z.shape[0]):
for b in range(Z.shape[1]): for b in range(Z.shape[1]):
ret += np.abs(np.dot(y, np.dot(Z[a, b], y.conjugate())) ret += np.abs(np.dot(y, np.dot(Z[a, b], y.conjugate()))
- (1.0 if a == b else 0.0)) ** 2 - (1.0 if a == b else 0.0))**2
return ret return ret
# use the minimization routine from scipy # use the minimization routine from scipy
@ -1372,8 +1345,7 @@ class SumkDFT(object):
# C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2) # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2)
if (ind1 < 0) and (ind2 >= 0): if (ind1 < 0) and (ind2 >= 0):
if conjugate: if conjugate:
self.deg_shells[ish][ind2][block1] = np.dot(T.conjugate().transpose(), self.deg_shells[ish][ind2][block1] = np.dot(T.conjugate().transpose(), v2[0].conjugate()), not v2[1]
v2[0].conjugate()), not v2[1]
else: else:
self.deg_shells[ish][ind2][block1] = np.dot(T.conjugate().transpose(), v2[0]), v2[1] self.deg_shells[ish][ind2][block1] = np.dot(T.conjugate().transpose(), v2[0]), v2[1]
# the first block is already present # the first block is already present
@ -1404,8 +1376,7 @@ class SumkDFT(object):
# a block was found, break out of the loop # a block was found, break out of the loop
break break
def calculate_diagonalization_matrix(self, prop_to_be_diagonal='eal', calc_in_solver_blocks=True, def calculate_diagonalization_matrix(self, prop_to_be_diagonal='eal', calc_in_solver_blocks=True, write_to_blockstructure = True, shells=None):
write_to_blockstructure=True, shells=None):
""" """
Calculates the diagonalisation matrix, and (optionally) stores it in the BlockStructure. Calculates the diagonalisation matrix, and (optionally) stores it in the BlockStructure.
@ -1436,13 +1407,13 @@ class SumkDFT(object):
if self.block_structure.transformation: if self.block_structure.transformation:
mpi.report( mpi.report(
"calculate_diagonalization_matrix: requires block_structure.transformation = None.") "calculate_diagonalization_matrix: requires block_structure.transformation = None.")
return 0 return 0
# Use all shells # Use all shells
if shells is None: if shells is None:
shells = range(self.n_corr_shells) shells = range(self.n_corr_shells)
elif max(shells) >= self.n_corr_shells: # Check if the shell indices are present elif max(shells) >= self.n_corr_shells: # Check if the shell indices are present
mpi.report("calculate_diagonalization_matrix: shells not correct.") mpi.report("calculate_diagonalization_matrix: shells not correct.")
return 0 return 0
@ -1465,7 +1436,7 @@ class SumkDFT(object):
prop[ish] = self.block_structure.convert_matrix(prop[ish], space_from='sumk', space_to='solver') prop[ish] = self.block_structure.convert_matrix(prop[ish], space_from='sumk', space_to='solver')
# Get diagonalisation matrix, if not already diagonal # Get diagonalisation matrix, if not already diagonal
for name in prop[ish]: for name in prop[ish]:
if np.sum(abs(prop[ish][name] - np.diag(np.diagonal(prop[ish][name])))) > 1e-13: if np.sum(abs(prop[ish][name]-np.diag(np.diagonal(prop[ish][name])))) > 1e-13:
trafo[name] = np.linalg.eigh(prop[ish][name])[1].conj().T trafo[name] = np.linalg.eigh(prop[ish][name])[1].conj().T
else: else:
trafo[name] = np.identity(np.shape(prop[ish][name])[0]) trafo[name] = np.identity(np.shape(prop[ish][name])[0])
@ -1503,7 +1474,7 @@ class SumkDFT(object):
ntoi = self.spin_names_to_ind[self.SO] ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO] spn = self.spin_block_names[self.SO]
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
dims = {sp: self.n_orbitals[ik, ntoi[sp]] for sp in spn} dims = {sp:self.n_orbitals[ik, ntoi[sp]] for sp in spn}
MMat = [np.zeros([dims[sp], dims[sp]], complex) for sp in spn] MMat = [np.zeros([dims[sp], dims[sp]], complex) for sp in spn]
for isp, sp in enumerate(spn): for isp, sp in enumerate(spn):
@ -1524,7 +1495,7 @@ class SumkDFT(object):
n_orb = self.n_orbitals[ik, ind] n_orb = self.n_orbitals[ik, ind]
projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb] projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb]
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()) projmat.transpose().conjugate())
# get data from nodes: # get data from nodes:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
@ -1541,14 +1512,14 @@ class SumkDFT(object):
for sp in dens_mat[icrsh]: for sp in dens_mat[icrsh]:
if self.rot_mat_time_inv[icrsh] == 1: if self.rot_mat_time_inv[icrsh] == 1:
dens_mat[icrsh][sp] = dens_mat[icrsh][sp].conjugate() dens_mat[icrsh][sp] = dens_mat[icrsh][sp].conjugate()
dens_mat[icrsh][sp] = np.dot( dens_mat[icrsh][sp] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh][sp]),
np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh][sp]), self.rot_mat[icrsh])
self.rot_mat[icrsh])
return dens_mat return dens_mat
def density_matrix(self, method='using_gf', mu=None, with_Sigma=True, with_dc=True, broadening=None, 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): transform_to_solver_blocks=True, show_warnings=True):
"""Calculate density matrices in one of two ways. """Calculate density matrices in one of two ways.
Parameters Parameters
@ -1588,8 +1559,7 @@ class SumkDFT(object):
transform_to_solver_blocks, show_warnings) transform_to_solver_blocks, show_warnings)
dens_mat = [G.density() for G in Gloc] dens_mat = [G.density() for G in Gloc]
elif method == "using_point_integration": elif method == "using_point_integration":
warn( warn("WARNING: density_matrix: method 'using_point_integration' is deprecated. Use 'density_matrix_using_point_integration' instead. All additionally provided arguments are ignored.")
"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() dens_mat = self.density_matrix_using_point_integration()
else: else:
raise ValueError("density_matrix: the method '%s' is not supported." % method) raise ValueError("density_matrix: the method '%s' is not supported." % method)
@ -1648,11 +1618,10 @@ class SumkDFT(object):
for ik in range(self.n_k): for ik in range(self.n_k):
n_orb = self.n_orbitals[ik, ind] n_orb = self.n_orbitals[ik, ind]
MMat = np.identity(n_orb, complex) MMat = np.identity(n_orb, complex)
MMat = self.hopping[ MMat = self.hopping[ik, ind, 0:n_orb, 0:n_orb] - (1 - 2 * isp) * self.h_field * MMat
ik, ind, 0:n_orb, 0:n_orb] - (1 - 2 * isp) * self.h_field * MMat
projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb] projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb]
self.Hsumk[icrsh][sp] += self.bz_weights[ik] * np.dot(np.dot(projmat, MMat), self.Hsumk[icrsh][sp] += self.bz_weights[ik] * np.dot(np.dot(projmat, MMat),
projmat.conjugate().transpose()) projmat.conjugate().transpose())
# symmetrisation: # symmetrisation:
if self.symm_op != 0: if self.symm_op != 0:
self.Hsumk = self.symmcorr.symmetrize(self.Hsumk) self.Hsumk = self.symmcorr.symmetrize(self.Hsumk)
@ -1664,9 +1633,8 @@ class SumkDFT(object):
if self.rot_mat_time_inv[icrsh] == 1: if self.rot_mat_time_inv[icrsh] == 1:
self.Hsumk[icrsh][sp] = self.Hsumk[ self.Hsumk[icrsh][sp] = self.Hsumk[
icrsh][sp].conjugate() icrsh][sp].conjugate()
self.Hsumk[icrsh][sp] = np.dot( self.Hsumk[icrsh][sp] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), self.Hsumk[icrsh][sp]),
np.dot(self.rot_mat[icrsh].conjugate().transpose(), self.Hsumk[icrsh][sp]), self.rot_mat[icrsh])
self.rot_mat[icrsh])
# add to matrix: # add to matrix:
for ish in range(self.n_inequiv_shells): for ish in range(self.n_inequiv_shells):
@ -1789,7 +1757,7 @@ class SumkDFT(object):
dim //= 2 dim //= 2
if use_dc_value is None: if use_dc_value is None:
# For legacy compatibility #For legacy compatibility
if use_dc_formula == 0: if use_dc_formula == 0:
mpi.report(f"Detected {use_dc_formula=}, changing to sFLL") mpi.report(f"Detected {use_dc_formula=}, changing to sFLL")
use_dc_formula = "sFLL" use_dc_formula = "sFLL"
@ -1801,8 +1769,7 @@ class SumkDFT(object):
use_dc_formula = "sAMF" use_dc_formula = "sAMF"
for sp in spn: for sp in spn:
DC_val, E_val = compute_DC_from_density(N_tot=Ncrtot, U=U_interact, J=J_hund, n_orbitals=dim, 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)
N_spin=Ncr[sp], method=use_dc_formula)
self.dc_imp[icrsh][sp] *= DC_val self.dc_imp[icrsh][sp] *= DC_val
self.dc_energ[icrsh] = E_val self.dc_energ[icrsh] = E_val
@ -1822,7 +1789,7 @@ class SumkDFT(object):
for sp in spn: for sp in spn:
T = self.block_structure.effective_transformation_sumk[icrsh][sp] T = self.block_structure.effective_transformation_sumk[icrsh][sp]
self.dc_imp[icrsh][sp] = np.dot(T.conjugate().transpose(), self.dc_imp[icrsh][sp] = np.dot(T.conjugate().transpose(),
np.dot(self.dc_imp[icrsh][sp], T)) np.dot(self.dc_imp[icrsh][sp], T))
def add_dc(self): def add_dc(self):
r""" r"""
@ -1844,8 +1811,7 @@ class SumkDFT(object):
for bname, gf in sigma_minus_dc[icrsh]: for bname, gf in sigma_minus_dc[icrsh]:
# Transform dc_imp to global coordinate system # Transform dc_imp to global coordinate system
if self.use_rotations: if self.use_rotations:
gf -= np.dot(self.rot_mat[icrsh], gf -= np.dot(self.rot_mat[icrsh], np.dot(self.dc_imp[icrsh][bname], self.rot_mat[icrsh].conjugate().transpose()))
np.dot(self.dc_imp[icrsh][bname], self.rot_mat[icrsh].conjugate().transpose()))
else: else:
gf -= self.dc_imp[icrsh][bname] gf -= self.dc_imp[icrsh][bname]
@ -1998,11 +1964,9 @@ class SumkDFT(object):
beta = 1 / broadening beta = 1 / broadening
if isinstance(self.mesh, MeshReFreq): if isinstance(self.mesh, MeshReFreq):
def tot_den(bgf): def tot_den(bgf): return bgf.total_density(beta)
return bgf.total_density(beta)
else: else:
def tot_den(bgf): def tot_den(bgf): return bgf.total_density()
return bgf.total_density()
dens = 0.0 dens = 0.0
ikarray = np.array(list(range(self.n_k))) ikarray = np.array(list(range(self.n_k)))
@ -2062,7 +2026,6 @@ class SumkDFT(object):
within specified precision. within specified precision.
""" """
def find_bounds(function, x_init, delta_x, max_loops=1000): def find_bounds(function, x_init, delta_x, max_loops=1000):
mpi.report("Finding bounds on chemical potential") mpi.report("Finding bounds on chemical potential")
x = x_init x = x_init
@ -2075,12 +2038,12 @@ class SumkDFT(object):
nbre_loop = 0 nbre_loop = 0
# abort the loop after maxiter is reached or when y1 and y2 have different sign # abort the loop after maxiter is reached or when y1 and y2 have different sign
while (nbre_loop <= max_loops) and (y2 * y1) > 0: while (nbre_loop <= max_loops) and (y2*y1) > 0:
nbre_loop += 1 nbre_loop += 1
x1 = x2 x1 = x2
y1 = y2 y1 = y2
x2 -= eps * delta_x x2 -= eps*delta_x
y2 = function(x2) y2 = function(x2)
if nbre_loop > (max_loops): if nbre_loop > (max_loops):
@ -2097,11 +2060,8 @@ class SumkDFT(object):
# previous implementation # previous implementation
def F_bisection(mu): def F_bisection(mu): return self.total_density(mu=mu, broadening=broadening, beta=beta).real
return self.total_density(mu=mu, broadening=broadening, beta=beta).real
density = self.density_required - self.charge_below density = self.density_required - self.charge_below
# using scipy.optimize # using scipy.optimize
def F_optimize(mu): def F_optimize(mu):
@ -2155,8 +2115,7 @@ class SumkDFT(object):
return self.chemical_potential return self.chemical_potential
def calc_density_correction(self, filename=None, dm_type=None, spinave=False, kpts_to_write=None, broadening=None, def calc_density_correction(self, filename=None, dm_type=None, spinave=False, kpts_to_write=None, broadening=None, beta=None):
beta=None):
r""" r"""
Calculates the charge density correction and stores it into a file. Calculates the charge density correction and stores it into a file.
@ -2196,7 +2155,7 @@ class SumkDFT(object):
the corresponing total charge `dens`. the corresponing total charge `dens`.
""" """
# automatically set dm_type if required #automatically set dm_type if required
if dm_type is None: if dm_type is None:
dm_type = self.dft_code dm_type = self.dft_code
@ -2281,15 +2240,13 @@ class SumkDFT(object):
G2 = np.sum(self.kpts_cart[ik, :] ** 2) G2 = np.sum(self.kpts_cart[ik, :] ** 2)
# Kerker mixing # Kerker mixing
mix_fac = self.charge_mixing_alpha * G2 / (G2 + self.charge_mixing_gamma ** 2) mix_fac = self.charge_mixing_alpha * G2 / (G2 + self.charge_mixing_gamma ** 2)
deltaN[bname][ik][diag_inds] = (1.0 - mix_fac) * self.deltaNOld[bname][ik][ deltaN[bname][ik][diag_inds] = (1.0 - mix_fac) * self.deltaNOld[bname][ik][diag_inds] + mix_fac * deltaN[bname][ik][diag_inds]
diag_inds] + mix_fac * deltaN[bname][ik][diag_inds]
dens[bname] -= self.bz_weights[ik] * dens_mat_dft[bname][ik].sum().real dens[bname] -= self.bz_weights[ik] * dens_mat_dft[bname][ik].sum().real
isp = ntoi[bname] isp = ntoi[bname]
b1, b2 = band_window[isp][ik, :2] b1, b2 = band_window[isp][ik, :2]
nb = b2 - b1 + 1 nb = b2 - b1 + 1
assert nb == self.n_orbitals[ik, ntoi[bname]], "Number of bands is inconsistent at ik = %s" % (ik) assert nb == self.n_orbitals[ik, ntoi[bname]], "Number of bands is inconsistent at ik = %s"%(ik)
band_en_correction += np.dot(deltaN[bname][ik], self.hopping[ik, isp, :nb, :nb]).trace().real * \ band_en_correction += np.dot(deltaN[bname][ik], self.hopping[ik, isp, :nb, :nb]).trace().real * self.bz_weights[ik]
self.bz_weights[ik]
# mpi reduce: # mpi reduce:
for bname in deltaN: for bname in deltaN:
@ -2299,6 +2256,8 @@ class SumkDFT(object):
self.deltaNOld = copy.copy(deltaN) self.deltaNOld = copy.copy(deltaN)
mpi.barrier() mpi.barrier()
band_en_correction = mpi.all_reduce(band_en_correction) band_en_correction = mpi.all_reduce(band_en_correction)
# now save to file: # now save to file:
@ -2325,10 +2284,8 @@ class SumkDFT(object):
f.write("%s\n" % self.n_orbitals[ik, 0]) f.write("%s\n" % self.n_orbitals[ik, 0])
for inu in range(self.n_orbitals[ik, 0]): for inu in range(self.n_orbitals[ik, 0]):
for imu in range(self.n_orbitals[ik, 0]): for imu in range(self.n_orbitals[ik, 0]):
valre = (deltaN['up'][ik][ valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0 valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
valim = (deltaN['up'][ik][
inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
f.write("%.14f %.14f " % (valre, valim)) f.write("%.14f %.14f " % (valre, valim))
f.write("\n") f.write("\n")
f.write("\n") f.write("\n")
@ -2347,8 +2304,7 @@ class SumkDFT(object):
fout.write("%s\n" % self.n_orbitals[ik, isp]) fout.write("%s\n" % self.n_orbitals[ik, isp])
for inu in range(self.n_orbitals[ik, isp]): for inu in range(self.n_orbitals[ik, isp]):
for imu in range(self.n_orbitals[ik, isp]): for imu in range(self.n_orbitals[ik, isp]):
fout.write("%.14f %.14f " % (deltaN[sp][ik][ fout.write("%.14f %.14f " % (deltaN[sp][ik][inu, imu].real, deltaN[sp][ik][inu, imu].imag))
inu, imu].real, deltaN[sp][ik][inu, imu].imag))
fout.write("\n") fout.write("\n")
fout.write("\n") fout.write("\n")
fout.close() fout.close()
@ -2396,60 +2352,55 @@ class SumkDFT(object):
f.write("\n") f.write("\n")
elif dm_type == 'elk': elif dm_type == 'elk':
# output each k-point density matrix for Elk # output each k-point density matrix for Elk
if mpi.is_master_node(): if mpi.is_master_node():
# read in misc data from .h5 file # read in misc data from .h5 file
things_to_read = ['band_window', 'vkl', 'nstsv'] things_to_read = ['band_window','vkl','nstsv']
self.subgroup_present, self.value_read = self.read_input_from_hdf( self.subgroup_present, self.value_read = self.read_input_from_hdf(
subgrp=self.misc_data, things_to_read=things_to_read) subgrp=self.misc_data, things_to_read=things_to_read)
# open file # open file
with open(filename, 'w') as f: with open(filename, 'w') as f:
# determine the number of spin blocks # determine the number of spin blocks
n_spin_blocks = self.SP + 1 - self.SO n_spin_blocks = self.SP + 1 - self.SO
nbmax = np.max(self.n_orbitals) nbmax = np.max(self.n_orbitals)
# output beta and mu in Hartrees # output beta and mu in Hartrees
beta = self.mesh.beta * self.energy_unit beta = self.mesh.beta * self.energy_unit
mu = self.chemical_potential / self.energy_unit mu = self.chemical_potential/self.energy_unit
# ouput n_k, nspin and max orbitals - a check # ouput n_k, nspin and max orbitals - a check
f.write(" %d %d %d %.14f %.14f ! nkpt, nspin, nstmax, beta, mu\n" % ( f.write(" %d %d %d %.14f %.14f ! nkpt, nspin, nstmax, beta, mu\n"%(self.n_k, n_spin_blocks, nbmax, beta, mu))
self.n_k, n_spin_blocks, nbmax, beta, mu))
for ik in range(self.n_k): for ik in range(self.n_k):
for ispn in range(n_spin_blocks): for ispn in range(n_spin_blocks):
# Determine the SO density matrix band indices from the spinor band indices #Determine the SO density matrix band indices from the spinor band indices
if (self.SO == 1): if(self.SO==1):
band0 = [self.band_window[0][ik, 0], self.band_window[1][ik, 0]] band0=[self.band_window[0][ik, 0],self.band_window[1][ik, 0]]
band1 = [self.band_window[0][ik, 1], self.band_window[1][ik, 1]] band1=[self.band_window[0][ik, 1],self.band_window[1][ik, 1]]
ib1 = int(min(band0)) ib1=int(min(band0))
ib2 = int(max(band1)) ib2=int(max(band1))
else: else:
# Determine the density matrix band indices from the spinor band indices #Determine the density matrix band indices from the spinor band indices
ib1 = self.band_window[ispn][ik, 0] ib1 = self.band_window[ispn][ik, 0]
ib2 = self.band_window[ispn][ik, 1] ib2 = self.band_window[ispn][ik, 1]
f.write(" %d %d %d %d ! ik, ispn, minist, maxist\n" % (ik + 1, ispn + 1, ib1, ib2)) f.write(" %d %d %d %d ! ik, ispn, minist, maxist\n"%(ik + 1, ispn + 1, ib1, ib2))
for inu in range(self.n_orbitals[ik, ispn]): for inu in range(self.n_orbitals[ik, ispn]):
for imu in range(self.n_orbitals[ik, ispn]): for imu in range(self.n_orbitals[ik, ispn]):
# output non-magnetic or spin-averaged density matrix #output non-magnetic or spin-averaged density matrix
if ((self.SP == 0) or (spinave)): if((self.SP==0) or (spinave)):
valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][ valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
inu, imu].real) / 2.0 valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][ else:
inu, imu].imag) / 2.0 valre = deltaN[spn[ispn]][ik][inu, imu].real
else: valim = deltaN[spn[ispn]][ik][inu, imu].imag
valre = deltaN[spn[ispn]][ik][inu, imu].real f.write(" %.14f %.14f"%(valre, valim))
valim = deltaN[spn[ispn]][ik][inu, imu].imag f.write("\n")
f.write(" %.14f %.14f" % (valre, valim))
f.write("\n")
elif dm_type == 'qe': elif dm_type == 'qe':
if self.SP == 0: if self.SP == 0:
mpi.report( mpi.report("SUMK calc_density_correction: WARNING! Averaging out spin-polarized correction in the density channel")
"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)
mpi.report(" %i -1 ! Number of k-points, default number of bands\n" % (self.n_k)) mpi.report(" %i -1 ! Number of k-points, default number of bands\n"%(self.n_k))
for ik in range(self.n_k): for ik in range(self.n_k):
ib1 = band_window[0][ik, 0] ib1 = band_window[0][ik, 0]
ib2 = band_window[0][ik, 1] ib2 = band_window[0][ik, 1]
@ -2458,7 +2409,7 @@ class SumkDFT(object):
valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0 valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0 valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
# write into delta_N # write into delta_N
delta_N[ik, inu, imu] = valre + 1j * valim delta_N[ik, inu, imu] = valre + 1j*valim
if mpi.is_master_node(): if mpi.is_master_node():
with HDFArchive(self.hdf_file, 'a') as ar: with HDFArchive(self.hdf_file, 'a') as ar:
if subgrp not in ar: if subgrp not in ar:
@ -2469,7 +2420,7 @@ class SumkDFT(object):
else: else:
raise NotImplementedError("Unknown density matrix type: '%s'" % (dm_type)) raise NotImplementedError("Unknown density matrix type: '%s'"%(dm_type))
res = deltaN, dens res = deltaN, dens
@ -2491,9 +2442,9 @@ class SumkDFT(object):
self.max_band_energy = max_band_energy self.max_band_energy = max_band_energy
return min_band_energy, max_band_energy return min_band_energy, max_band_energy
################ ################
# FIXME LEAVE UNDOCUMENTED # FIXME LEAVE UNDOCUMENTED
################ ################
def check_projectors(self): def check_projectors(self):
"""Calculated the density matrix from projectors (DM = P Pdagger) to check that it is correct and """Calculated the density matrix from projectors (DM = P Pdagger) to check that it is correct and
@ -2506,8 +2457,7 @@ class SumkDFT(object):
dim = self.corr_shells[icrsh]['dim'] dim = self.corr_shells[icrsh]['dim']
n_orb = self.n_orbitals[ik, 0] n_orb = self.n_orbitals[ik, 0]
projmat = self.proj_mat[ik, 0, icrsh, 0:dim, 0:n_orb] projmat = self.proj_mat[ik, 0, icrsh, 0:dim, 0:n_orb]
dens_mat[icrsh][ dens_mat[icrsh][:, :] += np.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik]
:, :] += np.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik]
if self.symm_op != 0: if self.symm_op != 0:
dens_mat = self.symmcorr.symmetrize(dens_mat) dens_mat = self.symmcorr.symmetrize(dens_mat)
@ -2517,8 +2467,7 @@ class SumkDFT(object):
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
if self.rot_mat_time_inv[icrsh] == 1: if self.rot_mat_time_inv[icrsh] == 1:
dens_mat[icrsh] = dens_mat[icrsh].conjugate() dens_mat[icrsh] = dens_mat[icrsh].conjugate()
dens_mat[icrsh] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh]), dens_mat[icrsh] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh]),self.rot_mat[icrsh])
self.rot_mat[icrsh])
return dens_mat return dens_mat
@ -2542,51 +2491,39 @@ class SumkDFT(object):
# after introducing the block_structure class # after introducing the block_structure class
def __get_gf_struct_sumk(self): def __get_gf_struct_sumk(self):
return self.block_structure.gf_struct_sumk return self.block_structure.gf_struct_sumk
def __set_gf_struct_sumk(self,value):
def __set_gf_struct_sumk(self, value):
self.block_structure.gf_struct_sumk = value self.block_structure.gf_struct_sumk = value
gf_struct_sumk = property(__get_gf_struct_sumk,__set_gf_struct_sumk)
gf_struct_sumk = property(__get_gf_struct_sumk, __set_gf_struct_sumk)
def __get_gf_struct_solver(self): def __get_gf_struct_solver(self):
return self.block_structure.gf_struct_solver return self.block_structure.gf_struct_solver
def __set_gf_struct_solver(self,value):
def __set_gf_struct_solver(self, value):
self.block_structure.gf_struct_solver = value self.block_structure.gf_struct_solver = value
gf_struct_solver = property(__get_gf_struct_solver,__set_gf_struct_solver)
gf_struct_solver = property(__get_gf_struct_solver, __set_gf_struct_solver)
def __get_solver_to_sumk(self): def __get_solver_to_sumk(self):
return self.block_structure.solver_to_sumk return self.block_structure.solver_to_sumk
def __set_solver_to_sumk(self,value):
def __set_solver_to_sumk(self, value):
self.block_structure.solver_to_sumk = value self.block_structure.solver_to_sumk = value
solver_to_sumk = property(__get_solver_to_sumk,__set_solver_to_sumk)
solver_to_sumk = property(__get_solver_to_sumk, __set_solver_to_sumk)
def __get_sumk_to_solver(self): def __get_sumk_to_solver(self):
return self.block_structure.sumk_to_solver return self.block_structure.sumk_to_solver
def __set_sumk_to_solver(self,value):
def __set_sumk_to_solver(self, value):
self.block_structure.sumk_to_solver = value self.block_structure.sumk_to_solver = value
sumk_to_solver = property(__get_sumk_to_solver,__set_sumk_to_solver)
sumk_to_solver = property(__get_sumk_to_solver, __set_sumk_to_solver)
def __get_solver_to_sumk_block(self): def __get_solver_to_sumk_block(self):
return self.block_structure.solver_to_sumk_block return self.block_structure.solver_to_sumk_block
def __set_solver_to_sumk_block(self,value):
def __set_solver_to_sumk_block(self, value):
self.block_structure.solver_to_sumk_block = value self.block_structure.solver_to_sumk_block = value
solver_to_sumk_block = property(__get_solver_to_sumk_block,__set_solver_to_sumk_block)
solver_to_sumk_block = property(__get_solver_to_sumk_block, __set_solver_to_sumk_block)
def __get_deg_shells(self): def __get_deg_shells(self):
return self.block_structure.deg_shells return self.block_structure.deg_shells
def __set_deg_shells(self,value):
def __set_deg_shells(self, value):
self.block_structure.deg_shells = value self.block_structure.deg_shells = value
deg_shells = property(__get_deg_shells,__set_deg_shells)
deg_shells = property(__get_deg_shells, __set_deg_shells)
@property @property
def gf_struct_solver_list(self): def gf_struct_solver_list(self):
@ -2606,16 +2543,12 @@ class SumkDFT(object):
def __get_corr_to_inequiv(self): def __get_corr_to_inequiv(self):
return self.block_structure.corr_to_inequiv return self.block_structure.corr_to_inequiv
def __set_corr_to_inequiv(self, value): def __set_corr_to_inequiv(self, value):
self.block_structure.corr_to_inequiv = value self.block_structure.corr_to_inequiv = value
corr_to_inequiv = property(__get_corr_to_inequiv, __set_corr_to_inequiv) corr_to_inequiv = property(__get_corr_to_inequiv, __set_corr_to_inequiv)
def __get_inequiv_to_corr(self): def __get_inequiv_to_corr(self):
return self.block_structure.inequiv_to_corr return self.block_structure.inequiv_to_corr
def __set_inequiv_to_corr(self, value): def __set_inequiv_to_corr(self, value):
self.block_structure.inequiv_to_corr = value self.block_structure.inequiv_to_corr = value
inequiv_to_corr = property(__get_inequiv_to_corr, __set_inequiv_to_corr) inequiv_to_corr = property(__get_inequiv_to_corr, __set_inequiv_to_corr)