sf-BSE

input/methods

@@ -7,13 +7,13 @@
 # drCCD rCCD lCCD pCCD
   F      F   F    F
 # CIS* CIS(D) CID CISD
-  T    F      F   F    
+  F    F      F   F    
 # RPA* RPAx* ppRPA 
   F   F    F    
 # G0F2 evGF2 G0F3 evGF3
   F    F     F    F
 # G0W0* evGW* qsGW 
-  F    F    F
+  T    F    F
 # G0T0 evGT qsGT 
   F    F    F
 # MCMP2

input/options

@@ -5,7 +5,7 @@
 # CC: maxSCF thresh   DIIS n_diis
       64     0.0000001  T    5
 # spin: TDA singlet triplet spin_conserved spin_flip 
-         F   T       T       T              T        
+         T   T       T       T              T        
 # 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 
@@ -13,6 +13,6 @@
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
-        T    F    T    F
+        T    T    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

src/CI/UCIS.f90

@@ -38,9 +38,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
 
   integer                       :: nS_ab,nS_ba,nS_sf
   double precision,allocatable  :: A_sf(:,:)
-  double precision,allocatable  :: Z_sf(:,:)
   double precision,allocatable  :: Omega_sf(:)
-  integer         ,allocatable  :: order(:)
 
 ! Hello world
 
@@ -80,9 +78,10 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
     endif
 
     call diagonalize_matrix(nS_sc,A_sc,Omega_sc)
+    A_sc(:,:) = transpose(A_sc)
     call print_excitation('UCIS        ',5,nS_sc,Omega_sc)
     call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
-                                               cHF,S,Omega_sc,transpose(A_sc),transpose(A_sc))
+                                               cHF,S,Omega_sc,A_sc,A_sc)
  
     if(dump_trans) then
       print*,'Spin-conserved CIS transition vectors'
@@ -108,7 +107,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
     nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1))
     nS_sf = nS_ab + nS_ba
 
-    allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf),Z_sf(nS_sf,nS_sf))
+    allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf))
     
     call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,lambda,eHF, & 
                                                ERI_aaaa,ERI_aabb,ERI_bbbb,A_sf)
@@ -120,13 +119,10 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
     endif
 
     call diagonalize_matrix(nS_sf,A_sf,Omega_sf)
-!   allocate(order(nS_sf))
+    A_sf(:,:) = transpose(A_sf)
-!   call diagonalize_general_matrix(nS_sf,A_sf,Omega_sf,Z_sf)
-!   call quick_sort(Omega_sf,order(:),nS_sf)
-
     call print_excitation('UCIS        ',6,nS_sf,Omega_sf)
     call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
-                                               cHF,S,Omega_sf,transpose(A_sf),transpose(A_sf))
+                                               cHF,S,Omega_sf,A_sf,A_sf)
  
     if(dump_trans) then
       print*,'Spin-flip CIS transition vectors'

src/LR/linear_response.f90

@@ -19,7 +19,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
   
 ! Local variables
 
-  integer                       :: ia
   double precision              :: trace_matrix
   double precision,allocatable  :: A(:,:)
   double precision,allocatable  :: B(:,:)
@@ -48,21 +47,26 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
 
 ! Tamm-Dancoff approximation
 
-  B = 0d0
+  if(TDA) then
-  if(.not. TDA) then
+ 
+    B(:,:)   = 0d0
+    XpY(:,:) = A(:,:)
+    XmY(:,:) = 0d0
+    call diagonalize_matrix(nS,XpY,Omega)
+    XpY(:,:) = transpose(XpY(:,:))
+
+  else
 
     call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
 
     if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,B)
 
-  end if
+    ! Build A + B and A - B matrices 
-
-! Build A + B and A - B matrices 
 
     ApB = A + B
     AmB = A - B
 
-! Diagonalize linear response matrix
+  ! Diagonalize linear response matrix
 
     call diagonalize_matrix(nS,AmB,Omega)
 
@@ -83,9 +87,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
     if(minval(Omega) < 0d0) & 
       call print_warning('You may have instabilities in linear response: negative excitations!!')
  
-! do ia=1,nS
+  ! do ia=1,nS
-!   if(Omega(ia) < 0d0) Omega(ia) = 0d0
+  !   if(Omega(ia) < 0d0) Omega(ia) = 0d0
-! end do
+  ! end do
 
     Omega = sqrt(Omega)
 
@@ -95,7 +99,9 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E
     XmY = matmul(transpose(Z),AmBIv)
     call DA(nS,1d0*sqrt(Omega),XmY)
 
