bug

examples/molecule.O3

@@ -1,6 +1,6 @@
-# nAt nEl nCore nRyd
+# nAt nEla nElb nCore nRyd
-  3   24  6     0 
+  3  12 12      6     0 
 # Znuc   x            y           z
-8. 0.0000  0.0000 0.0000
+O 0.0000  0.0000 0.0000
-8. 0.0000  2.05696689 1.26554959
+O 0.0000  2.05696689 1.26554959
-8. 0.0000 -2.05696689 1.26554959
+O 0.0000 -2.05696689 1.26554959

input/basis

@@ -1,9 +1,14 @@
-1   3
+1   2
 S   3   1.00
-     38.3600000              0.0238090        
+     18.7311370              0.03349460
-      5.7700000              0.1548910        
+      2.8253937              0.23472695
-      1.2400000              0.4699870        
+      0.6401217              0.81375733
 S   1   1.00
-      0.2976000              1.0000000        
+      0.1612778              1.0000000
-P   1   1.00
+2   2
-      1.2750000              1.0000000 
+S   3   1.00
+     18.7311370              0.03349460
+      2.8253937              0.23472695
+      0.6401217              0.81375733
+S   1   1.00
+      0.1612778              1.0000000

input/methods

@@ -9,8 +9,8 @@
 # GF2 GF3
   F   F
 # G0W0 evGW qsGW 
-  T    F    F
-# G0T0 evGT qsGT 
   F    F    F
+# G0T0 evGT qsGT 
+  T    F    F
 # MCMP2
   F

input/molecule

@@ -1,4 +1,5 @@
 # nAt nEla nElb nCore nRyd
-  1   1     1   0     0
+ 2   1     1   0     0
 # Znuc   x            y           z
-  He      0.0          0.0         0.0
+  H     0.           0.         0.
+  H     0.           0.         1.4

input/options

@@ -9,6 +9,6 @@
 # GF: maxSCF thresh  DIIS n_diis renormalization
       64     0.00001 T    10      3 
 # GW: maxSCF thresh  DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta
-      256    0.0000001  T  5      T      F     F   F   F   F   F         0.000
+      256    0.0000001  T  5      F      F     F   F   F   F   F         0.000
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

input/weight

@@ -1,9 +1,14 @@
-1   3
+1   2
 S   3   1.00
-     38.3600000              0.0238090        
+     18.7311370              0.03349460
-      5.7700000              0.1548910        
+      2.8253937              0.23472695
-      1.2400000              0.4699870        
+      0.6401217              0.81375733
 S   1   1.00
-      0.2976000              1.0000000        
+      0.1612778              1.0000000
-P   1   1.00
+2   2
-      1.2750000              1.0000000 
+S   3   1.00
+     18.7311370              0.03349460
+      2.8253937              0.23472695
+      0.6401217              0.81375733
+S   1   1.00
+      0.1612778              1.0000000

src/QuAcK/G0T0.f90

@@ -38,7 +38,7 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,
 
 ! Output variables
 
-  double precision              :: eG0T0(nBas)
+  double precision,intent(out)  :: eG0T0(nBas)
 
 ! Hello world
 
@@ -50,20 +50,20 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,
 
 ! Dimensions of the rr-RPA linear reponse matrices
 
-  nOOs = nO*(nO+1)/2
+  nOOs = nO*(nO + 1)/2
-  nVVs = nV*(nV+1)/2
+  nVVs = nV*(nV + 1)/2
 
-  nOOt = nO*(nO-1)/2
+  nOOt = nO*(nO - 1)/2
-  nVVt = nV*(nV-1)/2
+  nVVt = nV*(nV - 1)/2
 
 ! Memory allocation
 
   allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & 
            Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & 
-           rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & 
+           rho1s(nBas,nO-nC,nVVs),rho2s(nBas,nV-nR,nOOs), & 
            Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & 
-           Omega2t(nOOs),X2t(nVVs,nOOs),Y2t(nOOs,nOOs), & 
+           Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & 
-           rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), & 
+           rho1t(nBas,nO-nC,nVVt),rho2t(nBas,nV-nR,nOOt), & 
            SigT(nBas),Z(nBas))
 
 !----------------------------------------------
@@ -85,7 +85,7 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,
 
 ! Compute excitation densities for the T-matrix
 
-  call excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOOs,nVVs,ERI(:,:,:,:), & 
+  call excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & 
                                   X1s(:,:),Y1s(:,:),rho1s(:,:,:),                & 
                                   X2s(:,:),Y2s(:,:),rho2s(:,:,:))
 
