UBSE

input/options

@@ -9,10 +9,10 @@
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm 
       256    0.00001 T    5      T   0.0  3      
 # GW/GT: maxSCF thresh  DIIS n_diis lin eta        COHSEX SOSEX TDA_W G0W GW0 
-         256    0.00001   T    5     T  0.0         F      F     F     F   F  
+         256    0.00001   T    5     T  0.001   F      F     F     F   F  
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
-        F    T    F    F
+        T    F    F    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

src/QuAcK/G0W0.f90

@@ -71,23 +71,31 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,
 
 ! SOSEX correction
 
-  if(SOSEX) write(*,*) 'SOSEX correction activated!'
+  if(SOSEX) then 
-  write(*,*)
+    write(*,*) 'SOSEX correction activated!'
+    write(*,*)
+  end if
 
 ! COHSEX approximation
 
-  if(COHSEX) write(*,*) 'COHSEX approximation activated!'
+  if(COHSEX) then 
-  write(*,*)
+    write(*,*) 'COHSEX approximation activated!'
+    write(*,*)
+  end if
 
 ! TDA for W
 
-  if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
+  if(TDA_W) then 
-  write(*,*)
+    write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
+    write(*,*)
+  end if
 
 ! TDA 
 
-  if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!'
+  if(TDA) then 
-  write(*,*)
+    write(*,*) 'Tamm-Dancoff approximation activated!'
+    write(*,*)
+  end if
 
 ! Spin manifold 
 

src/QuAcK/QuAcK.f90

@@ -57,9 +57,9 @@ program QuAcK
   double precision,allocatable  :: ERI_MO(:,:,:,:)
   integer                       :: bra
   integer                       :: ket
-  double precision,allocatable  :: ERI_MO_aa(:,:,:,:)
+  double precision,allocatable  :: ERI_MO_aaaa(:,:,:,:)
-  double precision,allocatable  :: ERI_MO_ab(:,:,:,:)
+  double precision,allocatable  :: ERI_MO_aabb(:,:,:,:)
-  double precision,allocatable  :: ERI_MO_bb(:,:,:,:)
+  double precision,allocatable  :: ERI_MO_bbbb(:,:,:,:)
   double precision,allocatable  :: ERI_ERF_AO(:,:,:,:)
   double precision,allocatable  :: ERI_ERF_MO(:,:,:,:)
   double precision,allocatable  :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:)
@@ -330,25 +330,25 @@ program QuAcK
 
       ! Memory allocation
      
-      allocate(ERI_MO_aa(nBas,nBas,nBas,nBas),ERI_MO_ab(nBas,nBas,nBas,nBas),ERI_MO_bb(nBas,nBas,nBas,nBas))
+      allocate(ERI_MO_aaaa(nBas,nBas,nBas,nBas),ERI_MO_aabb(nBas,nBas,nBas,nBas),ERI_MO_bbbb(nBas,nBas,nBas,nBas))
      
       ! 4-index transform for (aa|aa) block
      
       bra = 1
       ket = 1
-      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aa)
+      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aaaa)
-      
-      ! 4-index transform for (bb|bb) block
-     
-      bra = 1
-      ket = 2
-      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_ab)
       
       ! 4-index transform for (aa|bb) block
      
+      bra = 1
+      ket = 2
+      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aabb)
+     
+      ! 4-index transform for (bb|bb) block
+     
       bra = 2
       ket = 2
-      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bb)
+      call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bbbb)
      
     else
 
@@ -382,7 +382,7 @@ program QuAcK
 
     if(unrestricted) then
 
-      call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,ENuc,EUHF,eHF,EcMP2)
+      call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ENuc,EUHF,eHF,EcMP2)
 
     else
 
@@ -749,7 +749,7 @@ program QuAcK
 
       call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
                  singlet_manifold,triplet_manifold,linGW,eta_GW,nBas,nC,nO,nV,nR,nS, & 
-                 ENuc,EUHF,Hc,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,PHF,cHF,eHF,eG0W0)
+                 ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,PHF,cHF,eHF,eG0W0)
     else
 
       call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,       & 

src/QuAcK/UG0W0.f90

@@ -1,6 +1,6 @@
 subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,  &
-                 singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,nO,nV,nR,nS, &
+                 spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS, &
-                 ENuc,EUHF,Hc,ERI_aa,ERI_ab,ERI_bb,PHF,cHF,eHF,eGW)
+                 ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,PHF,cHF,eHF,eGW)
 
 ! Perform unrestricted G0W0 calculation
 
