include B term in phBSE2@GW

2024-12-23 04:43:42 +01:00 · 2023-07-21 11:18:00 +02:00 · 2023-07-21 11:18:00 +02:00 · 6cc9802d00
commit 6cc9802d00
parent d7e059c330
8 changed files with 129 additions and 15 deletions
--- a/input/options
+++ b/input/options
@ -5,7 +5,7 @@
 # CC: maxSCF thresh   DIIS n_diis
      64     0.0000001  T    5
 # spin: TDA singlet triplet spin_conserved spin_flip 
-        F   T       T       T              T 
+        F   T       F       T              T 
 # GF: maxSCF thresh  DIIS n_diis lin eta renorm reg
      256    0.00001 T    5      T   0.0 0      F
 # GW: maxSCF thresh  DIIS n_diis lin eta COHSEX TDA_W reg
@ -15,4 +15,4 @@
 # ACFDT: AC Kx  XBS
         F  F   T
 # BSE: phBSE phBSE2 ppBSE dBSE dTDA
-       T     F      F     T    T     
+       T     T      F     T    T     
--- a/src/GW/GW_excitation_density.f90
+++ b/src/GW/GW_excitation_density.f90
@ -6,7 +6,11 @@ subroutine GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)

 ! Input variables

-  integer,intent(in)            :: nBas,nC,nO,nR,nS
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
  double precision,intent(in)   :: XpY(nS,nS)

--- a/src/GW/GW_phBSE2_dynamic_kernel_A.f90
+++ b/src/GW/GW_phBSE2_dynamic_kernel_A.f90
@ -1,4 +1,4 @@
-subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,ZA_dyn)
+subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,KA_dyn,ZA_dyn)

 ! Compute the dynamic part of the Bethe-Salpeter equation matrices

@ -7,7 +7,12 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,

 ! Input variables

-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
  double precision,intent(in)   :: eta
  double precision,intent(in)   :: W(nBas,nBas,nBas,nBas)
  double precision,intent(in)   :: eGW(nBas)
@ -21,14 +26,14 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,

 ! Output variables

-  double precision,intent(out)  :: A_dyn(nS,nS)
+  double precision,intent(out)  :: KA_dyn(nS,nS)
  double precision,intent(out)  :: ZA_dyn(nS,nS)

 ! Build dynamic A matrix

  jb = 0

-!$omp parallel do default(private) shared(A_dyn,ZA_dyn,OmBSE,W,num,dem,eGW,nO,nBas,eta,nC,nR)
+!$omp parallel do default(private) shared(KA_dyn,ZA_dyn,OmBSE,W,num,dem,eGW,nO,nBas,eta,nC,nR)
  do j=nC+1,nO
    do b=nO+1,nBas-nR
      jb = (b-nO) + (j-1)*(nBas-nO)
@ -44,13 +49,13 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,
              dem =  OmBSE - eGW(a) - eGW(c) + eGW(k) + eGW(j)
              num = 2d0*W(j,k,i,c)*W(a,c,b,k)

-              A_dyn(ia,jb) = A_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
+              KA_dyn(ia,jb) = KA_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
              ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2

              dem = OmBSE + eGW(i) - eGW(c) + eGW(k) - eGW(b)
              num = 2d0*W(j,c,i,k)*W(a,k,b,c)

-              A_dyn(ia,jb) = A_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
+              KA_dyn(ia,jb) = KA_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
              ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2

            enddo
@ -62,7 +67,7 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,
              dem = OmBSE + eGW(i) + eGW(j) - eGW(c) - eGW(d)
              num = 2d0*W(a,j,c,d)*W(c,d,i,b)

-              A_dyn(ia,jb) =  A_dyn(ia,jb) + num*dem/(dem**2 + eta**2)
+              KA_dyn(ia,jb) = KA_dyn(ia,jb) + num*dem/(dem**2 + eta**2)
              ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2

            enddo
@ -74,7 +79,7 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,
              dem = OmBSE - eGW(a) - eGW(b) + eGW(k) + eGW(l)
              num = 2d0*W(a,j,k,l)*W(k,l,i,b)

-              A_dyn(ia,jb) =  A_dyn(ia,jb) + num*dem/(dem**2 + eta**2)
+              KA_dyn(ia,jb) = KA_dyn(ia,jb) + num*dem/(dem**2 + eta**2)
              ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - num*(dem**2 - eta**2)/(dem**2 + eta**2)**2

            end do
