From 861f0f946119267250eb3763d56b0c36b714c589 Mon Sep 17 00:00:00 2001
From: eginer <giner.emmanuel@gmail.com>
Date: Sat, 29 Oct 2022 18:39:30 +0200
Subject: [PATCH] added a diagonalization by block, doing some cleaning

---
 src/tc_bi_ortho/tc_natorb.irp.f   | 177 ++++++++++++++++++++++++++++++
 src/tc_bi_ortho/tc_prop.irp.f     |   2 -
 src/tc_bi_ortho/test_natorb.irp.f |  51 +++++++++
 3 files changed, 228 insertions(+), 2 deletions(-)
 create mode 100644 src/tc_bi_ortho/tc_natorb.irp.f
 create mode 100644 src/tc_bi_ortho/test_natorb.irp.f

diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f
new file mode 100644
index 00000000..e1bac6e4
--- /dev/null
+++ b/src/tc_bi_ortho/tc_natorb.irp.f
@@ -0,0 +1,177 @@
+
+subroutine diagonalize_dm_per_fock_degen(fock_diag,mat_ref,n,thr_deg,leigvec,reigvec,eigval)
+ implicit none
+ integer, intent(in) :: n
+ double precision, intent(in) :: fock_diag(n),mat_ref(n,n),thr_deg
+ double precision, intent(out):: leigvec(n,n),reigvec(n,n),eigval(n)
+ 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"
+!
+! 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
+ double precision, allocatable :: leigvec_unsrtd(:,:),reigvec_unsrtd(:,:),eigval_unsrtd(:)
+ integer, allocatable :: list_degen(:,:),list_same_degen(:)
+ integer, allocatable :: iorder(:),list_degen_sorted(:)
+ integer :: n_degen_list,n_degen,size_mat,i,j,k,icount,m,index_degen
+ integer :: ii,jj,i_good,j_good,n_real
+ 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
+
+ allocate(list_degen(mo_num,0:mo_num))
+
+ ! obtain degeneracies 
+ call give_degen_full_list(fock_diag,mo_num,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)
+ integer :: icount_eigval
+ logical, allocatable :: is_ok(:)
+ 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)
+  print*,'diagonalizing for n_degen = ',n_degen
+  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
+  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*,'icount = ',icount
+  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&
+                     , 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
+ print*,'icount,n',icount,n
+ 
+ deallocate(iorder)
+ allocate(iorder(n))
+ print*,'before sort '
+ do i = 1, n
+  print*,'i,eigval(i) = ',i,eigval_unsrtd(i)
+  iorder(i) = i
+ enddo
+ call dsort(eigval_unsrtd,iorder,n)
+ print*,'after sort '
+ do i = 1, n
+  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
+end
+
+subroutine give_degen_full_list(A,n,thr,list_degen,n_degen_list)
+ implicit none
+ 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
+ END_DOC
+ double precision,intent(in) :: A(n)
+ double precision,intent(in)    :: thr
+ integer, intent(in) :: n
+ integer, intent(out) :: list_degen(n,0:n),n_degen_list
+ logical, allocatable :: is_ok(:)
+ allocate(is_ok(n))
+ integer :: i,j,icount,icheck
+ n_degen_list = 0
+ is_ok = .True.
+ 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
+!    print*,A(i),A(j)
+    is_ok(j) = .False.
+    icount += 1
+    list_degen(n_degen_list,icount) = j
+   endif
+  enddo
+  list_degen(n_degen_list,0) = icount
+ enddo
+ icheck = 0
+ do i = 1, n_degen_list
+  icheck += list_degen(i,0)
+ enddo
+ if(icheck.ne.n)then
+  print*,'pb ! :: icheck.ne.n'
+  print*,icheck,n
+  stop
+ endif
+end
diff --git a/src/tc_bi_ortho/tc_prop.irp.f b/src/tc_bi_ortho/tc_prop.irp.f
index daea11a0..bc755ada 100644
--- a/src/tc_bi_ortho/tc_prop.irp.f
+++ b/src/tc_bi_ortho/tc_prop.irp.f
@@ -66,8 +66,6 @@ BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_state
   do i = 1, mo_num
    write(*,'(100(F16.10,X))')-dm_tmp(:,i)
   enddo
-!    call non_hrmt_diag_split_degen( mo_num, dm_tmp&
-!    call non_hrmt_fock_mat( mo_num, dm_tmp&
    call non_hrmt_bieig( mo_num, dm_tmp&
                       , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& 
                       , n_real, natorb_tc_eigval )
diff --git a/src/tc_bi_ortho/test_natorb.irp.f b/src/tc_bi_ortho/test_natorb.irp.f
new file mode 100644
index 00000000..f99f3cea
--- /dev/null
+++ b/src/tc_bi_ortho/test_natorb.irp.f
@@ -0,0 +1,51 @@
+program test_natorb
+  implicit none
+  BEGIN_DOC
+! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. 
+  END_DOC
+  print *, 'Hello world'
+  my_grid_becke = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  read_wf = .True.
+  touch read_wf
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+ call routine
+! call test
+
+end
+
+subroutine routine
+ implicit none
+ double precision, allocatable :: fock_diag(:),eigval(:),leigvec(:,:),reigvec(:,:),mat_ref(:,:)
+ allocate(eigval(mo_num),leigvec(mo_num,mo_num),reigvec(mo_num,mo_num),fock_diag(mo_num),mat_ref(mo_num, mo_num))
+ double precision, allocatable :: eigval_ref(:),leigvec_ref(:,:),reigvec_ref(:,:)
+ allocate(eigval_ref(mo_num),leigvec_ref(mo_num,mo_num),reigvec_ref(mo_num,mo_num))
+ 
+ double precision :: thr_deg
+ integer :: i,n_real,j
+ print*,'fock_matrix'
+ do i = 1, mo_num
+  fock_diag(i) = Fock_matrix_mo(i,i)
+  print*,i,fock_diag(i)
+ enddo
+ thr_deg = 1.d-6
+ mat_ref = -one_e_dm_mo
+ print*,'diagonalization by block'
+ call diagonalize_dm_per_fock_degen(fock_diag,mat_ref,mo_num,thr_deg,leigvec,reigvec,eigval)
+ call non_hrmt_bieig( mo_num, mat_ref&
+                     , leigvec_ref, reigvec_ref& 
+                     , n_real, eigval_ref)
+ print*,'TEST ***********************************'
+ double precision :: accu_l, accu_r
+ do i = 1, mo_num
+  accu_l = 0.d0
+  accu_r = 0.d0
+  do j = 1, mo_num
+   accu_r += reigvec_ref(j,i) * reigvec(j,i)
+   accu_l += leigvec_ref(j,i) * leigvec(j,i)
+  enddo
+  print*,i
+  write(*,'(I3,X,100(F16.10,X))')i,eigval(i),eigval_ref(i),accu_l,accu_r
+ enddo
+end