-! Compute the RPA correlation energy
+  end if
+
+  ! Compute the RPA correlation energy
 
     EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A))
 

src/LR/print_excitation.f90

@@ -14,7 +14,7 @@ subroutine print_excitation(method,ispin,nS,Omega)
 ! Local variables
 
   character*14                       :: spin_manifold
-  integer,parameter                  :: maxS = 10
+  integer,parameter                  :: maxS = 20
   integer                            :: ia
 
   if(ispin == 1) spin_manifold = 'singlet'

src/LR/print_unrestricted_transition_vectors.f90

@@ -29,7 +29,7 @@ subroutine print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,n
 ! Local variables
 
   integer                       :: ia,jb,j,b
-  integer                       :: maxS = 10
+  integer                       :: maxS = 20
   double precision,parameter    :: thres_vec = 0.1d0
   double precision,allocatable  :: X(:)
   double precision,allocatable  :: Y(:)

src/LR/unrestricted_S2_expval.f90

@@ -62,7 +62,7 @@ subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S
 !-------------------------!
 
   S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0)
-  S2_gs    = S2_exact + nO(2) - sum(OO(:,:)**2)
+  S2_gs    = S2_exact + dble(nO(2)) - sum(OO(:,:)**2)
 
 !------------------------------------------!
 ! <S**2> for spin-conserved-excited states !
@@ -152,7 +152,7 @@ subroutine unrestricted_S2_expval(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,maxS,c,S
       Yat(:,:) = transpose(Ya(:,:))
       Ybt(:,:) = transpose(Yb(:,:))
 
-      S2(m) = S2_gs & 
+      S2(m) = S2_gs + dble(nO(2) - nO(1)) + 1d0 &
 
             + trace_matrix(nV(1),matmul(Xbt,matmul(OOt,matmul(OO,Xb)))) &
             - trace_matrix(nO(2),matmul(Xb,matmul(VO,matmul(VOt,Xbt)))) &

src/LR/unrestricted_linear_response.f90

@@ -33,7 +33,6 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
   
 ! Local variables
 
-  integer                       :: ia
   double precision,external     :: trace_matrix
   double precision,allocatable  :: A(:,:)
   double precision,allocatable  :: B(:,:)
@@ -65,8 +64,15 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
 
 ! Tamm-Dancoff approximation
 
-  B = 0d0
+  if(TDA) then
-  if(.not. TDA) then
+
+    B(:,:)   = 0d0
+    XpY(:,:) = A(:,:)
+    XmY(:,:) = 0d0
+    call diagonalize_matrix(nSt,XpY,Omega)
+    XpY(:,:) = transpose(XpY(:,:))
+
+  else
 
     call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
                                                ERI_aaaa,ERI_aabb,ERI_bbbb,B)
@@ -75,23 +81,21 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
       call unrestricted_Bethe_Salpeter_B_matrix(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda, & 
                                                 ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,B)
 
-  end if
+  ! Build A + B and A - B matrices 
-
-! Build A + B and A - B matrices 
 
     ApB = A + B
     AmB = A - B
 
-! Diagonalize linear response matrix
+  ! Diagonalize linear response matrix
 
    call diagonalize_matrix(nSt,AmB,Omega)
 
     if(minval(Omega) < 0d0) &
       call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
 
-! do ia=1,nSt
+  ! do ia=1,nSt
-!   if(Omega(ia) < 0d0) Omega(ia) = 0d0
+  !   if(Omega(ia) < 0d0) Omega(ia) = 0d0
-! end do
+  ! end do
 
     call ADAt(nSt,AmB,1d0*sqrt(Omega),AmBSq)
     call ADAt(nSt,AmB,1d0/sqrt(Omega),AmBIv)
@@ -103,9 +107,9 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
     if(minval(Omega) < 0d0) & 
       call print_warning('You may have instabilities in linear response: negative excitations!!')
  
-! do ia=1,nSt
+  ! do ia=1,nSt
-!   if(Omega(ia) < 0d0) Omega(ia) = 0d0
+  !   if(Omega(ia) < 0d0) Omega(ia) = 0d0
-! end do
+  ! end do
 
     Omega = sqrt(Omega)
  
@@ -115,6 +119,8 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
     XmY = matmul(transpose(Z),AmBIv)
     call DA(nSt,1d0*sqrt(Omega),XmY)
 
+  end if
+
 ! Compute the RPA correlation energy
 
   EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nSt,A))

src/MBPT/G0W0.f90

@@ -1,6 +1,6 @@
 subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & 
                 dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,            & 
-                nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI,dipole_int,PHF,cHF,eHF,eGW)
+                nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF,eGW)
 
 ! Perform G0W0 calculation
 
