corrected bug in spin flip

examples/molecule.Be

@@ -1,4 +1,4 @@
 # nAt nEla nElb nCore nRyd
-  1   2     2   0     0
+  1   3     1   0     0
 # Znuc   x            y           z
   Be      0.0          0.0         0.0

input/methods

@@ -7,9 +7,9 @@
 # drCCD rCCD lCCD pCCD
   F      F   F    F
 # CIS* CID CISD
-  F   F   F    
+  T   F   F    
 # RPA* RPAx* ppRPA 
-  F   T    F    
+  F   F    F    
 # G0F2 evGF2 G0F3 evGF3
   F    F     F    F
 # G0W0* evGW* qsGW 

src/CI/UCIS.f90

@@ -1,4 +1,5 @@
-subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int,eHF)
+subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab, &
+                dipole_int_aa,dipole_int_bb,eHF)
 
 ! Perform configuration interaction single calculation`
 
@@ -20,7 +21,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
   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,nspin)
+  double precision,intent(in)   :: dipole_int_bb(nBas,nBas,ncart,nspin)
 
 ! Local variables
 
@@ -103,7 +105,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E
 
     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_aa,nS_bb,nS_sf,lambda,eHF, & 
+    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,ERI_abab,A_sf)
 
     if(dump_matrix) then

src/MBPT/unrestricted_Bethe_Salpeter.f90

@@ -142,7 +142,7 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,
     !-------------------------------------------------
 
     if(dBSE) &
-      call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sf,nS_sc, &
+      call unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,nS_sc, &
                                                             eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb, &
                                                             OmRPA,rho_RPA,OmBSE_sf,XpY_BSE_sf,XmY_BSE_sf)
 

src/QuAcK/QuAcK.f90

@@ -617,7 +617,7 @@ program QuAcK
     if(unrestricted) then
 
       call UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_MO_aaaa,ERI_MO_aabb, & 
-                ERI_MO_bbbb,ERI_MO_abab,dipole_int,eHF)
+                ERI_MO_bbbb,ERI_MO_abab,dipole_int_aa,dipole_int_bb,eHF)
 
    else 
 
@@ -674,7 +674,7 @@ program QuAcK
     if(unrestricted) then
 
        call UdRPA(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
 

src/RPA/URPAx.f90

@@ -108,7 +108,7 @@ subroutine URPAx(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
 
     allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf))
 
-    call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, &
+    call unrestricted_linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, &
                                       ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPAx(ispin),Omega_sf,XpY_sf,XmY_sf)
     call print_excitation('URPAx  ',6,nS_sf,Omega_sf)
     call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &

src/RPA/UdRPA.f90

@@ -106,11 +106,11 @@ subroutine UdRPA(TDA,doACFDT,exchange_kernel,spin_conserved,spin_flip,eta,nBas,n
 
     allocate(Omega_sf(nS_sf),XpY_sf(nS_sf,nS_sf),XmY_sf(nS_sf,nS_sf))
 
-    call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sf,nS_sf,1d0,e, &
+    call unrestricted_linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,nS_sf,1d0,e, &
                                       ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,Omega_sf,rho_sf,EcRPA(ispin),Omega_sf,XpY_sf,XmY_sf)
     call print_excitation('URPA   ',6,nS_sf,Omega_sf)
-!   call print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
+    call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, &
-
+                                               Omega_sf,XpY_sf,XmY_sf)
 
   endif
 

src/RPA/print_excitation.f90

@@ -1,45 +0,0 @@
-subroutine print_excitation(method,ispin,nS,Omega)
-
-! Print excitation energies for a given spin manifold
-
-  implicit none
-  include 'parameters.h'
-
-! Input variables
-
-  character*12,intent(in)            :: method
-  integer,intent(in)                 :: ispin,nS
-  double precision,intent(in)        :: Omega(nS)
-
-! Local variables
-
-  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 = 'alpha-beta'
-  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,A14,A14,A14,A9)') 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) ','|'
-  write(*,*)'-------------------------------------------------------------'
-
-  do ia=1,min(nS,maxS)
-    write(*,'(1X,A1,1X,I5,1X,A1,1X,F23.6,1X,A1,1X,F23.6,1X,A1,1X)') & 
-      '|',ia,'|',Omega(ia),'|',Omega(ia)*HaToeV,'|'
-  enddo
-
-  write(*,*)'-------------------------------------------------------------'
-  write(*,*)
-
-end subroutine print_excitation
-
-

src/RPA/unrestricted_linear_response_A_matrix.f90

@@ -162,7 +162,8 @@ 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,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
-                                - (1d0 - delta_dRPA)*lambda*ERI_abab(b,j,i,a)
+                                - (1d0 - delta_dRPA)*lambda*ERI_abab(j,a,i,b)
+!           print*,'ia,jb,A(ia,jb) = ',i,a,j,b,ia,jb,A_lr(ia,jb)
 
           end  do
         end  do