pCCD and ppRPA

examples/molecule.LiF

@@ -2,4 +2,4 @@
  2   6     6   0     0
 # Znuc   x            y           z
   Li    0.           0.         0.
-  F     0.           0.         3.5
+  F     0.           0.         2.954767548731

examples/molecule.LiH

@@ -2,4 +2,4 @@
  2   2     2   0     0
 # Znuc   x            y           z
   Li    0.           0.         0.
-  H     0.           0.         3.5
+  H     0.           0.         3.051018991463

input/basis

@@ -1,24 +1,30 @@
-1   9
+1   6
-S   4
+S   8
-  1    234.0000000              0.0025870        
+  1  17880.0000000              0.0007380        
-  2     35.1600000              0.0195330        
+  2   2683.0000000              0.0056770        
-  3      7.9890000              0.0909980        
+  3    611.5000000              0.0288830        
-  4      2.2120000              0.2720500        
+  4    173.5000000              0.1085400        
+  5     56.6400000              0.2909070        
+  6     20.4200000              0.4483240        
+  7      7.8100000              0.2580260        
+  8      1.6530000              0.0150630        
+S   8
+  1  17880.0000000             -0.0001720        
+  2   2683.0000000             -0.0013570        
+  3    611.5000000             -0.0067370        
+  4    173.5000000             -0.0276630        
+  5     56.6400000             -0.0762080        
+  6     20.4200000             -0.1752270        
+  7      7.8100000             -0.1070380        
+  8      1.6530000              0.5670500        
 S   1
-  1      0.6669000              1.0000000        
+  1      0.4869000              1.0000000        
-S   1
+P   3
-  1      0.2089000              1.0000000        
+  1     28.3900000              0.0460870        
-S   1
+  2      6.2700000              0.2401810        
-  1      0.0513800              1.0000000        
+  3      1.6950000              0.5087440        
 P   1
-  1      3.0440000              1.0000000        
+  1      0.4317000              1.0000000        
-P   1
-  1      0.7580000              1.0000000        
-P   1
-  1      0.1993000              1.0000000        
 D   1
-  1      1.9650000              1.0000000        
+  1      2.2020000              1.0000000
-D   1
-  1      0.4592000              1.0000000
-
 

input/methods

@@ -3,7 +3,7 @@
 # MP2 MP3 MP2-F12
   T   F   F
 # CCD CCSD CCSD(T) drCCD rCCD lCCD pCCD
-  F   F    F       F      F   F    F
+  F   F    F       F      F   F    T
 # CIS RPA RPAx ppRPA ADC
   F   F   F    F     F
 # G0F2 evGF2 G0F3 evGF3
@@ -11,6 +11,6 @@
 # G0W0 evGW qsGW 
   F    F    F
 # G0T0 evGT qsGT 
-  T    F    F
+  F    F    F
 # MCMP2
   F

input/molecule

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

input/molecule.xyz

@@ -1,3 +1,3 @@
   1
 
- He     0.0000000000    0.0000000000    0.0000000000
+ Ne     0.0000000000    0.0000000000    0.0000000000

input/options

@@ -9,7 +9,7 @@
 # GF: maxSCF thresh  DIIS n_diis lin renorm
       256    0.00001 T    5      T   3 
 # GW: maxSCF thresh  DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta     
-      256    0.00001    T  5     F      F     T   F   F   F   T    0.000   
+      256    0.00001    T  5     F      F     F   F   F   F   F    0.000   
 # ACFDT: AC Kx  XBS
          T  F   T
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift

input/weight

@@ -1,24 +1,30 @@
-1   9
+1   6
-S   4
+S   8
-  1    234.0000000              0.0025870        
+  1  17880.0000000              0.0007380        
-  2     35.1600000              0.0195330        
+  2   2683.0000000              0.0056770        
-  3      7.9890000              0.0909980        
+  3    611.5000000              0.0288830        
-  4      2.2120000              0.2720500        
+  4    173.5000000              0.1085400        
+  5     56.6400000              0.2909070        
+  6     20.4200000              0.4483240        
+  7      7.8100000              0.2580260        
+  8      1.6530000              0.0150630        
+S   8
+  1  17880.0000000             -0.0001720        
+  2   2683.0000000             -0.0013570        
+  3    611.5000000             -0.0067370        
+  4    173.5000000             -0.0276630        
+  5     56.6400000             -0.0762080        
+  6     20.4200000             -0.1752270        
+  7      7.8100000             -0.1070380        
+  8      1.6530000              0.5670500        
 S   1
