diff --git a/python/triqs_dft_tools/sumk_dft.py b/python/triqs_dft_tools/sumk_dft.py index 72d8cd17..fa3a8c48 100644 --- a/python/triqs_dft_tools/sumk_dft.py +++ b/python/triqs_dft_tools/sumk_dft.py @@ -25,7 +25,7 @@ General SumK class and helper functions for combining ab-initio code and triqs """ from types import * -import numpy +import numpy as np import triqs.utility.dichotomy as dichotomy from triqs.gf import * 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 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. - mesh: MeshImFreq or MeshImFreq, optional. Frequency mesh of Sigma. + mesh: MeshImFreq or MeshReFreq, optional. Frequency mesh of Sigma. beta : real, optional Inverse temperature. Used to construct imaginary frequency if mesh is not given. n_iw : integer, optional @@ -105,13 +105,21 @@ class SumkDFT(object): self.h_field = h_field if mesh is None: - self.mesh = MeshImFreq(beta=beta, S='Fermion',n_max=n_iw) - elif isinstance(mesh, MeshImFreq) or isinstance(mesh, MeshReFreq): + self.mesh = MeshImFreq(beta=beta, S='Fermion', n_max=n_iw) + 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_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: raise ValueError('mesh must be a triqs mesh of type MeshImFreq or MeshReFreq') - self.block_structure = BlockStructure() # Read input from HDF: @@ -537,18 +545,32 @@ class SumkDFT(object): # Are we including Sigma? if with_Sigma: Sigma_imp = self.Sigma_imp - sigma_minus_dc = [s.copy() for s in Sigma_imp] if with_dc: sigma_minus_dc = self.add_dc() + else: + sigma_minus_dc = [s for s in Sigma_imp] if not mesh is None: 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(): - warn('lattice_gf called with Sigma and broadening > 0 (broadening = {}). You might want to explicitly set the broadening to 0.'.format(broadening)) + + if self.mesh != Sigma_imp[0].mesh: + 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: - 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: mesh = self.mesh + mesh_values = self.mesh_values # Set up G_latt if set_up_G_latt: @@ -567,28 +589,25 @@ class SumkDFT(object): block_list=glist(), make_copies=False) G_latt.zero() - if isinstance(mesh, MeshImFreq): - G_latt << iOmega_n - 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]): + idmat = [np.identity( + self.n_orbitals[ik, ntoi[sp]], np.complex_) for sp in spn] + # fill Glatt + for ibl, (block, gf) in enumerate(G_latt): 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 isinstance(mesh, MeshImFreq): + 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]) + 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: - 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) + if with_Sigma: + for icrsh in range(self.n_corr_shells): + gf -= self.upfold(ik, icrsh, block, + sigma_minus_dc[icrsh][block], gf) G_latt.invert() self.G_latt = G_latt @@ -649,7 +668,7 @@ class SumkDFT(object): if isinstance(self.mesh, MeshReFreq): if self.min_band_energy is None or self.max_band_energy is None: 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): 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): 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): if isinstance(self.mesh, MeshImFreq): G_latt = self.lattice_gf( @@ -766,11 +785,8 @@ class SumkDFT(object): 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 + for bname, gf in G_loc[icrsh]: + gf += self.downfold(ik, icrsh, bname, G_latt[bname], gf) # Collect data from mpi for icrsh in range(self.n_corr_shells): @@ -930,8 +946,8 @@ class SumkDFT(object): dm = {} for block, block_dim in self.gf_struct_solver[ish].items(): # get dm for the blocks: - dm[block] = numpy.zeros( - [block_dim, block_dim], numpy.complex_) + dm[block] = np.zeros( + [block_dim, block_dim], complex) for ind1 in range(block_dim): for ind2 in range(block_dim): block_sumk, ind1_sumk = self.solver_to_sumk[ @@ -992,7 +1008,7 @@ class SumkDFT(object): gf = [g_sh.copy() for g_sh in G] for ish in range(len(gf)): 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): def get_delta_from_mesh(mesh): w0 = None @@ -1002,15 +1018,15 @@ class SumkDFT(object): else: return w-w0 gf = [BlockGf(name_block_generator = [(name, GfReFreq( - window=(-numpy.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))), + 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))), n_points=len(block.mesh), indices=block.indices)) for name, block in g_sh], make_copies=False) for g_sh in G] for ish in range(len(gf)): for name, g in gf[ish]: 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: raise Exception("G must be a list of BlockGf of either GfImFreq, GfImTime, GfReFreq or GfReTime") return gf @@ -1061,7 +1077,7 @@ class SumkDFT(object): for ish in range(self.n_inequiv_shells)] # 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: # include all shells @@ -1160,7 +1176,7 @@ class SumkDFT(object): # helper function def null(A, eps=1e-15): """ 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_space = compress(null_mask, vh, axis=0) return null_space.conjugate().transpose() @@ -1215,9 +1231,9 @@ class SumkDFT(object): # one block to make it equal to the other (at least # for tau=0). - e1 = numpy.linalg.eigvalsh(gf1.data[0]) - e2 = numpy.linalg.eigvalsh(gf2.data[0]) - if numpy.any(abs(e1-e2) > threshold): continue + e1 = np.linalg.eigvalsh(gf1.data[0]) + e2 = np.linalg.eigvalsh(gf2.data[0]) + if np.any(abs(e1-e2) > threshold): continue for conjugate in [False,True]: if conjugate: @@ -1236,21 +1252,21 @@ class SumkDFT(object): # product to get a linear problem, which consists # 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) # now we get the intersection of the null spaces # of all values of tau 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 - M = numpy.dot(M, N) - N = numpy.dot(N, null(M, threshold)) - if numpy.size(N) == 0: + M = np.dot(M, N) + N = np.dot(N, null(M, threshold)) + if np.size(N) == 0: break # 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 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. 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 @@ -1278,16 +1294,16 @@ class SumkDFT(object): """ def chi2(y): # reinterpret y as complex number - y = y.view(numpy.complex_) + y = y.view(np.complex_) ret = 0.0 for a in range(Z.shape[0]): 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 return ret # 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 not res.success: continue @@ -1295,10 +1311,10 @@ class SumkDFT(object): if res.fun > threshold: continue # reinterpret the solution as a complex number - y = res.x.view(numpy.complex_) + y = res.x.view(np.complex_) # 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)): T += N[:, :, i] * y[i] @@ -1336,9 +1352,9 @@ class SumkDFT(object): # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2) if (ind1 < 0) and (ind2 >= 0): 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: - 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 # set v2 and C2 so that they are compatible with # C(T gf1 T^dagger) = gf2 @@ -1346,9 +1362,9 @@ class SumkDFT(object): # C1(v1^dagger G1 v1) = C2(v2^dagger G2 v2) elif (ind1 >= 0) and (ind2 < 0): 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: - 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 # we arbitrarily choose v1=eye and C1=False and # 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) elif (ind1 < 0) and (ind2 < 0): d = dict() - d[block1] = numpy.eye(*gf1.target_shape), False + d[block1] = np.eye(*gf1.target_shape), False if conjugate: d[block2] = T.conjugate(), True else: @@ -1418,7 +1434,7 @@ class SumkDFT(object): "calculate_diagonalization_matrix: Choices for prop_to_be_diagonal are 'eal' or 'dm'.") 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: trafo = {} @@ -1427,10 +1443,10 @@ class SumkDFT(object): prop[ish] = self.block_structure.convert_matrix(prop[ish], space_from='sumk', space_to='solver') # Get diagonalisation matrix, if not already diagonal for name in prop[ish]: - if numpy.sum(abs(prop[ish][name]-numpy.diag(numpy.diagonal(prop[ish][name])))) > 1e-13: - trafo[name] = numpy.linalg.eigh(prop[ish][name])[1].conj().T + 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 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! if calc_in_solver_blocks: 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)] 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_) + dens_mat[icrsh][sp] = np.zeros( + [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): if method == "using_gf": @@ -1485,7 +1501,7 @@ class SumkDFT(object): ntoi = self.spin_names_to_ind[self.SO] spn = self.spin_block_names[self.SO] dims = {sp:self.n_orbitals[ik, ntoi[sp]] for sp in spn} - MMat = [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): ind = ntoi[sp] @@ -1508,10 +1524,10 @@ class SumkDFT(object): 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]), + dens_mat[icrsh][sp] += np.dot(np.dot(projmat, MMat[isp]), projmat.transpose().conjugate()) elif method == "using_point_integration": - dens_mat[icrsh][sp] += self.bz_weights[ik] * 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()) # get data from nodes: @@ -1530,7 +1546,7 @@ class SumkDFT(object): 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]), + dens_mat[icrsh][sp] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh][sp]), self.rot_mat[icrsh]) return dens_mat @@ -1565,8 +1581,8 @@ class SumkDFT(object): 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] = np.identity( + self.corr_shells[self.inequiv_to_corr[ish]]['dim'], np.complex_) eff_atlevels[ish][sp] *= -self.chemical_potential eff_atlevels[ish][ 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): 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_) + self.Hsumk[icrsh][sp] = np.zeros( + [dim, dim], np.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 = np.identity(n_orb, np.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), + self.Hsumk[icrsh][sp] += self.bz_weights[ik] * np.dot(np.dot(projmat, MMat), projmat.conjugate().transpose()) # symmetrisation: if self.symm_op != 0: @@ -1603,7 +1619,7 @@ class SumkDFT(object): 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.Hsumk[icrsh][sp] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), self.Hsumk[icrsh][sp]), self.rot_mat[icrsh]) # add to matrix: @@ -1628,7 +1644,7 @@ class SumkDFT(object): 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_imp[icrsh][sp] = np.zeros([dim, dim], np.float_) self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)] def set_dc(self, dc_imp, dc_energ): @@ -1706,7 +1722,7 @@ class SumkDFT(object): Ncr[bl] += dens_mat[block].real.trace() Ncrtot = sum(Ncr.values()) 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: Ncr[sp] = Ncrtot / len(spn) # correction for SO: we have only one block in this case, but @@ -1771,8 +1787,8 @@ class SumkDFT(object): if transform: for sp in spn: T = self.block_structure.effective_transformation_sumk[icrsh][sp] - self.dc_imp[icrsh][sp] = numpy.dot(T.conjugate().transpose(), - numpy.dot(self.dc_imp[icrsh][sp], T)) + self.dc_imp[icrsh][sp] = np.dot(T.conjugate().transpose(), + np.dot(self.dc_imp[icrsh][sp], T)) def add_dc(self): r""" @@ -1798,12 +1814,10 @@ class SumkDFT(object): for bname, gf in sigma_minus_dc[icrsh]: # Transform dc_imp to global coordinate system if self.use_rotations: - dccont = numpy.dot(self.rot_mat[icrsh], numpy.dot(self.dc_imp[icrsh][ - bname], self.rot_mat[icrsh].conjugate().transpose())) + gf -= np.dot(self.rot_mat[icrsh], np.dot(self.dc_imp[icrsh][ + bname], self.rot_mat[icrsh].conjugate().transpose())) else: - dccont = self.dc_imp[icrsh][bname] - - sigma_minus_dc[icrsh][bname] -= dccont + gf -= self.dc_imp[icrsh][bname] return sigma_minus_dc @@ -1844,7 +1858,7 @@ class SumkDFT(object): v, C = degsh[key] else: # for backward compatibility, allow degsh to be a list - v = numpy.eye(*ss.target_shape) + v = np.eye(*ss.target_shape) C = False # 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) @@ -1858,7 +1872,7 @@ class SumkDFT(object): v, C = degsh[key] else: # for backward compatibility, allow degsh to be a list - v = numpy.eye(*ss.target_shape) + v = np.eye(*ss.target_shape) C = False if C: 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: mu = self.chemical_potential 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): G_latt = self.lattice_gf( 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 dens_mat_dft = {} 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: 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)] + deltaN[sp] = [np.zeros([self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[ + 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): G_latt_iw = self.lattice_gf( ik=ik, mu=self.chemical_potential, iw_or_w="iw") @@ -2075,11 +2089,11 @@ class SumkDFT(object): if dm_type in ['vasp','qe']: # In 'vasp'-mode subtract the DFT density matrix 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] 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 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] @@ -2088,7 +2102,7 @@ class SumkDFT(object): b1, b2 = band_window[isp][ik, :2] nb = b2 - b1 + 1 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: for bname in deltaN: @@ -2157,9 +2171,9 @@ class SumkDFT(object): fout.close() elif dm_type == 'vasp': if kpts_to_write is None: - kpts_to_write = numpy.arange(self.n_k) + kpts_to_write = np.arange(self.n_k) 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" @@ -2188,7 +2202,7 @@ class SumkDFT(object): with open(filename, 'w') as f: # determine the number of spin blocks 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 beta = G_latt_iw.mesh.beta * 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" 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)) for ik in range(self.n_k): ib1 = band_window[0][ik, 0] @@ -2255,10 +2269,10 @@ class SumkDFT(object): def calculate_min_max_band_energies(self): hop = self.hopping - diag_hop = numpy.zeros(hop.shape[:-1]) + diag_hop = np.zeros(hop.shape[:-1]) hop_slice = mpi.slice_array(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) min_band_energy = diag_hop.min().real max_band_energy = diag_hop.max().real @@ -2297,7 +2311,7 @@ class SumkDFT(object): 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_) + 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 ik in range(self.n_k): @@ -2306,7 +2320,7 @@ class SumkDFT(object): 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] + :, :] += np.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik] if self.symm_op != 0: dens_mat = self.symmcorr.symmetrize(dens_mat) @@ -2316,7 +2330,7 @@ class SumkDFT(object): 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]), + dens_mat[icrsh] = np.dot(np.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh]), self.rot_mat[icrsh]) return dens_mat