@@ -29,12 +29,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: ERHF
-  double precision,intent(in)   :: eHF(nBas)
-  double precision,intent(in)   :: cHF(nBas,nBas)
-  double precision,intent(in)   :: PHF(nBas,nBas)
-  double precision,intent(in)   :: Hc(nBas,nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
+  double precision,intent(in)   :: eHF(nBas)
 
 ! Local variables
 

src/MBPT/UG0W0.f90

@@ -1,5 +1,5 @@
 subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, &
-                 linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,          & 
+                 linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,          & 
                  dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eGW)
 
 ! Perform unrestricted G0W0 calculation
@@ -33,15 +33,15 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
   integer,intent(in)            :: nS(nspin)
   double precision,intent(in)   :: ENuc
   double precision,intent(in)   :: EUHF
-  double precision,intent(in)   :: eHF(nBas,nspin)
+  double precision,intent(in)   :: S(nBas,nBas)
-  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_aaaa(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: dipole_int_aa(nBas,nBas,ncart)
   double precision,intent(in)   :: dipole_int_bb(nBas,nBas,ncart)
+  double precision,intent(in)   :: eHF(nBas,nspin)
+  double precision,intent(in)   :: cHF(nBas,nBas,nspin)
+  double precision,intent(in)   :: PHF(nBas,nBas,nspin)
 
 ! Local variables
 
@@ -180,8 +180,8 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
 
   if(BSE) then
 
-    call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
+    call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, &
-                                     ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF,eGW,EcBSE)
+                                     ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE)
 
 !   if(exchange_kernel) then
 !

src/MBPT/evGW.f90

@@ -1,6 +1,6 @@
 subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,  & 
-                G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI, & 
+                G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI, & 
-                dipole_int,PHF,cHF,eHF,eG0W0)
+                dipole_int,eHF,eG0W0)
 
 ! Perform self-consistent eigenvalue-only GW calculation
 
@@ -32,10 +32,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
   double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
   double precision,intent(in)   :: eHF(nBas)
-  double precision,intent(in)   :: cHF(nBas,nBas)
-  double precision,intent(in)   :: PHF(nBas,nBas) 
   double precision,intent(in)   :: eG0W0(nBas)
-  double precision,intent(in)   :: Hc(nBas,nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: dipole_int(nBas,nBas,ncart)
 

src/MBPT/evUGW.f90

@@ -1,6 +1,6 @@
 subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,    & 
                 G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, &
-                ERHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0)
+                ERHF,S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eG0W0)
 
 ! Perform self-consistent eigenvalue-only GW calculation
 
@@ -39,9 +39,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
 
   double precision,intent(in)   :: eHF(nBas,nspin)
   double precision,intent(in)   :: cHF(nBas,nBas,nspin)
-  double precision,intent(in)   :: PHF(nBas,nBas,nspin) 
   double precision,intent(in)   :: eG0W0(nBas,nspin)
-  double precision,intent(in)   :: Hc(nBas,nBas)
+  double precision,intent(in)   :: S(nBas,nBas)
   double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
@@ -255,7 +254,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
   if(BSE) then
 
     call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, &
-                                     ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eGW,eGW,EcBSE)
+                                     S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE)
 
 !   if(exchange_kernel) then
 

src/MBPT/print_qsGW.f90

@@ -1,4 +1,4 @@
-subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigC,Z,EcRPA,EcGM,EqsGW)
+subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigC,Z,EcRPA,EcGM,EqsGW)
 
 
 ! Print one-electron energies and other stuff for qsGW
@@ -11,7 +11,7 @@ subroutine print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,Sig
   integer,intent(in)                 :: nBas,nO,nSCF
   double precision,intent(in)        :: ENuc,EcRPA,EcGM,Conv,thresh
   double precision,intent(in)        :: eHF(nBas),eGW(nBas),c(nBas),P(nBas,nBas) 
-  double precision,intent(in)        :: T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas)
+  double precision,intent(in)        :: T(nBas,nBas),V(nBas,nBas)
   double precision,intent(in)        :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas)
   double precision,intent(in)        :: Z(nBas),SigC(nBas,nBas)
 

src/MBPT/qsGW.f90

@@ -223,7 +223,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
     ! Print results
 
 !   call print_excitation('RPA   ',ispin,nS,Omega(:,ispin))
-    call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA,EcGM,EqsGW)
+    call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,J,K,F,SigCp,Z,EcRPA,EcGM,EqsGW)
 
     ! Increment
 

src/MBPT/unrestricted_Bethe_Salpeter.f90