@@ -20,8 +20,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
   logical,intent(in)            :: dBSE
   logical,intent(in)            :: dTDA
   logical,intent(in)            :: evDyn
-  logical,intent(in)            :: singlet_manifold
+  logical,intent(in)            :: spin_conserved
-  logical,intent(in)            :: triplet_manifold
+  logical,intent(in)            :: spin_flip
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
 
@@ -37,9 +37,9 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
   double precision,intent(in)   :: cHF(nBas,nBas,nspin)
   double precision,intent(in)   :: PHF(nBas,nBas,nspin)
   double precision,intent(in)   :: Hc(nBas,nBas,nspin)
-  double precision,intent(in)   :: ERI_aa(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: ERI_ab(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: ERI_bb(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
 
 ! Local variables
 
@@ -82,18 +82,24 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
 
 ! COHSEX approximation
 
-  if(COHSEX) write(*,*) 'COHSEX approximation activated!'
+  if(COHSEX) then 
-  write(*,*)
+    write(*,*) 'COHSEX approximation activated!'
+    write(*,*)
+  end if
 
 ! TDA for W
 
-  if(TDA_W) write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
+  if(TDA_W) then 
-  write(*,*)
+    write(*,*) 'Tamm-Dancoff approximation for dynamic screening!'
+    write(*,*)
+  end if
 
 ! TDA 
 
-  if(TDA) write(*,*) 'Tamm-Dancoff approximation activated!'
+  if(TDA) then 
-  write(*,*)
+    write(*,*) 'Tamm-Dancoff approximation activated!'
+    write(*,*)
+  end if
 
 ! Memory allocation
 
@@ -113,15 +119,15 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
   ispin = 1
 
   call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, &
-                                    eHF,ERI_aa,ERI_ab,ERI_bb,rho(:,:,:,ispin),EcRPA,Omega,XpY,XmY)
+                                    eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY)
 
-  if(print_W) call print_excitation('RPA@UHF',3,nSt,Omega)
+  if(print_W) call print_excitation('RPA@UHF',5,nSt,Omega)
 
 !----------------------!
 ! Excitation densities !
 !----------------------!
  
-  call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY,rho)
+  call unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
 
 !---------------------!
 ! Compute self-energy !
@@ -133,14 +139,12 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
 ! Compute renormalization factor !
 !--------------------------------!
 
-! call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eHF, & 
+  call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,eHF,Omega,rho,Z)
-!                             Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z(:))
 
 !-----------------------------------!
 ! Solve the quasi-particle equation !
 !-----------------------------------!
 
-  Z(:,:) = 1d0
   eGWlin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:)
 
   if(linearize) then 
@@ -163,14 +167,12 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
 
 ! Dump results
 
-  do ispin=1,nspin
+  call print_UG0W0(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGW,EcRPA)
-    call print_G0W0(nBas,nO(ispin),eHF(:,ispin),ENuc,EUHF,SigC(:,ispin),Z(:,ispin),eGW(:,ispin),EcRPA)
-  end do
 
 ! Compute the RPA correlation energy
 
-! call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & 
+  call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,1d0, &
-!                      rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+                                    eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY)
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
@@ -181,10 +183,11 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
 
 ! Perform BSE calculation
 
-! if(BSE) then
+  if(BSE) then
 
-!   call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, &
+    call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,  &
-!                       nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
+                                     nBas,nC,nO,nV,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb, & 
+                                     eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
 
 !   if(exchange_kernel) then
 !
@@ -239,6 +242,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
 
 !   end if
 
-! end if
+  end if
 
 end subroutine UG0W0

src/QuAcK/print_UG0W0.f90

@@ -1,44 +1,61 @@
-subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigmaC,Z,eGW,EcRPA)
+subroutine print_UG0W0(nBas,nO,e,ENuc,EHF,SigC,Z,eGW,EcRPA)
 
 ! Print one-electron energies and other stuff for G0W0
 
   implicit none
   include 'parameters.h'
 
-  integer,intent(in)                 :: nBas,nO
+  integer,intent(in)                 :: nBas
+  integer,intent(in)                 :: nO(nspin)
   double precision,intent(in)        :: ENuc
   double precision,intent(in)        :: EHF
   double precision,intent(in)        :: EcRPA
