From eb6e0bcc02f9cc58fe2748f4e098e8c59a256b2a Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Fri, 30 Aug 2024 19:10:54 +0200 Subject: [PATCH] added non-sym diag with biorthog condition --- src/LR/ppLR.f90 | 138 +++++---- src/make_ninja.py | 4 +- src/utils/non_sym_diag.f90 | 571 +++++++++++++++++++++++++++++++++++++ 3 files changed, 657 insertions(+), 56 deletions(-) create mode 100644 src/utils/non_sym_diag.f90 diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 2e54490..64d533c 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -1,93 +1,123 @@ -subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) -! Solve the pp-RPA linear eigenvalue problem +! --- + +subroutine ppLR(TDA, nOO, nVV, Bpp, Cpp, Dpp, Om1, X1, Y1, Om2, X2, Y2, EcRPA) + + ! + ! Solve the pp-RPA linear eigenvalue problem + ! + ! right eigen-problem: H R = R w + ! left eigen-problem: H.T L = L w + ! + ! where L.T R = 1 + ! + ! + ! (+C +B) + ! H = ( ) where C = C.T and D = D.T + ! (-B.T -D) + ! + ! (w1 0) (X1 X2) (+X1 +X2) + ! w = ( ), R = ( ) and L = ( ) + ! (0 w2) (Y1 Y2) (-Y1 -Y2) + ! + ! + ! the normalisation condition reduces to + ! + ! X1.T X2 - Y1.T Y2 = 0 + ! X1.T X1 - Y1.T Y1 = 1 + ! X2.T X2 - Y2.T Y2 = 1 + ! implicit none include 'parameters.h' -! Input variables - - logical,intent(in) :: TDA - integer,intent(in) :: nOO - integer,intent(in) :: nVV - double precision,intent(in) :: Bpp(nVV,nOO) - double precision,intent(in) :: Cpp(nVV,nVV) - double precision,intent(in) :: Dpp(nOO,nOO) + logical, intent(in) :: TDA + integer, intent(in) :: nOO, nVV + double precision, intent(in) :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO) + double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV) + double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO) + double precision, intent(out) :: EcRPA -! Local variables + logical :: imp_bio, verbose + integer :: i, j, N + double precision :: EcRPA1, EcRPA2 + double precision :: thr_d, thr_nd, thr_deg + double precision,allocatable :: M(:,:), Z(:,:), Om(:) - double precision :: trace_matrix - double precision :: EcRPA1 - double precision :: EcRPA2 - double precision,allocatable :: M(:,:) - double precision,allocatable :: Z(:,:) - double precision,allocatable :: Om(:) + double precision, external :: trace_matrix -! Output variables - double precision,intent(out) :: Om1(nVV) - double precision,intent(out) :: X1(nVV,nVV) - double precision,intent(out) :: Y1(nOO,nVV) - double precision,intent(out) :: Om2(nOO) - double precision,intent(out) :: X2(nVV,nOO) - double precision,intent(out) :: Y2(nOO,nOO) - double precision,intent(out) :: EcRPA -! Memory allocation + N = nOO + nVV - allocate(M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Om(nOO+nVV)) - -!-------------------------------------------------! -! Solve the p-p eigenproblem ! -!-------------------------------------------------! -! ! -! | C B | | X1 X2 | | w1 0 | | X1 X2 | ! -! | | | | = | | | | ! -! | -Bt -D | | Y1 Y2 | | 0 w2 | | Y1 Y2 | ! -! ! -!-------------------------------------------------! + allocate(M(N,N), Z(N,N), Om(N)) if(TDA) then X1(:,:) = +Cpp(:,:) Y1(:,:) = 0d0 - if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1) + if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1) X2(:,:) = 0d0 Y2(:,:) = -Dpp(:,:) - if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2) + if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2) else - ! Diagonal blocks - + ! Diagonal blocks M( 1:nVV , 1:nVV) = + Cpp(1:nVV,1:nVV) M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO) - ! Off-diagonal blocks - + ! Off-diagonal blocks M( 1:nVV ,nVV+1:nOO+nVV) = - Bpp(1:nVV,1:nOO) M(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp(1:nVV,1:nOO)) -! call matout(nOO,nOO,Dpp) + !! Diagonalize the p-p matrix + !if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z) + !! Split the various quantities in p-p and h-h parts + !call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2) - ! Diagonalize the p-p matrix - if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z) + thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 + thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 + thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not + imp_bio = .True. ! impose bi-orthogonality + verbose = .False. + call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) - ! Split the various quantities in p-p and h-h parts + do i = 1, nOO + Om2(i) = Om(i) + do j = 1, nVV + X2(j,i) = Z(j,i) + enddo + do j = 1, nOO + Y2(j,i) = Z(nVV+j,i) + enddo + enddo - call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) + do i = 1, nVV + Om1(i) = Om(nOO+i) + do j = 1, nVV + X1(j,i) = M(j,nOO+i) + enddo + do j = 1, nOO + Y1(j,i) = M(nVV+j,nOO+i) + enddo + enddo end if -! Compute the RPA correlation energy + ! Compute the RPA correlation energy + EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp)) + EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp) + EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp) - EcRPA = 0.5d0*( sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp) ) - EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp) - EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp) - - if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & + if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' + endif + + deallocate(M, Z, Om) end subroutine + + diff --git a/src/make_ninja.py b/src/make_ninja.py index 1e72639..a5e6e84 100755 --- a/src/make_ninja.py +++ b/src/make_ninja.py @@ -89,8 +89,8 @@ FIX_ORDER_OF_LIBS=-Wl,--start-group if sys.platform in ["linux", "linux2"]: # compiler = compile_gfortran_linux -# compiler = compile_ifort_linux - compiler = compile_olympe + compiler = compile_ifort_linux +# compiler = compile_olympe elif sys.platform == "darwin": compiler = compile_gfortran_mac else: diff --git a/src/utils/non_sym_diag.f90 b/src/utils/non_sym_diag.f90 new file mode 100644 index 0000000..04279ac --- /dev/null +++ b/src/utils/non_sym_diag.f90 @@ -0,0 +1,571 @@ + +! --- + +subroutine diagonalize_nonsym_matrix(N, A, L, e_re, thr_d, thr_nd, thr_deg, imp_bio, verbose) + + ! Diagonalize a non-symmetric matrix + ! + ! Output + ! right-eigenvectors are saved in A + ! left-eigenvectors are saved in L + ! eigenvalues are saved in e = e_re + i e_im + + implicit none + + integer, intent(in) :: N + logical, intent(in) :: imp_bio, verbose + double precision, intent(in) :: thr_d, thr_nd, thr_deg + double precision, intent(inout) :: A(N,N) + double precision, intent(out) :: e_re(N), L(N,N) + + integer :: i, j, ii + integer :: lwork, info + double precision :: accu_d, accu_nd + integer, allocatable :: iorder(:), deg_num(:) + double precision, allocatable :: Atmp(:,:), Ltmp(:,:), work(:), e_im(:) + double precision, allocatable :: S(:,:) + + if(verbose) then + print*, ' Starting a non-Hermitian diagonalization ...' + print*, ' Good Luck ;)' + print*, ' imp_bio = ', imp_bio + endif + + ! --- + ! diagonalize + + allocate(Atmp(N,N), e_im(N)) + Atmp(1:N,1:N) = A(1:N,1:N) + + allocate(work(1)) + lwork = -1 + call dgeev('V', 'V', N, Atmp, N, e_re, e_im, L, N, A, N, work, lwork, info) + if(info .gt. 0) then + print*,'dgeev failed !!', info + stop + endif + + lwork = max(int(work(1)), 1) + deallocate(work) + allocate(work(lwork)) + + call dgeev('V', 'V', N, Atmp, N, e_re, e_im, L, N, A, N, work, lwork, info) + if(info .ne. 0) then + print*,'dgeev failed !!', info + stop + endif + + deallocate(Atmp, WORK) + + + ! --- + ! check if eigenvalues are real + + i = 1 + ii = 0 + do while(i .le. N) + if(dabs(e_im(i)) .gt. 1.d-12) then + ii = ii + 1 + if(verbose) then + print*, ' Warning: complex eigenvalue !' + print*, i, e_re(i), e_im(i) + if(dabs(e_im(i)/e_re(i)) .lt. 1.d-6) then + print*, ' small enouph to be igored' + else + print*, ' IMAGINARY PART IS SIGNIFANT !!!' + endif + endif + endif + i = i + 1 + enddo + + if(verbose) then + if(ii .eq. 0) print*, ' congratulations :) eigenvalues are real-valued !!' + endif + + + ! --- + ! track & sort the real eigenvalues + + allocate(Atmp(N,N), Ltmp(N,N), iorder(N)) + + do i = 1, N + iorder(i) = i + enddo + call quick_sort(e_re, iorder, N) + + Atmp(:,:) = A(:,:) + Ltmp(:,:) = L(:,:) + do i = 1, N + do j = 1, N + A(j,i) = Atmp(j,iorder(i)) + L(j,i) = Ltmp(j,iorder(i)) + enddo + enddo + + deallocate(Atmp, Ltmp, iorder) + + + + + ! --- + ! check bi-orthog + + allocate(S(N,N)) + call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .false., verbose) + + if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(N))/dble(N) .lt. thr_d)) then + + if(verbose) then + print *, ' lapack vectors are normalized and bi-orthogonalized' + endif + + elseif((accu_nd .lt. thr_nd) .and. (dabs(accu_d - dble(N)) .gt. thr_d)) then + + if(verbose) then + print *, ' lapack vectors are not normalized but bi-orthogonalized' + endif + + call check_biorthog_binormalize(N, N, L, A, thr_d, thr_nd, .true.) + call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .true., verbose) + + else + + if(verbose) then + print *, ' lapack vectors are not normalized neither bi-orthogonalized' + endif + + allocate(deg_num(N)) + call reorder_degen_eigvec(N, thr_deg, deg_num, e_re, L, A) + call impose_biorthog_degen_eigvec(N, deg_num, e_re, L, A) + deallocate(deg_num) + + call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .false., verbose) + if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(N))/dble(N) .lt. thr_d)) then + if(verbose) then + print *, ' lapack vectors are now normalized and bi-orthogonalized' + endif + elseif((accu_nd .lt. thr_nd) .and. (dabs(accu_d - dble(N)) .gt. thr_d)) then + if(verbose) then + print *, ' lapack vectors are now not normalized but bi-orthogonalized' + endif + call check_biorthog_binormalize(N, N, L, A, thr_d, thr_nd, .true.) + call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .true., verbose) + else + if(verbose) then + print*, ' bi-orthogonalization failed !' + endif + if(imp_bio) then + print*, ' bi-orthogonalization failed !' + deallocate(S) + stop + endif + endif + + endif + + deallocate(S) + return + +end + +! --- + +subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ifnot, verbose) + + implicit none + + integer, intent(in) :: n, m + logical, intent(in) :: stop_ifnot, verbose + double precision, intent(in) :: Vl(n,m), Vr(n,m) + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(out) :: accu_d, accu_nd, S(m,m) + + integer :: i, j + double precision, allocatable :: SS(:,:) + + if(verbose) then + print *, ' check bi-orthogonality' + endif + + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + dabs(S(i,i)) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + + if(verbose) then + 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-orthogonals' + endif + endif + + ! --- + + if(stop_ifnot .and. ((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) + stop + endif + +end + +! --- + +subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot) + + implicit none + + integer, intent(in) :: n, m + logical, intent(in) :: stop_ifnot + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(inout) :: Vl(n,m), Vr(n,m) + + integer :: i, j + double precision :: accu_d, accu_nd, s_tmp + double precision, allocatable :: S(:,:) + + ! --- + + allocate(S(m,m)) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + + do i = 1, m + if(S(i,i) .lt. 0.d0) then + do j = 1, n + Vl(j,i) = -1.d0 * Vl(j,i) + enddo + S(i,i) = -S(i,i) + endif + enddo + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + S(i,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + + ! --- + + if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then + + do i = 1, m + if(S(i,i) <= 0.d0) then + print *, ' overap negative' + print *, i, S(i,i) + exit + endif + if(dabs(S(i,i) - 1.d0) .gt. thr_d) then + s_tmp = 1.d0 / dsqrt(S(i,i)) + do j = 1, n + Vl(j,i) = Vl(j,i) * s_tmp + Vr(j,i) = Vr(j,i) * s_tmp + enddo + endif + + enddo + + endif + + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + S(i,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + + deallocate(S) + + ! --- + + if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) ) then + print *, accu_nd, thr_nd + print *, dabs(accu_d-dble(m))/dble(m), thr_d + print *, ' biorthog_binormalize failed !' + stop + endif + +end + +! --- + +subroutine reorder_degen_eigvec(n, thr_deg, deg_num, e0, L0, R0) + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: thr_deg + double precision, intent(inout) :: e0(n), L0(n,n), R0(n,n) + integer, intent(out) :: deg_num(n) + + logical :: complex_root + integer :: i, j, k, m, ii, j_tmp + double precision :: ei, ej, de + double precision :: accu_d, accu_nd + double precision :: e0_tmp, L0_tmp(n), R0_tmp(n) + double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:) + + do i = 1, n + deg_num(i) = 1 + enddo + + do i = 1, n-1 + ei = e0(i) + + ! already considered in degen vectors + if(deg_num(i) .eq. 0) cycle + + ii = 0 + do j = i+1, n + ej = e0(j) + de = dabs(ei - ej) + + if(de .lt. thr_deg) then + ii = ii + 1 + + j_tmp = i + ii + + deg_num(j_tmp) = 0 + + e0_tmp = e0(j_tmp) + e0(j_tmp) = e0(j) + e0(j) = e0_tmp + + L0_tmp(1:n) = L0(1:n,j_tmp) + L0(1:n,j_tmp) = L0(1:n,j) + L0(1:n,j) = L0_tmp(1:n) + + R0_tmp(1:n) = R0(1:n,j_tmp) + R0(1:n,j_tmp) = R0(1:n,j) + R0(1:n,j) = R0_tmp(1:n) + endif + enddo + + deg_num(i) = ii + 1 + enddo + + ii = 0 + do i = 1, n + if(deg_num(i) .gt. 1) then + 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 + +end + +! --- + +subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0) + + implicit none + + integer, intent(in) :: n, deg_num(n) + double precision, intent(in) :: e0(n) + double precision, intent(inout) :: L0(n,n), R0(n,n) + + logical :: complex_root + integer :: i, j, k, m + double precision :: ei, ej, de + double precision :: accu_d, accu_nd + double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:) + + !do i = 1, n + ! if(deg_num(i) .gt. 1) then + ! print *, ' degen on', i, deg_num(i), e0(i) + ! endif + !enddo + + ! --- + + do i = 1, n + m = deg_num(i) + + if(m .gt. 1) then + + allocate(L(n,m), R(n,m), S(m,m)) + + do j = 1, m + L(1:n,j) = L0(1:n,i+j-1) + R(1:n,j) = R0(1:n,i+j-1) + enddo + + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + + accu_nd = 0.d0 + do j = 1, m + do k = 1, m + if(j==k) cycle + accu_nd = accu_nd + dabs(S(j,k)) + enddo + enddo + + if(accu_nd .lt. 1d-12) then + deallocate(S, L, R) + cycle + endif + + call impose_biorthog_svd(n, m, L, R) + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + accu_nd = 0.d0 + do j = 1, m + do k = 1, m + if(j==k) cycle + accu_nd = accu_nd + dabs(S(j,k)) + enddo + enddo + if(accu_nd .gt. 1d-12) then + print*, ' accu_nd =', accu_nd + print*, ' your strategy for degenerates orbitals failed !' + print*, m, 'deg on', i + stop + endif + + deallocate(S) + + ! --- + + do j = 1, m + L0(1:n,i+j-1) = L(1:n,j) + R0(1:n,i+j-1) = R(1:n,j) + enddo + + deallocate(L, R) + + endif + enddo + +end + +! --- + +subroutine impose_biorthog_svd(n, m, L, R) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(inout) :: L(n,m), R(n,m) + + integer :: i, j, num_linear_dependencies + double precision :: threshold + double precision, allocatable :: S(:,:), tmp(:,:) + double precision, allocatable :: U(:,:), V(:,:), Vt(:,:), D(:) + + allocate(S(m,m)) + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + + ! --- + + allocate(U(m,m), Vt(m,m), D(m)) + + call svd(S, m, U, m, D, Vt, m, m, m) + + deallocate(S) + + threshold = 1.d-6 + num_linear_dependencies = 0 + do i = 1, m + if(abs(D(i)) <= threshold) then + D(i) = 0.d0 + num_linear_dependencies = num_linear_dependencies + 1 + else + D(i) = 1.d0 / dsqrt(D(i)) + endif + enddo + if(num_linear_dependencies > 0) then + write(*,*) ' linear dependencies = ', num_linear_dependencies + write(*,*) ' m = ', m + stop + endif + + allocate(V(m,m)) + do i = 1, m + do j = 1, m + V(j,i) = Vt(i,j) + enddo + enddo + deallocate(Vt) + + ! --- + + ! R <-- R x V x D^{-0.5} + ! L <-- L x U x D^{-0.5} + + do i = 1, m + do j = 1, m + V(j,i) = V(j,i) * D(i) + U(j,i) = U(j,i) * D(i) + enddo + enddo + + allocate(tmp(n,m)) + tmp(:,:) = R(:,:) + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , tmp, size(tmp, 1), V, size(V, 1) & + , 0.d0, R, size(R, 1)) + + tmp(:,:) = L(:,:) + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , tmp, size(tmp, 1), U, size(U, 1) & + , 0.d0, L, size(L, 1)) + + deallocate(tmp, U, V, D) + +end + +! --- +