@@ -95,7 +95,7 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,
 
  ispin = 2
 
-  ! Compute linear response
+! Compute linear response
 
   call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR, & 
                           nOOt,nVVt,eHF(:),ERI(:,:,:,:),  & 
@@ -106,9 +106,9 @@ subroutine G0T0(BSE,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,ENuc,
   call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t(:))
   call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t(:))
 
-  ! Compute excitation densities for the T-matrix
+! Compute excitation densities for the T-matrix
 
-  call excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOOt,nVVt,ERI(:,:,:,:), & 
+  call excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & 
                                   X1t(:,:),Y1t(:,:),rho1t(:,:,:),                & 
                                   X2t(:,:),Y2t(:,:),rho2t(:,:,:))
 

src/QuAcK/MP2F12.f90

@@ -60,8 +60,8 @@ subroutine MP2F12(nBas,nC,nO,nV,ERI,F12,Yuk,FC,EHF,e,c)
       do a=1,nV
         do b=1,nV
           eps = eO(i) + eO(j) - eV(a) - eV(b)
-          EcMP2a = EcMP2a + ooCoo(i,j,a,b)/eps
+          EcMP2a = EcMP2a + ooCvv(i,j,a,b)/eps
-          EcMP2b = EcMP2b + ooCoo(i,j,b,a)/eps
+          EcMP2b = EcMP2b + ooCvv(i,j,b,a)/eps
         enddo
       enddo
     enddo

src/QuAcK/QuAcK.f90

@@ -534,8 +534,8 @@ program QuAcK
   if(doG0T0) then
     
     call cpu_time(start_G0T0)
-    call G0T0(BSE,singlet_manifold,triplet_manifold,eta, & 
+!   call G0T0(BSE,singlet_manifold,triplet_manifold,eta, & 
-              nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0)
+!             nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0)
     call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
     call cpu_time(end_G0T0)
   

src/QuAcK/excitation_density_Tmatrix.f90

@@ -1,4 +1,4 @@
-subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
+subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
 
 ! Compute excitation densities for T-matrix self-energy
 
@@ -7,12 +7,12 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1
 ! Input variables
 
   integer,intent(in)            :: ispin
-  integer,intent(in)            :: nBas,nC,nO,nR,nOO,nVV
+  integer,intent(in)            :: nBas,nC,nO,nV,nR,nOO,nVV
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
-  double precision,intent(out)  :: X1(nVV,nVV)
+  double precision,intent(in)   :: X1(nVV,nVV)
-  double precision,intent(out)  :: Y1(nOO,nVV)
+  double precision,intent(in)   :: Y1(nOO,nVV)
-  double precision,intent(out)  :: X2(nVV,nOO)
+  double precision,intent(in)   :: X2(nVV,nOO)
-  double precision,intent(out)  :: Y2(nOO,nOO)
+  double precision,intent(in)   :: Y2(nOO,nOO)
 
 ! Local variables
 
@@ -213,6 +213,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nR,nOO,nVV,ERI,X1,Y1,rho1
 
         end do
       end do
+
     end do
 
   end if

src/QuAcK/linear_response.f90

@@ -78,7 +78,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,EcRP
  
   Omega = sqrt(Omega)
   XpY = matmul(transpose(Z),AmBSq)
-  call DA(nS,1d0/sqrt(Omega),XpY)
+  call DA(nS,1d0/sqrt(abs(Omega)),XpY)
 
 ! print*,'X+Y'
 ! call matout(nS,nS,XpY)

src/QuAcK/linear_response_pp.f90

@@ -16,6 +16,7 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
   
 ! Local variables
 
+  integer                       :: p,q
   double precision              :: trace_matrix
   double precision,allocatable  :: B(:,:)
   double precision,allocatable  :: C(:,:)
@@ -63,13 +64,18 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
   M(    1:nVV    ,nVV+1:nOO+nVV) = -           B(1:nVV,1:nOO)
   M(nVV+1:nOO+nVV,    1:nVV)     = + transpose(B(1:nVV,1:nOO))
 
+! do p=1,nOO+nVV
+!   do q=1,nOO+nVV
+!     write(42,*) p,q,M(p,q)
+!   end do
+! end do
+
 ! print*, 'pp-RPA matrix'
 ! call matout(nOO+nVV,nOO+nVV,M(:,:))
 
 ! Diagonalize the p-h matrix
 