-  double precision,intent(in)        :: e(nBas),SigmaC(nBas),Z(nBas),eGW(nBas)
+  double precision,intent(in)        :: e(nBas,nspin)
+  double precision,intent(in)        :: SigC(nBas,nspin)
+  double precision,intent(in)        :: Z(nBas,nspin)
+  double precision,intent(in)        :: eGW(nBas,nspin)
 
-  integer                            :: x,HOMO,LUMO
+  integer                            :: p
+  double precision                   :: HOMO
+  double precision                   :: LUMO
   double precision                   :: Gap
 
 ! HOMO and LUMO
 
-  HOMO = nO
+  HOMO = max(eGW(nO(1),1),eGW(nO(2),2))
-  LUMO = HOMO + 1
+  LUMO = min(eGW(nO(1)+1,1),eGW(nO(2)+1,2))
-  Gap = eGW(LUMO)-eGW(HOMO)
+  Gap = LUMO - HOMO
 
 ! Dump results
 
-  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)'-------------------------------------------------------------------------------&
-  write(*,*)'  One-shot G0W0 calculation'
+              -------------------------------------------------'
-  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)'  Unrestricted one-shot G0W0 calculation (eV)'
-  write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
+  write(*,*)'-------------------------------------------------------------------------------& 
-            '|','#','|','e_HF (eV)','|','Sigma_c (eV)','|','Z','|','e_QP (eV)','|'
+              -------------------------------------------------'
-  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') &
+            '|','#','|','e_HF up','e_HF dw','|','Sig_c up','Sig_c dw','|', & 
+                        'Z up','Z dw','|','e_QP up','e_QP dw','|'
+  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)') &
+    write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') &
-    '|',x,'|',e(x)*HaToeV,'|',SigmaC(x)*HaToeV,'|',Z(x),'|',eGW(x)*HaToeV,'|'
+    '|',p,'|',e(p,1)*HaToeV,e(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', & 
+              Z(p,1),Z(p,2),'|',eGW(p,1)*HaToeV,eGW(p,2)*HaToeV,'|'
   enddo
 
-  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)'-------------------------------------------------------------------------------& 
-  write(*,'(2X,A30,F15.6)') 'G0W0 HOMO      energy (eV):',eGW(HOMO)*HaToeV
+              -------------------------------------------------'
-  write(*,'(2X,A30,F15.6)') 'G0W0 LUMO      energy (eV):',eGW(LUMO)*HaToeV
+  write(*,'(2X,A30,F15.6)') 'G0W0 HOMO      energy (eV):',HOMO*HaToeV
+  write(*,'(2X,A30,F15.6)') 'G0W0 LUMO      energy (eV):',LUMO*HaToeV
   write(*,'(2X,A30,F15.6)') 'G0W0 HOMO-LUMO gap    (eV):',Gap*HaToeV
-  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)'-------------------------------------------------------------------------------& 
+              -------------------------------------------------'
+  write(*,'(2X,A30,F15.6)') 'RPA@HF   total energy       =',ENuc + EHF + EcRPA
+  write(*,'(2X,A30,F15.6)') 'RPA@HF   correlation energy =',EcRPA
+  write(*,*)'-------------------------------------------------------------------------------& 
+              -------------------------------------------------'
   write(*,*)
 
 end subroutine print_UG0W0

src/QuAcK/print_excitation.f90

@@ -13,18 +13,20 @@ subroutine print_excitation(method,ispin,nS,Omega)
 
 ! Local variables
 
-  character*7                        :: spin_manifold
+  character*14                       :: spin_manifold
   integer,parameter                  :: maxS = 32
   integer                            :: ia
 
   if(ispin == 1) spin_manifold = 'singlet'
   if(ispin == 2) spin_manifold = 'triplet'
-  if(ispin == 3) spin_manifold = 'alp-bet'
+  if(ispin == 3) spin_manifold = 'alpha-beta'
-  if(ispin == 4) spin_manifold = 'alp-alp'
+  if(ispin == 4) spin_manifold = 'alpha-alpha'
+  if(ispin == 5) spin_manifold = 'spin-conserved'
+  if(ispin == 6) spin_manifold = 'spin-flip'
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------'
-  write(*,'(1X,A1,1X,A14,A14,A7,A9,A15)')'|',method,' calculation: ',spin_manifold,' manifold','         |'
+  write(*,'(1X,A1,1X,A14,A14,A14,A9,A13)')'|',method,' calculation: ',spin_manifold,' manifold','     |'
   write(*,*)'-------------------------------------------------------------'
   write(*,'(1X,A1,1X,A5,1X,A1,1X,A23,1X,A1,1X,A23,1X,A1,1X)') &
             '|','State','|',' Excitation energy (au) ','|',' Excitation energy (eV) ','|'

