LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-22 20:34:38 +01:00
Code Issues Releases Wiki Activity

improve performance of extract Gloc by factor 5! huge cleanup

Browse Source
This commit is contained in:
Alexander Hampel 2022-03-04 12:55:28 -05:00
parent a36f2fdc61
commit 66b3719c9b
1 changed files with 123 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

232
python/triqs_dft_tools/sumk_dft.py
View File

@ -25,7 +25,7 @@ General SumK class and helper functions for combining ab-initio code and triqs
""" """
from types import * from types import *
import numpy import numpy as np
import triqs.utility.dichotomy as dichotomy import triqs.utility.dichotomy as dichotomy
from triqs.gf import * from triqs.gf import *
import triqs.utility.mpi as mpi import triqs.utility.mpi as mpi
@ -57,7 +57,7 @@ class SumkDFT(object):
The value of magnetic field to add to the DFT Hamiltonian. The value of magnetic field to add to the DFT Hamiltonian.
The contribution -h_field*sigma is added to diagonal elements of the 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. It cannot be used with the spin-orbit coupling on; namely h_field is set to 0 if self.SO=True.
mesh: MeshImFreq or MeshImFreq, optional. Frequency mesh of Sigma. mesh: MeshImFreq or MeshReFreq, optional. Frequency mesh of Sigma.
beta : real, optional beta : real, optional
Inverse temperature. Used to construct imaginary frequency if mesh is not given. Inverse temperature. Used to construct imaginary frequency if mesh is not given.
n_iw : integer, optional n_iw : integer, optional
@ -105,13 +105,21 @@ class SumkDFT(object):
self.h_field = h_field self.h_field = h_field
if mesh is None: if mesh is None:
self.mesh = MeshImFreq(beta=beta, S='Fermion',n_max=n_iw) self.mesh = MeshImFreq(beta=beta, S='Fermion', n_max=n_iw)
elif isinstance(mesh, MeshImFreq) or isinstance(mesh, MeshReFreq): self.mesh_values = np.linspace(self.mesh(self.mesh.first_index()),
self.mesh(self.mesh.last_index()),
len(self.mesh))
elif isinstance(mesh, MeshImFreq):
self.mesh = mesh self.mesh = mesh
self.mesh_values = np.linspace(self.mesh(self.mesh.first_index()),
self.mesh(self.mesh.last_index()),
len(self.mesh))
elif isinstance(mesh, MeshReFreq):
self.mesh = mesh
self.mesh_values = np.linspace(self.mesh.omega_min, self.mesh.omega_max, len(self.mesh))
else: else:
raise ValueError('mesh must be a triqs mesh of type MeshImFreq or MeshReFreq') raise ValueError('mesh must be a triqs mesh of type MeshImFreq or MeshReFreq')
self.block_structure = BlockStructure() self.block_structure = BlockStructure()
# Read input from HDF: # Read input from HDF:
@ -537,18 +545,32 @@ class SumkDFT(object):
# Are we including Sigma? # Are we including Sigma?
if with_Sigma: if with_Sigma:
Sigma_imp = self.Sigma_imp Sigma_imp = self.Sigma_imp
sigma_minus_dc = [s.copy() for s in Sigma_imp]
if with_dc: if with_dc:
sigma_minus_dc = self.add_dc() sigma_minus_dc = self.add_dc()
else:
sigma_minus_dc = [s for s in Sigma_imp]
if not mesh is None: if not mesh is None:
warn('lattice_gf called with Sigma and given mesh. Mesh will be taken from Sigma.') warn('lattice_gf called with Sigma and given mesh. Mesh will be taken from Sigma.')
mesh = Sigma_imp[0].mesh
if isinstance(mesh, MeshReFreq) and broadening > 0 and mpi.is_master_node(): if self.mesh != Sigma_imp[0].mesh:
warn('lattice_gf called with Sigma and broadening > 0 (broadening = {}). You might want to explicitly set the broadening to 0.'.format(broadening)) mesh = Sigma_imp[0].mesh
if isinstance(mesh, MeshImFreq):
mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh))
else:
mesh_values = np.linspace(mesh.omega_min, mesh.omega_max, len(mesh))
else:
mesh = self.mesh
mesh_values = self.mesh_values
elif not mesh is None: elif not mesh is None:
assert isinstance(mesh, MreshReFreq) or isinstance(mesh, MeshImFreq), "mesh must be a triqs MeshReFreq or MeshImFreq" assert isinstance(mesh, MeshReFreq) or isinstance(mesh, MeshImFreq), "mesh must be a triqs MeshReFreq or MeshImFreq"
if isinstance(mesh, MeshImFreq):
mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh))
else:
mesh_values = np.linspace(mesh.omega_min, mesh.omega_max, len(mesh))
else: else:
mesh = self.mesh mesh = self.mesh
mesh_values = self.mesh_values
# Set up G_latt # Set up G_latt
if set_up_G_latt: if set_up_G_latt:
@ -567,28 +589,25 @@ class SumkDFT(object):
block_list=glist(), make_copies=False) block_list=glist(), make_copies=False)
G_latt.zero() G_latt.zero()
if isinstance(mesh, MeshImFreq): idmat = [np.identity(
G_latt << iOmega_n self.n_orbitals[ik, ntoi[sp]], np.complex_) for sp in spn]
else:
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]):
# fill Glatt
for ibl, (block, gf) in enumerate(G_latt):
ind = ntoi[spn[ibl]] ind = ntoi[spn[ibl]]
n_orb = self.n_orbitals[ik, ind] n_orb = self.n_orbitals[ik, ind]
M[ibl] = self.hopping[ik, ind, 0:n_orb, 0:n_orb] - \ if isinstance(mesh, MeshImFreq):
(idmat[ibl] * mu) - (idmat[ibl] * self.h_field * (1 - 2 * ibl)) gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl))
G_latt -= M - self.hopping[ik, ind, 0:n_orb, 0:n_orb])
else:
gf.data[:, :, :] = (idmat[ibl] *
(mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl) + 1j*broadening)
- self.hopping[ik, ind, 0:n_orb, 0:n_orb])
if with_Sigma: if with_Sigma:
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
for bname, gf in G_latt: gf -= self.upfold(ik, icrsh, block,
gf -= self.upfold(ik, icrsh, bname, sigma_minus_dc[icrsh][block], gf)
sigma_minus_dc[icrsh][bname], gf)
G_latt.invert() G_latt.invert()
self.G_latt = G_latt self.G_latt = G_latt
@ -649,7 +668,7 @@ class SumkDFT(object):
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 = numpy.array([i for i in self.mesh.values()]) mesh = np.array([i for i in self.mesh.values()])
if mesh[0] > (self.min_band_energy - self.chemical_potential) or mesh[-1] < (self.max_band_energy - self.chemical_potential): if mesh[0] > (self.min_band_energy - self.chemical_potential) or mesh[-1] < (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))
@ -754,7 +773,7 @@ class SumkDFT(object):
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
G_loc[icrsh].zero() # initialize to zero G_loc[icrsh].zero() # initialize to zero
ikarray = numpy.array(list(range(self.n_k))) ikarray = np.array(list(range(self.n_k)))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
if isinstance(self.mesh, MeshImFreq): if isinstance(self.mesh, MeshImFreq):
G_latt = self.lattice_gf( G_latt = self.lattice_gf(
@ -766,11 +785,8 @@ class SumkDFT(object):
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
# init temporary storage # init temporary storage
tmp = G_loc[icrsh].copy() for bname, gf in G_loc[icrsh]:
for bname, gf in tmp: gf += self.downfold(ik, icrsh, bname, G_latt[bname], gf)
tmp[bname] << self.downfold(
ik, icrsh, bname, G_latt[bname], gf)
G_loc[icrsh] += tmp
# Collect data from mpi # Collect data from mpi
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
@ -930,8 +946,8 @@ class SumkDFT(object):
dm = {} dm = {}
for block, block_dim in self.gf_struct_solver[ish].items(): for block, block_dim in self.gf_struct_solver[ish].items():
# get dm for the blocks: # get dm for the blocks:
dm[block] = numpy.zeros( dm[block] = np.zeros(
[block_dim, block_dim], numpy.complex_) [block_dim, block_dim], complex)
for ind1 in range(block_dim): for ind1 in range(block_dim):
for ind2 in range(block_dim): for ind2 in range(block_dim):
block_sumk, ind1_sumk = self.solver_to_sumk[ block_sumk, ind1_sumk = self.solver_to_sumk[
@ -992,7 +1008,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/numpy.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
@ -1002,15 +1018,15 @@ class SumkDFT(object):
else: else:
return w-w0 return w-w0
gf = [BlockGf(name_block_generator = [(name, GfReFreq( gf = [BlockGf(name_block_generator = [(name, GfReFreq(
window=(-numpy.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)),
numpy.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/numpy.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
@ -1061,7 +1077,7 @@ class SumkDFT(object):
for ish in range(self.n_inequiv_shells)] for ish in range(self.n_inequiv_shells)]
# the maximum value of each matrix element of each block and shell # the maximum value of each matrix element of each block and shell
max_gf = [{name:numpy.max(numpy.abs(g.data),0) for name, g in gf[ish]} for ish in range(self.n_inequiv_shells)] max_gf = [{name:np.max(np.abs(g.data),0) for name, g in gf[ish]} for ish in range(self.n_inequiv_shells)]
if include_shells is None: if include_shells is None:
# include all shells # include all shells
@ -1160,7 +1176,7 @@ class SumkDFT(object):
# helper function # helper function
def null(A, eps=1e-15): def null(A, eps=1e-15):
""" Calculate the null-space of matrix A """ """ Calculate the null-space of matrix A """
u, s, vh = numpy.linalg.svd(A) u, s, vh = np.linalg.svd(A)
null_mask = (s <= eps) null_mask = (s <= eps)
null_space = compress(null_mask, vh, axis=0) null_space = compress(null_mask, vh, axis=0)
return null_space.conjugate().transpose() return null_space.conjugate().transpose()
@ -1215,9 +1231,9 @@ class SumkDFT(object):
# one block to make it equal to the other (at least # one block to make it equal to the other (at least
# for tau=0). # for tau=0).
e1 = numpy.linalg.eigvalsh(gf1.data[0]) e1 = np.linalg.eigvalsh(gf1.data[0])
e2 = numpy.linalg.eigvalsh(gf2.data[0]) e2 = np.linalg.eigvalsh(gf2.data[0])
if numpy.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:
@ -1236,21 +1252,21 @@ 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 = numpy.kron(numpy.eye(*gf1.target_shape),gf2.data[0])-numpy.kron(gf1.data[0].transpose(),numpy.eye(*gf1.target_shape)) M = np.kron(np.eye(*gf1.target_shape),gf2.data[0])-np.kron(gf1.data[0].transpose(),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 = numpy.kron(numpy.eye(*gf1.target_shape),gf2.data[i])-numpy.kron(gf1.data[i].transpose(),numpy.eye(*gf1.target_shape)) M = np.kron(np.eye(*gf1.target_shape),gf2.data[i])-np.kron(gf1.data[i].transpose(),np.eye(*gf1.target_shape))
# transform M into current null space # transform M into current null space
M = numpy.dot(M, N) M = np.dot(M, N)
N = numpy.dot(N, null(M, threshold)) N = np.dot(N, null(M, threshold))
if numpy.size(N) == 0: if np.size(N) == 0:
break break
# no intersection of the null spaces -> no symmetry # no intersection of the null spaces -> no symmetry
if numpy.size(N) == 0: continue if np.size(N) == 0: continue
# reshape N: it then has the indices matrix, matrix, number of basis vectors of the null space # reshape N: it then has the indices matrix, matrix, number of basis vectors of the null space
N = N.reshape(gf1.target_shape[0], gf1.target_shape[1], -1).transpose([1, 0, 2]) N = N.reshape(gf1.target_shape[0], gf1.target_shape[1], -1).transpose([1, 0, 2])
@ -1268,7 +1284,7 @@ class SumkDFT(object):
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.
The object N[:,:,i] N[:,:,j] is a four-index object which we call Z. The object N[:,:,i] N[:,:,j] is a four-index object which we call Z.
""" """
Z = numpy.einsum('aci,bcj->abij', N, N.conjugate()) Z = np.einsum('aci,bcj->abij', N, N.conjugate())
""" """
function chi2 function chi2
@ -1278,16 +1294,16 @@ class SumkDFT(object):
""" """
def chi2(y): def chi2(y):
# reinterpret y as complex number # reinterpret y as complex number
y = y.view(numpy.complex_) y = y.view(np.complex_)
ret = 0.0 ret = 0.0
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 += numpy.abs(numpy.dot(y, numpy.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
res = minimize(chi2, numpy.ones(2 * N.shape[-1])) res = minimize(chi2, np.ones(2 * N.shape[-1]))
# if the minimization fails, there is probably no symmetry # if the minimization fails, there is probably no symmetry
if not res.success: continue if not res.success: continue
@ -1295,10 +1311,10 @@ class SumkDFT(object):
if res.fun > threshold: continue if res.fun > threshold: continue
# reinterpret the solution as a complex number # reinterpret the solution as a complex number
y = res.x.view(numpy.complex_) y = res.x.view(np.complex_)
# reconstruct the T matrix # reconstruct the T matrix
T = numpy.zeros(N.shape[:-1], dtype=numpy.complex_) T = np.zeros(N.shape[:-1], dtype=np.complex_)
for i in range(len(y)): for i in range(len(y)):
T += N[:, :, i] * y[i] T += N[:, :, i] * y[i]
@ -1336,9 +1352,9 @@ 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] = numpy.dot(T.conjugate().transpose(), v2[0].conjugate()), not v2[1] self.deg_shells[ish][ind2][block1] = np.dot(T.conjugate().transpose(), v2[0].conjugate()), not v2[1]
else: else:
self.deg_shells[ish][ind2][block1] = numpy.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
# set v2 and C2 so that they are compatible with # set v2 and C2 so that they are compatible with
# C(T gf1 T^dagger) = gf2 # C(T gf1 T^dagger) = gf2
@ -1346,9 +1362,9 @@ class SumkDFT(object):
# C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2) # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2)
elif (ind1 >= 0) and (ind2 < 0): elif (ind1 >= 0) and (ind2 < 0):
if conjugate: if conjugate:
self.deg_shells[ish][ind1][block2] = numpy.dot(T.conjugate(), v1[0].conjugate()), not v1[1] self.deg_shells[ish][ind1][block2] = np.dot(T.conjugate(), v1[0].conjugate()), not v1[1]
else: else:
self.deg_shells[ish][ind1][block2] = numpy.dot(T, v1[0]), v1[1] self.deg_shells[ish][ind1][block2] = np.dot(T, v1[0]), v1[1]
# the blocks are not already present # the blocks are not already present
# we arbitrarily choose v1=eye and C1=False and # we arbitrarily choose v1=eye and C1=False and
# set v2 and C2 so that they are compatible with # set v2 and C2 so that they are compatible with
@ -1357,7 +1373,7 @@ class SumkDFT(object):
# C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2) # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2)
elif (ind1 < 0) and (ind2 < 0): elif (ind1 < 0) and (ind2 < 0):
d = dict() d = dict()
d[block1] = numpy.eye(*gf1.target_shape), False d[block1] = np.eye(*gf1.target_shape), False
if conjugate: if conjugate:
d[block2] = T.conjugate(), True d[block2] = T.conjugate(), True
else: else:
@ -1418,7 +1434,7 @@ class SumkDFT(object):
"calculate_diagonalization_matrix: Choices for prop_to_be_diagonal are 'eal' or 'dm'.") "calculate_diagonalization_matrix: Choices for prop_to_be_diagonal are 'eal' or 'dm'.")
return 0 return 0
trans = [{block: numpy.eye(block_dim) for block, block_dim in gfs} for gfs in self.gf_struct_sumk] trans = [{block: np.eye(block_dim) for block, block_dim in gfs} for gfs in self.gf_struct_sumk]
for ish in shells: for ish in shells:
trafo = {} trafo = {}
@ -1427,10 +1443,10 @@ 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 numpy.sum(abs(prop[ish][name]-numpy.diag(numpy.diagonal(prop[ish][name])))) > 1e-13: if np.sum(abs(prop[ish][name]-np.diag(np.diagonal(prop[ish][name])))) > 1e-13:
trafo[name] = numpy.linalg.eigh(prop[ish][name])[1].conj().T trafo[name] = np.linalg.eigh(prop[ish][name])[1].conj().T
else: else:
trafo[name] = numpy.identity(numpy.shape(prop[ish][name])[0]) trafo[name] = np.identity(np.shape(prop[ish][name])[0])
# Transform back to sumk if necessay, blocks change in this step! # Transform back to sumk if necessay, blocks change in this step!
if calc_in_solver_blocks: if calc_in_solver_blocks:
trafo = self.block_structure.convert_matrix(trafo, space_from='solver', space_to='sumk') trafo = self.block_structure.convert_matrix(trafo, space_from='solver', space_to='sumk')
@ -1467,10 +1483,10 @@ class SumkDFT(object):
dens_mat = [{} for icrsh in range(self.n_corr_shells)] dens_mat = [{} for icrsh in range(self.n_corr_shells)]
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']]: for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
dens_mat[icrsh][sp] = numpy.zeros( dens_mat[icrsh][sp] = np.zeros(
[self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_) [self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], np.complex_)
ikarray = numpy.array(list(range(self.n_k))) ikarray = np.array(list(range(self.n_k)))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
if method == "using_gf": if method == "using_gf":
@ -1485,7 +1501,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]
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 = [numpy.zeros([dims[sp], dims[sp]], numpy.complex_) for sp in spn] MMat = [np.zeros([dims[sp], dims[sp]], np.complex_) for sp in spn]
for isp, sp in enumerate(spn): for isp, sp in enumerate(spn):
ind = ntoi[sp] ind = ntoi[sp]
@ -1508,10 +1524,10 @@ 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]
if method == "using_gf": if method == "using_gf":
dens_mat[icrsh][sp] += numpy.dot(numpy.dot(projmat, MMat[isp]), dens_mat[icrsh][sp] += np.dot(np.dot(projmat, MMat[isp]),
projmat.transpose().conjugate()) projmat.transpose().conjugate())
elif method == "using_point_integration": elif method == "using_point_integration":
dens_mat[icrsh][sp] += self.bz_weights[ik] * numpy.dot(numpy.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:
@ -1530,7 +1546,7 @@ 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] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh][sp]), dens_mat[icrsh][sp] = np.dot(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
@ -1565,8 +1581,8 @@ class SumkDFT(object):
eff_atlevels = [{} for ish in range(self.n_inequiv_shells)] eff_atlevels = [{} 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 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']]:
eff_atlevels[ish][sp] = numpy.identity( eff_atlevels[ish][sp] = np.identity(
self.corr_shells[self.inequiv_to_corr[ish]]['dim'], numpy.complex_) self.corr_shells[self.inequiv_to_corr[ish]]['dim'], np.complex_)
eff_atlevels[ish][sp] *= -self.chemical_potential eff_atlevels[ish][sp] *= -self.chemical_potential
eff_atlevels[ish][ eff_atlevels[ish][
sp] -= self.dc_imp[self.inequiv_to_corr[ish]][sp] sp] -= self.dc_imp[self.inequiv_to_corr[ish]][sp]
@ -1579,18 +1595,18 @@ class SumkDFT(object):
for icrsh in range(self.n_corr_shells): for icrsh in range(self.n_corr_shells):
dim = self.corr_shells[icrsh]['dim'] dim = self.corr_shells[icrsh]['dim']
for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]: for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]:
self.Hsumk[icrsh][sp] = numpy.zeros( self.Hsumk[icrsh][sp] = np.zeros(
[dim, dim], numpy.complex_) [dim, dim], np.complex_)
for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]): for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]):
ind = self.spin_names_to_ind[ ind = self.spin_names_to_ind[
self.corr_shells[icrsh]['SO']][sp] self.corr_shells[icrsh]['SO']][sp]
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 = numpy.identity(n_orb, numpy.complex_) MMat = np.identity(n_orb, np.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] * numpy.dot(numpy.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:
@ -1603,7 +1619,7 @@ 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] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), self.Hsumk[icrsh][sp]), self.Hsumk[icrsh][sp] = np.dot(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:
@ -1628,7 +1644,7 @@ class SumkDFT(object):
dim = self.corr_shells[icrsh]['dim'] dim = self.corr_shells[icrsh]['dim']
spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] spn = self.spin_block_names[self.corr_shells[icrsh]['SO']]
for sp in spn: for sp in spn:
self.dc_imp[icrsh][sp] = numpy.zeros([dim, dim], numpy.float_) self.dc_imp[icrsh][sp] = np.zeros([dim, dim], np.float_)
self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)] self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)]
def set_dc(self, dc_imp, dc_energ): def set_dc(self, dc_imp, dc_energ):
@ -1706,7 +1722,7 @@ class SumkDFT(object):
Ncr[bl] += dens_mat[block].real.trace() Ncr[bl] += dens_mat[block].real.trace()
Ncrtot = sum(Ncr.values()) Ncrtot = sum(Ncr.values())
for sp in spn: for sp in spn:
self.dc_imp[icrsh][sp] = numpy.identity(dim, numpy.float_) self.dc_imp[icrsh][sp] = np.identity(dim, np.float_)
if self.SP == 0: # average the densities if there is no SP: if self.SP == 0: # average the densities if there is no SP:
Ncr[sp] = Ncrtot / len(spn) Ncr[sp] = Ncrtot / len(spn)
# correction for SO: we have only one block in this case, but # correction for SO: we have only one block in this case, but
@ -1771,8 +1787,8 @@ class SumkDFT(object):
if transform: if transform:
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] = numpy.dot(T.conjugate().transpose(), self.dc_imp[icrsh][sp] = np.dot(T.conjugate().transpose(),
numpy.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"""
@ -1798,12 +1814,10 @@ 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:
dccont = numpy.dot(self.rot_mat[icrsh], numpy.dot(self.dc_imp[icrsh][ gf -= np.dot(self.rot_mat[icrsh], np.dot(self.dc_imp[icrsh][
bname], self.rot_mat[icrsh].conjugate().transpose())) bname], self.rot_mat[icrsh].conjugate().transpose()))
else: else:
dccont = self.dc_imp[icrsh][bname] gf -= self.dc_imp[icrsh][bname]
sigma_minus_dc[icrsh][bname] -= dccont
return sigma_minus_dc return sigma_minus_dc
@ -1844,7 +1858,7 @@ class SumkDFT(object):
v, C = degsh[key] v, C = degsh[key]
else: else:
# for backward compatibility, allow degsh to be a list # for backward compatibility, allow degsh to be a list
v = numpy.eye(*ss.target_shape) v = np.eye(*ss.target_shape)
C = False C = False
# the helper is in the basis where the blocks are all equal # the helper is in the basis where the blocks are all equal
helper.from_L_G_R(v.conjugate().transpose(), gf_to_symm[key], v) helper.from_L_G_R(v.conjugate().transpose(), gf_to_symm[key], v)
@ -1858,7 +1872,7 @@ class SumkDFT(object):
v, C = degsh[key] v, C = degsh[key]
else: else:
# for backward compatibility, allow degsh to be a list # for backward compatibility, allow degsh to be a list
v = numpy.eye(*ss.target_shape) v = np.eye(*ss.target_shape)
C = False C = False
if C: if C:
gf_to_symm[key].from_L_G_R(v, ss.transpose().copy(), v.conjugate().transpose()) gf_to_symm[key].from_L_G_R(v, ss.transpose().copy(), v.conjugate().transpose())
@ -1915,7 +1929,7 @@ class SumkDFT(object):
if mu is None: if mu is None:
mu = self.chemical_potential mu = self.chemical_potential
dens = 0.0 dens = 0.0
ikarray = numpy.array(list(range(self.n_k))) ikarray = np.array(list(range(self.n_k)))
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
G_latt = self.lattice_gf( G_latt = self.lattice_gf(
ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening) ik=ik, mu=mu, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening)
@ -2046,16 +2060,16 @@ class SumkDFT(object):
# Convert Fermi weights to a density matrix # Convert Fermi weights to a density matrix
dens_mat_dft = {} dens_mat_dft = {}
for sp in spn: for sp in spn:
dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(numpy.complex_) for ik in range(self.n_k)] dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(np.complex_) for ik in range(self.n_k)]
# Set up deltaN: # Set up deltaN:
deltaN = {} deltaN = {}
for sp in spn: for sp in spn:
deltaN[sp] = [numpy.zeros([self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[ deltaN[sp] = [np.zeros([self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[
ik, ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)] ik, ntoi[sp]]], np.complex_) for ik in range(self.n_k)]
ikarray = numpy.arange(self.n_k) ikarray = np.arange(self.n_k)
for ik in mpi.slice_array(ikarray): for ik in mpi.slice_array(ikarray):
G_latt_iw = self.lattice_gf( G_latt_iw = self.lattice_gf(
ik=ik, mu=self.chemical_potential, iw_or_w="iw") ik=ik, mu=self.chemical_potential, iw_or_w="iw")
@ -2075,11 +2089,11 @@ class SumkDFT(object):
if dm_type in ['vasp','qe']: if dm_type in ['vasp','qe']:
# In 'vasp'-mode subtract the DFT density matrix # In 'vasp'-mode subtract the DFT density matrix
nb = self.n_orbitals[ik, ntoi[bname]] nb = self.n_orbitals[ik, ntoi[bname]]
diag_inds = numpy.diag_indices(nb) diag_inds = np.diag_indices(nb)
deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik][:nb] deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik][:nb]
if self.charge_mixing and self.deltaNOld is not None: if self.charge_mixing and self.deltaNOld is not None:
G2 = numpy.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][diag_inds] + mix_fac * deltaN[bname][ik][diag_inds] deltaN[bname][ik][diag_inds] = (1.0 - mix_fac) * self.deltaNOld[bname][ik][diag_inds] + mix_fac * deltaN[bname][ik][diag_inds]
@ -2088,7 +2102,7 @@ class SumkDFT(object):
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 += numpy.dot(deltaN[bname][ik], self.hopping[ik, isp, :nb, :nb]).trace().real * self.bz_weights[ik] band_en_correction += np.dot(deltaN[bname][ik], self.hopping[ik, isp, :nb, :nb]).trace().real * self.bz_weights[ik]
# mpi reduce: # mpi reduce:
for bname in deltaN: for bname in deltaN:
@ -2157,9 +2171,9 @@ class SumkDFT(object):
fout.close() fout.close()
elif dm_type == 'vasp': elif dm_type == 'vasp':
if kpts_to_write is None: if kpts_to_write is None:
kpts_to_write = numpy.arange(self.n_k) kpts_to_write = np.arange(self.n_k)
else: else:
assert numpy.min(kpts_to_write) >= 0 and numpy.max(kpts_to_write) < self.n_k assert np.min(kpts_to_write) >= 0 and np.max(kpts_to_write) < self.n_k
assert self.SP == 0, "Spin-polarized density matrix is not implemented" assert self.SP == 0, "Spin-polarized density matrix is not implemented"
@ -2188,7 +2202,7 @@ class SumkDFT(object):
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 = numpy.max(self.n_orbitals) nbmax = np.max(self.n_orbitals)
# output beta and mu in Hartrees # output beta and mu in Hartrees
beta = G_latt_iw.mesh.beta * self.energy_unit beta = G_latt_iw.mesh.beta * self.energy_unit
mu = self.chemical_potential/self.energy_unit mu = self.chemical_potential/self.energy_unit
@ -2223,7 +2237,7 @@ class SumkDFT(object):
assert self.SP == 0, "Spin-polarized density matrix is not implemented" assert self.SP == 0, "Spin-polarized density matrix is not implemented"
subgrp = 'dft_update' subgrp = 'dft_update'
delta_N = numpy.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]
@ -2255,10 +2269,10 @@ class SumkDFT(object):
def calculate_min_max_band_energies(self): def calculate_min_max_band_energies(self):
hop = self.hopping hop = self.hopping
diag_hop = numpy.zeros(hop.shape[:-1]) diag_hop = np.zeros(hop.shape[:-1])
hop_slice = mpi.slice_array(hop) hop_slice = mpi.slice_array(hop)
diag_hop_slice = mpi.slice_array(diag_hop) diag_hop_slice = mpi.slice_array(diag_hop)
diag_hop_slice[:] = numpy.linalg.eigvalsh(hop_slice) diag_hop_slice[:] = np.linalg.eigvalsh(hop_slice)
diag_hop = mpi.all_reduce(mpi.world, diag_hop, lambda x, y: x + y) diag_hop = mpi.all_reduce(mpi.world, diag_hop, lambda x, y: x + y)
min_band_energy = diag_hop.min().real min_band_energy = diag_hop.min().real
max_band_energy = diag_hop.max().real max_band_energy = diag_hop.max().real
@ -2297,7 +2311,7 @@ class SumkDFT(object):
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
specifically that it matches DFT.""" specifically that it matches DFT."""
dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_) dens_mat = [np.zeros([self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], np.complex_)
for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells)]
for ik in range(self.n_k): for ik in range(self.n_k):
@ -2306,7 +2320,7 @@ class SumkDFT(object):
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][
:, :] += numpy.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)
@ -2316,7 +2330,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] = numpy.dot(numpy.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