-  1      0.6669000              1.0000000        
+  1      0.4869000              1.0000000        
-S   1
+P   3
-  1      0.2089000              1.0000000        
+  1     28.3900000              0.0460870        
-S   1
+  2      6.2700000              0.2401810        
-  1      0.0513800              1.0000000        
+  3      1.6950000              0.5087440        
 P   1
-  1      3.0440000              1.0000000        
+  1      0.4317000              1.0000000        
-P   1
-  1      0.7580000              1.0000000        
-P   1
-  1      0.1993000              1.0000000        
 D   1
-  1      1.9650000              1.0000000        
+  1      2.2020000              1.0000000
-D   1
-  1      0.4592000              1.0000000
-
 

src/QuAcK/CCD.f90

@@ -48,7 +48,7 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF)
   double precision,allocatable  :: r2(:,:,:,:)
   double precision,allocatable  :: t2(:,:,:,:)
 
-  integer                       :: n_diis
+  integer                       :: n_diis,i,j,a,b
   double precision              :: rcond
   double precision,allocatable  :: error_diis(:,:)
   double precision,allocatable  :: t_diis(:,:)

src/QuAcK/G0T0.f90

@@ -137,6 +137,9 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
 ! Compute T-matrix version of the self-energy 
 !----------------------------------------------
 
+  rho1s = 0d0
+  rho1t = 0d0
+
   call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & 
                                 Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & 
                                 Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & 

src/QuAcK/QuAcK.f90

@@ -367,6 +367,7 @@ program QuAcK
 
     call cpu_time(start_CCD)
     call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF)
+    call cpu_time(end_CCD)
 
     t_CCD = end_CCD - start_CCD
     write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds'
@@ -669,10 +670,10 @@ program QuAcK
   if(doG0T0) then
     
     call cpu_time(start_G0T0)
-!   call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,       &
-!             singlet_manifold,triplet_manifold,linGW,eta, &  
-!             nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF)
     call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
+ !  call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,       &
+ !            singlet_manifold,triplet_manifold,linGW,eta, &  
+ !            nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF)
 
     call cpu_time(end_G0T0)
   

src/QuAcK/excitation_density_Tmatrix.f90

@@ -46,9 +46,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do d=c,nBas-nR
               cd = cd + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
-!                          + ERI(p,i,c,d)*X1(cd,ab) 
+                           + ERI(p,i,c,d)*X1(cd,ab) 
-                           + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) &
+!                          + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) &
-                           /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
+!                          /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
 
             end do
           end do
@@ -58,9 +58,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do l=k,nO
               kl = kl + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
-!                          + ERI(p,i,k,l)*Y1(kl,ab)
+                           + ERI(p,i,k,l)*Y1(kl,ab)
-                           + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) &
+!                          + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) &
-                           /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
+!                          /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
 
             end do
           end do
@@ -76,9 +76,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do d=c,nBas-nR
               cd = cd + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-!                          + ERI(p,nO+a,c,d)*X2(cd,ij)
+                           + ERI(p,nO+a,c,d)*X2(cd,ij)
-                           + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) &
+!                          + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) &
-                           /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d)))
+!                          /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d)))
 
             end do
           end do
@@ -88,9 +88,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
            do l=k,nO
               kl = kl + 1
               rho2(p,a,ij) = rho2(p,a,ij) & 
-!                          + ERI(p,nO+a,k,l)*Y2(kl,ij) 
+                           + ERI(p,nO+a,k,l)*Y2(kl,ij) 
-                           + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) &
+!                          + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) &
-                           /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l)))
+!                          /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l)))
 
             end do
           end do
@@ -118,7 +118,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
               cd = cd + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
                            + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
-              print*,rho1(p,i,ab),ERI(p,i,c,d),X1(cd,ab)
 
             end do
           end do
@@ -129,7 +128,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
               kl = kl + 1
               rho1(p,i,ab) = rho1(p,i,ab) & 
                            + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) 
-              print*,rho1(p,i,ab),ERI(p,i,k,l),Y1(kl,ab)
             end do
           end do
 
@@ -222,16 +220,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
 
     end do
 
-!   do p=nC+1,nBas-nR
-!     do i=nC+1,nO
-!       do ab=1,nVV
-!         print*,p,i,ab,rho1(p,i,ab)
-!       end do
-!     end do
-!   end do
-
-! call matout(nVV,nOO,Y1)
-
   end if
 
 end subroutine excitation_density_Tmatrix

src/QuAcK/linear_response_C_pp.f90

@@ -25,7 +25,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
 ! Define the chemical potential
  
   eF = e(nO) + e(nO+1)
