ppBSE2 beyond TDA

src/GF/GF2_ppBSE2_dynamic_perturbation.f90

@@ -34,15 +34,15 @@ subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,n
 
 ! Local variables
 
-  integer                       :: ab,ij
+  integer                       :: ab,ij,kl
 
   integer                       :: maxOO = 10
   integer                       :: maxVV = 10
 
-  double precision,allocatable  :: Om1Dyn(:)
+  double precision,allocatable  :: Om1_dyn(:)
-  double precision,allocatable  :: Om2Dyn(:)
+  double precision,allocatable  :: Om2_dyn(:)
-  double precision,allocatable  :: Z1Dyn(:)
+  double precision,allocatable  :: Z1_dyn(:)
-  double precision,allocatable  :: Z2Dyn(:)
+  double precision,allocatable  :: Z2_dyn(:)
 
   double precision,allocatable  :: KB_dyn(:,:)
   double precision,allocatable  :: KC_dyn(:,:)
@@ -52,7 +52,7 @@ subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,n
 
 ! Memory allocation
 
-  allocate(Om1Dyn(maxVV),Om2Dyn(maxOO),Z1Dyn(maxVV),Z2Dyn(maxOO), & 
+  allocate(Om1_dyn(maxVV),Om2_dyn(maxOO),Z1_dyn(maxVV),Z2_dyn(maxOO), & 
            KB_dyn(nVV,nOO),KC_dyn(nVV,nVV),KD_dyn(nOO,nOO),   & 
            ZC_dyn(nVV,nVV),ZD_dyn(nOO,nOO))
 
@@ -70,18 +70,34 @@ subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,n
 
   do ab=1,min(nVV,maxVV)
 
-!   if(.not.dTDA) call GF2_ppBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,OmBSE(ab),KB_dyn)
+    if(dTDA) then
-    call GF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn)
 
+      call GF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn)
+  
+      Z1_dyn(ab)  = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab)))
+      Om1_dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab)))
 
-    Z1Dyn(ab)  = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab)))
+    else
-    Om1Dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab)))
 
-    Z1Dyn(ab)  = 1d0/(1d0 - Z1Dyn(ab))
+      call GF2_ppBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn)
-    Om1Dyn(ab) = Z1Dyn(ab)*Om1Dyn(ab)
+      call GF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn)
+      call GF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,Om1(ab),KD_dyn,ZD_dyn)
+
+      Z1_dyn(ab)  = dot_product(X1(:,ab),matmul(ZC_dyn,X1(:,ab))) &
+                  + dot_product(Y1(:,ab),matmul(ZD_dyn,Y1(:,ab)))
+
+      Om1_dyn(ab) = dot_product(X1(:,ab),matmul(KC_dyn - KC_sta,X1(:,ab))) &
+                  - dot_product(Y1(:,ab),matmul(KD_dyn - KD_sta,Y1(:,ab))) &
+                  + dot_product(X1(:,ab),matmul(KB_dyn - KB_sta,Y1(:,ab))) &
+                  - dot_product(Y1(:,ab),matmul(transpose(KB_dyn - KB_sta),X1(:,ab)))
+
+    end if
+  
+    Z1_dyn(ab)  = 1d0/(1d0 - Z1_dyn(ab))
+    Om1_dyn(ab) = Z1_dyn(ab)*Om1_dyn(ab)
 
     write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & 
-      ab,Om1(ab)*HaToeV,(Om1(ab)+Om1Dyn(ab))*HaToeV,Om1Dyn(ab)*HaToeV,Z1Dyn(ab)
+      ab,Om1(ab)*HaToeV,(Om1(ab)+Om1_dyn(ab))*HaToeV,Om1_dyn(ab)*HaToeV,Z1_dyn(ab)
 
   end do
   write(*,*) '---------------------------------------------------------------------------------------------------'
@@ -93,19 +109,38 @@ subroutine GF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,n
   write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)'
   write(*,*) '---------------------------------------------------------------------------------------------------'
 
-  do ij=1,min(nOO,maxOO)
+  kl = 0
+  do ij=max(1,nOO+1-maxOO),nOO
+    kl = kl + 1
 
-!   if(.not.dTDA) call GF2_ppBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,OmBSE(ab),KB_dyn)
+    if(dTDA) then
-    call GF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,Om2(ij),KD_dyn,ZD_dyn)
 
