mirror of https://github.com/triqs/dft_tools
1500 lines
66 KiB
Python
1500 lines
66 KiB
Python
|
|
##########################################################################
|
|
#
|
|
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
|
#
|
|
# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
|
|
#
|
|
# 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/>.
|
|
#
|
|
##########################################################################
|
|
|
|
from types import *
|
|
import numpy
|
|
import pytriqs.utility.dichotomy as dichotomy
|
|
from pytriqs.gf.local import *
|
|
import pytriqs.utility.mpi as mpi
|
|
from pytriqs.archive import *
|
|
from symmetry import *
|
|
from sets import Set
|
|
from itertools import product
|
|
|
|
|
|
class SumkDFT:
|
|
"""This class provides a general SumK method for combining ab-initio code and pytriqs."""
|
|
|
|
def __init__(self, hdf_file, h_field=0.0, use_dft_blocks=False,
|
|
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',
|
|
misc_data='dft_misc_input'):
|
|
r"""
|
|
Initialises the class from data previously stored into an hdf5 archive.
|
|
|
|
Parameters
|
|
----------
|
|
hdf_file : string
|
|
Name of hdf5 containing the data.
|
|
h_field : scalar, optional
|
|
The value of magnetic field to add to the DFT Hamiltonian.
|
|
The contribution -h_field*sigma is added to diagonal elements of the Hamiltonian.
|
|
It cannot be used with the spin-orbit coupling on; namely h_field is set to 0 if self.SO=True.
|
|
use_dft_blocks : boolean, optional
|
|
If True, the local Green's function matrix for each spin is divided into smaller blocks
|
|
with the block structure determined from the DFT density matrix of the corresponding correlated shell.
|
|
dft_data : string, optional
|
|
Name of hdf5 subgroup in which DFT data for projector and lattice Green's function construction are stored.
|
|
symmcorr_data : string, optional
|
|
Name of hdf5 subgroup in which DFT data on symmetries of correlated shells
|
|
(symmetry operations, permutaion matrices etc.) are stored.
|
|
parproj_data : string, optional
|
|
Name of hdf5 subgroup in which DFT data on non-normalized projectors for non-correlated
|
|
states (used in the partial density of states calculations) are stored.
|
|
symmpar_data : string, optional
|
|
Name of hdf5 subgroup in which DFT data on symmetries of the non-normalized projectors
|
|
are stored.
|
|
bands_data : string, optional
|
|
Name of hdf5 subgroup in which DFT data necessary for band-structure/k-resolved spectral
|
|
function calculations (projectors, DFT Hamiltonian for a chosen path in the Brillouin zone etc.)
|
|
are stored.
|
|
transp_data : string, optional
|
|
Name of hdf5 subgroup in which DFT data necessary for transport calculations are stored.
|
|
misc_data : string, optional
|
|
Name of hdf5 subgroup in which miscellaneous DFT data are stored.
|
|
"""
|
|
|
|
if not type(hdf_file) == StringType:
|
|
mpi.report("Give a string for the hdf5 filename to read the input!")
|
|
else:
|
|
self.hdf_file = hdf_file
|
|
self.dft_data = dft_data
|
|
self.symmcorr_data = symmcorr_data
|
|
self.parproj_data = parproj_data
|
|
self.symmpar_data = symmpar_data
|
|
self.bands_data = bands_data
|
|
self.transp_data = transp_data
|
|
self.misc_data = misc_data
|
|
self.h_field = h_field
|
|
|
|
# Read input from HDF:
|
|
things_to_read = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required',
|
|
'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat',
|
|
'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.value_read = self.read_input_from_hdf(
|
|
subgrp=self.dft_data, things_to_read=things_to_read)
|
|
if self.symm_op:
|
|
self.symmcorr = Symmetry(hdf_file, subgroup=self.symmcorr_data)
|
|
|
|
if self.SO and (abs(self.h_field) > 0.000001):
|
|
self.h_field = 0.0
|
|
mpi.report(
|
|
"For SO, the external magnetic field is not implemented, setting it to 0!")
|
|
|
|
self.spin_block_names = [['up', 'down'], ['ud']]
|
|
self.n_spin_blocks = [2, 1]
|
|
# Convert spin_block_names to indices -- if spin polarized,
|
|
# differentiate up and down blocks
|
|
self.spin_names_to_ind = [{}, {}]
|
|
for iso in range(2): # SO = 0 or 1
|
|
for isp in range(self.n_spin_blocks[iso]):
|
|
self.spin_names_to_ind[iso][
|
|
self.spin_block_names[iso][isp]] = isp * self.SP
|
|
|
|
# GF structure used for the local things in the k sums
|
|
# Most general form allowing for all hybridisation, i.e. largest
|
|
# blocks possible
|
|
self.gf_struct_sumk = [[(sp, range(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, range(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 <->
|
|
# gf_struct_solver
|
|
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)]
|
|
for ish in range(self.n_inequiv_shells):
|
|
for block, inner_list in self.gf_struct_sumk[self.inequiv_to_corr[ish]]:
|
|
self.solver_to_sumk_block[ish][block] = block
|
|
for inner in inner_list:
|
|
self.sumk_to_solver[ish][
|
|
(block, inner)] = (block, inner)
|
|
self.solver_to_sumk[ish][
|
|
(block, inner)] = (block, inner)
|
|
# assume no shells are degenerate
|
|
self.deg_shells = [[] for ish in range(self.n_inequiv_shells)]
|
|
|
|
self.chemical_potential = 0.0 # initialise mu
|
|
self.init_dc() # initialise the double counting
|
|
|
|
# Analyse the block structure and determine the smallest gf_struct
|
|
# blocks and maps, if desired
|
|
if use_dft_blocks:
|
|
self.analyse_block_structure()
|
|
|
|
################
|
|
# hdf5 FUNCTIONS
|
|
################
|
|
|
|
def read_input_from_hdf(self, subgrp, things_to_read):
|
|
r"""
|
|
Reads data from the HDF file. Prints a warning if a requested dataset is not found.
|
|
|
|
Parameters
|
|
----------
|
|
subgrp : string
|
|
Name of hdf5 file subgroup from which the data are to be read.
|
|
things_to_read : list of strings
|
|
List of datasets to be read from the hdf5 file.
|
|
|
|
Returns
|
|
-------
|
|
subgroup_present : boolean
|
|
Is the subgrp is present in hdf5 file?
|
|
value_read : boolean
|
|
Did the reading of requested datasets succeed?
|
|
|
|
"""
|
|
|
|
value_read = True
|
|
# initialise variables on all nodes to ensure mpi broadcast works at
|
|
# the end
|
|
for it in things_to_read:
|
|
setattr(self, it, 0)
|
|
subgroup_present = 0
|
|
|
|
if mpi.is_master_node():
|
|
ar = HDFArchive(self.hdf_file, 'r')
|
|
if subgrp in ar:
|
|
subgroup_present = True
|
|
# first read the necessary things:
|
|
for it in things_to_read:
|
|
if it in ar[subgrp]:
|
|
setattr(self, it, ar[subgrp][it])
|
|
else:
|
|
mpi.report("Loading %s failed!" % it)
|
|
value_read = False
|
|
else:
|
|
if (len(things_to_read) != 0):
|
|
mpi.report(
|
|
"Loading failed: No %s subgroup in hdf5!" % subgrp)
|
|
subgroup_present = False
|
|
value_read = False
|
|
del ar
|
|
# now do the broadcasting:
|
|
for it in things_to_read:
|
|
setattr(self, it, mpi.bcast(getattr(self, it)))
|
|
subgroup_present = mpi.bcast(subgroup_present)
|
|
value_read = mpi.bcast(value_read)
|
|
|
|
return subgroup_present, value_read
|
|
|
|
def save(self, things_to_save, subgrp='user_data'):
|
|
r"""
|
|
Saves data from a list into the HDF file. Prints a warning if a requested data is not found in SumkDFT object.
|
|
|
|
Parameters
|
|
----------
|
|
things_to_save : list of strings
|
|
List of datasets to be saved into the hdf5 file.
|
|
subgrp : string, optional
|
|
Name of hdf5 file subgroup in which the data are to be stored.
|
|
"""
|
|
|
|
if not (mpi.is_master_node()):
|
|
return # do nothing on nodes
|
|
ar = HDFArchive(self.hdf_file, 'a')
|
|
if not subgrp in ar:
|
|
ar.create_group(subgrp)
|
|
for it in things_to_save:
|
|
try:
|
|
ar[subgrp][it] = getattr(self, it)
|
|
except:
|
|
mpi.report("%s not found, and so not saved." % it)
|
|
del ar
|
|
|
|
def load(self, things_to_load, subgrp='user_data'):
|
|
r"""
|
|
Loads user data from the HDF file. Raises an exeption if a requested dataset is not found.
|
|
|
|
Parameters
|
|
----------
|
|
things_to_read : list of strings
|
|
List of datasets to be read from the hdf5 file.
|
|
subgrp : string, optional
|
|
Name of hdf5 file subgroup from which the data are to be read.
|
|
|
|
Returns
|
|
-------
|
|
list_to_return : list
|
|
A list containing data read from hdf5.
|
|
"""
|
|
|
|
if not (mpi.is_master_node()):
|
|
return # do nothing on nodes
|
|
ar = HDFArchive(self.hdf_file, 'r')
|
|
if not subgrp in ar:
|
|
mpi.report("Loading %s failed!" % subgrp)
|
|
list_to_return = []
|
|
for it in things_to_load:
|
|
try:
|
|
list_to_return.append(ar[subgrp][it])
|
|
except:
|
|
raise ValueError, "load: %s not found, and so not loaded." % it
|
|
del ar
|
|
return list_to_return
|
|
|
|
################
|
|
# CORE FUNCTIONS
|
|
################
|
|
|
|
def downfold(self, ik, ish, bname, gf_to_downfold, gf_inp, shells='corr', ir=None):
|
|
r"""
|
|
Downfolds a block of the Green's function for a given shell and k-point using the corresponding projector matrices.
|
|
|
|
Parameters
|
|
----------
|
|
ik : integer
|
|
k-point index for which the downfolding is to be done.
|
|
ish : integer
|
|
Shell index of GF to be downfolded.
|
|
|
|
- if shells='corr': ish labels all correlated shells (equivalent or not)
|
|
- if shells='all': ish labels only representative (inequivalent) non-correlated shells
|
|
|
|
bname : string
|
|
Block name of the target block of the lattice Green's function.
|
|
gf_to_downfold : Gf
|
|
Block of the Green's function that is to be downfolded.
|
|
gf_inp : Gf
|
|
FIXME
|
|
shells : string, optional
|
|
|
|
- if shells='corr': orthonormalized projectors for correlated shells are used for the downfolding.
|
|
- if shells='all': non-normalized projectors for all included shells are used for the downfolding.
|
|
|
|
ir : integer, optional
|
|
Index of equivalent site in the non-correlated shell 'ish', only used if shells='all'.
|
|
|
|
Returns
|
|
-------
|
|
gf_downfolded : Gf
|
|
Downfolded block of the lattice Green's function.
|
|
"""
|
|
|
|
gf_downfolded = gf_inp.copy()
|
|
# get spin index for proj. matrices
|
|
isp = self.spin_names_to_ind[self.SO][bname]
|
|
n_orb = self.n_orbitals[ik, isp]
|
|
if shells == 'corr':
|
|
dim = self.corr_shells[ish]['dim']
|
|
projmat = self.proj_mat[ik, isp, ish, 0:dim, 0:n_orb]
|
|
elif shells == 'all':
|
|
if ir is None:
|
|
raise ValueError, "downfold: provide ir if treating all shells."
|
|
dim = self.shells[ish]['dim']
|
|
projmat = self.proj_mat_all[ik, isp, ish, ir, 0:dim, 0:n_orb]
|
|
|
|
gf_downfolded.from_L_G_R(
|
|
projmat, gf_to_downfold, projmat.conjugate().transpose())
|
|
|
|
return gf_downfolded
|
|
|
|
def upfold(self, ik, ish, bname, gf_to_upfold, gf_inp, shells='corr', ir=None):
|
|
r"""
|
|
Upfolds a block of the Green's function for a given shell and k-point using the corresponding projector matrices.
|
|
|
|
Parameters
|
|
----------
|
|
ik : integer
|
|
k-point index for which the upfolding is to be done.
|
|
ish : integer
|
|
Shell index of GF to be upfolded.
|
|
|
|
- if shells='corr': ish labels all correlated shells (equivalent or not)
|
|
- if shells='all': ish labels only representative (inequivalent) non-correlated shells
|
|
|
|
bname : string
|
|
Block name of the target block of the lattice Green's function.
|
|
gf_to_upfold : Gf
|
|
Block of the Green's function that is to be upfolded.
|
|
gf_inp : Gf
|
|
FIXME
|
|
shells : string, optional
|
|
|
|
- if shells='corr': orthonormalized projectors for correlated shells are used for the upfolding.
|
|
- if shells='all': non-normalized projectors for all included shells are used for the upfolding.
|
|
|
|
ir : integer, optional
|
|
Index of equivalent site in the non-correlated shell 'ish', only used if shells='all'.
|
|
|
|
Returns
|
|
-------
|
|
gf_upfolded : Gf
|
|
Upfolded block of the lattice Green's function.
|
|
"""
|
|
|
|
gf_upfolded = gf_inp.copy()
|
|
# get spin index for proj. matrices
|
|
isp = self.spin_names_to_ind[self.SO][bname]
|
|
n_orb = self.n_orbitals[ik, isp]
|
|
if shells == 'corr':
|
|
dim = self.corr_shells[ish]['dim']
|
|
projmat = self.proj_mat[ik, isp, ish, 0:dim, 0:n_orb]
|
|
elif shells == 'all':
|
|
if ir is None:
|
|
raise ValueError, "upfold: provide ir if treating all shells."
|
|
dim = self.shells[ish]['dim']
|
|
projmat = self.proj_mat_all[ik, isp, ish, ir, 0:dim, 0:n_orb]
|
|
|
|
gf_upfolded.from_L_G_R(
|
|
projmat.conjugate().transpose(), gf_to_upfold, projmat)
|
|
|
|
return gf_upfolded
|
|
|
|
def rotloc(self, ish, gf_to_rotate, direction, shells='corr'):
|
|
r"""
|
|
Rotates a block of the local Green's function from the local frame to the global frame and vice versa.
|
|
|
|
Parameters
|
|
----------
|
|
ish : integer
|
|
Shell index of GF to be rotated.
|
|
|
|
- if shells='corr': ish labels all correlated shells (equivalent or not)
|
|
- if shells='all': ish labels only representative (inequivalent) non-correlated shells
|
|
|
|
gf_to_rotate : Gf
|
|
Block of the Green's function that is to be rotated.
|
|
direction : string
|
|
The direction of rotation can be either
|
|
|
|
- 'toLocal' : global -> local transformation,
|
|
- 'toGlobal' : local -> global transformation.
|
|
|
|
shells : string, optional
|
|
|
|
- if shells='corr': the rotation matrix for the correlated shell 'ish' is used,
|
|
- if shells='all': the rotation matrix for the generic (non-correlated) shell 'ish' is used.
|
|
|
|
Returns
|
|
-------
|
|
gf_rotated : Gf
|
|
Rotated block of the local Green's function.
|
|
"""
|
|
|
|
assert ((direction == 'toLocal') or (direction == 'toGlobal')
|
|
), "rotloc: Give direction 'toLocal' or 'toGlobal'."
|
|
gf_rotated = gf_to_rotate.copy()
|
|
if shells == 'corr':
|
|
rot_mat_time_inv = self.rot_mat_time_inv
|
|
rot_mat = self.rot_mat
|
|
elif shells == 'all':
|
|
rot_mat_time_inv = self.rot_mat_all_time_inv
|
|
rot_mat = self.rot_mat_all
|
|
|
|
if direction == 'toGlobal':
|
|
|
|
if (rot_mat_time_inv[ish] == 1) and self.SO:
|
|
gf_rotated << gf_rotated.transpose()
|
|
gf_rotated.from_L_G_R(rot_mat[ish].conjugate(
|
|
), gf_rotated, rot_mat[ish].transpose())
|
|
else:
|
|
gf_rotated.from_L_G_R(rot_mat[ish], gf_rotated, rot_mat[
|
|
ish].conjugate().transpose())
|
|
|
|
elif direction == 'toLocal':
|
|
|
|
if (rot_mat_time_inv[ish] == 1) and self.SO:
|
|
gf_rotated << gf_rotated.transpose()
|
|
gf_rotated.from_L_G_R(rot_mat[ish].transpose(
|
|
), gf_rotated, rot_mat[ish].conjugate())
|
|
else:
|
|
gf_rotated.from_L_G_R(rot_mat[ish].conjugate(
|
|
).transpose(), gf_rotated, rot_mat[ish])
|
|
|
|
return gf_rotated
|
|
|
|
def lattice_gf(self, ik, mu=None, iw_or_w="iw", beta=40, broadening=None, mesh=None, with_Sigma=True, with_dc=True):
|
|
r"""
|
|
Calculates the lattice Green function for a given k-point from the DFT Hamiltonian and the self energy.
|
|
|
|
Parameters
|
|
----------
|
|
ik : integer
|
|
k-point index.
|
|
mu : real, optional
|
|
Chemical potential for which the Green's function is to be calculated.
|
|
If not provided, self.chemical_potential is used for mu.
|
|
iw_or_w : string, optional
|
|
|
|
- `iw_or_w` = 'iw' for a imaginary-frequency self-energy
|
|
- `iw_or_w` = 'w' for a real-frequency self-energy
|
|
|
|
beta : real, optional
|
|
Inverse temperature.
|
|
broadening : real, 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'.
|
|
mesh : list, optional
|
|
Data defining mesh on which the real-axis GF will be calculated, given in the form
|
|
(om_min,om_max,n_points), where om_min is the minimum omega, om_max is the maximum omega and n_points is the number of points.
|
|
with_Sigma : boolean, optional
|
|
If True the GF will be calculated with the self-energy stored in self.Sigmaimp_(w/iw), for real/Matsubara GF, respectively.
|
|
In this case the mesh is taken from the self.Sigma_imp object.
|
|
If with_Sigma=True but self.Sigmaimp_(w/iw) is not present, with_Sigma is reset to False.
|
|
with_dc : boolean, optional
|
|
if True and with_Sigma=True, the dc correction is substracted from the self-energy before it is included into GF.
|
|
|
|
Returns
|
|
-------
|
|
G_latt : BlockGf
|
|
Lattice Green's function.
|
|
|
|
"""
|
|
if mu is None:
|
|
mu = self.chemical_potential
|
|
ntoi = self.spin_names_to_ind[self.SO]
|
|
spn = self.spin_block_names[self.SO]
|
|
if (iw_or_w != "iw") and (iw_or_w != "w"):
|
|
raise ValueError, "lattice_gf: Implemented only for Re/Im frequency functions."
|
|
if not hasattr(self, "Sigma_imp_" + iw_or_w):
|
|
with_Sigma = False
|
|
if broadening is None:
|
|
if mesh is None:
|
|
broadening = 0.01
|
|
else: # broadening = 2 * \Delta omega, where \Delta omega is the spacing of omega points
|
|
broadening = 2.0 * ((mesh[1] - mesh[0]) / (mesh[2] - 1))
|
|
|
|
# Are we including Sigma?
|
|
if with_Sigma:
|
|
Sigma_imp = getattr(self, "Sigma_imp_" + iw_or_w)
|
|
sigma_minus_dc = [s.copy() for s in Sigma_imp]
|
|
if with_dc:
|
|
sigma_minus_dc = self.add_dc(iw_or_w)
|
|
if iw_or_w == "iw":
|
|
# override beta if Sigma_iw is present
|
|
beta = Sigma_imp[0].mesh.beta
|
|
mesh = Sigma_imp[0].mesh
|
|
elif iw_or_w == "w":
|
|
mesh = Sigma_imp[0].mesh
|
|
else:
|
|
if iw_or_w == "iw":
|
|
if beta is None:
|
|
raise ValueError, "lattice_gf: Give the beta for the lattice GfReFreq."
|
|
# Default number of Matsubara frequencies
|
|
mesh = MeshImFreq(beta=beta, S='Fermion', n_max=1025)
|
|
elif iw_or_w == "w":
|
|
if mesh is None:
|
|
raise ValueError, "lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq."
|
|
mesh = MeshReFreq(mesh[0], mesh[1], mesh[2])
|
|
|
|
# Check if G_latt is present
|
|
set_up_G_latt = False # Assume not
|
|
if not hasattr(self, "G_latt_" + iw_or_w):
|
|
# Need to create G_latt_(i)w
|
|
set_up_G_latt = True
|
|
else: # Check that existing GF is consistent
|
|
G_latt = getattr(self, "G_latt_" + iw_or_w)
|
|
GFsize = [gf.N1 for bname, gf in G_latt]
|
|
unchangedsize = all([self.n_orbitals[ik, ntoi[spn[isp]]] == GFsize[
|
|
isp] for isp in range(self.n_spin_blocks[self.SO])])
|
|
if not unchangedsize:
|
|
set_up_G_latt = True
|
|
if (iw_or_w == "iw") and (self.G_latt_iw.mesh.beta != beta):
|
|
set_up_G_latt = True # additional check for ImFreq
|
|
|
|
# Set up G_latt
|
|
if set_up_G_latt:
|
|
block_structure = [
|
|
range(self.n_orbitals[ik, ntoi[sp]]) for sp in spn]
|
|
gf_struct = [(spn[isp], block_structure[isp])
|
|
for isp in range(self.n_spin_blocks[self.SO])]
|
|
block_ind_list = [block for block, inner in gf_struct]
|
|
if iw_or_w == "iw":
|
|
glist = lambda: [GfImFreq(indices=inner, mesh=mesh)
|
|
for block, inner in gf_struct]
|
|
elif iw_or_w == "w":
|
|
glist = lambda: [GfReFreq(indices=inner, mesh=mesh)
|
|
for block, inner in gf_struct]
|
|
G_latt = BlockGf(name_list=block_ind_list,
|
|
block_list=glist(), make_copies=False)
|
|
G_latt.zero()
|
|
|
|
if iw_or_w == "iw":
|
|
G_latt << iOmega_n
|
|
elif iw_or_w == "w":
|
|
G_latt << Omega + 1j * broadening
|
|
|
|
idmat = [numpy.identity(
|
|
self.n_orbitals[ik, ntoi[sp]], numpy.complex_) for sp in spn]
|
|
M = copy.deepcopy(idmat)
|
|
for ibl in range(self.n_spin_blocks[self.SO]):
|
|
ind = ntoi[spn[ibl]]
|
|
n_orb = self.n_orbitals[ik, ind]
|
|
M[ibl] = self.hopping[ik, ind, 0:n_orb, 0:n_orb] - \
|
|
(idmat[ibl] * mu) - (idmat[ibl] * self.h_field * (1 - 2 * ibl))
|
|
G_latt -= M
|
|
|
|
if with_Sigma:
|
|
for icrsh in range(self.n_corr_shells):
|
|
for bname, gf in G_latt:
|
|
gf -= self.upfold(ik, icrsh, bname,
|
|
sigma_minus_dc[icrsh][bname], gf)
|
|
|
|
G_latt.invert()
|
|
setattr(self, "G_latt_" + iw_or_w, G_latt)
|
|
|
|
return G_latt
|
|
|
|
def set_Sigma(self, Sigma_imp):
|
|
self.put_Sigma(Sigma_imp)
|
|
|
|
def put_Sigma(self, Sigma_imp):
|
|
r"""
|
|
Inserts the impurity self-energies into the sumk_dft class.
|
|
|
|
Parameters
|
|
----------
|
|
Sigma_imp : list of BlockGf (Green's function) objects
|
|
List containing impurity self-energy for all inequivalent correlated shells.
|
|
Self-energies for equivalent shells are then automatically set by this function.
|
|
The self-energies can be of the real or imaginary-frequency type.
|
|
"""
|
|
|
|
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!"
|
|
assert len(
|
|
Sigma_imp) == self.n_inequiv_shells, "put_Sigma: give exactly one Sigma for each inequivalent corr. shell!"
|
|
|
|
# init self.Sigma_imp_(i)w:
|
|
if all(type(gf) == GfImFreq for bname, gf in Sigma_imp[0]):
|
|
# Imaginary frequency Sigma:
|
|
self.Sigma_imp_iw = [BlockGf(name_block_generator=[(block, GfImFreq(indices=inner, mesh=Sigma_imp[0].mesh))
|
|
for block, inner in self.gf_struct_sumk[icrsh]], make_copies=False)
|
|
for icrsh in range(self.n_corr_shells)]
|
|
SK_Sigma_imp = self.Sigma_imp_iw
|
|
elif all(type(gf) == GfReFreq for bname, gf in Sigma_imp[0]):
|
|
# Real frequency Sigma:
|
|
self.Sigma_imp_w = [BlockGf(name_block_generator=[(block, GfReFreq(indices=inner, mesh=Sigma_imp[0].mesh))
|
|
for block, inner in self.gf_struct_sumk[icrsh]], make_copies=False)
|
|
for icrsh in range(self.n_corr_shells)]
|
|
SK_Sigma_imp = self.Sigma_imp_w
|
|
else:
|
|
raise ValueError, "put_Sigma: This type of Sigma is not handled."
|
|
|
|
# transform the CTQMC blocks to the full matrix:
|
|
for icrsh in range(self.n_corr_shells):
|
|
# ish is the index of the inequivalent shell corresponding to icrsh
|
|
ish = self.corr_to_inequiv[icrsh]
|
|
for block, inner in self.gf_struct_solver[ish].iteritems():
|
|
for ind1 in inner:
|
|
for ind2 in inner:
|
|
block_sumk, ind1_sumk = self.solver_to_sumk[
|
|
ish][(block, ind1)]
|
|
block_sumk, ind2_sumk = self.solver_to_sumk[
|
|
ish][(block, ind2)]
|
|
SK_Sigma_imp[icrsh][block_sumk][
|
|
ind1_sumk, ind2_sumk] << Sigma_imp[ish][block][ind1, ind2]
|
|
|
|
# rotation from local to global coordinate system:
|
|
if self.use_rotations:
|
|
for icrsh in range(self.n_corr_shells):
|
|
for bname, gf in SK_Sigma_imp[icrsh]:
|
|
gf << self.rotloc(icrsh, gf, direction='toGlobal')
|
|
|
|
def extract_G_loc(self, mu=None, iw_or_w='iw', with_Sigma=True, with_dc=True, broadening=None):
|
|
r"""
|
|
Extracts the local downfolded Green function by the Brillouin-zone integration of the lattice Green's function.
|
|
|
|
Parameters
|
|
----------
|
|
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.
|
|
|
|
Returns
|
|
-------
|
|
G_loc_inequiv : list of BlockGf (Green's function) objects
|
|
List of the local Green's functions for all inequivalent correlated shells,
|
|
rotated into the corresponding local frames.
|
|
|
|
"""
|
|
|
|
if mu is None:
|
|
mu = self.chemical_potential
|
|
|
|
if iw_or_w == "iw":
|
|
G_loc = [self.Sigma_imp_iw[icrsh].copy() for icrsh in range(
|
|
self.n_corr_shells)] # this list will be returned
|
|
beta = G_loc[0].mesh.beta
|
|
G_loc_inequiv = [BlockGf(name_block_generator=[(block, GfImFreq(indices=inner, mesh=G_loc[0].mesh)) for block, inner in self.gf_struct_solver[ish].iteritems()],
|
|
make_copies=False) for ish in range(self.n_inequiv_shells)]
|
|
elif iw_or_w == "w":
|
|
G_loc = [self.Sigma_imp_w[icrsh].copy() for icrsh in range(
|
|
self.n_corr_shells)] # this list will be returned
|
|
mesh = G_loc[0].mesh
|
|
G_loc_inequiv = [BlockGf(name_block_generator=[(block, GfReFreq(indices=inner, mesh=mesh)) for block, inner in self.gf_struct_solver[ish].iteritems()],
|
|
make_copies=False) for ish in range(self.n_inequiv_shells)]
|
|
|
|
for icrsh in range(self.n_corr_shells):
|
|
G_loc[icrsh].zero() # initialize to zero
|
|
|
|
ikarray = numpy.array(range(self.n_k))
|
|
for ik in mpi.slice_array(ikarray):
|
|
if iw_or_w == 'iw':
|
|
G_latt = self.lattice_gf(
|
|
ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, beta=beta)
|
|
elif iw_or_w == 'w':
|
|
G_latt = self.lattice_gf(
|
|
ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening, mesh=mesh)
|
|
G_latt *= self.bz_weights[ik]
|
|
|
|
for icrsh in range(self.n_corr_shells):
|
|
# init temporary storage
|
|
tmp = G_loc[icrsh].copy()
|
|
for bname, gf in tmp:
|
|
tmp[bname] << self.downfold(
|
|
ik, icrsh, bname, G_latt[bname], gf)
|
|
G_loc[icrsh] += tmp
|
|
|
|
# Collect data from mpi
|
|
for icrsh in range(self.n_corr_shells):
|
|
G_loc[icrsh] << mpi.all_reduce(
|
|
mpi.world, G_loc[icrsh], lambda x, y: x + y)
|
|
mpi.barrier()
|
|
|
|
# G_loc[:] is now the sum over k projected to the local orbitals.
|
|
# here comes the symmetrisation, if needed:
|
|
if self.symm_op != 0:
|
|
G_loc = self.symmcorr.symmetrize(G_loc)
|
|
|
|
# G_loc is rotated to the local coordinate system:
|
|
if self.use_rotations:
|
|
for icrsh in range(self.n_corr_shells):
|
|
for bname, gf in G_loc[icrsh]:
|
|
G_loc[icrsh][bname] << self.rotloc(
|
|
icrsh, gf, direction='toLocal')
|
|
|
|
# transform to CTQMC blocks:
|
|
for ish in range(self.n_inequiv_shells):
|
|
for block, inner in self.gf_struct_solver[ish].iteritems():
|
|
for ind1 in inner:
|
|
for ind2 in inner:
|
|
block_sumk, ind1_sumk = self.solver_to_sumk[
|
|
ish][(block, ind1)]
|
|
block_sumk, ind2_sumk = self.solver_to_sumk[
|
|
ish][(block, ind2)]
|
|
G_loc_inequiv[ish][block][ind1, ind2] << G_loc[
|
|
self.inequiv_to_corr[ish]][block_sumk][ind1_sumk, ind2_sumk]
|
|
|
|
# return only the inequivalent shells:
|
|
return G_loc_inequiv
|
|
|
|
def analyse_block_structure(self, threshold=0.00001, include_shells=None, dm=None):
|
|
r"""
|
|
Determines the block structure of local Green's functions by analysing the structure of
|
|
the corresponding density matrices. The resulting block structures for correlated shells
|
|
are stored in self.gf_struct_solver list.
|
|
|
|
Parameters
|
|
----------
|
|
threshold : real, optional
|
|
If the difference between density matrix elements is below threshold,
|
|
they are considered to be equal.
|
|
include_shells : list of integers, optional
|
|
List of correlated shells to be analysed.
|
|
If include_shells is not provided all correlated shells will be analysed.
|
|
dm : list of dict, optional
|
|
List of density matrices from which block stuctures are to be analysed.
|
|
Each density matrix is a dict {block names: 2d numpy arrays}.
|
|
If not provided, dm will be calculated from the DFT Hamiltonian by a simple-point BZ integration.
|
|
"""
|
|
|
|
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)]
|
|
|
|
if dm is None:
|
|
dm = self.density_matrix(method='using_point_integration')
|
|
dens_mat = [dm[self.inequiv_to_corr[ish]]
|
|
for ish in range(self.n_inequiv_shells)]
|
|
|
|
if include_shells is None:
|
|
include_shells = range(self.n_inequiv_shells)
|
|
for ish in include_shells:
|
|
|
|
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']
|
|
# gives an index list of entries larger that threshold
|
|
dmbool = (abs(dens_mat[ish][sp]) > threshold)
|
|
|
|
# Determine off-diagonal entries in upper triangular part of
|
|
# density matrix
|
|
offdiag = Set([])
|
|
for i in range(n_orb):
|
|
for j in range(i + 1, n_orb):
|
|
if dmbool[i, j]:
|
|
offdiag.add((i, j))
|
|
|
|
# Determine the number of non-hybridising blocks in the gf
|
|
blocs = [[i] for i in range(n_orb)]
|
|
while len(offdiag) != 0:
|
|
pair = offdiag.pop()
|
|
for b1, b2 in product(blocs, blocs):
|
|
if (pair[0] in b1) and (pair[1] in b2):
|
|
if blocs.index(b1) != blocs.index(b2): # In separate blocks?
|
|
# Merge two blocks
|
|
b1.extend(blocs.pop(blocs.index(b2)))
|
|
break # Move on to next pair in offdiag
|
|
|
|
# Set the gf_struct for the solver accordingly
|
|
num_blocs = len(blocs)
|
|
for i in range(num_blocs):
|
|
blocs[i].sort()
|
|
self.gf_struct_solver[ish].update(
|
|
[('%s_%s' % (sp, i), range(len(blocs[i])))])
|
|
|
|
# Construct sumk_to_solver taking (sumk_block, sumk_index) --> (solver_block, solver_inner)
|
|
# and solver_to_sumk taking (solver_block, solver_inner) -->
|
|
# (sumk_block, sumk_index)
|
|
for i in range(num_blocs):
|
|
for j in range(len(blocs[i])):
|
|
block_sumk = sp
|
|
inner_sumk = blocs[i][j]
|
|
block_solv = '%s_%s' % (sp, i)
|
|
inner_solv = j
|
|
self.sumk_to_solver[ish][(block_sumk, inner_sumk)] = (
|
|
block_solv, inner_solv)
|
|
self.solver_to_sumk[ish][(block_solv, inner_solv)] = (
|
|
block_sumk, inner_sumk)
|
|
self.solver_to_sumk_block[ish][block_solv] = block_sumk
|
|
|
|
# Now calculate degeneracies of orbitals
|
|
dm = {}
|
|
for block, inner in self.gf_struct_solver[ish].iteritems():
|
|
# get dm for the blocks:
|
|
dm[block] = numpy.zeros(
|
|
[len(inner), len(inner)], numpy.complex_)
|
|
for ind1 in inner:
|
|
for ind2 in inner:
|
|
block_sumk, ind1_sumk = self.solver_to_sumk[
|
|
ish][(block, ind1)]
|
|
block_sumk, ind2_sumk = self.solver_to_sumk[
|
|
ish][(block, ind2)]
|
|
dm[block][ind1, ind2] = dens_mat[ish][
|
|
block_sumk][ind1_sumk, ind2_sumk]
|
|
|
|
for block1 in self.gf_struct_solver[ish].iterkeys():
|
|
for block2 in self.gf_struct_solver[ish].iterkeys():
|
|
if dm[block1].shape == dm[block2].shape:
|
|
if ((abs(dm[block1] - dm[block2]) < threshold).all()) and (block1 != block2):
|
|
ind1 = -1
|
|
ind2 = -2
|
|
# check if it was already there:
|
|
for n, ind in enumerate(self.deg_shells[ish]):
|
|
if block1 in ind:
|
|
ind1 = n
|
|
if block2 in ind:
|
|
ind2 = n
|
|
if (ind1 < 0) and (ind2 >= 0):
|
|
self.deg_shells[ish][ind2].append(block1)
|
|
elif (ind1 >= 0) and (ind2 < 0):
|
|
self.deg_shells[ish][ind1].append(block2)
|
|
elif (ind1 < 0) and (ind2 < 0):
|
|
self.deg_shells[ish].append([block1, block2])
|
|
|
|
def density_matrix(self, method='using_gf', beta=40.0):
|
|
"""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.
|
|
|
|
Returns
|
|
-------
|
|
dens_mat : list of dicts
|
|
Density matrix for each spin in each correlated shell.
|
|
"""
|
|
dens_mat = [{} for icrsh in range(self.n_corr_shells)]
|
|
for icrsh in range(self.n_corr_shells):
|
|
for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
|
|
dens_mat[icrsh][sp] = numpy.zeros(
|
|
[self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_)
|
|
|
|
ikarray = numpy.array(range(self.n_k))
|
|
for ik in mpi.slice_array(ikarray):
|
|
|
|
if method == "using_gf":
|
|
|
|
G_latt_iw = self.lattice_gf(
|
|
ik=ik, mu=self.chemical_potential, iw_or_w="iw", beta=beta)
|
|
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]
|
|
unchangedsize = all(
|
|
[self.n_orbitals[ik, ntoi[sp]] == self.n_orbitals[0, ntoi[sp]] for sp in spn])
|
|
if unchangedsize:
|
|
dim = self.n_orbitals[0, ntoi[sp]]
|
|
else:
|
|
dim = self.n_orbitals[ik, ntoi[sp]]
|
|
MMat = [numpy.zeros([dim, dim], numpy.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
|
|
|
|
for icrsh in range(self.n_corr_shells):
|
|
for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
|
|
ind = self.spin_names_to_ind[
|
|
self.corr_shells[icrsh]['SO']][sp]
|
|
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] += numpy.dot(numpy.dot(projmat, MMat[isp]),
|
|
projmat.transpose().conjugate())
|
|
elif method == "using_point_integration":
|
|
dens_mat[icrsh][sp] += self.bz_weights[ik] * numpy.dot(numpy.dot(projmat, MMat[isp]),
|
|
projmat.transpose().conjugate())
|
|
|
|
# get data from nodes:
|
|
for icrsh in range(self.n_corr_shells):
|
|
for sp in dens_mat[icrsh]:
|
|
dens_mat[icrsh][sp] = mpi.all_reduce(
|
|
mpi.world, dens_mat[icrsh][sp], lambda x, y: x + y)
|
|
mpi.barrier()
|
|
|
|
if self.symm_op != 0:
|
|
dens_mat = self.symmcorr.symmetrize(dens_mat)
|
|
|
|
# Rotate to local coordinate system:
|
|
if self.use_rotations:
|
|
for icrsh in range(self.n_corr_shells):
|
|
for sp in dens_mat[icrsh]:
|
|
if self.rot_mat_time_inv[icrsh] == 1:
|
|
dens_mat[icrsh][sp] = dens_mat[icrsh][sp].conjugate()
|
|
dens_mat[icrsh][sp] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh][sp]),
|
|
self.rot_mat[icrsh])
|
|
|
|
return dens_mat
|
|
|
|
# For simple dft input, get crystal field splittings.
|
|
def eff_atomic_levels(self):
|
|
r"""
|
|
Calculates the effective local Hamiltonian required as an input for
|
|
the Hubbard I Solver.
|
|
The local Hamiltonian (effective atomic levels) is calculated by
|
|
projecting the on-site Bloch Hamiltonian:
|
|
|
|
.. math:: H^{loc}_{m m'} = \sum_{k} P_{m \nu}(k) H_{\nu\nu'}(k) P^{*}_{\nu' m'}(k),
|
|
|
|
where
|
|
|
|
.. math:: H_{\nu\nu'}(k) = [\epsilon_{\nu k} - h_{z} \sigma_{z}] \delta_{\nu\nu'}.
|
|
|
|
Parameters
|
|
----------
|
|
None
|
|
|
|
Returns
|
|
-------
|
|
eff_atlevels : gf_struct_sumk like
|
|
Effective local Hamiltonian :math:`H^{loc}_{m m'}` for each correlated shell.
|
|
|
|
"""
|
|
|
|
# define matrices for inequivalent shells:
|
|
eff_atlevels = [{} for ish in range(self.n_inequiv_shells)]
|
|
for ish in range(self.n_inequiv_shells):
|
|
for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]:
|
|
eff_atlevels[ish][sp] = numpy.identity(
|
|
self.corr_shells[self.inequiv_to_corr[ish]]['dim'], numpy.complex_)
|
|
eff_atlevels[ish][sp] *= -self.chemical_potential
|
|
eff_atlevels[ish][
|
|
sp] -= self.dc_imp[self.inequiv_to_corr[ish]][sp]
|
|
|
|
# sum over k:
|
|
if not hasattr(self, "Hsumk"):
|
|
# calculate the sum over k. Does not depend on mu, so do it only
|
|
# once:
|
|
self.Hsumk = [{} for icrsh in range(self.n_corr_shells)]
|
|
for icrsh in range(self.n_corr_shells):
|
|
dim = self.corr_shells[icrsh]['dim']
|
|
for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
|
|
self.Hsumk[icrsh][sp] = numpy.zeros(
|
|
[dim, dim], numpy.complex_)
|
|
for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
|
|
ind = self.spin_names_to_ind[
|
|
self.corr_shells[icrsh]['SO']][sp]
|
|
for ik in range(self.n_k):
|
|
n_orb = self.n_orbitals[ik, ind]
|
|
MMat = numpy.identity(n_orb, numpy.complex_)
|
|
MMat = self.hopping[
|
|
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]
|
|
self.Hsumk[icrsh][sp] += self.bz_weights[ik] * numpy.dot(numpy.dot(projmat, MMat),
|
|
projmat.conjugate().transpose())
|
|
# symmetrisation:
|
|
if self.symm_op != 0:
|
|
self.Hsumk = self.symmcorr.symmetrize(self.Hsumk)
|
|
|
|
# Rotate to local coordinate system:
|
|
if self.use_rotations:
|
|
for icrsh in range(self.n_corr_shells):
|
|
for sp in self.Hsumk[icrsh]:
|
|
if self.rot_mat_time_inv[icrsh] == 1:
|
|
self.Hsumk[icrsh][sp] = self.Hsumk[
|
|
icrsh][sp].conjugate()
|
|
self.Hsumk[icrsh][sp] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), self.Hsumk[icrsh][sp]),
|
|
self.rot_mat[icrsh])
|
|
|
|
# add to matrix:
|
|
for ish in range(self.n_inequiv_shells):
|
|
for sp in eff_atlevels[ish]:
|
|
eff_atlevels[ish][
|
|
sp] += self.Hsumk[self.inequiv_to_corr[ish]][sp]
|
|
|
|
return eff_atlevels
|
|
|
|
def init_dc(self):
|
|
r"""
|
|
Initializes the double counting terms.
|
|
|
|
Parameters
|
|
----------
|
|
None
|
|
|
|
"""
|
|
self.dc_imp = [{} for icrsh in range(self.n_corr_shells)]
|
|
for icrsh in range(self.n_corr_shells):
|
|
dim = self.corr_shells[icrsh]['dim']
|
|
spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
|
|
for sp in spn:
|
|
self.dc_imp[icrsh][sp] = numpy.zeros([dim, dim], numpy.float_)
|
|
self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)]
|
|
|
|
def set_dc(self, dc_imp, dc_energ):
|
|
r"""
|
|
Sets double counting corrections to given values.
|
|
|
|
Parameters
|
|
----------
|
|
dc_imp : gf_struct_sumk like
|
|
Double-counting self-energy term.
|
|
dc_energ : list of floats
|
|
Double-counting energy corrections for each correlated shell.
|
|
|
|
"""
|
|
|
|
self.dc_imp = dc_imp
|
|
self.dc_energ = dc_energ
|
|
|
|
def calc_dc(self, dens_mat, orb=0, U_interact=None, J_hund=None, use_dc_formula=0, use_dc_value=None):
|
|
r"""
|
|
Calculates and sets the double counting corrections.
|
|
|
|
If 'use_dc_value' is provided the double-counting term is uniformly initialized
|
|
with this constant and 'U_interact' and 'J_hund' are ignored.
|
|
|
|
If 'use_dc_value' is None the correction is evaluated according to
|
|
one of the following formulae:
|
|
|
|
* use_dc_formula = 0: fully-localised limit (FLL)
|
|
* use_dc_formula = 1: Held's formula, i.e. mean-field formula for the Kanamori
|
|
type of the interaction Hamiltonian
|
|
* use_dc_formula = 2: around mean-field (AMF)
|
|
|
|
Note that FLL and AMF formulae were derived assuming a full Slater-type interaction
|
|
term and should be thus used accordingly. For the Kanamori-type interaction
|
|
one should use formula 1.
|
|
|
|
The double-counting self-energy term is stored in `self.dc_imp` and the energy
|
|
correction in `self.dc_energ`.
|
|
|
|
Parameters
|
|
----------
|
|
dens_mat : gf_struct_solver like
|
|
Density matrix for the specified correlated shell.
|
|
orb : int, optional
|
|
Index of an inequivalent shell.
|
|
U_interact : float, optional
|
|
Value of interaction parameter `U`.
|
|
J_hund : float, optional
|
|
Value of interaction parameter `J`.
|
|
use_dc_formula : int, optional
|
|
Type of double-counting correction (see description).
|
|
use_dc_value : float, optional
|
|
Value of the double-counting correction. If specified
|
|
`U_interact`, `J_hund` and `use_dc_formula` are ignored.
|
|
|
|
"""
|
|
|
|
for icrsh in range(self.n_corr_shells):
|
|
|
|
# ish is the index of the inequivalent shell corresponding to icrsh
|
|
ish = self.corr_to_inequiv[icrsh]
|
|
if ish != orb:
|
|
continue # ignore this orbital
|
|
# *(1+self.corr_shells[icrsh]['SO'])
|
|
dim = self.corr_shells[icrsh]['dim']
|
|
spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
|
|
|
|
Ncr = {sp: 0.0 for sp in spn}
|
|
for block, inner in self.gf_struct_solver[ish].iteritems():
|
|
bl = self.solver_to_sumk_block[ish][block]
|
|
Ncr[bl] += dens_mat[block].real.trace()
|
|
Ncrtot = sum(Ncr.itervalues())
|
|
for sp in spn:
|
|
self.dc_imp[icrsh][sp] = numpy.identity(dim, numpy.float_)
|
|
if self.SP == 0: # average the densities if there is no SP:
|
|
Ncr[sp] = Ncrtot / len(spn)
|
|
# correction for SO: we have only one block in this case, but
|
|
# in DC we need N/2
|
|
elif self.SP == 1 and self.SO == 1:
|
|
Ncr[sp] = Ncrtot / 2.0
|
|
|
|
if use_dc_value is None:
|
|
|
|
if U_interact is None and J_hund is None:
|
|
raise ValueError, "set_dc: either provide U_interact and J_hund or set use_dc_value to dc value."
|
|
|
|
if use_dc_formula == 0: # FLL
|
|
|
|
self.dc_energ[icrsh] = U_interact / \
|
|
2.0 * Ncrtot * (Ncrtot - 1.0)
|
|
for sp in spn:
|
|
Uav = U_interact * (Ncrtot - 0.5) - \
|
|
J_hund * (Ncr[sp] - 0.5)
|
|
self.dc_imp[icrsh][sp] *= Uav
|
|
self.dc_energ[icrsh] -= J_hund / \
|
|
2.0 * (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 * 0.5 * 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: # use value provided for user to determine dc_energ and dc_imp
|
|
|
|
self.dc_energ[icrsh] = use_dc_value * Ncrtot
|
|
for sp in spn:
|
|
self.dc_imp[icrsh][sp] *= use_dc_value
|
|
|
|
mpi.report(
|
|
"DC for shell %(icrsh)i = %(use_dc_value)f" % locals())
|
|
mpi.report("DC energy = %s" % self.dc_energ[icrsh])
|
|
|
|
def add_dc(self, iw_or_w="iw"):
|
|
r"""
|
|
Subtracts the double counting term from the impurity self energy.
|
|
|
|
Parameters
|
|
----------
|
|
iw_or_w : string, optional
|
|
|
|
- `iw_or_w` = 'iw' for a imaginary-frequency self-energy
|
|
- `iw_or_w` = 'w' for a real-frequency self-energy
|
|
|
|
Returns
|
|
-------
|
|
sigma_minus_dc : gf_struct_sumk like
|
|
Self-energy with a subtracted double-counting term.
|
|
|
|
"""
|
|
|
|
# Be careful: Sigma_imp is already in the global coordinate system!!
|
|
sigma_minus_dc = [s.copy()
|
|
for s in getattr(self, "Sigma_imp_" + iw_or_w)]
|
|
for icrsh in range(self.n_corr_shells):
|
|
for bname, gf in sigma_minus_dc[icrsh]:
|
|
# Transform dc_imp to global coordinate system
|
|
dccont = numpy.dot(self.rot_mat[icrsh], numpy.dot(self.dc_imp[icrsh][
|
|
bname], self.rot_mat[icrsh].conjugate().transpose()))
|
|
sigma_minus_dc[icrsh][bname] -= dccont
|
|
|
|
return sigma_minus_dc
|
|
|
|
def symm_deg_gf(self, gf_to_symm, orb):
|
|
r"""
|
|
Averages a GF over degenerate shells.
|
|
|
|
Degenerate shells of an inequivalent correlated shell are defined by
|
|
`self.deg_shells`. This function enforces corresponding degeneracies
|
|
in the input GF.
|
|
|
|
Parameters
|
|
----------
|
|
gf_to_symm : gf_struct_solver like
|
|
Input GF.
|
|
orb : int
|
|
Index of an inequivalent shell.
|
|
|
|
"""
|
|
|
|
for degsh in self.deg_shells[orb]:
|
|
ss = gf_to_symm[degsh[0]].copy()
|
|
ss.zero()
|
|
n_deg = len(degsh)
|
|
for bl in degsh:
|
|
ss += gf_to_symm[bl] / (1.0 * n_deg)
|
|
for bl in degsh:
|
|
gf_to_symm[bl] << ss
|
|
|
|
def total_density(self, mu=None, iw_or_w="iw", with_Sigma=True, with_dc=True, broadening=None):
|
|
r"""
|
|
Calculates the total charge within the energy window for a given chemical potential.
|
|
The chemical potential is either given by parameter `mu` or, if it is not specified,
|
|
taken from `self.chemical_potential`.
|
|
|
|
The total charge is calculated from the trace of the GF in the Bloch basis.
|
|
By default, a full interacting GF is used. To use the non-interacting GF, set
|
|
parameter `with_Sigma = False`.
|
|
|
|
The number of bands within the energy windows generally depends on `k`. The trace is
|
|
therefore calculated separately for each `k`-point.
|
|
|
|
Since in general n_orbitals depends on k, the calculation is done in the following order:
|
|
|
|
.. math:: n_{tot} = \sum_{k} n(k),
|
|
|
|
with
|
|
|
|
.. math:: n(k) = Tr G_{\nu\nu'}(k, i\omega_{n}).
|
|
|
|
The calculation is done in the global coordinate system, if distinction is made between local/global.
|
|
|
|
Parameters
|
|
----------
|
|
mu : float, optional
|
|
Input chemical potential. If not specified, `self.chemical_potential` is used instead.
|
|
iw_or_w : string, optional
|
|
- `iw_or_w` = 'iw' for a imaginary-frequency self-energy
|
|
- `iw_or_w` = 'w' for a real-frequency self-energy
|
|
with_Sigma : boolean, optional
|
|
If `True` the full interacing GF is evaluated, otherwise the self-energy is not
|
|
included and the charge would correspond to a non-interacting system.
|
|
with_dc : boolean, optional
|
|
Whether or not to subtract the double-counting term from the self-energy.
|
|
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.
|
|
|
|
Returns
|
|
-------
|
|
dens : float
|
|
Total charge :math:`n_{tot}`.
|
|
|
|
"""
|
|
|
|
if mu is None:
|
|
mu = self.chemical_potential
|
|
dens = 0.0
|
|
ikarray = numpy.array(range(self.n_k))
|
|
for ik in mpi.slice_array(ikarray):
|
|
G_latt = self.lattice_gf(
|
|
ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening)
|
|
dens += self.bz_weights[ik] * G_latt.total_density()
|
|
# collect data from mpi:
|
|
dens = mpi.all_reduce(mpi.world, dens, lambda x, y: x + y)
|
|
mpi.barrier()
|
|
|
|
return dens
|
|
|
|
def set_mu(self, mu):
|
|
r"""
|
|
Sets a new chemical potential.
|
|
|
|
Parameters
|
|
----------
|
|
mu : float
|
|
New value of the chemical potential.
|
|
|
|
"""
|
|
self.chemical_potential = mu
|
|
|
|
def calc_mu(self, precision=0.01, iw_or_w='iw', broadening=None):
|
|
r"""
|
|
Searches for the chemical potential that gives the DFT total charge.
|
|
A simple bisection method is used.
|
|
|
|
Parameters
|
|
----------
|
|
precision : float, optional
|
|
A desired precision of the resulting total charge.
|
|
iw_or_w : string, optional
|
|
- `iw_or_w` = 'iw' for a imaginary-frequency self-energy
|
|
- `iw_or_w` = 'w' for a real-frequency self-energy
|
|
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.
|
|
|
|
Returns
|
|
-------
|
|
mu : float
|
|
Value of the chemical potential giving the DFT total charge
|
|
within specified precision.
|
|
|
|
"""
|
|
F = lambda mu: self.total_density(
|
|
mu=mu, iw_or_w=iw_or_w, broadening=broadening)
|
|
density = self.density_required - self.charge_below
|
|
|
|
self.chemical_potential = dichotomy.dichotomy(function=F,
|
|
x_init=self.chemical_potential, y_value=density,
|
|
precision_on_y=precision, delta_x=0.5, max_loops=100,
|
|
x_name="Chemical Potential", y_name="Total Density",
|
|
verbosity=3)[0]
|
|
|
|
return self.chemical_potential
|
|
|
|
def calc_density_correction(self, filename='dens_mat.dat'):
|
|
r"""
|
|
Calculates the charge density correction and stores it into a file.
|
|
|
|
The charge density correction is needed for charge-self-consistent DFT+DMFT calculations.
|
|
It represents a density matrix of the interacting system defined in Bloch basis
|
|
and it is calculated from the sum over Matsubara frequecies of the full GF,
|
|
|
|
..math:: N_{\nu\nu'}(k) = \sum_{i\omega_{n}} G_{\nu\nu'}(k, i\omega_{n})
|
|
|
|
The density matrix for every `k`-point is stored into a file.
|
|
|
|
Parameters
|
|
----------
|
|
filename : string
|
|
Name of the file to store the charge density correction.
|
|
|
|
Returns
|
|
-------
|
|
(deltaN, dens) : tuple
|
|
Returns a tuple containing the density matrix `deltaN` and
|
|
the corresponing total charge `dens`.
|
|
|
|
"""
|
|
|
|
assert type(
|
|
filename) == StringType, "calc_density_correction: filename has to be a string!"
|
|
|
|
ntoi = self.spin_names_to_ind[self.SO]
|
|
spn = self.spin_block_names[self.SO]
|
|
dens = {sp: 0.0 for sp in spn}
|
|
|
|
# Set up deltaN:
|
|
deltaN = {}
|
|
for sp in spn:
|
|
deltaN[sp] = [numpy.zeros([self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[
|
|
ik, ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)]
|
|
|
|
ikarray = numpy.array(range(self.n_k))
|
|
for ik in mpi.slice_array(ikarray):
|
|
G_latt_iw = self.lattice_gf(
|
|
ik=ik, mu=self.chemical_potential, iw_or_w="iw")
|
|
for bname, gf in G_latt_iw:
|
|
deltaN[bname][ik] = G_latt_iw[bname].density()
|
|
dens[bname] += self.bz_weights[ik] * \
|
|
G_latt_iw[bname].total_density()
|
|
|
|
# mpi reduce:
|
|
for bname in deltaN:
|
|
for ik in range(self.n_k):
|
|
deltaN[bname][ik] = mpi.all_reduce(
|
|
mpi.world, deltaN[bname][ik], lambda x, y: x + y)
|
|
dens[bname] = mpi.all_reduce(
|
|
mpi.world, dens[bname], lambda x, y: x + y)
|
|
mpi.barrier()
|
|
|
|
# now save to file:
|
|
if mpi.is_master_node():
|
|
if self.SP == 0:
|
|
f = open(filename, 'w')
|
|
else:
|
|
f = open(filename + 'up', 'w')
|
|
f1 = open(filename + 'dn', 'w')
|
|
# write chemical potential (in Rydberg):
|
|
f.write("%.14f\n" % (self.chemical_potential / self.energy_unit))
|
|
if self.SP != 0:
|
|
f1.write("%.14f\n" %
|
|
(self.chemical_potential / self.energy_unit))
|
|
# write beta in rydberg-1
|
|
f.write("%.14f\n" % (G_latt_iw.mesh.beta * self.energy_unit))
|
|
if self.SP != 0:
|
|
f1.write("%.14f\n" % (G_latt_iw.mesh.beta * self.energy_unit))
|
|
|
|
if self.SP == 0: # no spin-polarization
|
|
|
|
for ik in range(self.n_k):
|
|
f.write("%s\n" % self.n_orbitals[ik, 0])
|
|
for inu in range(self.n_orbitals[ik, 0]):
|
|
for imu in range(self.n_orbitals[ik, 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
|
|
f.write("%.14f %.14f " % (valre, valim))
|
|
f.write("\n")
|
|
f.write("\n")
|
|
f.close()
|
|
|
|
elif self.SP == 1: # with spin-polarization
|
|
|
|
# dict of filename: (spin index, block_name)
|
|
if self.SO == 0:
|
|
to_write = {f: (0, 'up'), f1: (1, 'down')}
|
|
if self.SO == 1:
|
|
to_write = {f: (0, 'ud'), f1: (0, 'ud')}
|
|
for fout in to_write.iterkeys():
|
|
isp, sp = to_write[fout]
|
|
for ik in range(self.n_k):
|
|
fout.write("%s\n" % self.n_orbitals[ik, isp])
|
|
for inu in range(self.n_orbitals[ik, isp]):
|
|
for imu in range(self.n_orbitals[ik, isp]):
|
|
fout.write("%.14f %.14f " % (deltaN[sp][ik][
|
|
inu, imu].real, deltaN[sp][ik][inu, imu].imag))
|
|
fout.write("\n")
|
|
fout.write("\n")
|
|
fout.close()
|
|
|
|
return deltaN, dens
|
|
|
|
################
|
|
# FIXME LEAVE UNDOCUMENTED
|
|
################
|
|
|
|
def calc_dc_for_density(self, orb, dc_init, dens_mat, density=None, precision=0.01):
|
|
"""Searches for DC in order to fulfill charge neutrality.
|
|
If density is given, then DC is set such that the LOCAL charge of orbital
|
|
orb coincides with the given density."""
|
|
|
|
def F(dc):
|
|
self.calc_dc(dens_mat=dens_mat, U_interact=0,
|
|
J_hund=0, orb=orb, use_dc_value=dc)
|
|
if dens_req is None:
|
|
return self.total_density(mu=mu)
|
|
else:
|
|
return self.extract_G_loc()[orb].total_density()
|
|
|
|
if density is None:
|
|
density = self.density_required - self.charge_below
|
|
|
|
dc = dichotomy.dichotomy(function=F,
|
|
x_init=dc_init, y_value=density,
|
|
precision_on_y=precision, delta_x=0.5,
|
|
max_loops=100, x_name="Double Counting", y_name="Total Density",
|
|
verbosity=3)[0]
|
|
|
|
return dc
|
|
|
|
def check_projectors(self):
|
|
"""Calculated the density matrix from projectors (DM = P Pdagger) to check that it is correct and
|
|
specifically that it matches DFT."""
|
|
dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_)
|
|
for icrsh in range(self.n_corr_shells)]
|
|
|
|
for ik in range(self.n_k):
|
|
for icrsh in range(self.n_corr_shells):
|
|
dim = self.corr_shells[icrsh]['dim']
|
|
n_orb = self.n_orbitals[ik, 0]
|
|
projmat = self.proj_mat[ik, 0, icrsh, 0:dim, 0:n_orb]
|
|
dens_mat[icrsh][
|
|
:, :] += numpy.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik]
|
|
|
|
if self.symm_op != 0:
|
|
dens_mat = self.symmcorr.symmetrize(dens_mat)
|
|
|
|
# Rotate to local coordinate system:
|
|
if self.use_rotations:
|
|
for icrsh in range(self.n_corr_shells):
|
|
if self.rot_mat_time_inv[icrsh] == 1:
|
|
dens_mat[icrsh] = dens_mat[icrsh].conjugate()
|
|
dens_mat[icrsh] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh]),
|
|
self.rot_mat[icrsh])
|
|
|
|
return dens_mat
|
|
|
|
def sorts_of_atoms(self, shells):
|
|
"""
|
|
Determine the number of inequivalent sorts.
|
|
"""
|
|
sortlst = [shells[i]['sort'] for i in range(len(shells))]
|
|
n_sorts = len(set(sortlst))
|
|
return n_sorts
|
|
|
|
def number_of_atoms(self, shells):
|
|
"""
|
|
Determine the number of inequivalent atoms.
|
|
"""
|
|
atomlst = [shells[i]['atom'] for i in range(len(shells))]
|
|
n_atoms = len(set(atomlst))
|
|
return n_atoms
|