BSE and spin

input/options

@@ -4,8 +4,8 @@
       
 # CC: maxSCF thresh   DIIS n_diis
       64     0.0000001  T    5
-# spin: singlet triplet spin_conserved spinf_flip TDA
+# spin: singlet triplet spin_conserved spin_flip TDA
-        T       T       T              F          F
+        T       T       T              T          F
 # 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 

src/QuAcK/QuAcK.f90

@@ -775,7 +775,7 @@ program QuAcK
 
       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,PHF,cHF,eHF,eG0W0)
+                 ENuc,EUHF,Hc,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,ERI_MO_abab,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,  &
                  spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS, &
-                 ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,PHF,cHF,eHF,eGW)
+                 ENuc,EUHF,Hc,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,PHF,cHF,eHF,eGW)
 
 ! Perform unrestricted G0W0 calculation
 
@@ -40,6 +40,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
   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)   :: ERI_abab(nBas,nBas,nBas,nBas)
 
 ! Local variables
 
@@ -111,7 +112,7 @@ 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,nS_aa,nS_bb,nS_sc,1d0, &
-                                    eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
+                                    eHF,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
 
   if(print_W) call print_excitation('RPA@UHF',5,nS_sc,Omega_sc)
 
@@ -164,7 +165,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev
 ! Compute the RPA correlation energy
 
   call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, &
-                                    eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
+                                    eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_sc,EcRPA(ispin),Omega_sc,XpY_sc,XmY_sc)
 
   write(*,*)
   write(*,*)'-------------------------------------------------------------------------------'
@@ -182,7 +183,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,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGW,EcRPA,EcBSE)
+                                     nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & 
+                                     eHF,eGW,EcRPA,EcBSE)
 
 !   if(exchange_kernel) then
 !

src/QuAcK/unrestricted_Bethe_Salpeter.f90

@@ -1,5 +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,eW,eGW,EcRPA,EcBSE)
+                                       nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, & 
+                                       eW,eGW,EcRPA,EcBSE)
 
 ! Compute the Bethe-Salpeter excitation energies
 
@@ -28,6 +29,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
   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)   :: ERI_abab(nBas,nBas,nBas,nBas)
 
 
 ! Local variables
@@ -71,16 +73,19 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
     ! Compute spin-conserved RPA screening 
 
     call unrestricted_linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, &
-                                      eW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA_sc,EcRPA(ispin),OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc)
+                                      eW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcRPA(ispin),        & 
+                                      OmRPA_sc,XpY_RPA_sc,XmY_RPA_sc)
 
-    call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA_sc,rho_RPA_sc)
+    call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb, & 
+                                         XpY_RPA_sc,rho_RPA_sc)
 
     ! Compute spin-conserved BSE excitation energies
 
     OmBSE_sc(:) = OmRPA_sc(:)
 
     call unrestricted_linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,1d0, &
-                                      eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,rho_RPA_sc,EcBSE(ispin),OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
+                                      eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho_RPA_sc,EcBSE(ispin),    & 
+                                      OmBSE_sc,XpY_BSE_sc,XmY_BSE_sc)
 
     call print_excitation('BSE@UG0W0',5,nS_sc,OmBSE_sc)
 

src/QuAcK/unrestricted_Bethe_Salpeter_A_matrix.f90

@@ -1,4 +1,6 @@
-subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,A_lr)
+subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
+                                                ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,     & 
+                                                Omega,rho,A_lr)
 
 ! Compute the extra term for Bethe-Salpeter equation for linear response in the unrestricted formalism
 
@@ -20,6 +22,7 @@ subroutine unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
   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)   :: ERI_abab(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: Omega(nSt)
   double precision,intent(in)   :: rho(nBas,nBas,nSt,nspin)
   

src/QuAcK/unrestricted_Bethe_Salpeter_B_matrix.f90

@@ -1,4 +1,6 @@
-subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega,rho,B_lr)
+subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
+                                                ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,     & 
+                                                Omega,rho,B_lr)
 
 ! Compute the extra term for Bethe-Salpeter equation for linear response 
 
@@ -20,6 +22,7 @@ subroutine unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt
   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)   :: ERI_abab(nBas,nBas,nBas,nBas)
   double precision,intent(in)   :: Omega(nSt)
   double precision,intent(in)   :: rho(nBas,nBas,nSt,nspin)
   

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_aaaa,ERI_aabb,ERI_bbbb,rho,EcRPA,Omega,XpY,XmY)
+                                        e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,rho,EcRPA,Omega,XpY,XmY)
 
 ! Compute linear response for unrestricted formalism
 
@@ -27,6 +27,7 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,
   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)   :: ERI_abab(nBas,nBas,nBas,nBas)
   
 ! Local variables
 
@@ -53,20 +54,24 @@ 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_aaaa,ERI_aabb,ERI_bbbb,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,ERI_abab,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)
+    call unrestricted_Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
+                                              ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,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_aaaa,ERI_aabb,ERI_bbbb,B)
+    call unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
+                                               ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,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)
+      call unrestricted_Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,lambda, & 
+                                                ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,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_aaaa,ERI_aabb,ERI_bbbb,A_lr)
+                                                 e,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_lr)
 
 ! Compute linear response
 