-    Z2Dyn(ij)  = dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))
+      call GF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,-Om2(ij),KD_dyn,ZD_dyn)
-    Om2Dyn(ij) = dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij)))
 
-    Z2Dyn(ij)  = 1d0/(1d0 - Z2Dyn(ij))
+      Z2_dyn(kl)  = dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))
-    Om2Dyn(ij) = Z2Dyn(ij)*Om2Dyn(ij)
+      Om2_dyn(kl) = dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij)))
+
+    else
+
+      call GF2_ppBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn)
+      call GF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,-Om2(ij),KC_dyn,ZC_dyn)
+      call GF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,-Om2(ij),KD_dyn,ZD_dyn)
+
+      Z2_dyn(kl)  = dot_product(X2(:,ij),matmul(ZC_dyn,X2(:,ij))) &
+                  + dot_product(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij)))
+
+      Om2_dyn(kl) = dot_product(X2(:,ij),matmul(KC_dyn - KC_sta,X2(:,ij))) &
+                  - dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij))) &
+                  + dot_product(X2(:,ij),matmul(KB_dyn - KB_sta,Y2(:,ij))) &
+                  - dot_product(Y2(:,ij),matmul(transpose(KB_dyn - KB_sta),X2(:,ij)))
+
+    end if
+
+    Z2_dyn(kl)  = 1d0/(1d0 - Z2_dyn(kl))
+    Om2_dyn(kl) = Z2_dyn(kl)*Om2_dyn(kl)
 
     write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & 
-      ij,Om2(ij)*HaToeV,(Om2(ij)+Om2Dyn(ij))*HaToeV,Om2Dyn(ij)*HaToeV,Z2Dyn(ij)
+      ij,Om2(ij)*HaToeV,(Om2(ij)+Om2_dyn(ij))*HaToeV,Om2_dyn(ij)*HaToeV,Z2_dyn(ij)
 
   end do
   write(*,*) '---------------------------------------------------------------------------------------------------'

src/GW/GW_ppBSE_dynamic_kernel_B.f90

@@ -1,4 +1,4 @@
-subroutine GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGW,Om,rho,OmBSE,KB_dyn)
+subroutine GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGW,Om,rho,KB_dyn)
 
 ! Compute the dynamic part of the Bethe-Salpeter equation matrices
 
@@ -21,7 +21,6 @@ subroutine GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd
   double precision,intent(in)   :: eGW(nBas)
   double precision,intent(in)   :: Om(nS)
   double precision,intent(in)   :: rho(nBas,nBas,nS)
-  double precision,intent(in)   :: OmBSE
   
 ! Local variables
 

src/GW/GW_ppBSE_dynamic_perturbation.f90

@@ -56,15 +56,10 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,
 
 ! Memory allocation
 
-  allocate(Om1_dyn(nOO),Om2_dyn(nOO),Z1_dyn(nOO),Z2_dyn(nOO), & 
+  allocate(Om1_dyn(maxVV),Om2_dyn(maxOO),Z1_dyn(maxVV),Z2_dyn(maxOO), & 
-           KB_dyn(nVV,nOO),KC_dyn(nVV,nVV),KD_dyn(nOO,nOO),   & 
+           KB_dyn(nVV,nOO),KC_dyn(nVV,nVV),KD_dyn(nOO,nOO),           & 
            ZC_dyn(nVV,nVV),ZD_dyn(nOO,nOO))
 
-  Om1_dyn(:) = 0d0
-  Om2_dyn(:) = 0d0
-  Z1_dyn(:) = 0d0
-  Z2_dyn(:) = 0d0
-  
   if(dTDA) then 
     write(*,*)
     write(*,*) '*** dynamical TDA activated ***'
@@ -88,7 +83,7 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,
 
     else 
 
-      call GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KB_dyn)
+      call GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn)
       call GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KC_dyn,ZC_dyn)
       call GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,Om1(ab),KD_dyn,ZD_dyn)
 
@@ -131,7 +126,7 @@ subroutine GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,
 
     else
 
-      call GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,-Om2(ij),KB_dyn)
+      call GW_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGW,OmRPA,rho_RPA,KB_dyn)
       call GW_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,eGW,OmRPA,rho_RPA,-Om2(ij),KC_dyn,ZC_dyn)
       call GW_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,eGW,OmRPA,rho_RPA,-Om2(ij),KD_dyn,ZD_dyn)