proper spin adaptation in GTpp

src/GT/RG0T0pp.f90

@@ -40,11 +40,10 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
 ! Local variables
 
-  logical                       :: print_T   = .false.
+  logical                       :: print_T   = .true.
   logical                       :: plot_self = .false.
 
   integer                       :: isp_T
-  integer                       :: iblock
   integer                       :: nOOs,nOOt
   integer                       :: nVVs,nVVt
   integer                       :: n_states, n_states_diag
@@ -90,11 +89,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 
 ! Dimensions of the pp-RPA linear reponse matrices
 
-  !nOOs = nO*(nO + 1)/2
+  nOOs = nO*(nO + 1)/2
-  !nVVs = nV*(nV + 1)/2
+  nVVs = nV*(nV + 1)/2
-
-  nOOs = nO*nO
-  nVVs = nV*nV
 
   nOOt = nO*(nO - 1)/2
   nVVt = nV*(nV - 1)/2
@@ -114,60 +110,54 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 !----------------------------------------------
 
   isp_T  = 1
-  !iblock = 1
-  iblock = 3
 
 ! Compute linear response
 
   allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
 
-  if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
+  if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
-                 call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
+                 call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp)
-                 call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
+                 call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp)
 
   call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
 
   deallocate(Bpp,Cpp,Dpp)
 
-  if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s)
+  if(print_T) call print_excitation_energies('ppRPA@RHF','2p (singlets)',nVVs,Om1s)
-  if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-beta)',nOOs,Om2s)
+  if(print_T) call print_excitation_energies('ppRPA@RHF','2h (singlets)',nOOs,Om2s)
 
 !----------------------------------------------
 ! alpha-alpha block
 !----------------------------------------------
 
   isp_T  = 2
-! iblock = 2
-  iblock = 4
 
 ! Compute linear response
 
   allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
 
-  if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
+  if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
-                 call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
+                 call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp)
-                 call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
+                 call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp)
 
   call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
 
   deallocate(Bpp,Cpp,Dpp)
 
-  if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t)
+  if(print_T) call print_excitation_energies('ppRPA@RHF','2p (triplets)',nVVt,Om1t)
-  if(print_T) call print_excitation_energies('ppRPA@RHF','2h (alpha-alpha)',nOOt,Om2t)
+  if(print_T) call print_excitation_energies('ppRPA@RHF','2h (triplets)',nOOt,Om2t)
 
 !----------------------------------------------
 ! Compute excitation densities
 !----------------------------------------------
 
-! iblock = 1
+  isp_T = 1
-  iblock = 3
 
-  call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
+  call RGTpp_excitation_density(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
 
-! iblock = 2
+  isp_T = 2
-  iblock = 4
 
-  call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
+  call RGTpp_excitation_density(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
 
 !----------------------------------------------
 ! Compute T-matrix version of the self-energy 
@@ -214,34 +204,30 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d
 ! Compute the ppRPA correlation energy
 
   isp_T  = 1
-! iblock = 1
-  iblock = 3
 
   allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs))
 
-  if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
+  if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp)
-                 call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
+                 call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp)
-                 call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
+                 call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp)
 
   call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T))
 
   deallocate(Bpp,Cpp,Dpp)
 
   isp_T  = 2
-! iblock = 2
-  iblock = 4
 
   allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt))
 
-  if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
+  if(.not.TDA_T) call ppLR_B(isp_T,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp)
-                 call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
+                 call ppLR_C(isp_T,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp)
-                 call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
+                 call ppLR_D(isp_T,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp)
 
   call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T))
 
   deallocate(Bpp,Cpp,Dpp)
 
-  EcRPA(1) = EcRPA(1) - EcRPA(2)
+  EcRPA(1) = 1d0*EcRPA(1)
   EcRPA(2) = 3d0*EcRPA(2)
 
   call print_RG0T0pp(nOrb,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA)

src/GT/RGTpp_excitation_density.f90

@@ -2,12 +2,6 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
 
 ! Compute excitation densities for T-matrix self-energy
 
-  ! TODO
-  ! debug DGEMM for nC != 0
-  !             and nR != 0
-
-
-
   implicit none
 
 ! Input variables
@@ -72,8 +66,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
                     cd = cd + 1
                     rho1(p,q,ab) = rho1(p,q,ab) & 
                                  + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/ & 
-                                   (1d0 + Kronecker_delta(c,d))
+                                   sqrt(1d0 + Kronecker_delta(c,d))/sqrt(2d0)
-!                                  sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d)))
                  end do
               end do
           
@@ -83,8 +76,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
                     kl = kl + 1
                     rho1(p,q,ab) = rho1(p,q,ab) & 
                                  + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/ & 
-                                   (1d0 + Kronecker_delta(k,l))
+                                   sqrt(1d0 + Kronecker_delta(k,l))/sqrt(2d0)
-!                                  sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l)))
                  end do
               end do
               
@@ -101,8 +93,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
                     cd = cd + 1
                     rho2(p,q,ij) = rho2(p,q,ij) &
                                  + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/ & 
-                                   (1d0 + Kronecker_delta(c,d))
+                                   sqrt(1d0 + Kronecker_delta(c,d))/sqrt(2d0)
-!                                  sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(c,d)))
                  end do
               end do
               
@@ -112,8 +103,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
                     kl = kl + 1
                     rho2(p,q,ij) = rho2(p,q,ij) &
                                  + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/ & 
-                                   (1d0 + Kronecker_delta(k,l))
+                                   sqrt(1d0 + Kronecker_delta(k,l))/sqrt(2d0)
-!                                  sqrt((1d0 + Kronecker_delta(p,q))*(1d0 + Kronecker_delta(k,l)))
                  end do
               end do
  
@@ -158,7 +148,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
                 cd = cd + 1
   
                 rho1(p,q,ab) = rho1(p,q,ab) & 
-                             + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) 
+                             + sqrt(3d0/2d0)*(ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) 
               end do ! d
             end do ! c
    
@@ -169,7 +159,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
                 kl = kl + 1
   
                 rho1(p,q,ab) = rho1(p,q,ab) & 
-                             + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) 
+                             + sqrt(3d0/2d0)*(ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) 
               end do ! l
             end do ! k 
           end do ! b
@@ -189,7 +179,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
                 cd = cd + 1
             
                 rho2(p,q,ij) = rho2(p,q,ij) & 
-                             + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) 
+                             + sqrt(3d0/2d0)*(ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) 
               end do ! d
             end do ! c
             
@@ -200,7 +190,7 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
                 kl = kl + 1
             
                 rho2(p,q,ij) = rho2(p,q,ij) & 
-                             + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) 
+                             + sqrt(3d0/2d0)*(ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) 
               end do ! l
             end do ! k
           end do ! j
@@ -259,6 +249,9 @@ subroutine RGTpp_excitation_density(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho
   
       deallocate(ERI_1, ERI_2)
   
+      rho1 = rho1*sqrt(3d0/2d0)
+      rho2 = rho2*sqrt(3d0/2d0)
+
     endif
   endif