From c3d257c7aceff4c8641fc751cf9996783c1ce96f Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 22 Jun 2023 11:58:32 +0200 Subject: [PATCH] added routines to rotate orbitals without touching core orbitals --- src/tc_bi_ortho/tc_natorb.irp.f | 13 +- src/tc_keywords/EZFIO.cfg | 6 + src/tc_scf/routines_rotates.irp.f | 16 +- src/utils/block_diag_degen_core.irp.f | 244 ++++++++++++++++++++++++++ 4 files changed, 271 insertions(+), 8 deletions(-) create mode 100644 src/utils/block_diag_degen_core.irp.f diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f index b7e5ae81..1b5a66f3 100644 --- a/src/tc_bi_ortho/tc_natorb.irp.f +++ b/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 & - , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval) + if(n_core_orb.ne.0)then +! 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 ) diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index a69f5bac..f984d53a 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/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 diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f index 755c35b9..588382b5 100644 --- a/src/tc_scf/routines_rotates.irp.f +++ b/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 - print *, ' Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...' + else if(max_angle_left_right .gt. thresh_lr_angle .and. max_angle_left_right .lt. 75.d0) then + 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 ...' diff --git a/src/utils/block_diag_degen_core.irp.f b/src/utils/block_diag_degen_core.irp.f new file mode 100644 index 00000000..5d46bd87 --- /dev/null +++ b/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 + +! --- +