diff --git a/src/GF/RGF2_ppBSE_static_kernel_C.f90 b/src/GF/RGF2_ppBSE_static_kernel_C.f90 index 340e55a..cf133ae 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_C.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_C.f90 @@ -26,9 +26,6 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, integer :: m integer :: a,b,c,d,e integer :: a0,aa,ab,cd - - double precision :: eta2 - double precision, allocatable :: Om_tmp(:,:) ! Output variables @@ -37,92 +34,32 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, ! Initialization KC_sta(:,:) = 0d0 -! eta2 = eta * eta - -! allocate(Om_tmp(nBas,nBas)) -! Om_tmp(:,:) = 0d0 - -! ! Compute the energy differences and denominator once and store them in a temporary array -! !$OMP PARALLEL DEFAULT(NONE) PRIVATE(m,e,dem) SHARED(nC,nO,nBas,nR, eta2, eGF, Om_tmp) -! !$OMP DO -! do m=nC+1,nO -! do e=nO+1,nBas-nR -! dem = eGF(m) - eGF(e) -! Om_tmp(m,e) = dem / (dem*dem + eta2) -! enddo -! enddo -! !$OMP END DO -! !$OMP END PARALLEL ! Second-order correlation kernel for the block C of the singlet manifold -! --- --- --- -! OpenMP implementation -! --- --- --- +! --- --- --- +! OpenMP implementation +! --- --- --- -! if(ispin == 1) then +if(ispin == 1) then -! a0 = nBas - nR - nO -! !$OMP PARALLEL DEFAULT(NONE) & -! !$OMP PRIVATE(a, b, aa, ab, c, d, cd, m, e, num) & -! !$OMP SHARED(nO, nBas, nR, nC, a0, ERI, Om_tmp, KC_sta, lambda) -! !$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 = aa + b + a0 = nBas - nR - nO + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(a, b, aa, ab, c, d, cd, m, e, num, dem) & + !$OMP SHARED(nO, nBas, nR, nC, a0, ERI, KC_sta, lambda, eGF, eta) + !$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 = aa + b -! 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 -! -! 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) - ERI(a,m,e,c)*ERI(b,e,m,d) + cd = 0 + do c=nO+1,nBas-nR + do d=c,nBas-nR + cd = cd + 1 -! KC_sta(ab,cd) = KC_sta(ab,cd) + num * Om_tmp(m,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) - ERI(b,m,e,c)*ERI(a,e,m,d) - -! KC_sta(ab,cd) = KC_sta(ab,cd) + num * Om_tmp(m,e) -! -! end do -! end do - -! KC_sta(ab,cd) = 2d0*lambda*KC_sta(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) -! -! end do -! end do - -! end do -! end do -! !$OMP END DO -! !$OMP END PARALLEL - -! end if - -! --- --- --- -! Naive implementation -! --- --- --- - - 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 + do m=nC+1,nO + do e=nO+1,nBas-nR dem = eGF(m) - eGF(e) num = 2d0*ERI(a,m,c,e)*ERI(e,b,m,d) - ERI(a,m,c,e)*ERI(e,b,d,m) & @@ -145,18 +82,73 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) - end do - end do - - KC_sta(ab,cd) = lambda*KC_sta(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) - + end do end do - end do - - end do + + KC_sta(ab,cd) = lambda*KC_sta(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + + end do + end do + end do - - end if + end do + !$OMP END DO + !$OMP END PARALLEL + +end if + +! --- --- --- +! Naive implementation +! --- --- --- + +! 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(e,b,m,d) - ERI(a,m,c,e)*ERI(e,b,d,m) & +! - ERI(a,m,e,c)*ERI(e,b,m,d) - ERI(a,m,e,c)*ERI(e,b,d,m) + +! KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + +! num = 2d0*ERI(a,e,c,m)*ERI(m,b,e,d) - ERI(a,e,c,m)*ERI(m,b,d,e) & +! - ERI(a,e,m,c)*ERI(m,b,e,d) - ERI(a,e,m,c)*ERI(m,b,d,e) + +! KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + +! num = 2d0*ERI(b,m,c,e)*ERI(e,a,m,d) - ERI(b,m,c,e)*ERI(e,a,d,m) & +! - ERI(b,m,e,c)*ERI(e,a,m,d) - ERI(b,m,e,c)*ERI(e,a,d,m) + +! KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + +! num = 2d0*ERI(b,e,c,m)*ERI(m,a,e,d) - ERI(b,e,c,m)*ERI(m,a,d,e) & +! - ERI(b,e,m,c)*ERI(m,a,e,d) - ERI(b,e,m,c)*ERI(m,a,d,e) + +! KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + +! end do +! end do + +! KC_sta(ab,cd) = lambda*KC_sta(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + +! end do +! end do + +! end do +! end do + +! end if ! Second-order correlation kernel for the block C of the triplet manifold @@ -164,68 +156,26 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, ! OpenMP implementation ! --- --- --- -! if(ispin == 2) then +if(ispin == 2) then -! a0 = nBas - nR - nO - 1 -! !$OMP PARALLEL DEFAULT(NONE) & -! !$OMP PRIVATE(a, b, aa, ab, c, d, cd, m, e, num) & -! !$OMP SHARED(nO, nBas, nR, nC, a0, ERI, Om_tmp, KC_sta) -! !$OMP DO -! do a = nO+1, nBas-nR -! aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO - 1 -! do b = a+1, nBas-nR -! ab = aa + b + a0 = nBas - nR - nO - 1 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(a, b, aa, ab, c, d, cd, m, e, num, dem) & + !$OMP SHARED(nO, nBas, nR, nC, a0, ERI, KC_sta, eGF, eta) + !$OMP DO + do a = nO+1, nBas-nR + aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO - 1 + do b = a+1, nBas-nR + ab = aa + b -! 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 -! -! 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) + ERI(a,m,e,c)*ERI(b,e,m,d) + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 -! KC_sta(ab,cd) = KC_sta(ab,cd) + 2d0 * num * Om_tmp(m,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) + ERI(b,m,e,c)*ERI(a,e,m,d) - -! KC_sta(ab,cd) = KC_sta(ab,cd) - 2d0 * num * Om_tmp(m,e) -! -! end do -! end do - -! end do -! end do - -! end do -! end do -! !$OMP END DO -! !$OMP END PARALLEL - -! end if - -! --- --- --- -! Naive implementation -! --- --- --- - - 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 - + do m=nC+1,nO + do e=nO+1,nBas-nR + dem = eGF(m) - eGF(e) num = 2d0*ERI(a,m,c,e)*ERI(e,b,m,d) - ERI(a,m,c,e)*ERI(e,b,d,m) & - ERI(a,m,e,c)*ERI(e,b,m,d) + ERI(a,m,e,c)*ERI(e,b,d,m) @@ -246,18 +196,71 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, - ERI(b,e,m,c)*ERI(m,a,e,d) + ERI(b,e,m,c)*ERI(m,a,d,e) KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2) - - end do - end do + + end do + end do - end do - end do + end do + end do - end do - end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +end if + +! --- --- --- +! Naive implementation +! --- --- --- + +! 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(e,b,m,d) - ERI(a,m,c,e)*ERI(e,b,d,m) & +! - ERI(a,m,e,c)*ERI(e,b,m,d) + ERI(a,m,e,c)*ERI(e,b,d,m) + +! KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + +! num = 2d0*ERI(a,e,c,m)*ERI(m,b,e,d) - ERI(a,e,c,m)*ERI(m,b,d,e) & +! - ERI(a,e,m,c)*ERI(m,b,e,d) + ERI(a,e,m,c)*ERI(m,b,d,e) + +! KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + +! num = 2d0*ERI(b,m,c,e)*ERI(e,a,m,d) - ERI(b,m,c,e)*ERI(e,a,d,m) & +! - ERI(b,m,e,c)*ERI(e,a,m,d) + ERI(b,m,e,c)*ERI(e,a,d,m) + +! KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2) + +! num = 2d0*ERI(b,e,c,m)*ERI(m,a,e,d) - ERI(b,e,c,m)*ERI(m,a,d,e) & +! - ERI(b,e,m,c)*ERI(m,a,e,d) + ERI(b,e,m,c)*ERI(m,a,d,e) + +! KC_sta(ab,cd) = KC_sta(ab,cd) - num*dem/(dem**2 + eta**2) + +! end do +! end do + +! end do +! end do + +! end do +! end do + +! end if - end if - if(ispin == 4) then ab = 0