-! eF = 0d0
+  eF = 0d0
 
 ! Build C matrix for the singlet manifold
 

src/QuAcK/linear_response_D_pp.f90

@@ -25,7 +25,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
 ! Define the chemical potential
 
   eF = e(nO) + e(nO+1)
-! eF = 0d0
+  eF = 0d0
  
 ! Build the D matrix for the singlet manifold
 

src/QuAcK/linear_response_pp.f90

@@ -1,5 +1,5 @@
 subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & 
-                              e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcppRPA)
+                              e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcRPA)
 
 ! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013)
 
@@ -22,6 +22,8 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
   integer                       :: ab,cd,ij,kl
   integer                       :: p,q,r,s
   double precision              :: trace_matrix
+  double precision              :: EcRPA1
+  double precision              :: EcRPA2
   double precision,allocatable  :: B(:,:)
   double precision,allocatable  :: C(:,:)
   double precision,allocatable  :: D(:,:)
@@ -37,7 +39,7 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
   double precision,intent(out)  :: Omega2(nOO)
   double precision,intent(out)  :: X2(nVV,nOO)
   double precision,intent(out)  :: Y2(nOO,nOO)
-  double precision,intent(out)  :: EcppRPA
+  double precision,intent(out)  :: EcRPA
 
 ! Memory allocation
 
@@ -117,10 +119,11 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
 
 ! print*, 'pp-RPA matrix'
 ! call matout(nOO+nVV,nOO+nVV,M(:,:))
+! write(*,*) 
 
 ! Diagonalize the p-h matrix
 
-  call diagonalize_general_matrix(nOO+nVV,M,Omega,Z)
+  if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z)
 
 ! write(*,*) 'pp-RPA excitation energies'
 ! call matout(nOO+nVV,1,Omega(:))
@@ -132,11 +135,11 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
 
 ! Compute the RPA correlation energy
 
-! EcppRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
+  EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
-! print*,+sum(Omega1(:)),- trace_matrix(nVV,C(:,:))
+  EcRPA1 = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
-! print*,-sum(Omega2(:)),- trace_matrix(nOO,D(:,:))
+  EcRPA2 = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:))
-  EcppRPA = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
+  if(abs(EcRPA - EcRPA1) > 1d-10 .or. abs(EcRPA - EcRPA2) > 1d-10) & 
-! EcppRPA = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:))
+    print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'
 
 ! write(*,*)'X1'
 ! call matout(nVV,nVV,X1)
@@ -147,12 +150,4 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
 ! write(*,*)'Y2'
 ! call matout(nOO,nOO,Y2)
 
-! 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/pCCD.f90

@@ -53,6 +53,8 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
   allocate(delta_OOVV(nO,nV))
 
+  delta_OOVV(:,:) = 0d0
+
   do i=nC+1,nO
     do a=1,nV-nR
       delta_OOVV(i,a) = 2d0*(eHF(nO+a) - eHF(i))
@@ -63,6 +65,11 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
 
   allocate(OOOO(nO,nO),OOVV(nO,nV),OVOV(nO,nV),OVVO(nO,nV),VVVV(nV,nV))
 
+  OOOO(:,:) = 0d0
+  OOVV(:,:) = 0d0
+  OVVO(:,:) = 0d0
+  VVVV(:,:) = 0d0
+
   do i=nC+1,nO
     do j=nC+1,nO
       OOOO(i,j) = ERI(i,i,j,j)
@@ -141,7 +148,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
                - 2d0*(2d0*OVOV(i,a) - OVVO(i,a) - OOVV(i,a)*t(i,a))*t(i,a)
 
         do j=nC+1,nO
-          r(i,a) = r(i,a) - 2d0*OOVV(j,a)*t(j,a)*t(i,a) + OOOO(j,i)*t(j,a) + y(i,j)*t(j,a)
+          r(i,a) = r(i,a) - 2d0*OOVV(j,a)*t(j,a)*t(i,a) + OOOO(j,i)*t(j,a) + y(i,j)*t(j,a) 
         end do 
 
         do b=1,nV-nR

src/QuAcK/print_excitation.f90

@@ -7,7 +7,7 @@ subroutine print_excitation(method,ispin,nS,Omega)
 
 ! Input variables
 
-  character*6,intent(in)             :: method
+  character*12,intent(in)             :: method
   integer,intent(in)                 :: ispin,nS
   double precision,intent(in)        :: Omega(nS)
 

src/QuAcK/self_energy_Tmatrix_diag.f90

