From 991e975f86a1dda781d38d467ece1538b1997516 Mon Sep 17 00:00:00 2001 From: pfloos Date: Mon, 24 Jul 2023 13:57:56 +0200 Subject: [PATCH] ppBSE@GF2 static kernels --- src/GF/GF2_phBSE2_static_kernel_A.f90 | 26 ++--- src/GF/GF2_phBSE2_static_kernel_B.f90 | 26 ++--- src/GF/GF2_ppBSE2_dynamical_kernel_D.f90 | 4 +- src/GF/GF2_ppBSE2_static_kernel_B.f90 | 136 ++++++++++++++++++++++ src/GF/GF2_ppBSE2_static_kernel_C.f90 | 137 ++++++++++++++++++++++ src/GF/GF2_ppBSE2_static_kernel_D.f90 | 138 +++++++++++++++++++++++ 6 files changed, 439 insertions(+), 28 deletions(-) create mode 100644 src/GF/GF2_ppBSE2_static_kernel_B.f90 create mode 100644 src/GF/GF2_ppBSE2_static_kernel_C.f90 create mode 100644 src/GF/GF2_ppBSE2_static_kernel_D.f90 diff --git a/src/GF/GF2_phBSE2_static_kernel_A.f90 b/src/GF/GF2_phBSE2_static_kernel_A.f90 index 932ceed..45016b1 100644 --- a/src/GF/GF2_phBSE2_static_kernel_A.f90 +++ b/src/GF/GF2_phBSE2_static_kernel_A.f90 @@ -1,4 +1,4 @@ -subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,A_sta) +subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KA_sta) ! Compute the resonant part of the static BSE2 matrix @@ -23,18 +23,18 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e ! Output variables - double precision,intent(out) :: A_sta(nS,nS) + double precision,intent(out) :: KA_sta(nS,nS) ! Initialization - A_sta(:,:) = 0d0 + KA_sta(:,:) = 0d0 ! Second-order correlation kernel for the block A of the singlet manifold if(ispin == 1) then jb = 0 -!$omp parallel do default(private) shared(A_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) +!$omp parallel do default(private) shared(KA_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) do j=nC+1,nO do b=nO+1,nBas-nR jb = (b-nO) + (j-1)*(nBas-nO) @@ -52,13 +52,13 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) & - ERI(j,k,c,i)*ERI(a,c,b,k) + 2d0*ERI(j,k,c,i)*ERI(a,c,k,b) - A_sta(ia,jb) = A_sta(ia,jb) - num*dem/(dem**2 + eta**2) + KA_sta(ia,jb) = KA_sta(ia,jb) - num*dem/(dem**2 + eta**2) dem = + (eGF(c) - eGF(k)) num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) & - ERI(j,c,k,i)*ERI(a,k,b,c) + 2d0*ERI(j,c,k,i)*ERI(a,k,c,b) - A_sta(ia,jb) = A_sta(ia,jb) + num*dem/(dem**2 + eta**2) + KA_sta(ia,jb) = KA_sta(ia,jb) + num*dem/(dem**2 + eta**2) end do end do @@ -70,7 +70,7 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e num = 2d0*ERI(a,j,c,d)*ERI(c,d,i,b) - ERI(a,j,c,d)*ERI(c,d,b,i) & - ERI(a,j,d,c)*ERI(c,d,i,b) + 2d0*ERI(a,j,d,c)*ERI(c,d,b,i) - A_sta(ia,jb) = A_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + KA_sta(ia,jb) = KA_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) end do end do @@ -82,7 +82,7 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e num = 2d0*ERI(a,j,k,l)*ERI(k,l,i,b) - ERI(a,j,k,l)*ERI(k,l,b,i) & - ERI(a,j,l,k)*ERI(k,l,i,b) + 2d0*ERI(a,j,l,k)*ERI(k,l,b,i) - A_sta(ia,jb) = A_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + KA_sta(ia,jb) = KA_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) end do end do @@ -101,7 +101,7 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e if(ispin == 2) then jb = 0 -!$omp parallel do default(private) shared(A_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) +!$omp parallel do default(private) shared(KA_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) do j=nC+1,nO do b=nO+1,nBas-nR jb = (b-nO) + (j-1)*(nBas-nO) @@ -117,12 +117,12 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e dem = - (eGF(c) - eGF(k)) num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) - ERI(j,k,c,i)*ERI(a,c,b,k) - A_sta(ia,jb) = A_sta(ia,jb) - num*dem/(dem**2 + eta**2) + KA_sta(ia,jb) = KA_sta(ia,jb) - num*dem/(dem**2 + eta**2) dem = + (eGF(c) - eGF(k)) num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) - ERI(j,c,k,i)*ERI(a,k,b,c) - A_sta(ia,jb) = A_sta(ia,jb) + num*dem/(dem**2 + eta**2) + KA_sta(ia,jb) = KA_sta(ia,jb) + num*dem/(dem**2 + eta**2) end do end do @@ -133,7 +133,7 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e dem = - (eGF(c) + eGF(d)) num = ERI(a,j,c,d)*ERI(c,d,b,i) + ERI(a,j,d,c)*ERI(c,d,i,b) - A_sta(ia,jb) = A_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + KA_sta(ia,jb) = KA_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) end do end do @@ -144,7 +144,7 @@ subroutine GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e dem = - (eGF(k) + eGF(l)) num = ERI(a,j,k,l)*ERI(k,l,b,i) + ERI(a,j,l,k)*ERI(k,l,i,b) - A_sta(ia,jb) = A_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + KA_sta(ia,jb) = KA_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) end do end do diff --git a/src/GF/GF2_phBSE2_static_kernel_B.f90 b/src/GF/GF2_phBSE2_static_kernel_B.f90 index fdf8a53..961aaa4 100644 --- a/src/GF/GF2_phBSE2_static_kernel_B.f90 +++ b/src/GF/GF2_phBSE2_static_kernel_B.f90 @@ -1,4 +1,4 @@ -subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_sta) +subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_sta) ! Compute the anti-resonant part of the static BSE2 matrix @@ -23,18 +23,18 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e ! Output variables - double precision,intent(out) :: B_sta(nS,nS) + double precision,intent(out) :: KB_sta(nS,nS) ! Initialization - B_sta(:,:) = 0d0 + KB_sta(:,:) = 0d0 ! Second-order correlation kernel for the block A of the singlet manifold if(ispin == 1) then jb = 0 -!$omp parallel do default(private) shared(B_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) +!$omp parallel do default(private) shared(KB_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) do j=nC+1,nO do b=nO+1,nBas-nR jb = (b-nO) + (j-1)*(nBas-nO) @@ -51,13 +51,13 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) & - ERI(b,k,c,i)*ERI(a,c,j,k) + 2d0*ERI(b,k,c,i)*ERI(a,c,k,j) - B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) dem = - eGF(c) + eGF(k) num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) & - ERI(b,c,k,i)*ERI(a,k,j,c) + 2d0*ERI(b,c,k,i)*ERI(a,k,c,j) - B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) end do end do @@ -69,7 +69,7 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e num = 2d0*ERI(a,b,c,d)*ERI(c,d,i,j) - ERI(a,b,c,d)*ERI(c,d,j,i) & - ERI(a,b,d,c)*ERI(c,d,i,j) + 2d0*ERI(a,b,d,c)*ERI(c,d,j,i) - B_sta(ia,jb) = B_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) end do end do @@ -81,7 +81,7 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e num = 2d0*ERI(a,b,k,l)*ERI(k,l,i,j) - ERI(a,b,k,l)*ERI(k,l,j,i) & - ERI(a,b,l,k)*ERI(k,l,i,j) + 2d0*ERI(a,b,l,k)*ERI(k,l,j,i) - B_sta(ia,jb) = B_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) end do end do @@ -100,7 +100,7 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e if(ispin == 2) then jb = 0 -!$omp parallel do default(private) shared(B_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) +!$omp parallel do default(private) shared(KB_sta,ERI,num,dem,eGF,nO,nBas,eta,nC,nR) do j=nC+1,nO do b=nO+1,nBas-nR jb = (b-nO) + (j-1)*(nBas-nO) @@ -116,12 +116,12 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e dem = + eGF(k) - eGF(c) num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) - ERI(b,k,c,i)*ERI(a,c,j,k) - B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) dem = - eGF(c) + eGF(k) num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) - ERI(b,c,k,i)*ERI(a,k,j,c) - B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) end do end do @@ -132,7 +132,7 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e dem = - eGF(c) - eGF(d) num = ERI(a,b,c,d)*ERI(c,d,j,i) + ERI(a,b,d,c)*ERI(c,d,i,j) - B_sta(ia,jb) = B_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) end do end do @@ -143,7 +143,7 @@ subroutine GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,e dem = + eGF(k) + eGF(l) num = ERI(a,b,k,l)*ERI(k,l,j,i) + ERI(a,b,l,k)*ERI(k,l,i,j) - B_sta(ia,jb) = B_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) end do end do diff --git a/src/GF/GF2_ppBSE2_dynamical_kernel_D.f90 b/src/GF/GF2_ppBSE2_dynamical_kernel_D.f90 index 40c0906..b40358b 100644 --- a/src/GF/GF2_ppBSE2_dynamical_kernel_D.f90 +++ b/src/GF/GF2_ppBSE2_dynamical_kernel_D.f90 @@ -29,8 +29,8 @@ subroutine GF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI ! Output variables - double precision,intent(out) :: KC_dyn(nOO,nOO) - double precision,intent(out) :: ZC_dyn(nOO,nOO) + double precision,intent(out) :: KD_dyn(nOO,nOO) + double precision,intent(out) :: ZD_dyn(nOO,nOO) ! Initialization diff --git a/src/GF/GF2_ppBSE2_static_kernel_B.f90 b/src/GF/GF2_ppBSE2_static_kernel_B.f90 new file mode 100644 index 0000000..f7acfcb --- /dev/null +++ b/src/GF/GF2_ppBSE2_static_kernel_B.f90 @@ -0,0 +1,136 @@ +subroutine GF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_sta) + +! Compute the resonant part of the dynamic BSE2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nVV + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + +! Local variables + + double precision :: dem,num + integer :: i,j,k,a,b,c + integer :: ab,ij + +! Output variables + + double precision,intent(out) :: KB_sta(nVV,nOO) + +! Initialization + + KB_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block B of the singlet manifold + + if(ispin == 1) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a,nBas-nR + ab = ab + 1 + + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = eGF(k) - eGF(c) + num = 2d0*ERI(a,k,i,c)*ERI(b,c,j,k) - ERI(a,k,i,c)*ERI(b,c,k,j) & + - ERI(a,k,c,i)*ERI(b,c,j,k) + 2d0*ERI(a,k,c,i)*ERI(b,c,k,j) + + KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(k) - eGF(c) + num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) & + - ERI(b,k,c,i)*ERI(a,c,j,k) + 2d0*ERI(b,k,c,i)*ERI(a,c,k,j) + + KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(k) - eGF(c) + num = 2d0*ERI(a,c,i,k)*ERI(b,k,j,c) - ERI(a,c,i,k)*ERI(b,k,c,j) & + - ERI(a,c,k,i)*ERI(b,k,j,c) + 2d0*ERI(a,c,k,i)*ERI(b,k,c,j) + + KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(k) - eGF(c) + num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) & + - ERI(b,c,k,i)*ERI(a,k,j,c) + 2d0*ERI(b,c,k,i)*ERI(a,k,c,j) + + KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +! Second-order correlation kernel for the block B of the triplet manifold + + if(ispin == 2) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = ab + 1 + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = eGF(k) - eGF(c) + num = 2d0*ERI(a,k,i,c)*ERI(b,c,j,k) - ERI(a,k,i,c)*ERI(b,c,k,j) - ERI(a,k,c,i)*ERI(b,c,j,k) + + KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(k) - eGF(c) + num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) - ERI(b,k,c,i)*ERI(a,c,j,k) + + KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(k) - eGF(c) + num = 2d0*ERI(a,c,i,k)*ERI(b,k,j,c) - ERI(a,c,i,k)*ERI(b,k,c,j) - ERI(a,c,k,i)*ERI(b,k,j,c) + + KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(k) - eGF(c) + num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) - ERI(b,c,k,i)*ERI(a,k,j,c) + + KB_sta(ia,jb) = KB_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +end subroutine diff --git a/src/GF/GF2_ppBSE2_static_kernel_C.f90 b/src/GF/GF2_ppBSE2_static_kernel_C.f90 new file mode 100644 index 0000000..d4c196f --- /dev/null +++ b/src/GF/GF2_ppBSE2_static_kernel_C.f90 @@ -0,0 +1,137 @@ +subroutine GF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,KC_sta) + +! Compute the resonant part of the static BSE2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nVV + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + +! Local variables + + double precision :: dem,num + integer :: m + integer :: a,b,c,d,e + integer :: ab,cd + +! Output variables + + double precision,intent(out) :: KC_sta(nVV,nVV) + +! Initialization + + KC_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block C of the singlet manifold + + if(ispin == 1) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a,nBas-nR + ab = ab + 1 + + cd = 0 + do c=nO+1,nBas-nR + do d=c,nBas-nR + cd = cd + 1 + + do m=nC+1,nO + do e=nO+1,nBas-nR + + dem = eGF(m) - eGF(e) + num = 2d0*ERI(a,m,c,e)*ERI(b,e,d,m) - ERI(a,m,c,e)*ERI(b,e,m,d) & + - ERI(a,m,e,c)*ERI(b,e,d,m) + 2d0*ERI(a,m,e,c)*ERI(b,e,m,d) + + KC_sta(ia,jb) = KC_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(m) - eGF(e) + num = 2d0*ERI(b,m,c,e)*ERI(a,e,d,m) - ERI(b,m,c,e)*ERI(a,e,m,d) & + - ERI(b,m,e,c)*ERI(a,e,d,m) + 2d0*ERI(b,m,e,c)*ERI(a,e,m,d) + + KC_sta(ia,jb) = KC_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(m) - eGF(e) + num = 2d0*ERI(a,e,c,m)*ERI(b,m,d,e) - ERI(a,e,c,m)*ERI(b,m,e,d) & + - ERI(a,e,m,c)*ERI(b,m,d,e) + 2d0*ERI(a,e,m,c)*ERI(b,m,e,d) + + KC_sta(ia,jb) = KC_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(m) - eGF(e) + num = 2d0*ERI(b,e,c,m)*ERI(a,m,d,e) - ERI(b,e,c,m)*ERI(a,m,e,d) & + - ERI(b,e,m,c)*ERI(a,m,d,e) + 2d0*ERI(b,e,c,m)*ERI(a,m,e,d) + + KC_sta(ia,jb) = KC_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +! Second-order correlation kernel for the block C of the triplet manifold + + if(ispin == 2) then + + ab = 0 + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = ab + 1 + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + + do m=nC+1,nO + do e=nO+1,nBas-nR + + dem = eGF(m) - eGF(e) + num = 2d0*ERI(a,m,c,e)*ERI(b,e,d,m) - ERI(a,m,c,e)*ERI(b,e,m,d) - ERI(a,m,e,c)*ERI(b,e,d,m) + + KC_sta(ia,jb) = KC_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(m) - eGF(e) + num = 2d0*ERI(b,m,c,e)*ERI(a,e,d,m) - ERI(b,m,c,e)*ERI(a,e,m,d) - ERI(b,m,e,c)*ERI(a,e,d,m) + + KC_sta(ia,jb) = KC_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(m) - eGF(e) + num = 2d0*ERI(a,e,c,m)*ERI(b,m,d,e) - ERI(a,e,c,m)*ERI(b,m,e,d) - ERI(a,e,m,c)*ERI(b,m,d,e) + + KC_sta(ia,jb) = KC_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = eGF(m) - eGF(e) + num = 2d0*ERI(b,e,c,m)*ERI(a,m,d,e) - ERI(b,e,c,m)*ERI(a,m,e,d) - ERI(b,e,m,c)*ERI(a,m,d,e) + + KC_sta(ia,jb) = KC_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +end subroutine diff --git a/src/GF/GF2_ppBSE2_static_kernel_D.f90 b/src/GF/GF2_ppBSE2_static_kernel_D.f90 new file mode 100644 index 0000000..2006e7d --- /dev/null +++ b/src/GF/GF2_ppBSE2_static_kernel_D.f90 @@ -0,0 +1,138 @@ +subroutine GF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,KD_sta) + +! Compute the resonant part of the static BSE2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOO + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: OmBSE + +! Local variables + + double precision :: dem,num + integer :: i,j,k,l,m + integer :: e + integer :: ij,kl + +! Output variables + + double precision,intent(out) :: KD_sta(nOO,nOO) + +! Initialization + + KD_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block D of the singlet manifold + + if(ispin == 1) then + + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 + + kl = 0 + do k=nC+1,nO + do l=k,nO + kl = kl + 1 + + do m=nC+1,nO + do e=nO+1,nBas-nR + + dem = - eGF(e) + eGF(m) + num = 2d0*ERI(i,e,k,m)*ERI(j,m,l,e) - ERI(i,e,k,m)*ERI(j,m,e,l) & + - ERI(i,e,m,k)*ERI(j,m,l,e) + 2d0*ERI(i,e,m,k)*ERI(j,m,e,l) + + KD_sta(ia,jb) = KD_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = - eGF(e) + eGF(m) + num = 2d0*ERI(j,e,k,m)*ERI(i,m,l,e) - ERI(j,e,k,m)*ERI(i,m,e,l) & + - ERI(j,e,m,k)*ERI(i,m,l,e) + 2d0*ERI(j,e,m,k)*ERI(i,m,e,l) + + KD_sta(ia,jb) = KD_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + dem = - eGF(e) + eGF(m) + num = 2d0*ERI(i,m,k,e)*ERI(j,e,l,m) - ERI(i,m,k,e)*ERI(j,e,m,l) & + - ERI(i,m,e,k)*ERI(j,e,l,m) + 2d0*ERI(i,m,e,k)*ERI(j,e,m,l) + + KD_sta(ia,jb) = KD_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = - eGF(e) + eGF(m) + num = 2d0*ERI(j,m,k,e)*ERI(i,e,l,m) - ERI(j,m,k,e)*ERI(i,e,m,l) & + - ERI(j,m,e,k)*ERI(i,e,l,m) + 2d0*ERI(j,m,e,k)*ERI(i,e,m,l) + + KD_sta(ia,jb) = KD_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +! Second-order correlation kernel for the block D of the triplet manifold + + if(ispin == 2) then + + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + + do m=nC+1,nO + do e=nO+1,nBas-nR + + dem = - eGF(e) + eGF(m) + num = 2d0*ERI(i,e,k,m)*ERI(j,m,l,e) - ERI(i,e,k,m)*ERI(j,m,e,l) - ERI(i,e,m,k)*ERI(j,m,l,e) + + KD_sta(ia,jb) = KD_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = - eGF(e) + eGF(m) + num = 2d0*ERI(j,e,k,m)*ERI(i,m,l,e) - ERI(j,e,k,m)*ERI(i,m,e,l) - ERI(j,e,m,k)*ERI(i,m,l,e) + + KD_sta(ia,jb) = KD_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + dem = - eGF(e) + eGF(m) + num = 2d0*ERI(i,m,k,e)*ERI(j,e,l,m) - ERI(i,m,k,e)*ERI(j,e,m,l) - ERI(i,m,e,k)*ERI(j,e,l,m) + + KD_sta(ia,jb) = KD_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + dem = - eGF(e) + eGF(m) + num = 2d0*ERI(j,m,k,e)*ERI(i,e,l,m) - ERI(j,m,k,e)*ERI(i,e,m,l) - ERI(j,m,e,k)*ERI(i,e,l,m) + + KD_sta(ia,jb) = KD_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +end subroutine