From 1e2f8ab44fc36b6696be4b6a1a67984cf88d194b Mon Sep 17 00:00:00 2001 From: eginer Date: Sat, 29 Oct 2022 20:07:33 +0200 Subject: [PATCH] adapted the tc_natorb with block diagonalization --- .../lapack_diag_non_hermit.irp.f | 19 +++++++++++++++---- src/tc_bi_ortho/tc_natorb.irp.f | 16 +++++++++++----- src/tc_scf/minimize_tc_angles.irp.f | 2 +- src/tools/NEED | 1 + src/utils/block_diag_degen.irp.f | 10 +++++----- 5 files changed, 33 insertions(+), 15 deletions(-) diff --git a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f index b92f12e5..d79e25ba 100644 --- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f +++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -930,11 +930,15 @@ subroutine check_EIGVEC(n, m, A, eigval, leigvec, reigvec, thr_diag, thr_norm, s tmp_abs = tmp_abs + tmp V_nrm = V_nrm + U_nrm - print *, j, tmp, U_nrm + write(*,'(I4,X,(100(F25.16,X)))')j,eigval(j), tmp, U_nrm enddo - tmp_rel = tmp_abs / tmp_nrm + if(tmp_abs.lt.10.d-10)then + tmp_rel = thr_diag/10.d0 + else + tmp_rel = tmp_abs / tmp_nrm + endif tmp_dif = dabs(V_nrm - dble(m)) if( stop_ifnot .and. ((tmp_rel .gt. thr_diag) .or. (tmp_dif .gt. thr_norm)) ) then @@ -969,13 +973,20 @@ subroutine check_EIGVEC(n, m, A, eigval, leigvec, reigvec, thr_diag, thr_norm, s tmp_abs = tmp_abs + tmp V_nrm = V_nrm + U_nrm - print *, j, tmp, U_nrm + write(*,'(I4,X,(100(F25.16,X)))')j,eigval(j), tmp, U_nrm enddo - tmp_rel = tmp_abs / tmp_nrm + if(tmp_abs.lt.10.d-10)then + tmp_rel = thr_diag/10.d0 + else + tmp_rel = tmp_abs / tmp_nrm + endif if( stop_ifnot .and. ((tmp_rel .gt. thr_diag) .or. (tmp_dif .gt. thr_norm)) ) then print *, ' error in left-eigenvectors' + print *, ' err tol = ',thr_diag, thr_norm + print *, '(tmp_rel .gt. thr_diag) = ',(tmp_rel .gt. thr_diag) + print *, '(tmp_dif .gt. thr_norm) = ',(tmp_dif .gt. thr_norm) print *, ' err estim = ', tmp_abs, tmp_rel print *, ' CR norm = ', V_nrm stop diff --git a/src/tc_bi_ortho/tc_natorb.irp.f b/src/tc_bi_ortho/tc_natorb.irp.f index 71fc6716..0e012b18 100644 --- a/src/tc_bi_ortho/tc_natorb.irp.f +++ b/src/tc_bi_ortho/tc_natorb.irp.f @@ -8,17 +8,23 @@ ! natorb_tc_leigvec_mo : LEFT eigenvectors of the ground state transition matrix (equivalent of natural orbitals) ! natorb_tc_eigval : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !! END_DOC - double precision, allocatable :: dm_tmp(:,:) + double precision, allocatable :: dm_tmp(:,:),fock_diag(:) + double precision :: thr_deg integer :: i,j,k,n_real - allocate( dm_tmp(mo_num,mo_num)) + allocate( dm_tmp(mo_num,mo_num),fock_diag(mo_num)) dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1) print*,'dm_tmp' do i = 1, mo_num + fock_diag(i) = fock_matrix_tc_mo_tot(i,i) write(*,'(100(F16.10,X))')-dm_tmp(:,i) enddo - call non_hrmt_bieig( mo_num, dm_tmp& - , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& - , n_real, natorb_tc_eigval ) + thr_deg = thr_degen_tc + call diag_mat_per_fock_degen(fock_diag,dm_tmp,mo_num,thr_deg,& + natorb_tc_leigvec_mo,natorb_tc_reigvec_mo,& + natorb_tc_eigval) +! call non_hrmt_bieig( mo_num, dm_tmp& +! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& +! , n_real, natorb_tc_eigval ) double precision :: accu accu = 0.d0 do i = 1, n_real diff --git a/src/tc_scf/minimize_tc_angles.irp.f b/src/tc_scf/minimize_tc_angles.irp.f index b52beec0..cb729eb2 100644 --- a/src/tc_scf/minimize_tc_angles.irp.f +++ b/src/tc_scf/minimize_tc_angles.irp.f @@ -6,7 +6,7 @@ program print_angles ! my_n_pt_r_grid = 10 ! small grid for quick debug ! my_n_pt_a_grid = 14 ! small grid for quick debug touch my_n_pt_r_grid my_n_pt_a_grid - call sort_by_tc_fock +! call sort_by_tc_fock call minimize_tc_orb_angles end diff --git a/src/tools/NEED b/src/tools/NEED index c07c9109..46507542 100644 --- a/src/tools/NEED +++ b/src/tools/NEED @@ -2,3 +2,4 @@ fci mo_two_e_erf_ints aux_quantities hartree_fock +dft_utils_in_r diff --git a/src/utils/block_diag_degen.irp.f b/src/utils/block_diag_degen.irp.f index 9884bf21..f6906e23 100644 --- a/src/utils/block_diag_degen.irp.f +++ b/src/utils/block_diag_degen.irp.f @@ -29,10 +29,10 @@ subroutine diag_mat_per_fock_degen(fock_diag,mat_ref,n,thr_deg,leigvec,reigvec,e reigvec_unsrtd = 0.d0 eigval_unsrtd = 0.d0 - allocate(list_degen(mo_num,0:mo_num)) + allocate(list_degen(n,0:n)) ! obtain degeneracies - call give_degen_full_list(fock_diag,mo_num,thr_deg,list_degen,n_degen_list) + call give_degen_full_list(fock_diag,n,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) @@ -102,9 +102,9 @@ subroutine diag_mat_per_fock_degen(fock_diag,mat_ref,n,thr_deg,leigvec,reigvec,e deallocate(mat_tmp,list_same_degen) deallocate(eigval_tmp,leigvec_tmp,reigvec_tmp) enddo - if(icount.ne.n)then - print*,'pb !! (icount.ne.n)' - print*,'icount,n',icount,n + if(icount_eigval.ne.n)then + print*,'pb !! (icount_eigval.ne.n)' + print*,'icount_eigval,n',icount_eigval,n stop endif