3
0
mirror of https://github.com/triqs/dft_tools synced 2024-06-22 05:02:20 +02:00

Fixed a bug in the orthogonalization routine

There was a very nasty bug in the preparation of the block matrix
'p_mat'. The point is that this matrix is created once for all k-points
with the band dimension being the maximum possible. However, only
a part of the matrix is used at every k-point but the orthogonalization
is done for the whole matrix. The problem was that if the number of
bands for a given k-point was smaller than that for the next k-point
them for the next k-point some part of 'p_mat' still contained data from
the previous step, which messed up the orthonormalization. Now, 'p_mat'
is set to zero at each step of the loop.
Also, property 'nion' was added to ProjectorShell since it is used
very often.
This commit is contained in:
Oleg E. Peil 2015-11-18 15:17:51 +01:00
parent f994b82704
commit c1b3000c00
2 changed files with 5 additions and 7 deletions

View File

@ -154,24 +154,26 @@ class ProjectorGroup:
ndim = i2_bl
p_mat = np.zeros((ndim, nb_max), dtype=np.complex128)
# Note that 'ns' and 'nk' are the same for all shells
for isp in xrange(ns):
for ik in xrange(nk):
nb = self.ib_win[ik, isp, 1] - self.ib_win[ik, isp, 0] + 1
# Combine all projectors of the group to one block projector
p_mat[:, :] = 0.0j # !!! Clean-up from the last k-point!
for ish in self.ishells:
shell = self.shells[ish]
blocks = bl_map[ish]['bmat_blocks']
nion = shell.nion # !!!
for ion in xrange(nion):
i1, i2 = blocks[ion]
p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
# Now orthogonalize the obtained block projector
p_orth, overl, eig = self.orthogonalize_projector_matrix(p_mat)
# print "ik = ", ik
# print overl.real
# Distribute back projectors in the same order
for ish in self.ishells:
shell = self.shells[ish]
blocks = bl_map[ish]['bmat_blocks']
nion = shell.nion # !!!
for ion in xrange(nion):
i1, i2 = blocks[ion]
shell.proj_win[ion, isp, ik, :nlm, :nb] = p_orth[i1:i2, :nb]

View File

@ -47,11 +47,7 @@ class ProjectorShell:
self.lm2 = (self.lorb+1)**2
self.ndim = self.extract_tmatrices(sh_pars)
# if self.tmatrix is None:
# self.ndim = self.lm2 - self.lm1
# else:
## TODO: generalize this to a tmatrix for every ion
# self.ndim = self.tmatrix.shape[0]
self.nion = len(self.ion_list)
# Pre-select a subset of projectors (this should be an array view => no memory is wasted)
# !!! This sucks but I have to change the order of 'ib' and 'ilm' indices here