diff --git a/src/GW/complex_RGW_self_energy_diag.f90 b/src/GW/complex_RGW_self_energy_diag.f90 index 1b8d61c..797e253 100644 --- a/src/GW/complex_RGW_self_energy_diag.f90 +++ b/src/GW/complex_RGW_self_energy_diag.f90 @@ -56,7 +56,8 @@ subroutine complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,O eps = Re_e(p) - Re_e(i) + real(Om(m)) eta_tilde = eta - Im_e(p) + Im_e(i) - aimag(Om(m)) num = 2d0*rho(p,i,m)**2 - tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),eta_tilde/(eps**2+eta_tilde**2),kind=8) + tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& + eta_tilde/(eps**2+eta_tilde**2),kind=8) Re_Sig(p) = Re_Sig(p) + real(tmp) Im_Sig(p) = Im_Sig(p) + aimag(tmp) tmp = num*cmplx(-(eps**2-eta_tilde**2)/(eps**2 + eta_tilde**2)**2,& @@ -77,8 +78,8 @@ subroutine complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_e,Im_e,O eps = Re_e(p) - Re_e(a) - real(Om(m)) eta_tilde = eta + Im_e(p) - Im_e(a) - aimag(Om(m)) num = 2d0*rho(p,a,m)**2 - tmp = num*cmplx(eps/(eps**2 + eta_tilde**2)**2,& - -eta_tilde/(eps**2 + eta_tilde**2)**2,kind=8) + tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),& + -eta_tilde/(eps**2 + eta_tilde**2),kind=8) Re_Sig(p) = Re_Sig(p) + real(tmp) Im_Sig(p) = Im_Sig(p) + aimag(tmp) tmp = num*cmplx(-(eps**2 - eta_tilde**2)/(eps**2 + eta_tilde**2)**2,& diff --git a/src/GW/complex_cRG0W0.f90 b/src/GW/complex_cRG0W0.f90 index c5c8e40..a97a1d1 100644 --- a/src/GW/complex_cRG0W0.f90 +++ b/src/GW/complex_cRG0W0.f90 @@ -119,7 +119,6 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) Re_Om(:) = real(Om(:)) Im_Om(:) = aimag(Om(:)) - call complex_vecout(nS,Om) !if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om) @@ -128,7 +127,6 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, !--------------------------! call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) - write(*,*) rho(1,1,1) !------------------------! ! Compute GW self-energy ! !------------------------! diff --git a/src/LR/complex_phRLR.f90 b/src/LR/complex_phRLR.f90 index ce17ef2..8b2261b 100644 --- a/src/LR/complex_phRLR.f90 +++ b/src/LR/complex_phRLR.f90 @@ -18,7 +18,6 @@ subroutine complex_phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) complex*16 :: t1, t2 complex*16,allocatable :: RPA_matrix(:,:) complex*16,allocatable :: Z(:,:) - complex*16,allocatable :: tmp(:,:) complex*16,allocatable :: OmOmminus(:) ! Output variables @@ -36,6 +35,7 @@ subroutine complex_phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) XpY(:,:) = Aph(:,:) call complex_diagonalize_matrix(nS,XpY,Om) + call complex_orthogonalize_matrix(nS,XpY) XpY(:,:) = transpose(XpY(:,:)) XmY(:,:) = XpY(:,:) @@ -46,19 +46,16 @@ subroutine complex_phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) RPA_matrix(1:nS,nS+1:2*nS) = Bph(:,:) RPA_matrix(nS+1:2*nS,1:nS) = -Bph(:,:) RPA_matrix(nS+1:2*nS,nS+1:2*nS) = -Aph(:,:) - call complex_diagonalize_matrix(2*nS,RPA_matrix,OmOmminus) - call complex_orthogonalize_matrix(2*nS,RPA_matrix) - !call complex_diagonalize_matrix_without_sort(2*nS,RPA_matrix,OmOmminus) - !call complex_sort_eigenvalues_RPA(2*nS,OmOmminus,RPA_matrix) - !call complex_normalize_RPA(nS,RPA_matrix) + call complex_diagonalize_matrix_without_sort(2*nS,RPA_matrix,OmOmminus) + call complex_sort_eigenvalues_RPA(2*nS,OmOmminus,RPA_matrix) + call complex_normalize_RPA(nS,RPA_matrix) Om(:) = OmOmminus(1:nS) if(maxval(abs(OmOmminus(1:nS)+OmOmminus(nS+1:2*nS))) > 1e-12) & call print_warning('We dont find a Om and -Om structure as solution of the RPA. There might be a problem somewhere.') if(minval(abs(Om(:))) < 0d0) & call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') - XpY(:,:) = RPA_matrix(1:nS,1:nS) + RPA_matrix(1:nS,nS+1:2*nS) - XmY(:,:) = RPA_matrix(1:nS,1:nS) - RPA_matrix(1:nS,nS+1:2*nS) - call complex_matout(nS,nS,XpY) + XpY(:,:) = transpose(RPA_matrix(1:nS,1:nS) + RPA_matrix(nS+1:2*nS,1:nS)) + XmY(:,:) = transpose(RPA_matrix(1:nS,1:nS) - RPA_matrix(nS+1:2*nS,1:nS)) deallocate(RPA_matrix,OmOmminus) end if diff --git a/src/LR/phRLR.f90 b/src/LR/phRLR.f90 index 4eb74bc..8e9bd85 100644 --- a/src/LR/phRLR.f90 +++ b/src/LR/phRLR.f90 @@ -69,13 +69,12 @@ subroutine phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) call print_warning('You may have instabilities in linear response: negative excitations!!') Om = sqrt(Om) - call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBSq,size(AmBSq,1),0d0,XpY,size(XpY,1)) call DA(nS,1d0/dsqrt(Om),XpY) call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBIv,size(AmBIv,1),0d0,XmY,size(XmY,1)) call DA(nS,1d0*dsqrt(Om),XmY) - + ! XpY = matmul(transpose(Z),AmBSq) ! call DA(nS,1d0/sqrt(Om),XpY) diff --git a/src/utils/complex_orthogonalization.f90 b/src/utils/complex_orthogonalization.f90 index 7f8a551..883ffe1 100644 --- a/src/utils/complex_orthogonalization.f90 +++ b/src/utils/complex_orthogonalization.f90 @@ -32,15 +32,12 @@ subroutine complex_orthonormalize(N,vectors,A) complex*16, intent(inout) :: A(N, N) ! Local variables - integer :: i, j complex*16,allocatable :: L(:,:),Linv(:,:),tmp(:,:) - complex*16 :: proj - complex*16 :: norm ! Copy the input matrix to a temporary matrix allocate(L(N,N),Linv(N,N),tmp(N,N)) tmp = matmul(transpose(vectors),A) - L = matmul(L,vectors) + L = matmul(tmp,vectors) call complex_cholesky_decomp(N,L) call complex_inverse_matrix(N,L,Linv) vectors = matmul(vectors,transpose(Linv)) @@ -62,8 +59,8 @@ subroutine complex_normalize_RPA(nS,XYYX) allocate(A(2*nS,2*nS)) A(:,:) = cmplx(0d0,0d0,kind=8) do i=1,nS - A(i,i) = 1 - A(i+nS,i+nS) = -1 + A(i,i) = cmplx(1d0,0d0,kind=8) + A(i+nS,i+nS) = cmplx(-1d0,0d0,kind=8) end do call complex_orthonormalize(2*nS,XYYX,A) deallocate(A) @@ -142,6 +139,9 @@ subroutine complex_cholesky_decomp(n,A) end do if (i > j) then + if(abs(A(j,j))<1e-8) then + call print_warning('Diagonalelement in Cholesky Element is smaller than 1e-8.') + end if A(i, j) = s / A(j, j) ! Compute lower triangular elements else A(i, i) = sqrt(s)