fully complex G0W0 consistent with real one

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,&
+        tmp = num*cmplx(eps/(eps**2 + eta_tilde**2),&
-                -eta_tilde/(eps**2 + eta_tilde**2)**2,kind=8)
+                -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,&

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 !
 !------------------------!

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_diagonalize_matrix_without_sort(2*nS,RPA_matrix,OmOmminus)
-    call complex_orthogonalize_matrix(2*nS,RPA_matrix)
+    call complex_sort_eigenvalues_RPA(2*nS,OmOmminus,RPA_matrix)
-    !call complex_diagonalize_matrix_without_sort(2*nS,RPA_matrix,OmOmminus)
+    call complex_normalize_RPA(nS,RPA_matrix)
-    !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) 
+    XpY(:,:) = transpose(RPA_matrix(1:nS,1:nS) + RPA_matrix(nS+1:2*nS,1:nS)) 
-    XmY(:,:) = RPA_matrix(1:nS,1:nS) - RPA_matrix(1:nS,nS+1:2*nS) 
+    XmY(:,:) = transpose(RPA_matrix(1:nS,1:nS) - RPA_matrix(nS+1:2*nS,1:nS))
-    call complex_matout(nS,nS,XpY)
     deallocate(RPA_matrix,OmOmminus) 
   end if
 

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)
 

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,i) = cmplx(1d0,0d0,kind=8)
-    A(i+nS,i+nS) = -1
+    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)