-  Z(:,:) = M(:,:)
+  call diagonalize_general_matrix(nOO+nVV,M,Omega,Z)
-  call diagonalize_general_matrix(nOO+nVV,M(:,:),Omega(:),Z(:,:))
 
 ! write(*,*) 'pp-RPA excitation energies'
 ! call matout(nOO+nVV,1,Omega(:))
@@ -96,4 +102,10 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1
 
 ! print*,'Ec(pp-RPA) = ',EcppRPA
 
+! print*,'Eigenvalues'
+! call matout(nOO+nVV,1,Omega)
+
+! print*,'Eigenvectors'
+! call matout(nOO+nVV,nOO+nVV,matmul(transpose(Z),Z))
+
 end subroutine linear_response_pp

src/QuAcK/print_G0T0.f90

@@ -9,9 +9,12 @@ subroutine print_G0T0(nBas,nO,e,ENuc,ERHF,SigT,Z,eGW,EcRPA)
   double precision,intent(in)        :: ENuc
   double precision,intent(in)        :: ERHF
   double precision,intent(in)        :: EcRPA(nspin)
-  double precision,intent(in)        :: e(nBas),SigT(nBas),Z(nBas),eGW(nBas)
+  double precision,intent(in)        :: e(nBas)
+  double precision,intent(in)        :: SigT(nBas)
+  double precision,intent(in)        :: Z(nBas)
+  double precision,intent(in)        :: eGW(nBas)
 
-  integer                            :: x,HOMO,LUMO
+  integer                            :: p,HOMO,LUMO
   double precision                   :: Gap
 
 ! HOMO and LUMO
@@ -29,9 +32,9 @@ subroutine print_G0T0(nBas,nO,e,ENuc,ERHF,SigT,Z,eGW,EcRPA)
             '|','#','|','e_HF (eV)','|','Sigma_T (eV)','|','Z','|','e_QP (eV)','|'
   write(*,*)'-------------------------------------------------------------------------------'
 
-  do x=1,nBas
+  do p=1,nBas
     write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
-    '|',x,'|',e(x)*HaToeV,'|',SigT(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|'
+    '|',p,'|',e(p)*HaToeV,'|',SigT(p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|'
   enddo
 
   write(*,*)'-------------------------------------------------------------------------------'

src/QuAcK/soG0T0.f90

@@ -36,9 +36,6 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
   double precision,allocatable  :: seHF(:)
   double precision,allocatable  :: sERI(:,:,:,:)
 
-! Output variables
-
-
 ! Hello world
 
   write(*,*)
@@ -82,45 +79,43 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 
 ! Compute linear response
 
-  call linear_response_pp(ispin,.false.,nBas2,nC2,nO2,nV2,nR2, & 
+  call linear_response_pp(ispin,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV, & 
-                          nOO,nVV,seHF(:),sERI(:,:,:,:),  & 
+                          seHF,sERI,Omega1,X1,Y1,Omega2,X2,Y2,         & 
-                          Omega1(:),X1(:,:),Y1(:,:),   & 
-                          Omega2(:),X2(:,:),Y2(:,:),   & 
                           EcRPA)
 
-  call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:))
+! call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:))
-  call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:))
+! call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:))
 
 ! Compute excitation densities for the T-matrix
 
-  call excitation_density_Tmatrix(ispin,nBas2,nC2,nO2,nR2,nOO,nVV,sERI(:,:,:,:), & 
+! call excitation_density_Tmatrix(ispin,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), & 
-                                  X1(:,:),Y1(:,:),rho1(:,:,:),             & 
+!                                 X1(:,:),Y1(:,:),rho1(:,:,:),             & 
-                                  X2(:,:),Y2(:,:),rho2(:,:,:))
+!                                 X2(:,:),Y2(:,:),rho2(:,:,:))
 
 !----------------------------------------------
 ! Compute T-matrix version of the self-energy 
 !----------------------------------------------
 
-  call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & 
+! call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & 
-                                Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), & 
+!                               Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), & 
-                                SigT(:))
+!                               SigT(:))
 
 ! Compute renormalization factor for T-matrix self-energy
 
-  call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & 
+! call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & 
-                                      Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), & 
+!                                     Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), & 
-                                      Z(:))
+!                                     Z(:))
 
 !----------------------------------------------
 ! Solve the quasi-particle equation
 !----------------------------------------------
 