@@ -23,6 +23,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
   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)   :: ERI_abab(nBas,nBas,nBas,nBas) 
   
 ! Local variables
 
@@ -46,7 +47,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
 
   if(ispin == 1) then 
 
-    ! alpha-alpha block
+    ! aaaa block
 
     ia = 0
     do i=nC(1)+1,nO(1)
@@ -65,7 +66,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
       end  do
     end  do
 
-    ! alpha-beta block
+    ! aabb block
 
     ia = 0
     do i=nC(1)+1,nO(1)
@@ -83,7 +84,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
       end  do
     end  do
 
-    ! beta-alpha block
+    ! bbaa block
 
     ia = 0
     do i=nC(2)+1,nO(2)
@@ -101,7 +102,7 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
       end  do
     end  do
 
-    ! beta-beta block
+    ! bbbb block
 
     ia = 0
     do i=nC(2)+1,nO(2)
@@ -128,7 +129,45 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
 
   if(ispin == 2) then 
 
-    print*,'spin-flip transition NYI'
+    A_lr(:,:) = 0d0
+
+    ! abab block
+
+    ia = 0
+    do i=nC(1)+1,nO(1)
+      do a=nO(2)+1,nBas-nR(2)
+        ia = ia + 1
+        jb = 0
+        do j=nC(1)+1,nO(1)
+          do b=nO(2)+1,nBas-nR(2)
+            jb = jb + 1
+ 
+            A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                        + lambda*ERI_abab(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a)
+
+          end  do
+        end  do
+      end  do
+    end  do
+
+    ! baba block
+
+    ia = 0
+    do i=nC(2)+1,nO(2)
+      do a=nO(1)+1,nBas-nR(1)
+        ia = ia + 1
+        jb = 0
+        do j=nC(2)+1,nO(2)
+          do b=nO(1)+1,nBas-nR(1)
+            jb = jb + 1
+ 
+            A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                                + lambda*ERI_abab(b,i,j,a) - (1d0 - delta_dRPA)*lambda*ERI_abab(b,i,a,j)
+
+          end  do
+        end  do
+      end  do
+    end  do
 
   end if
 

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_aaaa,ERI_aabb,ERI_bbbb,B_lr)
+                                                 ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,B_lr)
 
 ! Compute linear response
 
@@ -22,6 +22,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
   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)   :: ERI_abab(nBas,nBas,nBas,nBas) 
   
 ! Local variables
 
@@ -40,12 +41,12 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
   if(dRPA) delta_dRPA = 1d0
 
 !-----------------------------------------------
-! Build A matrix for spin-conserving transitions
+! Build B matrix for spin-conserving transitions
 !-----------------------------------------------
 
   if(ispin == 1) then 
 
-    ! alpha-alpha block
+    ! aaaa block
 
     ia = 0
     do i=nC(1)+1,nO(1)
@@ -63,7 +64,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
       end  do
     end  do
 
-    ! alpha-beta block
+    ! aabb block
 
     ia = 0
     do i=nC(1)+1,nO(1)
@@ -81,7 +82,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
       end  do
     end  do
 
-    ! beta-alpha block
+    ! bbaa block
 
     ia = 0
     do i=nC(2)+1,nO(2)
@@ -99,7 +100,7 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
       end  do
     end  do
 
-    ! beta-beta block
+    ! bbbb block
 
     ia = 0
     do i=nC(2)+1,nO(2)
@@ -120,12 +121,48 @@ subroutine unrestricted_linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa
   end if
 
 !-----------------------------------------------
-! Build A matrix for spin-flip transitions
+! Build B matrix for spin-flip transitions
 !-----------------------------------------------
 
   if(ispin == 2) then
 
-    print*,'spin-flip transition NYI'
+    B_lr(:,:) = 0d0
+
+    ! abab block
+
+    ia = 0
+    do i=nC(1)+1,nO(1)
+      do a=nO(2)+1,nBas-nR(2)
+        ia = ia + 1
+        jb = 0
+        do j=nC(2)+1,nO(2)
+          do b=nO(2)+1,nBas-nR(2)
+            jb = jb + 1
+ 
+            B_lr(ia,jb) = lambda*ERI_abab(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI_abab(i,j,b,a)
+
+          end  do
+        end  do
+      end  do
+    end  do
+
+    ! bbbb block
+
+    ia = 0
+    do i=nC(2)+1,nO(2)
+      do a=nO(1)+1,nBas-nR(1)
+        ia = ia + 1
+        jb = 0
+        do j=nC(2)+1,nO(2)
+          do b=nO(1)+1,nBas-nR(1)
+            jb = jb + 1
+ 
+            B_lr(nSa+ia,nSa+jb) = lambda*ERI_abab(j,i,b,a) - (1d0 - delta_dRPA)*lambda*ERI_abab(j,i,a,b)
+
+          end  do
+        end  do
+      end  do
+    end  do
   
   end if