added routines to rotate orbitals without touching core orbitals

src/tc_bi_ortho/tc_natorb.irp.f

@@ -23,7 +23,7 @@
 
   dm_tmp(1:mo_num,1:mo_num) = -tc_transition_matrix_mo(1:mo_num,1:mo_num,1,1)
 
-  print *, ' dm_tmp'
+  print *, ' Transition density matrix '
   do i = 1, mo_num
     fock_diag(i) = fock_matrix_tc_mo_tot(i,i)
     write(*, '(100(F16.10,X))') -dm_tmp(:,i)
@@ -32,8 +32,15 @@
   thr_d   = 1.d-6
   thr_nd  = 1.d-6
   thr_deg = 1.d-3
-  call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & 
+  if(n_core_orb.ne.0)then
-                              , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
+!   print*,'core orbitals'
+!   pause
+   call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg & 
+                                    , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
+  else
+   call diag_mat_per_fock_degen( fock_diag, dm_tmp, mo_num, thr_d, thr_nd, thr_deg & 
+                               , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)
+  endif
 !   call non_hrmt_bieig( mo_num, dm_tmp&
 !                      , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& 
 !                      , mo_num, natorb_tc_eigval )

src/tc_keywords/EZFIO.cfg

@@ -208,6 +208,12 @@ doc: Threshold to determine if diagonal elements of the bi-orthogonal condition
 interface: ezfio,provider,ocaml
 default: 1.e-6
 
+[thresh_lr_angle]
+type: double precision
+doc: Maximum value of the angle between the couple of left and right orbital for the rotations
+interface: ezfio,provider,ocaml
+default: 20.0
+
 [thresh_biorthog_nondiag]
 type: Threshold
 doc: Threshold to determine if non-diagonal elements of L.T x R are close enouph to 0

src/tc_scf/routines_rotates.irp.f

@@ -140,7 +140,11 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
  ! compute the overlap between the left and rescaled right
   call build_s_matrix(ao_num, mo_num, mo_r_coef_new, mo_r_coef_new, ao_overlap, s_mat)
 ! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
-  call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list)
+  if(n_core_orb.ne.0)then
+   call give_degen_full_listcore(fock_diag, mo_num, list_core, n_core_orb, thr_deg, list_degen, n_degen_list)
+  else
+   call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list)
+  endif
   print *, ' fock_matrix_mo'
   do i = 1, mo_num
     print *, i, fock_diag(i), angle_left_right(i)
@@ -152,6 +156,8 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
 !  n_degen = ilast - ifirst +1
 
     n_degen = list_degen(i,0)
+    if(n_degen .ge. 1000)n_degen = 1 ! convention for core orbitals
+     
     if(n_degen .eq. 1) cycle
 
     allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen))
@@ -279,7 +285,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
   allocate(new_angles(mo_num))
   new_angles(1:mo_num) = dabs(angle_left_right(1:mo_num))
   max_angle = maxval(new_angles)
-  good_angles = max_angle.lt.45.d0
+  good_angles = max_angle.lt.thresh_lr_angle
   print *, ' max_angle = ', max_angle
   deallocate(new_angles)
 
@@ -397,11 +403,11 @@ subroutine print_energy_and_mos(good_angles)
   print *, ' TC SCF energy gradient = ', grad_non_hermit
   print *, ' Max angle Left/right   = ', max_angle_left_right
 
-  if(max_angle_left_right .lt. 45.d0) then
+  if(max_angle_left_right .lt. thresh_lr_angle) then
     print *, ' Maximum angle BELOW 45 degrees, everthing is OK !'
     good_angles = .true.
-  else if(max_angle_left_right .gt. 45.d0 .and. max_angle_left_right .lt. 75.d0) then
+  else if(max_angle_left_right .gt. thresh_lr_angle .and. max_angle_left_right .lt. 75.d0) then
-    print *, ' Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...'
+    print *, ' Maximum angle between thresh_lr_angle and 75 degrees, this is not the best for TC-CI calculations ...'
     good_angles = .false.
   else if(max_angle_left_right .gt. 75.d0) then
     print *, ' Maximum angle between ABOVE 75 degrees, YOU WILL CERTAINLY FIND TROUBLES IN TC-CI calculations ...'

src/utils/block_diag_degen_core.irp.f

@@ -0,0 +1,244 @@
+
+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
+
+! ---
+