-  eG0T0(:) = seHF(:) + Z(:)*SigT(:)
+! eG0T0(:) = seHF(:) + Z(:)*SigT(:)
 
 !----------------------------------------------
 ! Dump results
 !----------------------------------------------
 
-  call print_G0T0(nBas2,nO2,seHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA)
+! call print_G0T0(nBas2,nO2,seHF(:),ENuc,ERHF,SigT(:),Z(:),eG0T0(:),EcRPA)
 
 end subroutine soG0T0

src/QuAcK/sort_ppRPA.f90

@@ -28,7 +28,12 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 ! Initializatiom
 
   Omega1(:) = 0d0
+  X1(:,:)   = 0d0
+  Y1(:,:)   = 0d0
+
   Omega2(:) = 0d0
+  X2(:,:)   = 0d0
+  Y2(:,:)   = 0d0
 
   ab = 0
   ij = 0
@@ -56,12 +61,20 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
   if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
   if(maxval(Omega2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
 
-! write(*,*) 'pp-RPA positive excitation energies'
+  write(*,*) 'pp-RPA positive excitation energies'
-! call matout(nVV,1,Omega1(:))
+  call matout(nVV,1,Omega1(:))
-! write(*,*)
+  write(*,*)
+
+  write(*,*) 'pp-RPA negative excitation energies'
+  call matout(nOO,1,Omega2(:))
+  write(*,*)
+
+! write(*,*) 'X1.X1 - Y1.Y1'
+! call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1))
+! write(*,*) 'X2.X2 - Y2.Y2'
+! call matout(nOO,nOO,matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2))
+! write(*,*) 'X1.X2 - Y1.Y2'
+! call matout(nVV,nOO,matmul(transpose(X1),X2) - matmul(transpose(Y1),Y2))
 
-! write(*,*) 'pp-RPA negative excitation energies'
-! call matout(nOO,1,Omega2(:))
-! write(*,*)
 
 end subroutine sort_ppRPA

src/utils/wrap_lapack.f90

@@ -30,7 +30,7 @@
 !
 !end subroutine diagonalize_matrix_lowest
 
-subroutine diagonalize_general_matrix(N,A,e,X)
+subroutine diagonalize_general_matrix(N,A,WR,VR)
 
 ! Diagonalize a non-symmetric square matrix
 
@@ -38,27 +38,55 @@ subroutine diagonalize_general_matrix(N,A,e,X)
 
 ! Input variables
 
+  integer :: i,j,k
   integer,intent(in)            :: N
   double precision,intent(inout):: A(N,N)
-  double precision,intent(out)  :: X(N,N)
+  double precision,intent(out)  :: VR(N,N)
-  double precision,intent(out)  :: e(N)
+  double precision,intent(out)  :: WR(N)
 
 ! Local variables
 
   integer                       :: lwork,info
-  double precision,allocatable  :: work(:),WI(:),VL(:,:)
+  double precision,allocatable  :: work(:),WI(:),VL(:,:),tmp1(:,:),tmp2(:,:)
 
 ! Memory allocation
 
   lwork = 4*N
-  allocate(WI(N),VL(N,N),work(lwork))
+  allocate(WI(N),VL(N,N),work(lwork),tmp1(N,N),tmp2(N,N))
-
+  tmp1 = A
-  call dgeev('N','V',N,A,N,e,WI,VL,N,X,N,work,lwork,info)
+  call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
 
   if(info /= 0) then 
-    print*,'Problem in diagonalize_matrix (dgeev)!!'
+    print*,'Problem in diagonalize_general_matrix (dgeev)!!'
   endif
 
+  call matout(N,1,WI)
+
+  tmp2 = 0d0
+  do i=1,N
+    do j=1,N
+      do k=1,N
+        tmp2(i,j) = tmp2(i,j) + vl(k,i)*tmp1(k,j)
+      end do
+    end do
+  end do
+
+print*,'tmp2'
+call matout(N,N,tmp2)
+
+  tmp1 = 0d0
+  do i=1,N
+    do j=1,N
+      tmp1(i,j) = wr(i)*vl(i,j)
+    end do
+  end do
+
+print*,'tmp1'
+call matout(N,N,tmp1)
+
+  print*,'coucou'
+  print*,maxval(tmp1-tmp2),minval(tmp1-tmp2)
+
 end subroutine diagonalize_general_matrix
 
 subroutine diagonalize_matrix(N,A,e)