3
0
mirror of https://github.com/triqs/dft_tools synced 2024-12-22 12:23:41 +01:00

fixed analyze_block_structure in sumk

was buggy when the number of off-diagonal elements was larger than the number of orbitals
This commit is contained in:
pdelange 2015-06-26 18:47:20 +02:00 committed by Priyanka Seth
parent cb792604b1
commit 6eef3bd172

View File

@ -27,6 +27,7 @@ from pytriqs.gf.local import *
import pytriqs.utility.mpi as mpi import pytriqs.utility.mpi as mpi
from pytriqs.archive import * from pytriqs.archive import *
from symmetry import * from symmetry import *
from sets import Set
class SumkDFT: class SumkDFT:
"""This class provides a general SumK method for combining ab-initio code and pytriqs.""" """This class provides a general SumK method for combining ab-initio code and pytriqs."""
@ -655,29 +656,27 @@ class SumkDFT:
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']]:
n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim']
dmbool = (abs(dens_mat[ish][sp]) > threshold) # gives an index list of entries larger that threshold dmbool = (abs(dens_mat[ish][sp]) > threshold) # gives an index list of entries larger that threshold
# Determine off-diagonal entries in upper triangular part of density matrix # Determine off-diagonal entries in upper triangular part of density matrix
offdiag = [] offdiag = Set([])
for i in range(len(dmbool)): for i in range(n_orb):
for j in range(i+1,len(dmbool)): for j in range(i+1,n_orb):
if dmbool[i,j]: offdiag.append([i,j]) if dmbool[i,j]: offdiag.add((i,j))
# Determine the number of non-hybridising blocks in the gf # Determine the number of non-hybridising blocks in the gf
num_blocs = len(dmbool) blocs = [ [i] for i in range(n_orb) ]
blocs = [ [i] for i in range(num_blocs) ] while len(offdiag) != 0:
for i in range(len(offdiag)): pair = offdiag.pop()
for j in range(len(blocs[offdiag[i][1]])): blocs[offdiag[i][0]].append(blocs[offdiag[i][1]][j]) for b1,b2 in product(blocs,blocs):
del blocs[offdiag[i][1]] if (pair[0] in b1) and (pair[1] in b2):
for j in range(i+1,len(offdiag)): if blocs.index(b1) != blocs.index(b2): # In separate blocks?
if offdiag[j][0] == offdiag[i][1]: offdiag[j][0] = offdiag[i][0] b1.extend(blocs.pop(blocs.index(b2))) # Merge two blocks
if offdiag[j][1] == offdiag[i][1]: offdiag[j][1] = offdiag[i][0] break # Move on to next pair in offdiag
if offdiag[j][0] > offdiag[i][1]: offdiag[j][0] -= 1
if offdiag[j][1] > offdiag[i][1]: offdiag[j][1] -= 1
offdiag[j].sort()
num_blocs -= 1
# Set the gf_struct for the solver accordingly # Set the gf_struct for the solver accordingly
num_blocs = len(blocs)
for i in range(num_blocs): for i in range(num_blocs):
blocs[i].sort() blocs[i].sort()
self.gf_struct_solver[ish].update( [('%s_%s'%(sp,i),range(len(blocs[i])))] ) self.gf_struct_solver[ish].update( [('%s_%s'%(sp,i),range(len(blocs[i])))] )