--- a/src/GW/GW_phBSE2_dynamic_kernel_B.f90
+++ b/src/GW/GW_phBSE2_dynamic_kernel_B.f90
@ -0,0 +1,88 @@
+subroutine GW_phBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,eGW,W,KB_dyn)
+
+! Compute the dynamic part of the Bethe-Salpeter equation matrices
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: W(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: eGW(nBas)
+  
+! Local variables
+
+  double precision              :: num,dem
+  double precision              :: eps
+  integer                       :: i,j,a,b,ia,jb,kc,k,c,l,d
+
+! Output variables
+
+  double precision,intent(out)  :: KB_dyn(nS,nS)
+
+! Build dynamic B matrix
+
+  jb = 0
+
+!$omp parallel do default(private) shared(KB_dyn,W,num,dem,eGW,nO,nBas,eta,nC,nR)
+  do j=nC+1,nO
+    do b=nO+1,nBas-nR
+      jb = (b-nO) + (j-1)*(nBas-nO)
+
+      ia = 0
+      do i=nC+1,nO
+        do a=nO+1,nBas-nR
+          ia = (a-nO) + (i-1)*(nBas-nO)
+
+          do k=nC+1,nO
+            do c=nO+1,nBas-nR
+
+              dem =        - eGW(a) - eGW(c) + eGW(k) + eGW(j)
+              num = 2d0*W(b,k,i,c)*W(a,c,j,k)
+
+              KB_dyn(ia,jb) = KB_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
+
+              dem =       + eGW(i) - eGW(c) + eGW(k) - eGW(b)
+              num = 2d0*W(b,c,i,k)*W(a,k,j,c)
+
+              KB_dyn(ia,jb) = KB_dyn(ia,jb) - num*dem/(dem**2 + eta**2)
+
+            enddo
+          enddo
+
+          do c=nO+1,nBas-nR
+            do d=nO+1,nBas-nR 
+
+              dem =       + eGW(i) + eGW(j) - eGW(c) - eGW(d)
+              num = 2d0*W(a,b,c,d)*W(c,d,i,j)
+
+              KB_dyn(ia,jb) =  KB_dyn(ia,jb) + num*dem/(dem**2 + eta**2)
+
+            enddo
+          enddo
+
+          do k=nC+1,nO
+            do l=nC+1,nO
+
+              dem =       - eGW(a) - eGW(b) + eGW(k) + eGW(l)
+              num = 2d0*W(a,b,k,l)*W(k,l,i,j)
+
+              KB_dyn(ia,jb) =  KB_dyn(ia,jb) + num*dem/(dem**2 + eta**2)
+
+            end do
+          end do 
+
+        enddo
+      enddo
+    enddo
+  enddo
+!$omp end parallel do
+
+end subroutine 
--- a/src/GW/GW_phBSE_dynamic_kernel_A.f90
+++ b/src/GW/GW_phBSE_dynamic_kernel_A.f90
@ -7,7 +7,12 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh

 ! Input variables

-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
  double precision,intent(in)   :: eta
  double precision,intent(in)   :: lambda
  double precision,intent(in)   :: eGW(nBas)
--- a/src/GW/GW_phBSE_dynamic_kernel_B.f90
+++ b/src/GW/GW_phBSE_dynamic_kernel_B.f90
@ -7,7 +7,12 @@ subroutine GW_phBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh

 ! Input variables

-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
  double precision,intent(in)   :: eta
  double precision,intent(in)   :: lambda
  double precision,intent(in)   :: eGW(nBas)
--- a/src/GW/GW_phBSE_dynamic_perturbation.f90
+++ b/src/GW/GW_phBSE_dynamic_perturbation.f90
@ -65,7 +65,7 @@ subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,e
  if(dophBSE2) then 

    write(*,*)
-    write(*,*) '*** Second-order BSE static kernel activated! ***'
+    write(*,*) '*** Second-order BSE dynamic kernel activated! ***'
    write(*,*)

    allocate(W(nBas,nBas,nBas,nBas))
@ -108,6 +108,8 @@ subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,e
      call GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,-OmBSE(ia),KAm_dyn,ZAm_dyn)
      call GW_phBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,KB_dyn)

+      if(dophBSE2) call GW_phBSE2_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,eGW,W,KB_dyn)
+
      ! Renormalization factor of the resonant and anti-resonant parts

      ZDyn(ia)  = dot_product(X,matmul(ZAp_dyn,X)) &
--- a/src/GW/GW_phBSE_static_kernel_A.f90
+++ b/src/GW/GW_phBSE_static_kernel_A.f90
@ -7,7 +7,12 @@ subroutine GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,KA

 ! Input variables

-  integer,intent(in)            :: nBas,nC,nO,nV,nR,nS
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
  double precision,intent(in)   :: eta
  double precision,intent(in)   :: lambda
  double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)