GW_ppBSE_static_kernel_C for ispin = 1

src/GW/GW_ppBSE.f90

@@ -66,6 +66,11 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,
 
   double precision,intent(out)  :: EcBSE(nspin)
 
+  double precision :: t0, t1
+  double precision :: tt0, tt1
+
+  call wall_time(t0)
+
 !---------------------------------
 ! Compute (singlet) RPA screening 
 !---------------------------------
@@ -76,11 +81,25 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,
   allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), &
            Aph(nS,nS),Bph(nS,nS))
  
+  call wall_time(tt0)
   call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_A =',tt1-tt0,' seconds'
 
+  call wall_time(tt0)
+  if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_B =',tt1-tt0,' seconds'
+
+  call wall_time(tt0)
   call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR =',tt1-tt0,' seconds'
+
+  call wall_time(tt0)
   call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_excitation_density =',tt1-tt0,' seconds'
 
   deallocate(XpY_RPA,XmY_RPA,Aph,Bph)
 
@@ -108,32 +127,63 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,
 
     ! Compute BSE excitation energies
 
-    if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+    call wall_time(tt0)
     call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
-                 call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_C =',tt1-tt0,' seconds'
 
-    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+    call wall_time(tt0)
+    call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_D =',tt1-tt0,' seconds'
+
+    call wall_time(tt0)
+    if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_B =',tt1-tt0,' seconds'
+
+    call wall_time(tt0)
     call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_C =',tt1-tt0,' seconds'
+
+    call wall_time(tt0)
     call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_D =',tt1-tt0,' seconds'
+
+    call wall_time(tt0)
+    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_B =',tt1-tt0,' seconds'
 
     Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
     Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
     Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
 
+    call wall_time(tt0)
     call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR =',tt1-tt0,' seconds'
 
+    call wall_time(tt0)
     call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_transition_vectors =',tt1-tt0,' seconds'
 
     !----------------------------------------------------!
     ! Compute the dynamical screening at the ppBSE level !
     !----------------------------------------------------!
 
+    call wall_time(tt0)
     if(dBSE) &
         call GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
                                            Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
-    write(*,*) "Deallocate not done"
+
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_dynamic_perturbation =',tt1-tt0,' seconds'
+
     deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
-    write(*,*) "Deallocate done"
   end if
 
 !-------------------
@@ -160,33 +210,66 @@ subroutine GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,
 
     ! Compute BSE excitation energies
 
-    if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+    call wall_time(tt0)
     call GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_C =',tt1-tt0,' seconds'
+
+    call wall_time(tt0)
     call GW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_D =',tt1-tt0,' seconds'
 
+    call wall_time(tt0)
+    if(.not.TDA) call GW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_static_kernel_B =',tt1-tt0,' seconds'
 
-    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+    call wall_time(tt0)
     call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_C =',tt1-tt0,' seconds'
+
+    call wall_time(tt0)
     call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_D =',tt1-tt0,' seconds'
+
+    call wall_time(tt0)
+    if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_B =',tt1-tt0,' seconds'
 
     Bpp(:,:) = Bpp(:,:) + KB_sta(:,:)
     Cpp(:,:) = Cpp(:,:) + KC_sta(:,:)
     Dpp(:,:) = Dpp(:,:) + KD_sta(:,:)
 
+    call wall_time(tt0)
     call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin))
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR =',tt1-tt0,' seconds'
 
+    call wall_time(tt0)
     call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppLR_transition_vectors =',tt1-tt0,' seconds'
 
     !----------------------------------------------------!
     ! Compute the dynamical screening at the ppBSE level !
     !----------------------------------------------------!
 
+    call wall_time(tt0)
     if(dBSE) &
         call GW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, &
                                            Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta)
+    call wall_time(tt1)
+    write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE_dynamic_perturbation =',tt1-tt0,' seconds'
 
     deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta)
 
   end if
 
