dipole and f OK

input/methods

@@ -9,11 +9,11 @@
 # CIS* CID CISD
   F   F   F    
 # RPA* RPAx* ppRPA 
-  F   T    F    
+  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

@@ -13,6 +13,6 @@
 # ACFDT: AC Kx  XBS
          F  F   T
 # BSE: BSE dBSE dTDA evDyn
-        F    T    T    F
+        T    F    T    F
 # MCMP2: nMC    nEq    nWalk dt  nPrint iSeed doDrift
          1000000 100000 10    0.3 10000 1234  T

src/QuAcK/QuAcK.f90

@@ -58,7 +58,9 @@ program QuAcK
   double precision,allocatable  :: Hc(:,:)
   double precision,allocatable  :: H(:,:)
   double precision,allocatable  :: X(:,:)
-  double precision,allocatable  :: dipole_int(:,:,:,:)
+  double precision,allocatable  :: dipole_int(:,:,:)
+  double precision,allocatable  :: dipole_int_aa(:,:,:)
+  double precision,allocatable  :: dipole_int_bb(:,:,:)
   double precision,allocatable  :: ERI_AO(:,:,:,:)
   double precision,allocatable  :: ERI_MO(:,:,:,:)
   integer                       :: ixyz
@@ -233,8 +235,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), &
+           S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas))
-           dipole_int(nBas,nBas,ncart,nspin),ERI_AO(nBas,nBas,nBas,nBas))
 
 ! Read integrals
 
@@ -341,11 +342,13 @@ program QuAcK
 
       ! Read and transform dipole-related integrals
     
-      call read_dipole_integrals(nBas,dipole_int)
+      allocate(dipole_int_aa(nBas,nBas,ncart),dipole_int_bb(nBas,nBas,ncart))
+
+      call read_dipole_integrals(nBas,dipole_int_aa)
+      call read_dipole_integrals(nBas,dipole_int_bb)
       do ixyz=1,ncart
-        do ispin=1,nspin
+          call AOtoMO_transform(nBas,cHF(:,:,1),dipole_int_aa(:,:,ixyz))
-          call AOtoMO_transform(nBas,cHF(:,:,ispin),dipole_int(:,:,ixyz,ispin))
+          call AOtoMO_transform(nBas,cHF(:,:,2),dipole_int_bb(:,:,ixyz))
-        end do 
       end do 
 
       ! Memory allocation
@@ -399,10 +402,10 @@ program QuAcK
  
       ! Read and transform dipole-related integrals
     
-      ispin = 1
+      allocate(dipole_int(nBas,nBas,ncart))
       call read_dipole_integrals(nBas,dipole_int)
       do ixyz=1,ncart
-        call AOtoMO_transform(nBas,cHF,dipole_int(:,:,ixyz,ispin))
+        call AOtoMO_transform(nBas,cHF,dipole_int(:,:,ixyz))
       end do 
 
       ! 4-index transform 
@@ -696,7 +699,7 @@ program QuAcK
     if(unrestricted) then
 
        call URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,0d0,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
-                  ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,eHF)
+                  ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF)
 
     else 
 
@@ -822,9 +825,9 @@ program QuAcK
     call cpu_time(start_G0W0)
     if(unrestricted) then 
 
-      call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn, &
+      call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,   & 
-                 spin_conserved,spin_flip,linGW,eta_GW,nBas,nC,nO,nV,nR,nS, & 
+                 linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab, & 
-                 ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,PHF,cHF,eHF,eG0W0)
+                 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, &
@@ -849,9 +852,10 @@ program QuAcK
     call cpu_time(start_evGW)
     if(unrestricted) then 
 
-      call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
+      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,  &
+                G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,    &
-                ERHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int,PHF,cHF,eHF,eG0W0)
+                EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb, & 
+                PHF,cHF,eHF,eG0W0)
 
     else
 

src/QuAcK/UG0W0.f90

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

src/QuAcK/URPAx.f90

@@ -1,5 +1,5 @@
 subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & 
-                 ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,e)
+                 ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e)
 
 ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
 
@@ -28,7 +28,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
   double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_abab(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart,nspin)
+  double precision,intent(in)   :: dipole_int_aa(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int_bb(nBas,nBas,ncart)
 
 ! Local variables
 
@@ -78,7 +79,8 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
     call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & 
                                       ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPAx(ispin),Omega_sc,XpY_sc,XmY_sc)
     call print_excitation('URPAx  ',5,nS_sc,Omega_sc)