@@ -72,7 +72,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do i=nC+1,nO
       do cd=1,nVVt
         eps = e(p) + e(i) - Omega1t(cd)
-        SigT(p) = SigT(p) + rho1t(p,i,cd)*rho1t(p,i,cd)/eps
+ !      SigT(p) = SigT(p) + rho1t(p,i,cd)*rho1t(p,i,cd)/eps
       enddo
     enddo
   enddo
@@ -83,7 +83,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
     do a=1,nV-nR
       do kl=1,nOOt
         eps = e(p) + e(nO+a) - Omega2t(kl)
-        SigT(p) = SigT(p) + rho2t(p,a,kl)*rho2t(p,a,kl)/eps
+ !      SigT(p) = SigT(p) + rho2t(p,a,kl)*rho2t(p,a,kl)/eps
       enddo
     enddo
   enddo

src/QuAcK/soG0T0.f90

@@ -89,31 +89,28 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
 ! Compute excitation densities for the T-matrix
 
   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(:), & 
-                                Omega1(:),rho1(:,:,:), &
+                                Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:),  &
-                                Omega2(:),rho2(:,:,:), &
                                 SigT(:))
 
 ! Compute renormalization factor for T-matrix self-energy
 
   call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & 
-                                      Omega1(:),rho1(:,:,:), &
+                                      Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:),  &
-                                      Omega2(:),rho2(:,:,:), &
                                       Z(:))
 
 !----------------------------------------------
 ! Solve the quasi-particle equation
 !----------------------------------------------
 
-  eG0T0(:) = seHF(:) + SigT(:)
+! eG0T0(:) = seHF(:) + SigT(:)
-! eG0T0(:) = seHF(:) + Z(:)*SigT(:)
+  eG0T0(:) = seHF(:) + Z(:)*SigT(:)
 
 !----------------------------------------------
 ! Dump results

src/QuAcK/sort_ppRPA.f90

@@ -74,16 +74,12 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
       ab = ab + 1
       Omega1(ab) = Omega(pq)
       Z1(1:nOO+nVV,ab) = Z(1:nOO+nVV,pq)
-!     X1(1:nVV,ab) = Z(    1:    nVV,pq)
-!     Y1(1:nOO,ab) = Z(nVV+1:nOO+nVV,pq)
 
     else
 
       ij = ij + 1
       Omega2(ij) = Omega(pq)
       Z2(1:nOO+nVV,ij) = Z(1:nOO+nVV,pq)
-!     X2(1:nVV,ij) = Z(    1:    nVV,pq)
-!     Y2(1:nOO,ij) = Z(nVV+1:nOO+nVV,pq)
 
     end if
 
@@ -115,11 +111,13 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 
   end if
 
-  X1(1:nVV,1:nVV) = Z1(    1:    nVV,1:nVV)
+! Define submatrices X1, Y1, X2, & Y2
-  Y1(1:nOO,1:nVV) = Z1(nVV+1:nOO+nVV,1:nVV)
 
-  X2(1:nVV,1:nOO) = Z2(    1:    nVV,1:nOO)
+  X1(1:nVV,1:nVV) = + Z1(    1:    nVV,1:nVV)
-  Y2(1:nOO,1:nOO) = Z2(nVV+1:nOO+nVV,1:nOO)
+  Y1(1:nOO,1:nVV) = - Z1(nVV+1:nOO+nVV,1:nVV)
+
+  X2(1:nVV,1:nOO) = + Z2(    1:    nVV,1:nOO)
+  Y2(1:nOO,1:nOO) = - Z2(nVV+1:nOO+nVV,1:nOO)
 
 ! write(*,*) 'Z1t.M.Z1'
 ! call matout(nVV,nVV,matmul(matmul(transpose(Z1),M),Z1))

src/eDFT/LIM_RKS.f90

@@ -73,7 +73,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
 !------------------------------------------------------------------------
 
   write(*,'(A40)')           '*************************************************'
-  write(*,'(A40)')           '            ZERO-WEIGHT CALCULATION              '
+  write(*,'(A40)')           '            EQUI-WEIGHT CALCULATION              '
   write(*,'(A40)')           '*************************************************'
 
   wLIM(1:nEns) = 1d0/dble(nEns)

src/eDFT/quadrature_grid.f90

@@ -66,7 +66,7 @@ subroutine quadrature_grid(nRad,nAng,nGrid,root,weight)
     enddo
   enddo
  
-! Compute values of the basis functions (and the its gradient if required) at each grid point
+! Print grids
 
 20 format(T2,58('-'))
 30 format(T20,'Radial Quadrature',/, &