ppBSE upfolded written

src/GW/RGW_phBSE.f90

@@ -1,5 +1,5 @@
 subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, & 
-                     nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE)
+                     nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE)
 
 ! Compute the Bethe-Salpeter excitation energies
 
@@ -18,16 +18,16 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
   logical,intent(in)            :: triplet
 
   double precision,intent(in)   :: eta
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
-  double precision,intent(in)   :: eW(nBas)
+  double precision,intent(in)   :: eW(nOrb)
-  double precision,intent(in)   :: eGW(nBas)
+  double precision,intent(in)   :: eGW(nOrb)
-  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
-  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int(nOrb,nOrb,ncart)
 
 ! Local variables
 
@@ -61,10 +61,19 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
 
 ! Memory allocation
 
-  allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
+  allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), &
            Aph(nS,nS),Bph(nS,nS),KA_sta(nS,nS),KB_sta(nS,nS), &
            OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS))
 
+!-----!
+! TDA !
+!-----!
+
+  if(TDA) then
+    write(*,*) 'Tamm-Dancoff approximation activated in phBSE!'
+    write(*,*)
+  end if
+
 ! Initialization
 
   EcBSE(:) = 0d0
@@ -76,23 +85,14 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
   isp_W = 1
   EcRPA = 0d0
 
-                 call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
+                 call phLR_A(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
   call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
-  call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
 
-  call RGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
+  call RGW_phBSE_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
-  call RGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+  call RGW_phBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
-
-!-----!
-! TDA !
-!-----!
-
-  if(TDA) then
-    write(*,*) 'Tamm-Dancoff approximation activated in phBSE!'
-    write(*,*)
-  end if
 
 !-------------------
 ! Singlet manifold
@@ -104,23 +104,23 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
 
     ! Compute BSE excitation energies
 
-                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+                 call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+    if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
     ! Second-order BSE static kernel
   
     if(dophBSE2) then 
 
-      allocate(W(nBas,nBas,nBas,nBas))
+      allocate(W(nOrb,nOrb,nOrb,nOrb))
 
       write(*,*) 
       write(*,*) '*** Second-order BSE static kernel activated! ***'
       write(*,*) 
 
-      call RGW_phBSE_static_kernel(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W)
+      call RGW_phBSE_static_kernel(eta,nOrb,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W)
-      call RGW_phBSE2_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta)
+      call RGW_phBSE2_static_kernel_A(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta)
 
-      if(.not.TDA) call RGW_phBSE2_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KB_sta)
+      if(.not.TDA) call RGW_phBSE2_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,1d0,eW,W,KB_sta)
 
       deallocate(W)
 
@@ -132,22 +132,22 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
     call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
     call print_excitation_energies('phBSE@GW@RHF','singlet',nS,OmBSE)
-    call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
+    call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
-
-    !--------------------!
-    ! Cumulant expansion !
-    !--------------------!
-
-    ! call RGWC(.false.,nBas,nC,nO,nR,nS,OmBSE,rho_RPA,eGW)
 
     !----------------------------------------------------!
     ! Compute the dynamical screening at the phBSE level !
     !----------------------------------------------------!
 
     if(dBSE) &
-        call RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
+        call RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nOrb,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
                                            OmBSE,XpY_BSE,XmY_BSE,KA_sta,KB_sta)
 
+    !----------------!
+    ! Upfolded phBSE !
+    !----------------!
+
+    call RGW_phBSE_upfolded(ispin,nOrb,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW)
+
   end if
 
 !-------------------
@@ -160,8 +160,8 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
 
     ! Compute BSE excitation energies
 
-                 call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
+                 call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-    if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+    if(.not.TDA) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
 
                  Aph(:,:) = Aph(:,:) + KA_sta(:,:)
     if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:)
@@ -169,16 +169,22 @@ subroutine RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triple
     call phLR(TDA,nS,Aph,Bph,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE)
 
     call print_excitation_energies('phBSE@GW@RHF','triplet',nS,OmBSE)