src/QuAcK/renormalization_factor.f90

@@ -36,41 +36,9 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
     Z(:) = 1d0
     return
  
-  end if
+! SOSEX correction
  
-! Occupied part of the correlation self-energy
+  elseif(SOSEX) then
-
-  do x=nC+1,nBas-nR
-    do i=nC+1,nO
-      jb = 0
-      do j=nC+1,nO
-        do b=nO+1,nBas-nR
-          jb = jb + 1
-          eps = e(x) - e(i) + Omega(jb) 
-          Z(x) = Z(x)  - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2
-        end do
-      end do
-    end do
-  end do
-
-! Virtual part of the correlation self-energy
-
-  do x=nC+1,nBas-nR
-    do a=nO+1,nBas-nR
-      jb = 0
-      do j=nC+1,nO
-        do b=nO+1,nBas-nR
-          jb = jb + 1
-          eps = e(x) - e(a) - Omega(jb) 
-          Z(x) = Z(x)  - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2
-        end do
-      end do
-    end do
-  end do
-
-  ! SOSEX correction
-
-  if(SOSEX) then
 
     ! Occupied part of the correlation self-energy
 
@@ -102,7 +70,39 @@ subroutine renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r
       end do
     end do
 
-  endif
+  else
+
+    ! Occupied part of the correlation self-energy
+
+    do x=nC+1,nBas-nR
+      do i=nC+1,nO
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+            eps = e(x) - e(i) + Omega(jb) 
+            Z(x) = Z(x)  - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2
+          end do
+        end do
+      end do
+    end do
+
+    ! Virtual part of the correlation self-energy
+
+    do x=nC+1,nBas-nR
+      do a=nO+1,nBas-nR
+        jb = 0
+        do j=nC+1,nO
+          do b=nO+1,nBas-nR
+            jb = jb + 1
+            eps = e(x) - e(a) - Omega(jb) 
+            Z(x) = Z(x)  - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2
+          end do
+        end do
+      end do
+    end do
+
+  end if
 
 ! Compute renormalization factor from derivative of SigC
  

src/QuAcK/unrestricted_excitation_density.f90

@@ -1,4 +1,4 @@
-subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_ab,ERI_bb,XpY,rho)
+subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY,rho)
 
 ! Compute excitation densities for unrestricted reference
 
@@ -14,9 +14,9 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
   integer,intent(in)            :: nSa
   integer,intent(in)            :: nSb
   integer,intent(in)            :: nSt
-  double precision,intent(in)   :: ERI_aa(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: ERI_ab(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: ERI_bb(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: XpY(nSt,nSt)
 
 ! Local variables
@@ -45,7 +45,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
           do b=nO(1)+1,nBas-nR(1)
             jb = jb + 1
 
-            rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aa(p,j,q,b)*XpY(ia,jb)
+            rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aaaa(p,j,q,b)*XpY(ia,jb)
 
           enddo
         enddo
@@ -58,7 +58,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
           do b=nO(2)+1,nBas-nR(2)
             jb = jb + 1
 
-            rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_ab(p,j,q,b)*XpY(ia,jb)
+            rho(p,q,ia,1) = rho(p,q,ia,1) + ERI_aabb(p,j,q,b)*XpY(ia,jb)
 
           enddo
         enddo
@@ -81,7 +81,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
           do b=nO(1)+1,nBas-nR(1)
             jb = jb + 1
 
-            rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_ab(j,p,b,q)*XpY(ia,jb)
+            rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_aabb(j,p,b,q)*XpY(ia,jb)
 
           enddo
         enddo
@@ -94,7 +94,7 @@ subroutine unrestricted_excitation_density(nBas,nC,nO,nR,nSa,nSb,nSt,ERI_aa,ERI_
           do b=nO(2)+1,nBas-nR(2)
             jb = jb + 1
 
-            rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bb(p,j,q,b)*XpY(ia,jb)
+            rho(p,q,ia,2) = rho(p,q,ia,2) + ERI_bbbb(p,j,q,b)*XpY(ia,jb)
 
           enddo
         enddo

src/QuAcK/unrestricted_linear_response.f90

@@ -1,5 +1,5 @@
 subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
-                                        e,ERI_aa,ERI_ab,ERI_bb,rho,EcRPA,Omega,XpY,XmY)
+                                        e,ERI_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY)
 
 ! Compute linear response for unrestricted formalism
 
@@ -24,9 +24,9 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
   double precision,intent(in)   :: lambda
   double precision,intent(in)   :: e(nBas,nspin)
   double precision,intent(in)   :: rho(nBas,nBas,nSt,nspin)
-  double precision,intent(in)   :: ERI_aa(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: ERI_ab(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: ERI_bb(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
   
 ! Local variables
 
@@ -53,18 +53,20 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
 
 ! Build A and B matrices 
 
-  call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aa,ERI_ab,ERI_bb,A)
+  call unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,A)
 
-! if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A)
+  if(BSE) & 
+    call unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A)
 
 ! Tamm-Dancoff approximation
 
   B = 0d0
   if(.not. TDA) then
 
-    call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aa,ERI_ab,ERI_bb,B)
+    call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,B)
 
