From 6235c2015d98c2ed1f89eeca13555cff0e7c8785 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Fri, 22 Dec 2023 20:15:58 +0100 Subject: [PATCH] added non-sym diag for tc-rpa --- .../dav_ext_rout_nonsym_B1space.irp.f | 2 +- src/hartree_fock/print_scf_int.irp.f | 114 +++++++++++ .../lapack_diag_non_hermit.irp.f | 41 ++-- src/tc_bi_ortho/drpa_matrix.irp.f | 116 ----------- src/tc_bi_ortho/tc_effect_int.irp.f | 39 ---- src/tc_bi_ortho/tc_rpa.irp.f | 181 ------------------ src/utils/util.irp.f | 40 ++++ 7 files changed, 182 insertions(+), 351 deletions(-) create mode 100644 src/hartree_fock/print_scf_int.irp.f delete mode 100644 src/tc_bi_ortho/drpa_matrix.irp.f delete mode 100644 src/tc_bi_ortho/tc_effect_int.irp.f delete mode 100644 src/tc_bi_ortho/tc_rpa.irp.f diff --git a/src/dav_general_mat/dav_ext_rout_nonsym_B1space.irp.f b/src/dav_general_mat/dav_ext_rout_nonsym_B1space.irp.f index 4b7b9cc9..d89aaadb 100644 --- a/src/dav_general_mat/dav_ext_rout_nonsym_B1space.irp.f +++ b/src/dav_general_mat/dav_ext_rout_nonsym_B1space.irp.f @@ -346,7 +346,7 @@ subroutine davidson_general_ext_rout_nonsym_b1space(u_in, H_jj, energies, sze, N endif if(i_omax(l) .ne. l) then - print *, ' !!! WARNONG !!!' + print *, ' !!! WARNING !!!' print *, ' index of state', l, i_omax(l) endif enddo diff --git a/src/hartree_fock/print_scf_int.irp.f b/src/hartree_fock/print_scf_int.irp.f new file mode 100644 index 00000000..ee7590f6 --- /dev/null +++ b/src/hartree_fock/print_scf_int.irp.f @@ -0,0 +1,114 @@ + +program print_scf_int + + call main() + +end + +subroutine main() + + implicit none + integer :: i, j, k, l + + print *, " Hcore:" + do j = 1, ao_num + do i = 1, ao_num + print *, i, j, ao_one_e_integrals(i,j) + enddo + enddo + + print *, " P:" + do j = 1, ao_num + do i = 1, ao_num + print *, i, j, SCF_density_matrix_ao_alpha(i,j) + enddo + enddo + + + double precision :: integ, density_a, density_b, density + double precision :: J_scf(ao_num, ao_num) + double precision :: K_scf(ao_num, ao_num) + + + double precision, external :: get_ao_two_e_integral + PROVIDE ao_integrals_map + + print *, " J:" + !do j = 1, ao_num + ! do l = 1, ao_num + ! do i = 1, ao_num + ! do k = 1, ao_num + ! ! < 1:k, 2:l | 1:i, 2:j > + ! print *, '< k l | i j >', k, l, i, j + ! print *, get_ao_two_e_integral(i, j, k, l, ao_integrals_map) + ! enddo + ! enddo + ! enddo + !enddo + + !do k = 1, ao_num + ! do i = 1, ao_num + ! do j = 1, ao_num + ! do l = 1, ao_num + ! ! ( 1:k, 1:i | 2:l, 2:j ) + ! print *, '(k i | l j)', k, i, l, j + ! print *, get_ao_two_e_integral(l, j, k, i, ao_integrals_map) + ! enddo + ! enddo + ! print *, '' + ! enddo + !enddo + + J_scf = 0.d0 + K_scf = 0.d0 + do i = 1, ao_num + do k = 1, ao_num + do j = 1, ao_num + do l = 1, ao_num + + density_a = SCF_density_matrix_ao_alpha(l,j) + density_b = SCF_density_matrix_ao_beta (l,j) + density = density_a + density_b + + integ = get_ao_two_e_integral(l, j, k, i, ao_integrals_map) + J_scf(k,i) += density * integ + integ = get_ao_two_e_integral(l, i, k, j, ao_integrals_map) + K_scf(k,i) -= density_a * integ + enddo + enddo + enddo + enddo + + print *, 'J x P' + do i = 1, ao_num + do k = 1, ao_num + print *, k, i, J_scf(k,i) + enddo + enddo + + print *, '' + print *, 'K x P' + do i = 1, ao_num + do k = 1, ao_num + print *, k, i, K_scf(k,i) + enddo + enddo + + print *, '' + print *, 'F in AO' + do i = 1, ao_num + do k = 1, ao_num + print *, k, i, Fock_matrix_ao(k,i) + enddo + enddo + + print *, '' + print *, 'F in MO' + do i = 1, ao_num + do k = 1, ao_num + print *, k, i, 2.d0 * Fock_matrix_mo_alpha(k,i) + enddo + enddo + +end + 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 09fcee24..1144f29f 100644 --- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f +++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -1883,8 +1883,13 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ enddo accu_nd = dsqrt(accu_nd) / dble(m) - print *, ' accu_nd = ', accu_nd - print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + if((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d) then + print *, ' non bi-orthogonal vectors !' + print *, ' accu_nd = ', accu_nd + print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + else + print *, ' vectors are bi-orthogonaly' + endif ! --- @@ -1994,10 +1999,13 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0) ii = ii + 1 endif enddo + if(ii .eq. 0) then print*, ' WARNING: bi-orthogonality is lost but there is no degeneracies' print*, ' rotations may change energy' + stop endif + print *, ii, ' type of degeneracies' ! --- @@ -2018,17 +2026,18 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0) call dgemm( 'T', 'N', m, m, n, 1.d0 & , L, size(L, 1), R, size(R, 1) & , 0.d0, S, size(S, 1) ) - print*,'Overlap matrix ' - accu_nd = 0.D0 + + print*, 'Overlap matrix ' + accu_nd = 0.d0 do j = 1, m - write(*,'(100(F16.10,X))') S(1:m,j) - do k = 1, m - if(j==k)cycle - accu_nd += dabs(S(j,k)) - enddo + write(*,'(100(F16.10,X))') S(1:m,j) + do k = 1, m + if(j==k) cycle + accu_nd += dabs(S(j,k)) + enddo enddo print*,'accu_nd = ',accu_nd -! if(accu_nd .gt.1.d-10)then +! if(accu_nd .gt.1.d-10) then ! stop ! endif do j = 1, m @@ -2036,13 +2045,15 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0) R0(1:n,i+j-1) = R(1:n,j) enddo - deallocate(L, R,S) + deallocate(L, R, S) endif enddo end subroutine reorder_degen_eigvec +! --- + subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0) implicit none @@ -2108,8 +2119,10 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0) ! --- -! call impose_orthog_svd(n, m, L) call impose_orthog_svd(n, m, R) + L(:,:) = R(:,:) + + !call impose_orthog_svd(n, m, L) !call impose_orthog_GramSchmidt(n, m, L) !call impose_orthog_GramSchmidt(n, m, R) @@ -2128,8 +2141,8 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0) !call bi_ortho_s_inv_half(m, L, R, S_inv_half) !deallocate(S, S_inv_half) - call impose_biorthog_svd(n, m, L, R) -! call impose_biorthog_inverse(n, m, L, R) + !call impose_biorthog_svd(n, m, L, R) + !call impose_biorthog_inverse(n, m, L, R) !call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R) diff --git a/src/tc_bi_ortho/drpa_matrix.irp.f b/src/tc_bi_ortho/drpa_matrix.irp.f deleted file mode 100644 index 56891ca2..00000000 --- a/src/tc_bi_ortho/drpa_matrix.irp.f +++ /dev/null @@ -1,116 +0,0 @@ - -BEGIN_PROVIDER [double precision, M_RPA, (2*nS_exc, 2*nS_exc)] - - BEGIN_DOC - ! - ! full matrix for direct RPA calculation - ! with the TC-Hamiltonian - ! - END_DOC - - implicit none - integer :: ia, i, a, jb, j, b - double precision :: e(mo_num) - double precision, external :: Kronecker_delta - - PROVIDE mo_tc_effec2e_int - PROVIDE Fock_matrix_tc_diag_mo_tot - - e(1:mo_num) = Fock_matrix_tc_diag_mo_tot(1:mo_num) - - - ! --- --- --- - ! block A - - ia = 0 - do i = nC_orb+1, nO_orb - do a = nO_orb+1, mo_num-nR_orb - ia = ia + 1 - - jb = 0 - do j = nC_orb+1, nO_orb - do b = nO_orb+1, mo_num-nR_orb - jb = jb + 1 - - M_RPA(ia,jb) = (e(a) - e(i)) * Kronecker_delta(i,j) * Kronecker_delta(a,b) + 2.d0 * mo_tc_effec2e_int(a,j,i,b) - enddo - enddo - enddo - enddo - - ! - ! --- --- --- - - - ! --- --- --- - ! block B - - ia = 0 - do i = nC_orb+1, nO_orb - do a = nO_orb+1, mo_num-nR_orb - ia = ia + 1 - - jb = nS_exc - do j = nC_orb+1, nO_orb - do b = nO_orb+1, mo_num-nR_orb - jb = jb + 1 - - M_RPA(ia,jb) = 2.d0 * mo_tc_effec2e_int(a,b,i,j) - enddo - enddo - enddo - enddo - - ! - ! --- --- --- - - - ! --- --- --- - ! block C - - ia = nS_exc - do i = nC_orb+1, nO_orb - do a = nO_orb+1, mo_num-nR_orb - ia = ia + 1 - - jb = 0 - do j = nC_orb+1, nO_orb - do b = nO_orb+1, mo_num-nR_orb - jb = jb + 1 - - M_RPA(ia,jb) = 2.d0 * mo_tc_effec2e_int(i,j,a,b) - enddo - enddo - enddo - enddo - - ! - ! --- --- --- - - - ! --- --- --- - ! block D - - ia = nS_exc - do i = nC_orb+1, nO_orb - do a = nO_orb+1, mo_num-nR_orb - ia = ia + 1 - - jb = nS_exc - do j = nC_orb+1, nO_orb - do b = nO_orb+1, mo_num-nR_orb - jb = jb + 1 - - M_RPA(ia,jb) = (e(a) - e(i)) * Kronecker_delta(i,j) * Kronecker_delta(a,b) + 2.d0 * mo_tc_effec2e_int(i,b,a,j) - enddo - enddo - enddo - enddo - - ! - ! --- --- --- - - -END_PROVIDER - - diff --git a/src/tc_bi_ortho/tc_effect_int.irp.f b/src/tc_bi_ortho/tc_effect_int.irp.f deleted file mode 100644 index 48a786d2..00000000 --- a/src/tc_bi_ortho/tc_effect_int.irp.f +++ /dev/null @@ -1,39 +0,0 @@ - - -BEGIN_PROVIDER [double precision, mo_tc_effec2e_int, (mo_num, mo_num, mo_num, mo_num)] - - BEGIN_DOC - ! - ! mo_tc_effec2e_int(p,q,s,t) = < p q| V(12) | s t > + \sum_i < p q i | L(123)| s t i > - ! - ! the potential V(12) contains ALL TWO-E CONTRIBUTION OF THE TC-HAMILTONIAN - ! - END_DOC - - implicit none - integer :: i, j, k, l, ii - double precision :: integral - - PROVIDE mo_bi_ortho_tc_two_e_chemist - - do j = 1, mo_num - do i = 1, mo_num - do l = 1, mo_num - do k = 1, mo_num - mo_tc_effec2e_int(k,l,i,j) = mo_bi_ortho_tc_two_e_chemist(k,i,l,j) - - do ii = 1, elec_alpha_num - call give_integrals_3_body_bi_ort(k, l, ii, i, j, ii, integral) - mo_tc_effec2e_int(k,l,i,j) -= 2.d0 * integral - enddo - enddo - enddo - enddo - enddo - - FREE mo_bi_ortho_tc_two_e_chemist - -END_PROVIDER - -! --- - diff --git a/src/tc_bi_ortho/tc_rpa.irp.f b/src/tc_bi_ortho/tc_rpa.irp.f deleted file mode 100644 index c9818a1d..00000000 --- a/src/tc_bi_ortho/tc_rpa.irp.f +++ /dev/null @@ -1,181 +0,0 @@ -program tc_rpa - - BEGIN_DOC - ! - ! - ! - END_DOC - - my_grid_becke = .True. - PROVIDE tc_grid1_a tc_grid1_r - my_n_pt_r_grid = tc_grid1_r - my_n_pt_a_grid = tc_grid1_a - touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - - if(j1b_type .ge. 100) then - my_extra_grid_becke = .True. - PROVIDE tc_grid2_a tc_grid2_r - my_n_pt_r_extra_grid = tc_grid2_r - my_n_pt_a_extra_grid = tc_grid2_a - touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid - - call write_int(6, my_n_pt_r_extra_grid, 'radial internal grid over') - call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over') - endif - - - call main() - -end - -! --- - -subroutine main() - - implicit none - integer :: i, j, n - integer :: n_good, n_real_eigv - double precision :: thr_cpx, thr_d, thr_nd - double precision :: accu_d, accu_nd - integer, allocatable :: list_good(:), iorder(:) - double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:) - double precision, allocatable :: Omega_p(:), Reigvec_p(:,:), Leigvec_p(:,:) - double precision, allocatable :: Omega_m(:), Reigvec_m(:,:), Leigvec_m(:,:) - double precision, allocatable :: S(:,:) - - PROVIDE M_RPA - - print *, ' ' - print *, ' Computing left/right eigenvectors for TC-RPA ...' - print *, ' ' - - - n = 2 * nS_exc - - thr_cpx = 1d-7 - thr_d = 1d-07 - thr_nd = 1d-07 - - - allocate(WR(n), WI(n), VL(n,n), VR(n,n)) - call lapack_diag_non_sym(n, M_RPA, WR, WI, VL, VR) - FREE M_RPA - - print *, ' excitation energies:' - do i = 1, nS_exc - write(*, '(I3, X, 1000(F16.10,X))') i, WR(i), WI(i) - if(dabs(WI(i)) .gt. thr_cpx) then - print *, ' WARNING ! IMAGINARY EIGENVALUES !!!' - write(*, '(1000(F16.10,X))') WR(i), WI(i+1) - endif - enddo - - print *, ' ' - print *, ' desexcitation energies:' - do i = nS_exc+1, n - write(*, '(I3, X, 1000(F16.10,X))') i, WR(i), WI(i) - if(dabs(WI(i)) .gt. thr_cpx) then - print *, ' WARNING ! IMAGINARY EIGENVALUES !!!' - write(*, '(1000(F16.10,X))') WR(i), WI(i+1) - endif - enddo - - - ! track & sort the real eigenvalues - - n_good = 0 - do i = 1, nS_exc - if(dabs(WI(i)) .lt. thr_cpx) then - if(dabs(WI(nS_exc+i)) .lt. thr_cpx) then - n_good += 1 - endif - endif - enddo - n_real_eigv = n_good - - print *, ' ' - print *, ' nb of real eigenvalues = ', n_real_eigv - print *, ' total nb of eigenvalues = ', nS_exc - - allocate(Omega_p(n_real_eigv), Reigvec_p(n,n_real_eigv), Leigvec_p(n,n_real_eigv)) - allocate(Omega_m(n_real_eigv), Reigvec_m(n,n_real_eigv), Leigvec_m(n,n_real_eigv)) - - n_good = 0 - do i = 1, nS_exc - if(dabs(WI(i)) .lt. thr_cpx) then - if(dabs(WI(nS_exc+i)) .lt. thr_cpx) then - n_good += 1 - - Omega_p(n_good) = WR(i) - do j = 1, n - Reigvec_p(j,n_good) = VR(j,n_good) - Leigvec_p(j,n_good) = VL(j,n_good) - enddo - - Omega_m(n_good) = WR(nS_exc+i) - do j = 1, n - Reigvec_m(j,n_good) = VR(j,nS_exc+n_good) - Leigvec_m(j,n_good) = VL(j,nS_exc+n_good) - enddo - endif - endif - enddo - - deallocate(WR, WI, VL, VR) - - - ! check bi-orthogonality - - ! first block - - allocate(S(n_real_eigv,n_real_eigv)) - - call check_biorthog(n, n_real_eigv, Leigvec_p, Reigvec_p, accu_d, accu_nd, S, thr_d, thr_nd, .false.) - print *, ' accu_d = ', accu_d - print *, ' accu_nd = ', accu_nd - - if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d)) then - print *, ' RPA first-block eigenvectors are normalized and bi-orthogonalized' - else - print *, ' RPA first-block eigenvectors are neither normalized nor bi-orthogonalized' - - call reorder_degen_eigvec(n, Omega_p, Leigvec_p, Reigvec_p) - call impose_biorthog_degen_eigvec(n, Omega_p, Leigvec_p, Reigvec_p) - - call check_biorthog(n, n_real_eigv, Leigvec_p, Reigvec_p, accu_d, accu_nd, S, thr_d, thr_nd, .false.) - if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .gt. thr_d) ) then - call check_biorthog_binormalize(n, n_real_eigv, Leigvec_p, Reigvec_p, thr_d, thr_nd, .true.) - endif - call check_biorthog(n, n_real_eigv, Leigvec_p, Reigvec_p, accu_d, accu_nd, S, thr_d, thr_nd, .true.) - endif - - - ! second block - - call check_biorthog(n, n_real_eigv, Leigvec_m, Reigvec_m, accu_d, accu_nd, S, thr_d, thr_nd, .false.) - print *, ' accu_d = ', accu_d - print *, ' accu_nd = ', accu_nd - - if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d)) then - print *, ' RPA first-block eigenvectors are normalized and bi-orthogonalized' - else - print *, ' RPA first-block eigenvectors are neither normalized nor bi-orthogonalized' - - call reorder_degen_eigvec(n, Omega_m, Leigvec_m, Reigvec_m) - call impose_biorthog_degen_eigvec(n, Omega_m, Leigvec_m, Reigvec_m) - - call check_biorthog(n, n_real_eigv, Leigvec_m, Reigvec_m, accu_d, accu_nd, S, thr_d, thr_nd, .false.) - if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv)) .gt. thr_d) ) then - call check_biorthog_binormalize(n, n_real_eigv, Leigvec_m, Reigvec_m, thr_d, thr_nd, .true.) - endif - call check_biorthog(n, n_real_eigv, Leigvec_m, Reigvec_m, accu_d, accu_nd, S, thr_d, thr_nd, .true.) - endif - - deallocate(S) - - return - -end - -! --- - diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f index 785d6539..97cbde67 100644 --- a/src/utils/util.irp.f +++ b/src/utils/util.irp.f @@ -600,3 +600,43 @@ end function Kronecker_delta ! --- +subroutine diagonalize_sym_matrix(N, A, e) + + BEGIN_DOC + ! + ! Diagonalize a symmetric matrix + ! + END_DOC + + implicit none + + integer, intent(in) :: N + double precision, intent(inout) :: A(N,N) + double precision, intent(out) :: e(N) + + integer :: lwork, info + double precision, allocatable :: work(:) + + allocate(work(1)) + + lwork = -1 + call dsyev('V', 'U', N, A, N, e, work, lwork, info) + lwork = int(work(1)) + + deallocate(work) + + allocate(work(lwork)) + + call dsyev('V', 'U', N, A, N, e, work, lwork, info) + deallocate(work) + + if(info /= 0) then + print*,'Problem in diagonalize_sym_matrix (dsyev)!!' + endif + +end subroutine diagonalize_sym_matrix + +! --- + + +