+  call wall_time(t1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_ppBSE =',t1-t0,' seconds'
+
 end subroutine 

src/GW/GW_ppBSE_static_kernel_C.f90

@@ -26,44 +26,94 @@ subroutine GW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI
   double precision,external     :: Kronecker_delta
   double precision              :: chi
   double precision              :: eps
+  double precision              :: tmp_ab, lambda4, eta2
   integer                       :: a,b,c,d,ab,cd,m
+  integer                       :: a0, aa
+
+  double precision, allocatable :: Om_tmp(:)
 
 ! Output variables
 
   double precision,intent(out)  :: KC(nVV,nVV)
 
-! Initialization
-
-  KC(:,:) = 0d0
-
 !---------------!
 ! Singlet block !
 !---------------!
 
   if(ispin == 1) then
 
-    ab = 0
+    a0 = nBas - nR - nO
+    lambda4 = 4.d0 * lambda
+    eta2 = eta * eta
+
+    allocate(Om_tmp(nS))
+
+    !$OMP PARALLEL DEFAULT(NONE) PRIVATE(m) SHARED(nS, eta2, Om, Om_tmp)
+    !$OMP DO
+    do m = 1, nS
+      Om_tmp(m) = Om(m) / (Om(m)*Om(m) + eta2)
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+    !$OMP PARALLEL DEFAULT(NONE)                              &
+    !$OMP          PRIVATE(a, b, aa, ab, c, d, cd, m, tmp_ab) &
+    !$OMP          SHARED(nO, nBas, nR, nS, a0, lambda4, Om_tmp, rho, KC)
+    !$OMP DO
     do a = nO+1, nBas-nR
+      aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO
       do b = a, nBas-nR
-        ab = ab + 1
+        ab = aa + b
+
+        tmp_ab = lambda4
+        if(a .eq. b) then
+          tmp_ab = 0.7071067811865475d0 * lambda4
+        endif
+
         cd = 0
         do c = nO + 1, nBas-nR
           do d = c, nBas-nR
             cd = cd + 1
 
-              chi = 0d0
+            KC(ab,cd) = 0d0
             do m = 1, nS
-                eps = Om(m)**2 + eta**2
+              KC(ab,cd) = KC(ab,cd) - rho(a,c,m) * rho(b,d,m) * Om_tmp(m) &
-                chi = chi - rho(a,c,m)*rho(b,d,m)*Om(m)/eps &
+                                    - rho(a,d,m) * rho(b,c,m) * Om_tmp(m)
-                          - rho(a,d,m)*rho(b,c,m)*Om(m)/eps
             end do
 
-              KC(ab,cd) = 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
+            KC(ab,cd) = tmp_ab * KC(ab,cd)
+            if(c .eq. d) then
+              KC(ab,cd) = 0.7071067811865475d0 * KC(ab,cd)
+            endif
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
 
-          end do
+!    ab = 0
-        end do
+!    do a=nO+1,nBas-nR
-      end do
+!      do b=a,nBas-nR
-    end do
+!        ab = ab + 1
+!        cd = 0
+!        do c=nO+1,nBas-nR
+!          do d=c,nBas-nR
+!            cd = cd + 1
+!
+!              chi = 0d0
+!              do m=1,nS
+!                eps = Om(m)**2 + eta**2
+!                chi = chi - rho(a,c,m)*rho(b,d,m)*Om(m)/eps &
+!                          - rho(a,d,m)*rho(b,c,m)*Om(m)/eps
+!              end do
+!
+!              KC(ab,cd) = 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d)))
+!
+!          end do
+!        end do
+!      end do
+!    end do
 
   end if
 

src/GW/RG0W0.f90

@@ -59,6 +59,11 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
   double precision,allocatable  :: eGWlin(:)
   double precision,allocatable  :: eGW(:)
 
+  double precision :: t0, t1
+  double precision :: tt0, tt1
+
+  call wall_time(t0)
+
 ! Output variables
 
 ! Hello world
@@ -101,26 +106,48 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 ! Compute screening !
 !-------------------!
 
+  call wall_time(tt0)
   call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_A =',tt1-tt0,' seconds'
+
+  call wall_time(tt0)
   if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_B =',tt1-tt0,' seconds'
 
+  call wall_time(tt0)
   call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR =',tt1-tt0,' seconds'
 
+  call wall_time(tt0)
   if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for print_excitation_energies =',tt1-tt0,' seconds'
 
 !--------------------------!
 ! Compute spectral weights !
 !--------------------------!
 
+  call wall_time(tt0)
   call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_excitation_density =',tt1-tt0,' seconds'
 
 !------------------------!
 ! Compute GW self-energy !
 !------------------------!
 
+  call wall_time(tt0)
   if(regularize) call GW_regularization(nBas,nC,nO,nV,nR,nS,eHF,Om,rho)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_regularization =',tt1-tt0,' seconds'
 
+  call wall_time(tt0)
   call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for GW_self_energy_diag =',tt1-tt0,' seconds'
 
 !-----------------------------------!
 ! Solve the quasi-particle equation !
@@ -128,6 +155,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
   ! Linearized or graphical solution?
 
+  call wall_time(tt0)
   eGWlin(:) = eHF(:) + Z(:)*SigC(:)
 
   if(linearize) then 
@@ -145,6 +173,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
     call GW_QP_graph(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,eGWlin,eHF,eGW,Z)
 
   end if
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for QP =',tt1-tt0,' seconds'
 
 ! Plot self-energy, renormalization factor, and spectral function
 
@@ -158,19 +188,33 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
 ! Compute the RPA correlation energy
 
+  call wall_time(tt0)
   call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
-  if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_A =',tt1-tt0,' seconds'
 
+  call wall_time(tt0)
+  if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR_B =',tt1-tt0,' seconds'
+
+  call wall_time(tt0)
   call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phLR =',tt1-tt0,' seconds'
 
 !--------------!
 ! Dump results !
 !--------------!
 
+  call wall_time(tt0)
   call print_RG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM)
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for print_RG0W0 =',tt1-tt0,' seconds'
 
 ! Perform BSE calculation
 
+  call wall_time(tt0)
   if(dophBSE) then
 
     call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
@@ -221,7 +265,10 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
     end if
 
   end if
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for phBSE =',tt1-tt0,' seconds'
 
+  call wall_time(tt0)
   if(doppBSE) then
 
     call GW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE)
@@ -238,6 +285,8 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
     write(*,*)
 
   end if
+  call wall_time(tt1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for ppBSE =',tt1-tt0,' seconds'
 
 ! end if
   
@@ -251,4 +300,7 @@ subroutine RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA
 
   end if
 
+  call wall_time(t1)
+  write(*,'(A65,1X,F9.3,A8)') 'Wall time for RG0W0 =',t1-t0,' seconds'
+
 end subroutine