-    call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
+    call phLR_transition_vectors(.false.,nOrb,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)
 
     !-------------------------------------------------
     ! Compute the dynamical screening at the BSE level
     !-------------------------------------------------
 
     if(dBSE) &
-        call RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
+        call RGW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nOrb,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
                                            OmBSE,XpY_BSE,XmY_BSE,KA_sta,KB_sta)
 
+    !----------------!
+    ! Upfolded phBSE !
+    !----------------!
+
+    call RGW_phBSE_upfolded(ispin,nOrb,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW)
+
   end if
 
 ! Scale properly correlation energy if exchange is included in interaction kernel

src/GW/RGW_phBSE_upfolded.f90

@@ -1,12 +1,13 @@
-subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
+subroutine RGW_phBSE_upfolded(ispin,nBas,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
 
-! Upfolded phBSE@GW (TDA singlets only!)
+! Upfolded phBSE@GW (TDA only!)
 
   implicit none
   include 'parameters.h'
 
 ! Input variables
 
+  integer,intent(in)            :: ispin
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
@@ -14,10 +15,9 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
-  double precision,intent(in)   :: ENuc
-  double precision,intent(in)   :: ERHF
   double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
-  double precision,intent(in)   :: eHF(nOrb)
+  double precision,intent(in)   :: rho(nOrb,nOrb,nS)
+  double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: eGW(nOrb)
 
 ! Local variables
@@ -25,17 +25,17 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
   integer                       :: s
   integer                       :: i,j,k,l
   integer                       :: a,b,c,d
-  integer                       :: ia,jb,kc,iajb,kcld
+  integer                       :: m,ia,jb,iam,kcm
   integer,parameter             :: maxH = 20
   double precision              :: tmp
 
-  integer                       :: n1h1p,n2h2p,nH
+  integer                       :: n1h,n1p,n1h1p,n2h2p,nH
   double precision,external     :: Kronecker_delta
 
   integer,allocatable           :: order(:)
   double precision,allocatable  :: H(:,:)
   double precision,allocatable  :: X(:,:)
-  double precision,allocatable  :: Om(:)
+  double precision,allocatable  :: OmBSE(:)
   double precision,allocatable  :: Z(:)
 
 ! Output variables
@@ -55,13 +55,16 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 
 ! Dimension of the supermatrix
 
-  n1h1p = nO*nV
+  n1h = nO
-  n2h2p = nO*nO*nV*nV
+  n1p = nV
+
+  n1h1p = n1h*n1p
+  n2h2p = n1h1p*n1h1p
      nH = n1h1p + n2h2p + n2h2p
 
 ! Memory allocation
 
-  allocate(order(nH),H(nH,nH),X(nH,nH),Om(nH),Z(nH))
+  allocate(order(nH),H(nH,nH),X(nH,nH),OmBSE(nH),Z(nH))
 
 ! Initialization
 
@@ -79,9 +82,11 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 !                           !
 !---------------------------!
 
-  !---------!
+  !----------------------!
-  ! Block A !
+  ! Block A for singlets !
-  !---------!
+  !----------------------!
+
+  if(ispin == 1) then
 
     ia = 0
     do i=nC+1,nO
@@ -102,6 +107,35 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
       end do
     end do
 
+  end if
+
+  !----------------------!
+  ! Block A for triplets !
+  !----------------------!
+
+  if(ispin == 2) then
+
+    ia = 0
+    do i=nC+1,nO
+      do a=nO+1,nOrb-nR
+      ia = ia + 1
+ 
+      jb = 0
+      do j=nC+1,nO
+        do b=nO+1,nOrb-nR
+        jb = jb + 1
+ 
+          H(ia,jb) = (eGW(a) - eGW(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                   - ERI(i,b,j,a)
+ 
+          end do
+        end do
+ 
+      end do
+    end do
+
+  end if
+
   !------------------!
   ! Blocks Jph & Kph !
   !------------------!
@@ -111,56 +145,45 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
     do a=nO+1,nOrb-nR
       ia = ia + 1
 
+      kcm = 0
       do k=nC+1,nO
-        do c=nO+1,nOrn-nR
+        do c=nO+1,nOrb-nR
           do m=1,nS
             kcm = kcm + 1
 
-            tmp = sqrt(2d0)*Kronecker_delta(k,j)*ERI(b,a,c,i)
+            ! Jph
-            H(n1h1p+iajb,kc              ) = +tmp
-            H(kc        ,n1h1p+n2h2p+iajb) = -tmp
 
-            tmp = sqrt(2d0)*Kronecker_delta(b,c)*ERI(a,k,i,j)
+            tmp = - sqrt(2d0)*Kronecker_delta(a,c)*rho(i,k,m)
-            H(n1h1p+n2h2p+iajb,kc        ) = +tmp
+            H(ia             ,n1h1p+kcm) = tmp
-            H(kc              ,n1h1p+iajb) = -tmp
+            H(n1h1p+n2h2p+kcm,ia)        = tmp
+
+            ! Kph
+           
+            tmp = sqrt(2d0)*Kronecker_delta(i,k)*rho(a,c,m)
+            H(ia       ,n1h1p+n2h2p+kcm) = tmp
+            H(n1h1p+kcm,ia)              = tmp
 
           end do
         end do
+      end do
 
     end do
   end do
-    end do
-  end do
 
   !-------------!
   ! Block 2h2p !
   !-------------!
 
-  iajb = 0
+  iam = 0
   do i=nC+1,nO
     do a=nO+1,nOrb-nR
-      do j=nC+1,nO
+      do m=1,nS
-        do b=nO+1,nOrb-nR
+        iam = iam + 1
-          iajb = iajb + 1
 
-          kcld = 0
+        tmp = Om(m) + eGW(a) - eGW(i) 
-          do k=nC+1,nO
+        H(n1h1p      +iam,n1h1p      +iam) = tmp
-            do c=nO+1,nOrb-nR
+        H(n1h1p+n2h2p+iam,n1h1p+n2h2p+iam) = tmp
-              do l=nC+1,nO
-                do d=nO+1,nOrb-nR
-                kcld = kcld + 1
 
-                tmp = ((eHF(a) + eGW(b) - eHF(i) - eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & 
-                    + 2d0*ERI(a,k,i,c))*Kronecker_delta(j,l)*Kronecker_delta(b,d)
-                H(n1h1p      +iajb,n1h1p      +kcld) = tmp
-                H(n1h1p+n2h2p+iajb,n1h1p+n2h2p+kcld) = tmp
-
-                end do
-              end do
-            end do
-          end do
-
-        end do
       end do
     end do
   end do
@@ -169,15 +192,13 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 ! Diagonalize supermatrix !
 !-------------------------!
 
-! call matout(nH,nH,H)
+  call diagonalize_general_matrix(nH,H,OmBSE,X)
-
-  call diagonalize_general_matrix(nH,H,Om,X)
 
   do s=1,nH
     order(s) = s
   end do
 
-  call quick_sort(Om,order,nH)
+  call quick_sort(OmBSE,order,nH)
   call set_order(X,order,nH,nH)
 
 !-----------------!
@@ -196,16 +217,16 @@ subroutine RGW_phBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 !--------------!
 
   write(*,*)'-------------------------------------------'
-  write(*,*)'   unfolded BSE excitation energies (eV)   '
+  write(*,*)'  Upfolded phBSE excitation energies (eV)  '
   write(*,*)'-------------------------------------------'
   write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
-            '|','#','|','Omega (eV)','|','Z','|'
+            '|','#','|','OmBSE (eV)','|','Z','|'
   write(*,*)'-------------------------------------------'
 
   do s=1,min(nH,maxH)
     if(Z(s) > 1d-7) &
       write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
-      '|',s,'|',Om(s)*HaToeV,'|',Z(s),'|'
+      '|',s,'|',OmBSE(s)*HaToeV,'|',Z(s),'|'
   end do
 
   write(*,*)'-------------------------------------------'

src/GW/RGW_ppBSE.f90

@@ -76,12 +76,24 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
 
   double precision,intent(out)  :: EcBSE(nspin)
 
+!-----!
+! TDA !
+!-----!
+
+  if(TDA) then
+    write(*,*) 'Tamm-Dancoff approximation activated in ppBSE!'
+    write(*,*)
+  end if
+
+! Initialization
+
+  EcBSE(:) = 0d0
+
 !---------------------------------
 ! Compute (singlet) RPA screening 
 !---------------------------------
 
   isp_W = 1
-  EcRPA = 0d0
 
   allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), &
            Aph(nS,nS),Bph(nS,nS))
@@ -108,7 +120,6 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
     write(*,*) 
 
     ispin = 1
-    EcBSE(ispin) = 0d0
 
     nOO = nO*(nO+1)/2
     nVV = nV*(nV+1)/2
@@ -259,7 +270,6 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
     write(*,*) 
 
     ispin = 2
-    EcBSE(ispin) = 0d0
 
     nOO = nO*(nO-1)/2
     nVV = nV*(nV-1)/2
@@ -354,6 +364,13 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS
       call RGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
                                           Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
 
+ 
+    !----------------!
+    ! Upfolded ppBSE !
+    !----------------!
+
+    call RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho_RPA,OmRPA,eGW)
+
     deallocate(KB_sta,KC_sta,KD_sta)
     deallocate(Om1,X1,Y1,Om2,X2,Y2)
 

src/GW/RGW_ppBSE_upfolded.f90

@@ -1,23 +1,22 @@
-subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
+subroutine RGW_ppBSE_upfolded(ispin,nOrb,nC,nO,nV,nR,nS,ERI,rho,Om,eGW)
 
-! Upfolded ppBSE@GW (TDA triplets only!)
+! Upfolded ppBSE@GW (TDA only!)
 
   implicit none
   include 'parameters.h'
 
 ! Input variables
 
-  integer,intent(in)            :: nBas
+  integer,intent(in)            :: ispin
   integer,intent(in)            :: nOrb
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
   integer,intent(in)            :: nV
   integer,intent(in)            :: nR
   integer,intent(in)            :: nS
-  double precision,intent(in)   :: ENuc
-  double precision,intent(in)   :: ERHF
   double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
-  double precision,intent(in)   :: eHF(nOrb)
+  double precision,intent(in)   :: rho(nOrb,nOrb,nS)
+  double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: eGW(nOrb)
 
 ! Local variables
@@ -25,17 +24,17 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
   integer                       :: s
   integer                       :: i,j,k,l
   integer                       :: a,b,c,d
-  integer                       :: ia,jb,kc,iajb,kcld
+  integer                       :: m,ij,kl,ijm
   integer,parameter             :: maxH = 20
-  double precision              :: tmp
+  double precision              :: tmp,tmp1,tmp2,tmp3,tmp4
 
-  integer                       :: n1h1p,n2h2p,nH
+  integer                       :: n1h,n1p,n2h,n2p,n1h1p,n3h1p,n3p1h,n2h2p,nH
   double precision,external     :: Kronecker_delta
 
   integer,allocatable           :: order(:)
   double precision,allocatable  :: H(:,:)
   double precision,allocatable  :: X(:,:)
-  double precision,allocatable  :: Om(:)
+  double precision,allocatable  :: OmBSE(:)
   double precision,allocatable  :: Z(:)
 
 ! Output variables
@@ -55,13 +54,33 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 
 ! Dimension of the supermatrix
 
-  n1h1p = nO*nV
+  n1h = nO
-  n2h2p = nO*nO*nV*nV
+  n1p = nV
-     nH = n1h1p + n2h2p + n2h2p
+
+  if(ispin == 1) then
+
+    n2h = nO*(nO+1)/2
+    n2p = nV*(nV+1)/2
+
+  end if
+
+  if(ispin == 2) then
+
+    n2h = nO*(nO-1)/2
+    n2p = nV*(nV-1)/2
+
+  end if
+
+  n1h1p = n1h*n1p
+
+  n3h1p = n2h*n1h1p
+  n3p1h = n2p*n1h1p
+
+  nH = n1h1p + 4*n3h1p 
 
 ! Memory allocation
 
-  allocate(order(nH),H(nH,nH),X(nH,nH),Om(nH),Z(nH))
+  allocate(order(nH),H(nH,nH),X(nH,nH),OmBSE(nH),Z(nH))
 
 ! Initialization
 
@@ -83,9 +102,38 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 !                                        !
 !----------------------------------------!
 
-  !---------!
+  !----------------------!
-  ! Block D !
+  ! Block D for singlets !
-  !---------!
+  !----------------------!
+
+  if(ispin == 1) then
+
+    ij = 0
+    do i=nC+1,nO
+      do j=i,nO
+      ij = ij + 1
+ 
+      kl = 0
+      do k=nC+1,nO
+        do l=k,nO
+        kl = kl + 1
+ 
+            H(ij,kl) = (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+                     - (ERI(i,j,k,l) + ERI(i,j,l,k))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
+ 
+          end do
+        end do
+ 
+      end do
+    end do
+
+  end if
+
+  !----------------------!
+  ! Block D for triplets !
+  !----------------------!
+
+  if(ispin == 2) then
 
     ij = 0
     do i=nC+1,nO
@@ -97,8 +145,8 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
         do l=k+1,nO
         kl = kl + 1
  
-          H(ij,kl) = - (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
+            H(ij,kl) = (eGW(i) + eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
-                   + (ERI(i,j,k,l) - ERI(i,j,l,k))
+                     - (ERI(i,j,k,l) - ERI(i,j,l,k))
  
           end do
         end do
@@ -106,8 +154,10 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
       end do
     end do
 
+  end if
+
   !----------------!
-  ! Blocks M1      !
+  ! Blocks M1 & M2 !
   !----------------!
 
   ijm = 0
@@ -121,21 +171,22 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
           do l=k+1,nO
             kl = kl + 1
 
-            tmp = sqrt(2d0)*Kronecker_delta(k,j)*M(i,k,m)
+            tmp1 = sqrt(2d0)*Kronecker_delta(j,l)*rho(i,k,m)
-            H(n2h+ijm,kl     ) = +tmp
+            tmp2 = sqrt(2d0)*Kronecker_delta(j,l)*rho(i,l,m)
-            H(kl     ,n2h+ijm) = -tmp
+            tmp3 = sqrt(2d0)*Kronecker_delta(i,l)*rho(j,k,m)
+            tmp4 = sqrt(2d0)*Kronecker_delta(i,k)*rho(j,l,m)
 
-            tmp = sqrt(2d0)*Kronecker_delta(k,j)*M(i,k,m)
+            H(n2h+0*n3h1p+ijm,kl             ) = -tmp1
-            H(n2h+1*n3h1p+ijm,kl     ) = +tmp
+            H(kl             ,n2h+0*n3h1p+ijm) = +tmp3
-            H(kl     ,n2h+n3h1p+ijm) = -tmp
            
-            tmp = sqrt(2d0)*Kronecker_delta(b,c)*M(j,k,m)
+            H(n2h+1*n3h1p+ijm,kl             ) = -tmp1
-            H(n2h+2*n3h1p+iajb,kc        ) = +tmp
+            H(kl             ,n2h+1*n3h1p+ijm) = +tmp4
-            H(kc              ,n1h1p+iajb) = -tmp
            
-            tmp = sqrt(2d0)*Kronecker_delta(b,c)*M(j,k,m)
+            H(n2h+2*n3h1p+ijm,kl             ) = -tmp2
-            H(n2h+3*n2h2p+iajb,kc        ) = +tmp
+            H(kl             ,n2h+2*n3h1p+ijm) = +tmp3
-            H(kc              ,n1h1p+iajb) = -tmp
+
+            H(n2h+3*n3h1p+ijm,kl             ) = -tmp2
+            H(kl             ,n2h+3*n3h1p+ijm) = +tmp4
 
           end do
         end do
@@ -143,54 +194,39 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
       end do
     end do
   end do
-  end do
 
-  !-------------!
+  !------------!
-  ! Block 2h2p !
+  ! Block 3h1p !
-  !-------------!
+  !------------!
 
-  iajb = 0
+  ijm = 0
   do i=nC+1,nO
-    do a=nO+1,nOrb-nR
+    do j=i+1,nO
-      do j=nC+1,nO
+      do m=1,nS
-        do b=nO+1,nOrb-nR
+        ijm = ijm + 1
-          iajb = iajb + 1
 
-          kcld = 0
+        tmp = - eGW(i) - eGW(j) + Om(m)
-          do k=nC+1,nO
-            do c=nO+1,nOrb-nR
-              do l=nC+1,nO
-                do d=nO+1,nOrb-nR
-                kcld = kcld + 1
 
-                tmp = ((eHF(a) + eGW(b) - eHF(i) - eGW(j))*Kronecker_delta(i,k)*Kronecker_delta(a,c) & 
+        H(n2h+0*n3h1p+ijm,n2h+0*n3h1p+ijm) = tmp
-                    + 2d0*ERI(a,k,i,c))*Kronecker_delta(j,l)*Kronecker_delta(b,d)
+        H(n2h+1*n3h1p+ijm,n2h+1*n3h1p+ijm) = tmp
-                H(n1h1p      +iajb,n1h1p      +kcld) = tmp
+        H(n2h+2*n3h1p+ijm,n2h+2*n3h1p+ijm) = tmp
-                H(n1h1p+n2h2p+iajb,n1h1p+n2h2p+kcld) = tmp
+        H(n2h+3*n3h1p+ijm,n2h+3*n3h1p+ijm) = tmp
 
       end do
     end do
   end do
-          end do
-
-        end do
-      end do
-    end do
-  end do
 
 !-------------------------!
 ! Diagonalize supermatrix !
 !-------------------------!
 
-! call matout(nH,nH,H)
+  call diagonalize_general_matrix(nH,H,OmBSE,X)
-
-  call diagonalize_general_matrix(nH,H,Om,X)
 
   do s=1,nH
     order(s) = s
   end do
 
-  call quick_sort(Om,order,nH)
+  call quick_sort(OmBSE,order,nH)
   call set_order(X,order,nH,nH)
 
 !-----------------!
@@ -199,8 +235,8 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 
   Z(:) = 0d0
   do s=1,nH
-    do ia=1,n1h1p
+    do ij=1,n2h
-      Z(s) = Z(s) + X(ia,s)**2
+      Z(s) = Z(s) + X(ij,s)**2
     end do
   end do
 
@@ -209,16 +245,16 @@ subroutine RGW_ppBSE_upfolded(nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW)
 !--------------!
 
   write(*,*)'-------------------------------------------'
-  write(*,*)'   unfolded BSE excitation energies (eV)   '
+  write(*,*)'  Upfolded ppBSE excitation energies (eV)  '
   write(*,*)'-------------------------------------------'
   write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X)') &
-            '|','#','|','Omega (eV)','|','Z','|'
+            '|','#','|','OmBSE (eV)','|','Z','|'
   write(*,*)'-------------------------------------------'
 
   do s=1,min(nH,maxH)
     if(Z(s) > 1d-7) &
       write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
-      '|',s,'|',Om(s)*HaToeV,'|',Z(s),'|'
+      '|',s,'|',OmBSE(s)*HaToeV,'|',Z(s),'|'
   end do
 
   write(*,*)'-------------------------------------------'