From 796e05ea646b392955ec14e0c61c93b344dde1ab Mon Sep 17 00:00:00 2001 From: Manuel Date: Mon, 12 Aug 2019 19:23:27 -0400 Subject: [PATCH] Rewrite calculate_diagonalization_matrix --- python/block_structure.py | 2 +- python/sumk_dft.py | 81 ++++++++++++++++++--------------------- 2 files changed, 39 insertions(+), 44 deletions(-) diff --git a/python/block_structure.py b/python/block_structure.py index 5739482c..4b40a69a 100644 --- a/python/block_structure.py +++ b/python/block_structure.py @@ -90,7 +90,7 @@ class BlockStructure(object): inequivalent correlated shell is given transformation : list of numpy.array or list of dict a list with entries for each ``ish`` giving transformation matrices - that are used on the Green's function in ``sumk`` space when before + that are used on the Green's function in ``sumk`` space before converting to the ``solver`` space Up to the change in block structure, diff --git a/python/sumk_dft.py b/python/sumk_dft.py index ef7f7dd8..bdbccddd 100644 --- a/python/sumk_dft.py +++ b/python/sumk_dft.py @@ -1330,7 +1330,7 @@ class SumkDFT(object): # a block was found, break out of the loop break - def calculate_diagonalization_matrix(self, prop_to_be_diagonal='eal', calc_in_solver_blocks=False, write_to_blockstructure = True, ish=0): + def calculate_diagonalization_matrix(self, prop_to_be_diagonal='eal', calc_in_solver_blocks=False, write_to_blockstructure = True, shells=None): """ Calculates the diagonalisation matrix, and (optionally) stores it in the BlockStructure. @@ -1349,62 +1349,57 @@ class SumkDFT(object): write_to_blockstructure : bool, optional Whether the diagonalization matrix shall be written to the BlockStructure directly. - ish : int, optional - Number of the correlated shell to be diagonalized. + shells : list of int, optional + Indices of correlated shells to be diagonalized. + None: all shells Returns ------- trafo : dict The transformation matrix for each spin-block in the correlated shell - """ - trafo = {} - - if prop_to_be_diagonal == 'eal': - prop = self.eff_atomic_levels()[ish] - elif prop_to_be_diagonal == 'dm': - prop = self.density_matrix(method='using_point_integration')[ish] - else: - mpi.report( - "calculate_diagonalization_matrix: not a valid quantitiy to be diagonal. Choices are 'eal' or 'dm'.") + # Use all shells + if shells is None: + shells = range(self.n_corr_shells) + elif max(shells) >= self.n_corr_shells: # Check if the shell indices are present + mpi.report("calculate_diagonalization_matrix: shells not correct.") return 0 - if calc_in_solver_blocks: - trafo_tmp = self.block_structure.transformation - self.block_structure.transformation = None - - prop_solver = self.block_structure.convert_matrix(prop, space_from='sumk', space_to='solver') - t= {} - for name in prop_solver: - t[name] = numpy.linalg.eigh(prop_solver[name])[1].conjugate().transpose() - trafo = self.block_structure.convert_matrix(t, space_from='solver', space_to='sumk') - #self.T = numpy.dot(self.T.transpose().conjugate(), - # self.w).conjugate().transpose() - self.block_structure.transformation = trafo_tmp + if prop_to_be_diagonal == 'eal': + prop = [self.eff_atomic_levels()[self.corr_to_inequiv[ish]] + for ish in range(self.n_corr_shells)] + elif prop_to_be_diagonal == 'dm': + prop = self.density_matrix(method='using_point_integration') else: - for name in prop: - t = numpy.linalg.eigh(prop[name])[1].conjugate().transpose() - trafo[name] = t - # calculate new Transformation matrix - #self.T = numpy.dot(self.T.transpose().conjugate(), - # self.w).conjugate().transpose() + mpi.report( + "calculate_diagonalization_matrix: Choices for prop_to_be_diagonal are 'eal' or 'dm'.") + return 0 - # measure for the 'unity' of the transformation: - #wsqr = sum(abs(self.w.diagonal())**2) / self.w.diagonal().size - #return wsqr + trans = [{block: numpy.eye(len(indices)) for block, indices in gfs} for gfs in self.gf_struct_sumk] + + for ish in shells: + trafo = {} + # Transform to solver basis if desired, blocks of prop change in this step! + if calc_in_solver_blocks: + 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 abs(numpy.sum(prop[ish][name]-numpy.diag(numpy.diagonal(prop[ish][name])))) > 1e-13: + trafo[name] = numpy.linalg.eigh(prop[ish][name])[1].conj().T + else: + trafo[name] = numpy.identity(numpy.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') + trans[ish] = trafo + + # Write to block_structure object if write_to_blockstructure: - if self.block_structure.transformation == None: - self.block_structure.transformation = [{} 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']]: - self.block_structure.transformation[icrsh][sp] = numpy.eye(self.corr_shells[icrsh]['dim'], dtype=numpy.complex_) + self.block_structure.transformation = trans - - self.block_structure.transformation[ish] = trafo - - return trafo + return trans def density_matrix(self, method='using_gf', beta=40.0):