@@ -1,6 +1,6 @@
 subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,  & 
-                                       nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb, & 
+                                       nBas,nC,nO,nV,nR,nS,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & 
-                                       dipole_int_aa,dipole_int_bb,eW,eGW,EcBSE)
+                                       dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE)
 
 ! Compute the Bethe-Salpeter excitation energies
 
@@ -24,6 +24,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
   integer,intent(in)            :: nV(nspin)
   integer,intent(in)            :: nR(nspin)
   integer,intent(in)            :: nS(nspin)
+  double precision,intent(in)   :: S(nBas,nBas)
+  double precision,intent(in)   :: cW(nBas,nBas,nspin)
   double precision,intent(in)   :: eW(nBas,nspin)
   double precision,intent(in)   :: eGW(nBas,nspin)
   double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
@@ -96,9 +98,9 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
     call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, &
                                       eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin),       & 
                                       OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
-    call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc)
+    call print_excitation('BSE@UGW  ',5,nS_sc,OmBSE_sc)
-!   call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
+    call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
-!                                              OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
+                                               cW,S,OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
 
     !-------------------------------------------------
     ! Compute the dynamical screening at the BSE level
@@ -136,7 +138,9 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
                                       eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcBSE(ispin),       & 
                                       OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
 
-    call print_excitation('BSE@UG0W0',6,nS_sf,OmBSE_sf)
+    call print_excitation('BSE@UGW  ',6,nS_sf,OmBSE_sf)
+    call print_unrestricted_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
+                                               cW,S,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
 
     !-------------------------------------------------
     ! Compute the dynamical screening at the BSE level

src/QuAcK/QuAcK.f90

@@ -56,7 +56,6 @@ program QuAcK
   double precision,allocatable  :: T(:,:)
   double precision,allocatable  :: V(:,:)
   double precision,allocatable  :: Hc(:,:)
-  double precision,allocatable  :: H(:,:)
   double precision,allocatable  :: X(:,:)
   double precision,allocatable  :: dipole_int(:,:,:)
   double precision,allocatable  :: dipole_int_aa(:,:,:)
@@ -64,7 +63,6 @@ program QuAcK
   double precision,allocatable  :: ERI_AO(:,:,:,:)
   double precision,allocatable  :: ERI_MO(:,:,:,:)
   integer                       :: ixyz
-  integer                       :: ispin
   integer                       :: bra1,bra2
   integer                       :: ket1,ket2
   double precision,allocatable  :: ERI_MO_aaaa(:,:,:,:)
@@ -234,7 +232,7 @@ program QuAcK
 ! Memory allocation for one- and two-electron integrals
 
   allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), &
-           S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), &
+           S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas), &
            dipole_int(nBas,nBas,ncart))
 
 ! Read integrals
@@ -322,10 +320,6 @@ program QuAcK
 ! AO to MO integral transform for post-HF methods
 !------------------------------------------------------------------------
 
-! Compute Hartree Hamiltonian in the MO basis
-
-  call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO,H)
-
   call cpu_time(start_AOtoMO)
 
   write(*,*)
@@ -811,12 +805,12 @@ program QuAcK
     if(unrestricted) then 
 
       call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,   & 
-                 linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & 
+                 linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & 
                  dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0)
     else
 
       call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, &
-                linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,dipole_int,PHF,cHF,eHF,eG0W0)
+                linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF,eG0W0)
 
     end if
 
@@ -839,14 +833,14 @@ program QuAcK
 
       call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,   &
                 G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,    &
-                EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & 
+                EUHF,S,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & 
-                PHF,cHF,eHF,eG0W0)
+                cHF,eHF,eG0W0)
 
     else
 
       call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,       &
                 BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, &
-                nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,dipole_int,PHF,cHF,eHF,eG0W0)
+                nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int,eHF,eG0W0)
     end if
     call cpu_time(end_evGW)
 
@@ -1016,7 +1010,7 @@ program QuAcK
     call cpu_time(start_G0W0)
     call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,       & 
               dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & 
-              nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,dipole_int,PHF,cHF,eHF,eG0W0)
+              nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int,eHF,eG0W0)
     call cpu_time(end_G0W0)
   
     t_G0W0 = end_G0W0 - start_G0W0

src/QuAcK/read_options.f90

@@ -71,8 +71,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
   
 ! Local variables
 
-  character(len=1)              :: answer1,answer2,answer3,answer4,answer5, & 
+  character(len=1)              :: answer1,answer2,answer3,answer4,answer5,answer6,answer7
-                                   answer6,answer7,answer8,answer9,answer10
 
 ! Open file with method specification