-    call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_sc,dipole_int,Omega_sc,XpY_sc,XmY_sc)
+    call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & 
+                                               Omega_sc,XpY_sc,XmY_sc)
 
     deallocate(Omega_sc,XpY_sc,XmY_sc)
 

src/QuAcK/UdRPA.f90

@@ -1,5 +1,5 @@
 subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & 
-                 ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,e)
+                 ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,e)
 
 ! Perform random phase approximation calculation with exchange (aka TDHF) in the unrestricted formalism
 
@@ -28,7 +28,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
   double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_abab(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart,nspin)
+  double precision,intent(in)   :: dipole_int_aa(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int_bb(nBas,nBas,ncart)
 
 ! Local variables
 
@@ -78,7 +79,8 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
     call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0,e, & 
                                       ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sc,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
     call print_excitation('URPA   ',5,nS_sc,Omega_sc)
-!   call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+    call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, &
+                                               Omega_sc,XpY_sc,XmY_sc)
 
 
   endif

src/QuAcK/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,ERI_abab,dipole_int,PHF,cHF,eHF,eG0W0)
+                ERHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,eG0W0)
 
 ! Perform self-consistent eigenvalue-only GW calculation
 
@@ -46,7 +46,8 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE
   double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_abab(nBas,nBas,nBas,nBas)
-  double precision,intent(in)   :: dipole_int(nBas,nBas,ncart,nspin)
+  double precision,intent(in)   :: dipole_int_aa(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int_bb(nBas,nBas,ncart)
 
 ! Local variables
 
@@ -255,7 +256,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,ERI_abab,eGW,eGW,EcBSE)
+                                     ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb,eGW,eGW,EcBSE)
 
 !   if(exchange_kernel) then
 

src/QuAcK/print_transition_vectors.f90

@@ -21,6 +21,7 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int,
 
 ! Local variables
 
+  logical                       :: debug = .false.
   integer                       :: ia,jb,i,j,a,b
   integer                       :: ixyz
   integer,parameter             :: maxS = 10
@@ -54,21 +55,25 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int,
     end do
     f(:,:) = sqrt(2d0)*f(:,:)
  
-    write(*,*) '------------------------'
+    if(debug) then
-    write(*,*) ' Dipole moments (X Y Z) '
+
-    write(*,*) '------------------------'
+      write(*,*) '------------------------'
-    call matout(nS,ncart,f)
+      write(*,*) ' Dipole moments (X Y Z) '
-    write(*,*)
+      write(*,*) '------------------------'
- 
+      call matout(nS,ncart,f)
-    do ia=1,nS
+      write(*,*)
-      os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2)
+  
-    end do
+      do ia=1,nS
-    
+        os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2)
-    write(*,*) '----------------------'
+      end do
-    write(*,*) ' Oscillator strengths '
+      
-    write(*,*) '----------------------'
+      write(*,*) '----------------------'
-    call matout(nS,1,os)
+      write(*,*) ' Oscillator strengths '
-    write(*,*)
+      write(*,*) '----------------------'
+      call matout(nS,1,os)
+      write(*,*)
+
+    end if
 
   end if
   

src/QuAcK/print_unrestricted_transition_vectors.f90

@@ -1,4 +1,5 @@
-subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,nSt,dipole_int,Omega,XpY,XmY)
+subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,dipole_int_aa,dipole_int_bb, & 
+                                                 Omega,XpY,XmY)
 
 ! Print transition vectors for linear response calculation
 
@@ -14,14 +15,18 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
   integer,intent(in)            :: nV(nspin)
   integer,intent(in)            :: nR(nspin)
   integer,intent(in)            :: nS(nspin)
+  integer,intent(in)            :: nSa
+  integer,intent(in)            :: nSb
   integer,intent(in)            :: nSt
-  double precision              :: dipole_int(nBas,nBas,ncart,nspin)
+  double precision              :: dipole_int_aa(nBas,nBas,ncart)
+  double precision              :: dipole_int_bb(nBas,nBas,ncart)
   double precision,intent(in)   :: Omega(nSt)
   double precision,intent(in)   :: XpY(nSt,nSt)
   double precision,intent(in)   :: XmY(nSt,nSt)
 
 ! Local variables
 
+  logical                       :: debug = .false.
   integer                       :: ia,jb,i,j,a,b
   integer                       :: ixyz
   integer                       :: ispin
@@ -43,35 +48,47 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
   f(:,:) = 0d0
   if(spin_allowed) then
 
