correct a few bugs in ppBSE

input/options

@@ -1,5 +1,5 @@
 # HF: maxSCF thresh   DIIS n_diis guess_type ortho_type mix_guess level_shift stability
-      512    0.0000001  T    5     2          1          F         0.0         F
+      512    0.0000001  T    5     2          1          F         1.0         F
 # MP: 
       
 # CC: maxSCF thresh   DIIS n_diis

mol/formaldehyde.xyz

@@ -1,5 +1,5 @@
 4
-Formaldehyde,^1A_1,CC3,aug-cc-pVTZ
+
 C    0.00000000    0.00000000   -0.60298508
 O    0.00000000    0.00000000    0.60539399
 H    0.00000000    0.93467313   -1.18217476

mol/h2.xyz

@@ -1,4 +1,4 @@
 2
 
 H 0. 0. 0.
-H 0. 0. 1.0
+H 0. 0. 0.7

src/GW/Bethe_Salpeter_pp.f90

@@ -94,8 +94,8 @@ subroutine Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,E
     call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, &
                                 Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin))
 
-    call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1)
+!   call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1)
-    call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2)
+!   call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2)
 
     call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2)
 
@@ -128,8 +128,8 @@ subroutine Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,E
     call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, &
                                 Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin))
 
-    call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1)
+!   call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1)
-    call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2)
+!   call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2)
 
     call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2)
 

src/GW/static_screening_WB_pp.f90

@@ -24,6 +24,7 @@ subroutine static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E
 
 ! Local variables
 
+  double precision,external     :: Kronecker_delta
   double precision              :: chi
   double precision              :: eps
   integer                       :: a,b,i,j,ab,ij,m
@@ -57,7 +58,7 @@ subroutine static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E
               chi = chi + rho(a,j,m)*rho(b,i,m)*Omega(m)/eps
             enddo
            
-            WB(ab,ij) = WB(ab,ij) - 4d0*lambda*chi
+            WB(ab,ij) = WB(ab,ij) + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
 
           end do
         end do

src/GW/static_screening_WC_pp.f90

@@ -24,6 +24,7 @@ subroutine static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E
 
 ! Local variables
 
+  double precision,external     :: Kronecker_delta
   double precision              :: chi
   double precision              :: eps
   integer                       :: a,b,c,d,ab,cd,m
@@ -57,8 +58,7 @@ subroutine static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E
               chi = chi + rho(a,d,m)*rho(b,c,m)*Omega(m)/eps
             enddo
             
-            WC(ab,cd) = WC(ab,cd) - 4d0*lambda*chi
+            WC(ab,cd) = WC(ab,cd) + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
-
 
           end do
         end do

src/GW/static_screening_WD_pp.f90

@@ -24,6 +24,7 @@ subroutine static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E
 
 ! Local variables
 
+  double precision,external     :: Kronecker_delta
   double precision              :: chi
   double precision              :: eps
   integer                       :: i,j,k,l,ij,kl,m
@@ -57,7 +58,7 @@ subroutine static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E
               chi = chi + rho(i,l,m)*rho(j,k,m)*Omega(m)/eps
             enddo
  
-            WD(ij,kl) = WD(ij,kl) - 4d0*lambda*chi
+            WD(ij,kl) = WD(ij,kl) + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l)))
 
           end do
         end do

src/LR/print_transition_vectors_pp.f90

@@ -27,24 +27,24 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip
 
   integer                       :: a,b,c,d,ab,cd
   integer                       :: i,j,k,l,ij,kl
-  integer                       :: maxOO
+  integer                       :: maxOO = 10
   integer                       :: maxVV = 10
   double precision              :: S2
   double precision,parameter    :: thres_vec = 0.1d0
-  double precision,allocatable  :: osOO(:)
+  double precision,allocatable  :: os1(:)
-  double precision,allocatable  :: osVV(:)
+  double precision,allocatable  :: os2(:)
 
 ! Memory allocation
 
   maxOO = min(nOO,maxOO)
   maxVV = min(nVV,maxVV)
 
-  allocate(osOO(maxOO),osVV(maxVV))
+  allocate(os1(nVV),os2(nOO))
 
 ! Compute oscillator strengths
 
-  osOO(:) = 0d0
+  os1(:) = 0d0
-  osVV(:) = 0d0
+  os2(:) = 0d0
 
 ! if(spin_allowed) call oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY,XmY,os)
 
@@ -63,8 +63,8 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip
     end if
 
     print*,'-------------------------------------------------------------'
-    write(*,'(A20,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
+    write(*,'(A20,I3,A2,F10.4,A3,A6,F6.4,A11,F6.4)') &
-            ' p-p excitation n. ',ab,': ',Omega1(ab)*HaToeV,' eV','  f = ',osVV(ab),'  <S**2> = ',S2
+            ' p-p excitation n. ',ab,': ',Omega1(ab)*HaToeV,' eV','  f = ',os1(ab),'  <S**2> = ',S2
     print*,'-------------------------------------------------------------'
 
    if(spin_allowed) then
@@ -111,14 +111,14 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip
 
 ! Thomas-Reiche-Kuhn sum rule
 
-  write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for p-p sector = ',sum(osVV(:))
+  write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for p-p sector = ',sum(os1(:))
   write(*,*)
 
 !-----------------------------------------------!
 ! Print details about excitations for hh sector !
 !-----------------------------------------------!
 
-  do ij=1,maxOO
+  do ij=nOO,nOO-maxOO+1,-1
 
     ! <S**2> values
 
@@ -129,8 +129,8 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip
     end if
 
     print*,'-------------------------------------------------------------'
-    write(*,'(A20,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') &
+    write(*,'(A20,I3,A2,F10.4,A3,A6,F6.4,A11,F6.4)') &
-            ' h-h excitation n. ',ij,': ',Omega2(ij)*HaToeV,' eV','  f = ',osOO(ij),'  <S**2> = ',S2
+            ' h-h excitation n. ',ij,': ',Omega2(ij)*HaToeV,' eV','  f = ',os2(ij),'  <S**2> = ',S2
     print*,'-------------------------------------------------------------'
 
    if(spin_allowed) then
@@ -177,7 +177,7 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip
 
 ! Thomas-Reiche-Kuhn sum rule
 
-  write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(osOO(:))
+  write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(os2(:))
   write(*,*)
 
 end subroutine print_transition_vectors_pp

src/RPA/ppRPA.f90

@@ -69,8 +69,8 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo
     call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & 
                             Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA(ispin))
 
-    call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1)
+!   call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1)
-    call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2)
+!   call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2)
 
     call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2)
 
@@ -92,8 +92,8 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo
     call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, &
                             Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA(ispin))
 
-    call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1)
+!   call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1)
-    call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2)
+!   call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2)
 
     call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2)