From 9fc4b6d63bbfa3f91d29a7a8f2c5452cb357bed9 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 28 Oct 2023 21:53:04 +0200 Subject: [PATCH] v0 of tc-dRPA --- .../lapack_diag_non_hermit.irp.f | 12 +- src/tc_bi_ortho/ORBITALS.irp.f | 38 ++++ 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 | 19 ++ 6 files changed, 403 insertions(+), 2 deletions(-) create mode 100644 src/tc_bi_ortho/ORBITALS.irp.f create mode 100644 src/tc_bi_ortho/drpa_matrix.irp.f create mode 100644 src/tc_bi_ortho/tc_effect_int.irp.f create mode 100644 src/tc_bi_ortho/tc_rpa.irp.f 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 836bf707..09fcee24 100644 --- a/src/non_hermit_dav/lapack_diag_non_hermit.irp.f +++ b/src/non_hermit_dav/lapack_diag_non_hermit.irp.f @@ -1944,6 +1944,7 @@ subroutine check_orthog(n, m, V, accu_d, accu_nd, S) end subroutine check_orthog ! --- + subroutine reorder_degen_eigvec(n, e0, L0, R0) implicit none @@ -1953,7 +1954,7 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0) double precision, intent(inout) :: L0(n,n), R0(n,n) logical :: complex_root - integer :: i, j, k, m + integer :: i, j, k, m, ii double precision :: ei, ej, de, de_thr double precision :: accu_d, accu_nd integer, allocatable :: deg_num(:) @@ -1986,11 +1987,18 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0) enddo enddo + ii = 0 do i = 1, n if(deg_num(i) .gt. 1) then print *, ' degen on', i, deg_num(i), e0(i) + 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' + endif + print *, ii, ' type of degeneracies' ! --- @@ -2013,7 +2021,7 @@ subroutine reorder_degen_eigvec(n, e0, L0, R0) print*,'Overlap matrix ' accu_nd = 0.D0 do j = 1, m - write(*,'(100(F16.10,X))')S(1:m,j) + write(*,'(100(F16.10,X))') S(1:m,j) do k = 1, m if(j==k)cycle accu_nd += dabs(S(j,k)) diff --git a/src/tc_bi_ortho/ORBITALS.irp.f b/src/tc_bi_ortho/ORBITALS.irp.f new file mode 100644 index 00000000..fdc4758d --- /dev/null +++ b/src/tc_bi_ortho/ORBITALS.irp.f @@ -0,0 +1,38 @@ + +! --- + + BEGIN_PROVIDER [integer, nC_orb] +&BEGIN_PROVIDER [integer, nO_orb] +&BEGIN_PROVIDER [integer, nV_orb] +&BEGIN_PROVIDER [integer, nR_orb] +&BEGIN_PROVIDER [integer, nS_exc] + + BEGIN_DOC + ! + ! nC_orb = number of core orbitals + ! nO_orb = number of occupied orbitals + ! nV_orb = number of virtual orbitals + ! nR_orb = number of Rydberg orbitals + ! nS_exc = number of single excitation + ! + END_DOC + + implicit none + + nC_orb = 0 + nO_orb = elec_beta_num - nC_orb + nV_orb = mo_num - (nC_orb + nO_orb) + nR_orb = 0 + nS_exc = (nO_orb-nC_orb) * (nV_orb-nR_orb) + + print *, ' nC_orb = ', nC_orb + print *, ' nO_orb = ', nO_orb + print *, ' nV_orb = ', nV_orb + print *, ' nR_orb = ', nR_orb + print *, ' nS_exc = ', nS_exc + +END_PROVIDER + +! --- + + diff --git a/src/tc_bi_ortho/drpa_matrix.irp.f b/src/tc_bi_ortho/drpa_matrix.irp.f new file mode 100644 index 00000000..56891ca2 --- /dev/null +++ b/src/tc_bi_ortho/drpa_matrix.irp.f @@ -0,0 +1,116 @@ + +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 new file mode 100644 index 00000000..48a786d2 --- /dev/null +++ b/src/tc_bi_ortho/tc_effect_int.irp.f @@ -0,0 +1,39 @@ + + +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 new file mode 100644 index 00000000..c9818a1d --- /dev/null +++ b/src/tc_bi_ortho/tc_rpa.irp.f @@ -0,0 +1,181 @@ +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 ebb13781..785d6539 100644 --- a/src/utils/util.irp.f +++ b/src/utils/util.irp.f @@ -579,5 +579,24 @@ logical function is_same_spin(sigma_1, sigma_2) end function is_same_spin ! --- + +function Kronecker_delta(i, j) result(delta) + BEGIN_DOC + ! Kronecker Delta + END_DOC + + implicit none + integer, intent(in) :: i, j + double precision :: delta + + if(i == j) then + delta = 1.d0 + else + delta = 0.d0 + endif + +end function Kronecker_delta + +! ---