-    do ispin=1,nspin
+    do ia=1,nSt
-      do ia=1,nSt
+      do ixyz=1,ncart
-        do ixyz=1,ncart
+
-          jb = 0
+        jb = 0
-          do j=nC(ispin)+1,nO(ispin)
+        do j=nC(1)+1,nO(1)
-            do b=nO(ispin)+1,nBas-nR(ispin)
+          do b=nO(1)+1,nBas-nR(1)
-              jb = jb + 1
+            jb = jb + 1
-              f(ia,ixyz) = f(ia,ixyz) + dipole_int(j,b,ixyz,ispin)*XpY(ia,jb)
+            f(ia,ixyz) = f(ia,ixyz) + dipole_int_aa(j,b,ixyz)*XpY(ia,jb)
-            end do
           end do
         end do
+
+        jb = 0
+        do j=nC(2)+1,nO(2)
+          do b=nO(2)+1,nBas-nR(2)
+            jb = jb + 1
+            f(ia,ixyz) = f(ia,ixyz) + dipole_int_bb(j,b,ixyz)*XpY(ia,nSa+jb)
+          end do
+        end do
+
       end do
     end do
  
-    write(*,*) '----------------'
+    if(debug) then
-    write(*,*) ' Dipole moments '
+
-    write(*,*) '----------------'
+      write(*,*) '----------------'
-    call matout(nSt,ncart,f(:,:))
+      write(*,*) ' Dipole moments '
-    write(*,*)
+      write(*,*) '----------------'
- 
+      call matout(nSt,ncart,f(:,:))
-    do ia=1,nSt
+      write(*,*)
-      os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2)
+  
-    end do
+      do ia=1,nSt
-    
+        os(ia) = 2d0/3d0*Omega(ia)*sum(f(ia,:)**2)
-    write(*,*) '----------------------'
+      end do
-    write(*,*) ' Oscillator strengths '
+      
-    write(*,*) '----------------------'
+      write(*,*) '----------------------'
-    call matout(nSt,1,os(:))
+      write(*,*) ' Oscillator strengths '
-    write(*,*)
+      write(*,*) '----------------------'
+      call matout(nSt,1,os(:))
+      write(*,*)
+
+    end if
 
   end if
   
@@ -92,7 +109,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
     do j=nC(1)+1,nO(1)
       do b=nO(1)+1,nBas-nR(1)
         jb = jb + 1
-        if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' -> ',b,' = ',X(jb)/sqrt(2d0)
+        if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb)/sqrt(2d0)
       end do
     end do
  
@@ -100,10 +117,9 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
     do j=nC(1)+1,nO(1)
       do b=nO(1)+1,nBas-nR(1)
         jb = jb + 1
-        if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' <- ',b,' = ',Y(jb)/sqrt(2d0)
+        if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb)/sqrt(2d0)
       end do
     end do
-   write(*,*)
 
     ! Spin-down transitions
 
@@ -111,7 +127,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
     do j=nC(2)+1,nO(2)
       do b=nO(2)+1,nBas-nR(2)
         jb = jb + 1
-        if(abs(X(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' -> ',b,' = ',X(jb)/sqrt(2d0)
+        if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(jb)/sqrt(2d0)
       end do
     end do
  
@@ -119,7 +135,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n
     do j=nC(2)+1,nO(2)
       do b=nO(2)+1,nBas-nR(2)
         jb = jb + 1
-        if(abs(Y(jb)) > thres_vec) write(*,'(I3,A4,I3,A3,F10.6)') j,' <- ',b,' = ',Y(jb)/sqrt(2d0)
+        if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(jb)/sqrt(2d0)
       end do
     end do
    write(*,*)

src/QuAcK/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,ERI_abab, & 
-                                       eW,eGW,EcBSE)
+                                       dipole_int_aa,dipole_int_bb,eW,eGW,EcBSE)
 
 ! Compute the Bethe-Salpeter excitation energies
 
@@ -30,7 +30,8 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
   double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: ERI_abab(nBas,nBas,nBas,nBas)
-
+  double precision,intent(in)   :: dipole_int_aa(nBas,nBas,ncart)
+  double precision,intent(in)   :: dipole_int_bb(nBas,nBas,ncart)
 
 ! Local variables
 
@@ -96,8 +97,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,ERI_abab,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_unrestricted_transition_vectors(.true.,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)
 
     !-------------------------------------------------
     ! Compute the dynamical screening at the BSE level