-!   if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B)
+    if(BSE) &
+      call unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B)
 
   end if
 

src/QuAcK/unrestricted_linear_response_A_matrix.f90

@@ -1,5 +1,5 @@
 subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
-                                                 e,ERI_aa,ERI_ab,ERI_bb,A_lr)
+                                                 e,ERI_aaaa,ERI_aabb,ERI_bbbb,A_lr)
 
 ! Compute linear response
 
@@ -20,9 +20,9 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
   integer,intent(in)            :: nSt
   double precision,intent(in)   :: lambda
   double precision,intent(in)   :: e(nBas,nspin)
-  double precision,intent(in)   :: ERI_aa(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas) 
-  double precision,intent(in)   :: ERI_ab(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas) 
-  double precision,intent(in)   :: ERI_bb(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas) 
   
 ! Local variables
 
@@ -58,7 +58,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
             jb = jb + 1
  
             A_lr(ia,jb) = (e(a,1) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
-                        + lambda*ERI_aa(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_aa(i,b,j,a)
+                        + lambda*ERI_aaaa(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,b,j,a)
 
           end  do
         end  do
@@ -76,7 +76,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
           do b=nO(2)+1,nBas-nR(2)
             jb = jb + 1
  
-            A_lr(ia,nSa+jb) = lambda*ERI_ab(i,b,a,j) 
+            A_lr(ia,nSa+jb) = lambda*ERI_aabb(i,b,a,j) 
 
           end  do
         end  do
@@ -94,7 +94,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
           do b=nO(1)+1,nBas-nR(1)
             jb = jb + 1
  
-            A_lr(nSa+ia,jb) = lambda*ERI_ab(b,i,j,a)
+            A_lr(nSa+ia,jb) = lambda*ERI_aabb(b,i,j,a)
 
           end  do
         end  do
@@ -113,7 +113,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
             jb = jb + 1
  
             A_lr(nSa+ia,nSa+jb) = (e(a,2) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
-                                + lambda*ERI_bb(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_bb(i,b,j,a)
+                                + lambda*ERI_bbbb(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,b,j,a)
 
           end  do
         end  do
@@ -122,5 +122,15 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
 
   end if
 
+!-----------------------------------------------
+! Build A matrix for spin-flip transitions
+!-----------------------------------------------
+
+  if(ispin == 2) then 
+
+    print*,'spin-flip transition NYI'
+
+  end if
+
 
 end subroutine unrestricted_linear_response_A_matrix

src/QuAcK/unrestricted_linear_response_B_matrix.f90

@@ -1,5 +1,5 @@
 subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
-                                                 ERI_aa,ERI_ab,ERI_bb,B_lr)
+                                                 ERI_aaaa,ERI_aabb,ERI_bbbb,B_lr)
 
 ! Compute linear response
 
@@ -19,9 +19,9 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
   integer,intent(in)            :: nSb
   integer,intent(in)            :: nSt
   double precision,intent(in)   :: lambda
-  double precision,intent(in)   :: ERI_aa(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas) 
-  double precision,intent(in)   :: ERI_ab(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas) 
-  double precision,intent(in)   :: ERI_bb(nBas,nBas,nBas,nBas) 
+  double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas) 
   
 ! Local variables
 
@@ -56,7 +56,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
           do b=nO(1)+1,nBas-nR(1)
             jb = jb + 1
  
-            B_lr(ia,jb) = lambda*ERI_aa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aa(i,j,b,a)
+            B_lr(ia,jb) = lambda*ERI_aaaa(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_aaaa(i,j,b,a)
 
           end  do
         end  do
@@ -74,7 +74,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
           do b=nO(2)+1,nBas-nR(2)
             jb = jb + 1
  
-            B_lr(ia,nSa+jb) = lambda*ERI_ab(i,j,a,b) 
+            B_lr(ia,nSa+jb) = lambda*ERI_aabb(i,j,a,b) 
 
           end  do
         end  do
@@ -92,7 +92,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
           do b=nO(1)+1,nBas-nR(1)
             jb = jb + 1
  
-            B_lr(nSa+ia,jb) = lambda*ERI_ab(j,i,b,a)
+            B_lr(nSa+ia,jb) = lambda*ERI_aabb(j,i,b,a)
 
           end  do
         end  do
@@ -110,7 +110,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
           do b=nO(2)+1,nBas-nR(2)
             jb = jb + 1
  
-            B_lr(nSa+ia,nSa+jb) = lambda*ERI_bb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bb(i,j,b,a)
+            B_lr(nSa+ia,nSa+jb) = lambda*ERI_bbbb(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_bbbb(i,j,b,a)
 
           end  do
         end  do
@@ -119,5 +119,15 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
 
   end if
 
+!-----------------------------------------------
+! Build A matrix for spin-flip transitions
+!-----------------------------------------------
+
+  if(ispin == 2) then
+
+    print*,'spin-flip transition NYI'
+  
+  end if
+
 
 end subroutine unrestricted_linear_response_B_matrix

src/QuAcK/unrestricted_renormalization_factor.f90

@@ -0,0 +1,93 @@
+subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,e,Omega,rho,Z)
+
+! Compute the renormalization factor in the unrestricted formalism
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  double precision,intent(in)   :: eta
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nSa
+  integer,intent(in)            :: nSb
+  integer,intent(in)            :: nSt
+  double precision,intent(in)   :: e(nBas,nspin)
+  double precision,intent(in)   :: Omega(nSt)
+  double precision,intent(in)   :: rho(nBas,nBas,nSt,nspin)
+
+! Local variables
+
+  integer                       :: i,j,a,b,p,q,jb
+  double precision              :: eps
+
+! Output variables
+
+  double precision,intent(out)  :: Z(nBas,nspin)
+
+! Initialize 
+
+  Z(:,:) = 0d0
+
+!--------------!
+! Spin-up part !
+!--------------!
+
+  ! Occupied part of the renormalization factor
+
+  do p=nC(1)+1,nBas-nR(1)
+    do i=nC(1)+1,nO(1)
+      do jb=1,nSt
+        eps = e(p,1) - e(i,1) + Omega(jb)
+        Z(p,1) = Z(p,1) + rho(p,i,jb,1)**2*(eps/(eps**2 + eta**2))**2
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC(1)+1,nBas-nR(1)
+    do a=nO(1)+1,nBas-nR(1)
+      do jb=1,nSt
+        eps = e(p,1) - e(a,1) - Omega(jb)
+        Z(p,1) = Z(p,1) + rho(p,a,jb,1)**2*(eps/(eps**2 + eta**2))**2
+      end do
+    end do
+  end do
+
+!----------------!
+! Spin-down part !
+!----------------!
+
+  ! Occupied part of the correlation self-energy
+
+  do p=nC(2)+1,nBas-nR(2)
+    do i=nC(2)+1,nO(2)
+      do jb=1,nSt
+        eps = e(p,2) - e(i,2) + Omega(jb)
+        Z(p,2) = Z(p,2) + rho(p,i,jb,2)**2*(eps/(eps**2 + eta**2))**2
+      end do
+    end do
+  end do
+
+  ! Virtual part of the correlation self-energy
+
+  do p=nC(2)+1,nBas-nR(2)
+    do a=nO(2)+1,nBas-nR(2)
+      do jb=1,nSt
+        eps = e(p,2) - e(a,2) - Omega(jb)
+        Z(p,2) = Z(p,2) + rho(p,a,jb,2)**2*(eps/(eps**2 + eta**2))**2
+      end do
+    end do
+  end do
+
+! Final rescaling
+
+  Z(:,:) = 1d0/(1d0 + Z(:,:))
+
+
+end subroutine unrestricted_renormalization_factor