245 lines
6.9 KiB
Fortran
245 lines
6.9 KiB
Fortran
|
|
subroutine diag_mat_per_fock_degen_core(fock_diag, mat_ref, listcore,ncore, n, thr_d, thr_nd, thr_deg, leigvec, reigvec, eigval)
|
|
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! subroutine that diagonalizes a matrix mat_ref BY BLOCK
|
|
!
|
|
! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag"
|
|
!
|
|
! the elements of listcore are untouched
|
|
!
|
|
! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together
|
|
!
|
|
! : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together
|
|
!
|
|
! : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together
|
|
!
|
|
! etc... the advantage is to guarentee no spurious mixing because of numerical problems.
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
integer, intent(in) :: n,ncore, listcore(ncore)
|
|
double precision, intent(in) :: fock_diag(n), mat_ref(n,n), thr_d, thr_nd, thr_deg
|
|
double precision, intent(out) :: leigvec(n,n), reigvec(n,n), eigval(n)
|
|
|
|
integer :: n_degen_list, n_degen,size_mat, i, j, k, icount, m, index_degen
|
|
integer :: ii, jj, i_good, j_good, n_real
|
|
integer :: icount_eigval
|
|
logical, allocatable :: is_ok(:)
|
|
integer, allocatable :: list_degen(:,:), list_same_degen(:)
|
|
integer, allocatable :: iorder(:), list_degen_sorted(:)
|
|
double precision, allocatable :: leigvec_unsrtd(:,:), reigvec_unsrtd(:,:), eigval_unsrtd(:)
|
|
double precision, allocatable :: mat_tmp(:,:), eigval_tmp(:), leigvec_tmp(:,:), reigvec_tmp(:,:)
|
|
|
|
allocate(leigvec_unsrtd(n,n), reigvec_unsrtd(n,n), eigval_unsrtd(n))
|
|
leigvec_unsrtd = 0.d0
|
|
reigvec_unsrtd = 0.d0
|
|
eigval_unsrtd = 0.d0
|
|
|
|
! obtain degeneracies
|
|
allocate(list_degen(n,0:n))
|
|
call give_degen_full_listcore(fock_diag, n, listcore, ncore, thr_deg, list_degen, n_degen_list)
|
|
|
|
allocate(iorder(n_degen_list), list_degen_sorted(n_degen_list))
|
|
do i = 1, n_degen_list
|
|
n_degen = list_degen(i,0)
|
|
list_degen_sorted(i) = n_degen
|
|
iorder(i) = i
|
|
enddo
|
|
|
|
! sort by number of degeneracies
|
|
call isort(list_degen_sorted, iorder, n_degen_list)
|
|
|
|
allocate(is_ok(n_degen_list))
|
|
is_ok = .True.
|
|
icount_eigval = 0
|
|
|
|
! loop over degeneracies
|
|
do i = 1, n_degen_list
|
|
if(.not.is_ok(i)) cycle
|
|
|
|
is_ok(i) = .False.
|
|
n_degen = list_degen_sorted(i)
|
|
|
|
|
|
if(n_degen.ge.1000)then
|
|
print*,'core orbital '
|
|
else
|
|
print *, ' diagonalizing for n_degen = ', n_degen
|
|
endif
|
|
|
|
k = 1
|
|
|
|
! group all the entries having the same degeneracies
|
|
!! do while (list_degen_sorted(i+k)==n_degen)
|
|
do m = i+1, n_degen_list
|
|
if(list_degen_sorted(m)==n_degen) then
|
|
is_ok(i+k) = .False.
|
|
k += 1
|
|
endif
|
|
enddo
|
|
|
|
print *, ' number of identical degeneracies = ', k
|
|
if(n_degen.ge.1000)then
|
|
n_degen = 1
|
|
endif
|
|
size_mat = k*n_degen
|
|
print *, ' size_mat = ', size_mat
|
|
allocate(mat_tmp(size_mat,size_mat), list_same_degen(size_mat))
|
|
allocate(eigval_tmp(size_mat), leigvec_tmp(size_mat,size_mat), reigvec_tmp(size_mat,size_mat))
|
|
! group all the elements sharing the same degeneracy
|
|
icount = 0
|
|
do j = 1, k ! jth set of degeneracy
|
|
index_degen = iorder(i+j-1)
|
|
do m = 1, n_degen
|
|
icount += 1
|
|
list_same_degen(icount) = list_degen(index_degen,m)
|
|
enddo
|
|
enddo
|
|
|
|
print *, ' list of elements '
|
|
do icount = 1, size_mat
|
|
print *, icount, list_same_degen(icount)
|
|
enddo
|
|
|
|
! you copy subset of matrix elements having all the same degeneracy in mat_tmp
|
|
do ii = 1, size_mat
|
|
i_good = list_same_degen(ii)
|
|
do jj = 1, size_mat
|
|
j_good = list_same_degen(jj)
|
|
mat_tmp(jj,ii) = mat_ref(j_good,i_good)
|
|
enddo
|
|
enddo
|
|
|
|
call non_hrmt_bieig( size_mat, mat_tmp, thr_d, thr_nd &
|
|
, leigvec_tmp, reigvec_tmp &
|
|
, n_real, eigval_tmp )
|
|
|
|
do ii = 1, size_mat
|
|
icount_eigval += 1
|
|
eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues
|
|
do jj = 1, size_mat ! copy the eigenvectors
|
|
j_good = list_same_degen(jj)
|
|
leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii)
|
|
reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii)
|
|
enddo
|
|
enddo
|
|
|
|
deallocate(mat_tmp, list_same_degen)
|
|
deallocate(eigval_tmp, leigvec_tmp, reigvec_tmp)
|
|
enddo
|
|
|
|
if(icount_eigval .ne. n) then
|
|
print *, ' pb !! (icount_eigval.ne.n)'
|
|
print *, ' icount_eigval,n', icount_eigval, n
|
|
stop
|
|
endif
|
|
|
|
deallocate(iorder)
|
|
allocate(iorder(n))
|
|
do i = 1, n
|
|
iorder(i) = i
|
|
enddo
|
|
call dsort(eigval_unsrtd, iorder, n)
|
|
|
|
do i = 1, n
|
|
print*,'sorted eigenvalues '
|
|
i_good = iorder(i)
|
|
eigval(i) = eigval_unsrtd(i)
|
|
print*,'i,eigval(i) = ',i,eigval(i)
|
|
do j = 1, n
|
|
leigvec(j,i) = leigvec_unsrtd(j,i_good)
|
|
reigvec(j,i) = reigvec_unsrtd(j,i_good)
|
|
enddo
|
|
enddo
|
|
|
|
deallocate(leigvec_unsrtd, reigvec_unsrtd, eigval_unsrtd)
|
|
deallocate(list_degen)
|
|
deallocate(iorder, list_degen_sorted)
|
|
deallocate(is_ok)
|
|
|
|
end
|
|
|
|
! ---
|
|
|
|
subroutine give_degen_full_listcore(A, n, listcore, ncore, thr, list_degen, n_degen_list)
|
|
|
|
BEGIN_DOC
|
|
! you enter with an array A(n) and spits out all the elements degenerated up to thr
|
|
!
|
|
! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL
|
|
!
|
|
! list_degen(i,0) = number of degenerate entries
|
|
!
|
|
! list_degen(i,1) = index of the first degenerate entry
|
|
!
|
|
! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries
|
|
!
|
|
! if list_degen(i,0) == 1 it means that there is no degeneracy for that element
|
|
!
|
|
! if list_degen(i,0) >= 1000 it means that it is core orbitals
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
double precision, intent(in) :: A(n)
|
|
double precision, intent(in) :: thr
|
|
integer, intent(in) :: n,ncore, listcore(ncore)
|
|
integer, intent(out) :: list_degen(n,0:n), n_degen_list
|
|
integer :: i, j, icount, icheck,k
|
|
logical, allocatable :: is_ok(:)
|
|
|
|
|
|
allocate(is_ok(n))
|
|
n_degen_list = 0
|
|
is_ok = .True.
|
|
! you first exclude the "core" orbitals
|
|
do i = 1, ncore
|
|
j=listcore(i)
|
|
is_ok(j) = .False.
|
|
enddo
|
|
do i = 1, n
|
|
if(.not.is_ok(i)) cycle
|
|
n_degen_list +=1
|
|
is_ok(i) = .False.
|
|
list_degen(n_degen_list,1) = i
|
|
icount = 1
|
|
do j = i+1, n
|
|
if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then
|
|
is_ok(j) = .False.
|
|
icount += 1
|
|
list_degen(n_degen_list,icount) = j
|
|
endif
|
|
enddo
|
|
|
|
list_degen(n_degen_list,0) = icount
|
|
enddo
|
|
! you set all the core orbitals as separate entities
|
|
icheck = 0
|
|
do i = 1, n_degen_list
|
|
icheck += list_degen(i,0)
|
|
enddo
|
|
if(icheck.ne.(n-ncore))then
|
|
print *, ' pb ! :: icheck.ne.n-ncore'
|
|
print *, icheck, n-ncore
|
|
stop
|
|
endif
|
|
k=1000
|
|
do i = 1, ncore
|
|
n_degen_list+= 1
|
|
j=listcore(i)
|
|
list_degen(n_degen_list,1) = i
|
|
list_degen(n_degen_list,0) = k
|
|
k+=1
|
|
enddo
|
|
|
|
|
|
|
|
end
|
|
|
|
! ---
|
|
|