From 3111666d1583f70c25f35e789ebc431c66f81ecc Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Fri, 6 Sep 2024 11:26:47 +0200 Subject: [PATCH 1/9] openMP implementation of the VV block for ppBSE@GF2 --- mol/GW20/Ar.xyz | 3 - mol/GW20/BF.xyz | 4 - mol/GW20/BH3.xyz | 6 - mol/GW20/BN.xyz | 4 - mol/GW20/BeO.xyz | 4 - mol/GW20/CH4.xyz | 7 -- mol/GW20/CO.xyz | 4 - mol/GW20/Cu2.xyz | 4 - mol/GW20/F2.xyz | 4 - mol/GW20/H2.xyz | 4 - mol/GW20/H2O.xyz | 5 - mol/GW20/HCl.xyz | 4 - mol/GW20/HF.xyz | 4 - mol/GW20/He.xyz | 3 - mol/GW20/Li2.xyz | 4 - mol/GW20/LiF.xyz | 4 - mol/GW20/LiH.xyz | 4 - mol/GW20/N2.xyz | 4 - mol/GW20/NH3.xyz | 6 - mol/GW20/Ne.xyz | 3 - mol/GW20/O3.xyz | 5 - mol/GW20/SH2.xyz | 5 - src/GF/RGF2_ppBSE2_static_kernel_C.f90 | 162 ++++++++++++++++++++++--- src/GF/ufRG0F02.f90 | 2 +- src/GW/RGW_ppBSE.f90 | 4 +- src/LR/ppLR_transition_vectors.f90 | 8 +- src/make_ninja.py | 6 +- 27 files changed, 153 insertions(+), 124 deletions(-) delete mode 100644 mol/GW20/Ar.xyz delete mode 100644 mol/GW20/BF.xyz delete mode 100644 mol/GW20/BH3.xyz delete mode 100644 mol/GW20/BN.xyz delete mode 100644 mol/GW20/BeO.xyz delete mode 100644 mol/GW20/CH4.xyz delete mode 100644 mol/GW20/CO.xyz delete mode 100644 mol/GW20/Cu2.xyz delete mode 100644 mol/GW20/F2.xyz delete mode 100644 mol/GW20/H2.xyz delete mode 100644 mol/GW20/H2O.xyz delete mode 100644 mol/GW20/HCl.xyz delete mode 100644 mol/GW20/HF.xyz delete mode 100644 mol/GW20/He.xyz delete mode 100644 mol/GW20/Li2.xyz delete mode 100644 mol/GW20/LiF.xyz delete mode 100644 mol/GW20/LiH.xyz delete mode 100644 mol/GW20/N2.xyz delete mode 100644 mol/GW20/NH3.xyz delete mode 100644 mol/GW20/Ne.xyz delete mode 100644 mol/GW20/O3.xyz delete mode 100644 mol/GW20/SH2.xyz diff --git a/mol/GW20/Ar.xyz b/mol/GW20/Ar.xyz deleted file mode 100644 index cc7c90f..0000000 --- a/mol/GW20/Ar.xyz +++ /dev/null @@ -1,3 +0,0 @@ -1 -Argon; atom; s -Ar 0.0 0.0 0.0 diff --git a/mol/GW20/BF.xyz b/mol/GW20/BF.xyz deleted file mode 100644 index 6c5f116..0000000 --- a/mol/GW20/BF.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Fluoroborane; experimental structure from HCP92; s -B 0.0000 0.0000 0.0000 -F 0.0000 0.0000 1.2626 diff --git a/mol/GW20/BH3.xyz b/mol/GW20/BH3.xyz deleted file mode 100644 index 5a56085..0000000 --- a/mol/GW20/BH3.xyz +++ /dev/null @@ -1,6 +0,0 @@ -4 -Borane; experimental structure from HCP92; s -B 0.0000 0.0000 0.0000 -H 0.0000 0.0000 1.19 -H 0.0000 1.0306 -0.595 -H 0.0000 -1.0306 -0.595 diff --git a/mol/GW20/BN.xyz b/mol/GW20/BN.xyz deleted file mode 100644 index 47257d3..0000000 --- a/mol/GW20/BN.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Boron nitride; experimental structure from HCP92; s -B 0.0000 0.0000 0.0000 -N 0.0000 0.0000 1.281 diff --git a/mol/GW20/BeO.xyz b/mol/GW20/BeO.xyz deleted file mode 100644 index b443ac2..0000000 --- a/mol/GW20/BeO.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Beryllium monoxide; experimental structure from HCP92; s -Be 0.0000 0.0000 0.0000 -O 0.0000 0.0000 1.3308 diff --git a/mol/GW20/CH4.xyz b/mol/GW20/CH4.xyz deleted file mode 100644 index 2f8dc8c..0000000 --- a/mol/GW20/CH4.xyz +++ /dev/null @@ -1,7 +0,0 @@ -5 -Methane; experimental structure from HCP92; m -C 0.0000 0.0000 0.0000 -H 0.6276 -0.6275 0.6276 -H -0.6276 0.6276 0.6276 -H -0.6276 -0.6276 -0.6276 -H 0.6276 0.6276 -0.6276 diff --git a/mol/GW20/CO.xyz b/mol/GW20/CO.xyz deleted file mode 100644 index 0485aa3..0000000 --- a/mol/GW20/CO.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Carbon monoxide; experimental structure from HCP92; s -C 0.0000 0.0000 0.0000 -O 0.0000 0.0000 1.283 diff --git a/mol/GW20/Cu2.xyz b/mol/GW20/Cu2.xyz deleted file mode 100644 index 4a12699..0000000 --- a/mol/GW20/Cu2.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Copper dimer; experimental structure form HCP92; s -Cu 0.0 0.0 0.0 -Cu 0.0 0.0 2.2197 diff --git a/mol/GW20/F2.xyz b/mol/GW20/F2.xyz deleted file mode 100644 index 3316685..0000000 --- a/mol/GW20/F2.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Fluorine; experimental structure from HCP92; s -F 0.0000 0.0000 0.0000 -F 0.0000 0.0000 1.4119 diff --git a/mol/GW20/H2.xyz b/mol/GW20/H2.xyz deleted file mode 100644 index 55105ed..0000000 --- a/mol/GW20/H2.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Hydrogen; experimental structure from HCP92; s -H 0.0000 0.0000 0.0000 -H 0.0000 0.0000 0.74144 diff --git a/mol/GW20/H2O.xyz b/mol/GW20/H2O.xyz deleted file mode 100644 index 4e1148e..0000000 --- a/mol/GW20/H2O.xyz +++ /dev/null @@ -1,5 +0,0 @@ -3 -Water; experimental structure from HCP92; s -O 0.0000 0.0000 0.0000 -H 0.7571 0.0000 0.5861 -H -0.7571 0.0000 0.5861 diff --git a/mol/GW20/HCl.xyz b/mol/GW20/HCl.xyz deleted file mode 100644 index ef7c213..0000000 --- a/mol/GW20/HCl.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Hydrogen chloride; experimental structure from HCP92; s -H 0.0000 0.0000 0.0000 -Cl 0.0000 0.0000 1.2746 diff --git a/mol/GW20/HF.xyz b/mol/GW20/HF.xyz deleted file mode 100644 index 2b500c4..0000000 --- a/mol/GW20/HF.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Hydrogen fluoride; experimental structure from HCP92; s -H 0.0000 0.0000 0.0000 -F 0.0000 0.0000 0.9169 diff --git a/mol/GW20/He.xyz b/mol/GW20/He.xyz deleted file mode 100644 index 4feac49..0000000 --- a/mol/GW20/He.xyz +++ /dev/null @@ -1,3 +0,0 @@ -1 -Helium; atom; s -He 0.0 0.0 0.0 diff --git a/mol/GW20/Li2.xyz b/mol/GW20/Li2.xyz deleted file mode 100644 index ac00597..0000000 --- a/mol/GW20/Li2.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Lithium dimer; experimental structure from HCP92; s -Li 0.0000 0.0000 0.0000 -Li 0.0000 0.0000 2.6729 diff --git a/mol/GW20/LiF.xyz b/mol/GW20/LiF.xyz deleted file mode 100644 index 10e6fe1..0000000 --- a/mol/GW20/LiF.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Lithium fluoride; experimental structure from HCP92; s -Li 0.0000 0.0000 0.0000 -F 0.0000 0.0000 1.5639 diff --git a/mol/GW20/LiH.xyz b/mol/GW20/LiH.xyz deleted file mode 100644 index 6c56774..0000000 --- a/mol/GW20/LiH.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Lithium hydride; experimental structure from HCP92; s -Li 0.0000 0.0000 0.0000 -H 0.0000 0.0000 1.5949 diff --git a/mol/GW20/N2.xyz b/mol/GW20/N2.xyz deleted file mode 100644 index 9498bee..0000000 --- a/mol/GW20/N2.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 -Nitrogen; experimental structure from HCP92; s -N 0.0000 0.0000 0.0000 -N 0.0000 0.0000 1.0977 diff --git a/mol/GW20/NH3.xyz b/mol/GW20/NH3.xyz deleted file mode 100644 index 595067d..0000000 --- a/mol/GW20/NH3.xyz +++ /dev/null @@ -1,6 +0,0 @@ -4 -Amonia; experimental structure from HCP92; s -N 0.0000 0.0000 0.0000 -H 0.0000 -0.9377 -0.3816 -H 0.8121 0.4689 -0.3816 -H -0.8121 0.4689 -0.3816 diff --git a/mol/GW20/Ne.xyz b/mol/GW20/Ne.xyz deleted file mode 100644 index f148885..0000000 --- a/mol/GW20/Ne.xyz +++ /dev/null @@ -1,3 +0,0 @@ -1 -Neon; atom; s -Ne 0.0 0.0 0.0 diff --git a/mol/GW20/O3.xyz b/mol/GW20/O3.xyz deleted file mode 100644 index 17a9481..0000000 --- a/mol/GW20/O3.xyz +++ /dev/null @@ -1,5 +0,0 @@ -3 -Ozon; experimental structure from HCP92; s -O 0.0000 0.0000 0.0000 -O 1.0869 0.0000 0.6600 -O -1.0869 0.0000 0.6600 diff --git a/mol/GW20/SH2.xyz b/mol/GW20/SH2.xyz deleted file mode 100644 index f4069dc..0000000 --- a/mol/GW20/SH2.xyz +++ /dev/null @@ -1,5 +0,0 @@ -3 -Hydrogen sulfide; experimental structure from HCP92; s -S 0.0000 0.0000 0.0000 -H 0.9617 0.0000 0.9268 -H -0.9617 0.0000 0.9268 diff --git a/src/GF/RGF2_ppBSE2_static_kernel_C.f90 b/src/GF/RGF2_ppBSE2_static_kernel_C.f90 index 80ac15f..7ffeadb 100644 --- a/src/GF/RGF2_ppBSE2_static_kernel_C.f90 +++ b/src/GF/RGF2_ppBSE2_static_kernel_C.f90 @@ -25,8 +25,11 @@ subroutine RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI double precision :: dem,num integer :: m integer :: a,b,c,d,e - integer :: ab,cd + integer :: a0,aa,ab,cd + double precision :: eta2 + double precision, allocatable :: Om_tmp(:,:) + ! Output variables double precision,intent(out) :: KC_sta(nVV,nVV) @@ -34,15 +37,39 @@ subroutine RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI ! Initialization KC_sta(:,:) = 0d0 + eta2 = eta * eta + allocate(Om_tmp(nO,nV)) + + ! 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 +! --- --- --- + if(ispin == 1) then - ab = 0 + 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 = ab + 1 + ab = aa + b cd = 0 do c=nO+1,nBas-nR @@ -52,39 +79,91 @@ subroutine RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI 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) - ERI(a,m,e,c)*ERI(b,e,m,d) - KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + KC_sta(ab,cd) = KC_sta(ab,cd) + num * Om_tmp(m,e) - 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) - ERI(b,m,e,c)*ERI(a,e,m,d) - KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + KC_sta(ab,cd) = KC_sta(ab,cd) + num * Om_tmp(m,e) end do - 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))) + 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 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 + +! 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) - ERI(a,m,e,c)*ERI(b,e,m,d) + +! KC_sta(ab,cd) = KC_sta(ab,cd) + 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) - ERI(b,m,e,c)*ERI(a,e,m,d) + +! KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + +! 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 + +! end if + ! Second-order correlation kernel for the block C of the triplet manifold +! --- --- --- +! OpenMP implementation +! --- --- --- + if(ispin == 2) then - ab = 0 - do a=nO+1,nBas-nR - do b=a+1,nBas-nR - ab = ab + 1 + 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 cd = 0 do c=nO+1,nBas-nR @@ -94,17 +173,15 @@ subroutine RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI 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) + ERI(a,m,e,c)*ERI(b,e,m,d) - KC_sta(ab,cd) = KC_sta(ab,cd) + 2d0*num*dem/(dem**2 + eta**2) + KC_sta(ab,cd) = KC_sta(ab,cd) + 2d0 * num * Om_tmp(m,e) - 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) + ERI(b,m,e,c)*ERI(a,e,m,d) - KC_sta(ab,cd) = KC_sta(ab,cd) - 2d0*num*dem/(dem**2 + eta**2) + KC_sta(ab,cd) = KC_sta(ab,cd) - 2d0 * num * Om_tmp(m,e) end do end do @@ -114,7 +191,54 @@ subroutine RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI 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(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) + +! KC_sta(ab,cd) = KC_sta(ab,cd) + 2d0*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) + ERI(b,m,e,c)*ERI(a,e,m,d) + +! KC_sta(ab,cd) = KC_sta(ab,cd) - 2d0*num*dem/(dem**2 + eta**2) + +! end do +! end do + +! end do +! end do + +! end do +! end do + + ! end if + + + deallocate(Om_tmp) + end subroutine diff --git a/src/GF/ufRG0F02.f90 b/src/GF/ufRG0F02.f90 index c5ea84f..5719022 100644 --- a/src/GF/ufRG0F02.f90 +++ b/src/GF/ufRG0F02.f90 @@ -71,7 +71,7 @@ subroutine ufRG0F02(dotest,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Main loop over orbitals ! !-------------------------! - do p=nO-1,nO + do p=nO-3,nO H(:,:) = 0d0 Reigv(:,:) = 0d0 diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90 index 0df6c35..7963b9a 100644 --- a/src/GW/RGW_ppBSE.f90 +++ b/src/GW/RGW_ppBSE.f90 @@ -61,7 +61,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS double precision,allocatable :: KB_sta(:,:) double precision,allocatable :: KC_sta(:,:) double precision,allocatable :: KD_sta(:,:) - + ! Output variables double precision,intent(out) :: EcBSE(nspin) @@ -114,7 +114,7 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS call RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) call RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) - + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) diff --git a/src/LR/ppLR_transition_vectors.f90 b/src/LR/ppLR_transition_vectors.f90 index 2c29aa7..d4155fd 100644 --- a/src/LR/ppLR_transition_vectors.f90 +++ b/src/LR/ppLR_transition_vectors.f90 @@ -117,8 +117,8 @@ subroutine ppLR_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dipole_ ! Thomas-Reiche-Kuhn sum rule - if(nVV > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for p-p sector = ',sum(os1(:)) - write(*,*) +! if(nVV > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for p-p sector = ',sum(os1(:)) +! write(*,*) !-----------------------------------------------! ! Print details about excitations for hh sector ! @@ -188,7 +188,7 @@ subroutine ppLR_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dipole_ ! Thomas-Reiche-Kuhn sum rule - if(nOO > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(os2(:)) - write(*,*) +! if(nOO > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(os2(:)) +! write(*,*) end subroutine diff --git a/src/make_ninja.py b/src/make_ninja.py index a5e6e84..c97e979 100755 --- a/src/make_ninja.py +++ b/src/make_ninja.py @@ -88,9 +88,9 @@ FIX_ORDER_OF_LIBS=-Wl,--start-group """ if sys.platform in ["linux", "linux2"]: -# compiler = compile_gfortran_linux - compiler = compile_ifort_linux -# compiler = compile_olympe +# compiler = compile_gfortran_linux +# compiler = compile_ifort_linux + compiler = compile_olympe elif sys.platform == "darwin": compiler = compile_gfortran_mac else: From 95a7372e2fab608706aca92c634937b7099e1a13 Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Fri, 6 Sep 2024 11:50:43 +0200 Subject: [PATCH 2/9] rename BSE files in GF2 --- src/GF/GG0F2.f90 | 8 +++---- src/GF/{GGF2_phBSE2.f90 => GGF2_phBSE.f90} | 12 +++++----- src/GF/RG0F2.f90 | 8 +++---- src/GF/{RGF2_phBSE2.f90 => RGF2_phBSE.f90} | 22 +++++++++---------- ..._A.f90 => RGF2_phBSE_dynamic_kernel_A.f90} | 4 ++-- ..._B.f90 => RGF2_phBSE_dynamic_kernel_B.f90} | 4 ++-- ...90 => RGF2_phBSE_dynamic_perturbation.f90} | 8 +++---- ...l_A.f90 => RGF2_phBSE_static_kernel_A.f90} | 4 ++-- ...l_B.f90 => RGF2_phBSE_static_kernel_B.f90} | 4 ++-- src/GF/{RGF2_ppBSE2.f90 => RGF2_ppBSE.f90} | 18 +++++++-------- ..._B.f90 => RGF2_ppBSE_dynamic_kernel_B.f90} | 4 ++-- ..._C.f90 => RGF2_ppBSE_dynamic_kernel_C.f90} | 4 ++-- ..._D.f90 => RGF2_ppBSE_dynamic_kernel_D.f90} | 4 ++-- ...90 => RGF2_ppBSE_dynamic_perturbation.f90} | 22 +++++++++---------- ...l_B.f90 => RGF2_ppBSE_static_kernel_B.f90} | 4 ++-- ...l_C.f90 => RGF2_ppBSE_static_kernel_C.f90} | 4 ++-- ...l_D.f90 => RGF2_ppBSE_static_kernel_D.f90} | 4 ++-- src/GF/UG0F2.f90 | 2 +- src/GF/evGGF2.f90 | 8 +++---- src/GF/evRGF2.f90 | 8 +++---- src/GF/evUGF2.f90 | 2 +- src/GF/qsRGF2.f90 | 8 +++---- src/GF/qsUGF2.f90 | 2 +- 23 files changed, 84 insertions(+), 84 deletions(-) rename src/GF/{GGF2_phBSE2.f90 => GGF2_phBSE.f90} (78%) rename src/GF/{RGF2_phBSE2.f90 => RGF2_phBSE.f90} (73%) rename src/GF/{RGF2_phBSE2_dynamic_kernel_A.f90 => RGF2_phBSE_dynamic_kernel_A.f90} (97%) rename src/GF/{RGF2_phBSE2_dynamic_kernel_B.f90 => RGF2_phBSE_dynamic_kernel_B.f90} (96%) rename src/GF/{RGF2_phBSE2_dynamic_perturbation.f90 => RGF2_phBSE_dynamic_perturbation.f90} (89%) rename src/GF/{RGF2_phBSE2_static_kernel_A.f90 => RGF2_phBSE_static_kernel_A.f90} (96%) rename src/GF/{RGF2_phBSE2_static_kernel_B.f90 => RGF2_phBSE_static_kernel_B.f90} (96%) rename src/GF/{RGF2_ppBSE2.f90 => RGF2_ppBSE.f90} (82%) rename src/GF/{RGF2_ppBSE2_dynamic_kernel_B.f90 => RGF2_ppBSE_dynamic_kernel_B.f90} (96%) rename src/GF/{RGF2_ppBSE2_dynamic_kernel_C.f90 => RGF2_ppBSE_dynamic_kernel_C.f90} (96%) rename src/GF/{RGF2_ppBSE2_dynamic_kernel_D.f90 => RGF2_ppBSE_dynamic_kernel_D.f90} (96%) rename src/GF/{RGF2_ppBSE2_dynamic_perturbation.f90 => RGF2_ppBSE_dynamic_perturbation.f90} (81%) rename src/GF/{RGF2_ppBSE2_static_kernel_B.f90 => RGF2_ppBSE_static_kernel_B.f90} (95%) rename src/GF/{RGF2_ppBSE2_static_kernel_C.f90 => RGF2_ppBSE_static_kernel_C.f90} (97%) rename src/GF/{RGF2_ppBSE2_static_kernel_D.f90 => RGF2_ppBSE_static_kernel_D.f90} (95%) diff --git a/src/GF/GG0F2.f90 b/src/GF/GG0F2.f90 index 2e3ee7b..1ae1a66 100644 --- a/src/GF/GG0F2.f90 +++ b/src/GF/GG0F2.f90 @@ -86,11 +86,11 @@ subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, call GMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec) call print_RG0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec) -! Perform BSE2 calculation +! Perform BSE@GF2 calculation if(dophBSE) then - call GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + call GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -101,11 +101,11 @@ subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, end if -! Perform ppBSE2 calculation +! Perform ppBSE@GF2 calculation ! if(doppBSE) then ! -! call GGF2_ppBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) +! call GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) ! write(*,*) ! write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GF/GGF2_phBSE2.f90 b/src/GF/GGF2_phBSE.f90 similarity index 78% rename from src/GF/GGF2_phBSE2.f90 rename to src/GF/GGF2_phBSE.f90 index 307c9f9..f8a75b3 100644 --- a/src/GF/GGF2_phBSE2.f90 +++ b/src/GF/GGF2_phBSE.f90 @@ -1,4 +1,4 @@ -subroutine GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) +subroutine GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) ! Compute the second-order Bethe-Salpeter excitation energies @@ -51,21 +51,21 @@ subroutine GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF, ! Compute static kernel - call RGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) - if(.not.TDA) call RGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) + call RGF2_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) + if(.not.TDA) call RGF2_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) - ! Compute phBSE2@GF2 excitation energies + ! Compute phBSE@GF2 excitation energies call phLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY) - call print_excitation_energies('phBSE2@GGF2','spinorbital',nS,OmBSE) + call print_excitation_energies('phBSE@GGF2','spinorbital',nS,OmBSE) call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) ! Compute dynamic correction for BSE via perturbation theory ! if(dBSE) & -! call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) +! call GF2_phBSE_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) end subroutine diff --git a/src/GF/RG0F2.f90 b/src/GF/RG0F2.f90 index 0f3646e..546a7d0 100644 --- a/src/GF/RG0F2.f90 +++ b/src/GF/RG0F2.f90 @@ -87,11 +87,11 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, call RMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec) call print_RG0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec) -! Perform BSE2 calculation +! Perform BSE@GF2 calculation if(dophBSE) then - call RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -104,11 +104,11 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, end if -! Perform ppBSE2 calculation +! Perform ppBSE@GF2 calculation if(doppBSE) then - call RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) + call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) EcBSE(2) = 3d0*EcBSE(2) diff --git a/src/GF/RGF2_phBSE2.f90 b/src/GF/RGF2_phBSE.f90 similarity index 73% rename from src/GF/RGF2_phBSE2.f90 rename to src/GF/RGF2_phBSE.f90 index 9c3ddef..0e424b5 100644 --- a/src/GF/RGF2_phBSE2.f90 +++ b/src/GF/RGF2_phBSE.f90 @@ -1,4 +1,4 @@ -subroutine RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) +subroutine RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) ! Compute the second-order Bethe-Salpeter excitation energies @@ -59,22 +59,22 @@ subroutine RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI ! Compute static kernel - call RGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) - if(.not.TDA) call RGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) + call RGF2_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) + if(.not.TDA) call RGF2_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) - ! Compute phBSE2@GF2 excitation energies + ! Compute phBSE@GF2 excitation energies call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY) - call print_excitation_energies('phBSE2@GF2','singlet',nS,OmBSE) + call print_excitation_energies('phBSE@GF2','singlet',nS,OmBSE) call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) ! Compute dynamic correction for BSE via perturbation theory if(dBSE) & - call RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) + call RGF2_phBSE_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) end if @@ -92,22 +92,22 @@ subroutine RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI ! Compute static kernel - call RGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) - if(.not.TDA) call RGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) + call RGF2_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) + if(.not.TDA) call RGF2_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) - ! Compute phBSE2@GF2 excitation energies + ! Compute phBSE@GF2 excitation energies call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY) - call print_excitation_energies('phBSE2@GF2','triplet',nS,OmBSE) + call print_excitation_energies('phBSE@GF2','triplet',nS,OmBSE) call phLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) ! Compute dynamic correction for BSE via perturbation theory if(dBSE) & - call RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) + call RGF2_phBSE_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) end if diff --git a/src/GF/RGF2_phBSE2_dynamic_kernel_A.f90 b/src/GF/RGF2_phBSE_dynamic_kernel_A.f90 similarity index 97% rename from src/GF/RGF2_phBSE2_dynamic_kernel_A.f90 rename to src/GF/RGF2_phBSE_dynamic_kernel_A.f90 index cadd452..3aa5a3a 100644 --- a/src/GF/RGF2_phBSE2_dynamic_kernel_A.f90 +++ b/src/GF/RGF2_phBSE_dynamic_kernel_A.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,OmBSE,KA_dyn,ZA_dyn) +subroutine RGF2_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,OmBSE,KA_dyn,ZA_dyn) -! Compute the resonant part of the dynamic BSE2 matrix +! Compute the resonant part of the dynamic BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/RGF2_phBSE2_dynamic_kernel_B.f90 b/src/GF/RGF2_phBSE_dynamic_kernel_B.f90 similarity index 96% rename from src/GF/RGF2_phBSE2_dynamic_kernel_B.f90 rename to src/GF/RGF2_phBSE_dynamic_kernel_B.f90 index 26779ce..55cb206 100644 --- a/src/GF/RGF2_phBSE2_dynamic_kernel_B.f90 +++ b/src/GF/RGF2_phBSE_dynamic_kernel_B.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_dyn) +subroutine RGF2_phBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_dyn) -! Compute the anti-resonant part of the dynamic BSE2 matrix +! Compute the anti-resonant part of the dynamic BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/RGF2_phBSE2_dynamic_perturbation.f90 b/src/GF/RGF2_phBSE_dynamic_perturbation.f90 similarity index 89% rename from src/GF/RGF2_phBSE2_dynamic_perturbation.f90 rename to src/GF/RGF2_phBSE_dynamic_perturbation.f90 index 486d9d3..5587eda 100644 --- a/src/GF/RGF2_phBSE2_dynamic_perturbation.f90 +++ b/src/GF/RGF2_phBSE_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) +subroutine RGF2_phBSE_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -76,7 +76,7 @@ subroutine RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,E ! Resonant part of the BSE correction for dynamical TDA - call RGF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),KAp_dyn,ZAp_dyn) + call RGF2_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),KAp_dyn,ZAp_dyn) if(dTDA) then @@ -87,9 +87,9 @@ subroutine RGF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,E ! Second part of the resonant and anti-resonant part of the BSE correction (frequency independent) - call RGF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),KAm_dyn,ZAm_dyn) + call RGF2_phBSE_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),KAm_dyn,ZAm_dyn) - call RGF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_dyn) + call RGF2_phBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & + dot_product(Y,matmul(ZAm_dyn,Y)) diff --git a/src/GF/RGF2_phBSE2_static_kernel_A.f90 b/src/GF/RGF2_phBSE_static_kernel_A.f90 similarity index 96% rename from src/GF/RGF2_phBSE2_static_kernel_A.f90 rename to src/GF/RGF2_phBSE_static_kernel_A.f90 index c5b1e85..91de55d 100644 --- a/src/GF/RGF2_phBSE2_static_kernel_A.f90 +++ b/src/GF/RGF2_phBSE_static_kernel_A.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KA_sta) +subroutine RGF2_phBSE_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 +! Compute the resonant part of the static BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/RGF2_phBSE2_static_kernel_B.f90 b/src/GF/RGF2_phBSE_static_kernel_B.f90 similarity index 96% rename from src/GF/RGF2_phBSE2_static_kernel_B.f90 rename to src/GF/RGF2_phBSE_static_kernel_B.f90 index 406d568..6a4a36f 100644 --- a/src/GF/RGF2_phBSE2_static_kernel_B.f90 +++ b/src/GF/RGF2_phBSE_static_kernel_B.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_sta) +subroutine RGF2_phBSE_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 +! Compute the anti-resonant part of the static BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/RGF2_ppBSE2.f90 b/src/GF/RGF2_ppBSE.f90 similarity index 82% rename from src/GF/RGF2_ppBSE2.f90 rename to src/GF/RGF2_ppBSE.f90 index 8fb476c..fe626c7 100644 --- a/src/GF/RGF2_ppBSE2.f90 +++ b/src/GF/RGF2_ppBSE.f90 @@ -1,4 +1,4 @@ -subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) +subroutine RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) ! Compute the Bethe-Salpeter excitation energies at the pp level @@ -74,9 +74,9 @@ subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,di ! Compute BSE excitation energies - if(.not.TDA) call RGF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta) - call RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta) - call RGF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta) + if(.not.TDA) call RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta) + call RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta) + call RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp) @@ -95,7 +95,7 @@ subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,di !----------------------------------------------------! if(dBSE) & - call RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & + call RGF2_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) @@ -126,9 +126,9 @@ subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,di ! Compute BSE excitation energies - if(.not.TDA) call RGF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta) - call RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta) - call RGF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta) + if(.not.TDA) call RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta) + call RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta) + call RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta) if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) @@ -148,7 +148,7 @@ subroutine RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,di !----------------------------------------------------! if(dBSE) & - call RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & + call RGF2_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) diff --git a/src/GF/RGF2_ppBSE2_dynamic_kernel_B.f90 b/src/GF/RGF2_ppBSE_dynamic_kernel_B.f90 similarity index 96% rename from src/GF/RGF2_ppBSE2_dynamic_kernel_B.f90 rename to src/GF/RGF2_ppBSE_dynamic_kernel_B.f90 index 7f4cb9c..53f1bba 100644 --- a/src/GF/RGF2_ppBSE2_dynamic_kernel_B.f90 +++ b/src/GF/RGF2_ppBSE_dynamic_kernel_B.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_ppBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_dyn) +subroutine RGF2_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_dyn) -! Compute the resonant part of the dynamic BSE2 matrix +! Compute the resonant part of the dynamic BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/RGF2_ppBSE2_dynamic_kernel_C.f90 b/src/GF/RGF2_ppBSE_dynamic_kernel_C.f90 similarity index 96% rename from src/GF/RGF2_ppBSE2_dynamic_kernel_C.f90 rename to src/GF/RGF2_ppBSE_dynamic_kernel_C.f90 index 985f990..a6052ff 100644 --- a/src/GF/RGF2_ppBSE2_dynamic_kernel_C.f90 +++ b/src/GF/RGF2_ppBSE_dynamic_kernel_C.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,OmBSE,KC_dyn,ZC_dyn) +subroutine RGF2_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,OmBSE,KC_dyn,ZC_dyn) -! Compute the resonant part of the dynamic BSE2 matrix +! Compute the resonant part of the dynamic BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/RGF2_ppBSE2_dynamic_kernel_D.f90 b/src/GF/RGF2_ppBSE_dynamic_kernel_D.f90 similarity index 96% rename from src/GF/RGF2_ppBSE2_dynamic_kernel_D.f90 rename to src/GF/RGF2_ppBSE_dynamic_kernel_D.f90 index 1e1f47a..6d47784 100644 --- a/src/GF/RGF2_ppBSE2_dynamic_kernel_D.f90 +++ b/src/GF/RGF2_ppBSE_dynamic_kernel_D.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,OmBSE,KD_dyn,ZD_dyn) +subroutine RGF2_ppBSE_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,OmBSE,KD_dyn,ZD_dyn) -! Compute the resonant part of the dynamic BSE2 matrix +! Compute the resonant part of the dynamic BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/RGF2_ppBSE2_dynamic_perturbation.f90 b/src/GF/RGF2_ppBSE_dynamic_perturbation.f90 similarity index 81% rename from src/GF/RGF2_ppBSE2_dynamic_perturbation.f90 rename to src/GF/RGF2_ppBSE_dynamic_perturbation.f90 index 75f9296..8920a5a 100644 --- a/src/GF/RGF2_ppBSE2_dynamic_perturbation.f90 +++ b/src/GF/RGF2_ppBSE_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & +subroutine RGF2_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) ! Compute dynamical effects via perturbation theory for BSE @@ -63,7 +63,7 @@ subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO, end if write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) ' First-order dynamical correction to static ppBSE2 double electron attachment energies ' + write(*,*) ' First-order dynamical correction to static ppBSE@GF2 double electron attachment energies ' write(*,*) '---------------------------------------------------------------------------------------------------' write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' @@ -72,16 +72,16 @@ subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO, if(dTDA) then - call RGF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn) + call RGF2_ppBSE_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))) else - call RGF2_ppBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn) - call RGF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn) - call RGF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,Om1(ab),KD_dyn,ZD_dyn) + call RGF2_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn) + call RGF2_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,Om1(ab),KC_dyn,ZC_dyn) + call RGF2_ppBSE_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))) @@ -104,7 +104,7 @@ subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO, write(*,*) write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) ' First-order dynamical correction to static ppBSE2 double electron detachment energies ' + write(*,*) ' First-order dynamical correction to static ppBSE@GF2 double electron detachment energies ' write(*,*) '---------------------------------------------------------------------------------------------------' write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' @@ -115,16 +115,16 @@ subroutine RGF2_ppBSE2_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO, if(dTDA) then - call RGF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,-Om2(ij),KD_dyn,ZD_dyn) + call RGF2_ppBSE_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(Y2(:,ij),matmul(ZD_dyn,Y2(:,ij))) Om2_dyn(kl) = dot_product(Y2(:,ij),matmul(KD_dyn - KD_sta,Y2(:,ij))) else - call RGF2_ppBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn) - call RGF2_ppBSE2_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,-Om2(ij),KC_dyn,ZC_dyn) - call RGF2_ppBSE2_dynamic_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,-Om2(ij),KD_dyn,ZD_dyn) + call RGF2_ppBSE_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_dyn) + call RGF2_ppBSE_dynamic_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,-Om2(ij),KC_dyn,ZC_dyn) + call RGF2_ppBSE_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))) diff --git a/src/GF/RGF2_ppBSE2_static_kernel_B.f90 b/src/GF/RGF2_ppBSE_static_kernel_B.f90 similarity index 95% rename from src/GF/RGF2_ppBSE2_static_kernel_B.f90 rename to src/GF/RGF2_ppBSE_static_kernel_B.f90 index 1797062..f8525a5 100644 --- a/src/GF/RGF2_ppBSE2_static_kernel_B.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_B.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_ppBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_sta) +subroutine RGF2_ppBSE_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 +! Compute the resonant part of the dynamic BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/RGF2_ppBSE2_static_kernel_C.f90 b/src/GF/RGF2_ppBSE_static_kernel_C.f90 similarity index 97% rename from src/GF/RGF2_ppBSE2_static_kernel_C.f90 rename to src/GF/RGF2_ppBSE_static_kernel_C.f90 index 7ffeadb..d26dcfb 100644 --- a/src/GF/RGF2_ppBSE2_static_kernel_C.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_C.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_ppBSE2_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,KC_sta) +subroutine RGF2_ppBSE_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 +! Compute the resonant part of the static BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/RGF2_ppBSE2_static_kernel_D.f90 b/src/GF/RGF2_ppBSE_static_kernel_D.f90 similarity index 95% rename from src/GF/RGF2_ppBSE2_static_kernel_D.f90 rename to src/GF/RGF2_ppBSE_static_kernel_D.f90 index 6f52272..ab02ade 100644 --- a/src/GF/RGF2_ppBSE2_static_kernel_D.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_D.f90 @@ -1,6 +1,6 @@ -subroutine RGF2_ppBSE2_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,KD_sta) +subroutine RGF2_ppBSE_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 +! Compute the resonant part of the static BSE@GF2 matrix implicit none include 'parameters.h' diff --git a/src/GF/UG0F2.f90 b/src/GF/UG0F2.f90 index 3874808..aaff64e 100644 --- a/src/GF/UG0F2.f90 +++ b/src/GF/UG0F2.f90 @@ -124,7 +124,7 @@ subroutine UG0F2(dotest,BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta if(BSE) then - print*,'!!! BSE2 NYI for UG0F2 !!!' + print*,'!!! BSE@GF2 NYI for UG0F2 !!!' end if diff --git a/src/GF/evGGF2.f90 b/src/GF/evGGF2.f90 index 986a130..2631d30 100644 --- a/src/GF/evGGF2.f90 +++ b/src/GF/evGGF2.f90 @@ -142,11 +142,11 @@ subroutine evGGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & end if -! Perform BSE2 calculation +! Perform BSE@GF2 calculation if(dophBSE) then - call GGF2_phBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + call GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -157,11 +157,11 @@ subroutine evGGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & end if -! Perform ppBSE2 calculation +! Perform ppBSE@GF2 calculation ! if(doppBSE) then -! call GGF2_ppBSE2(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) +! call GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) ! write(*,*) ! write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GF/evRGF2.f90 b/src/GF/evRGF2.f90 index 41851a1..4af5c8b 100644 --- a/src/GF/evRGF2.f90 +++ b/src/GF/evRGF2.f90 @@ -145,11 +145,11 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si end if -! Perform BSE2 calculation +! Perform BSE@GF2 calculation if(dophBSE) then - call RGF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) + call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -162,11 +162,11 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si end if -! Perform ppBSE2 calculation +! Perform ppBSE@GF2 calculation if(doppBSE) then - call RGF2_ppBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) + call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GF/evUGF2.f90 b/src/GF/evUGF2.f90 index 1735ffa..fb9ed09 100644 --- a/src/GF/evUGF2.f90 +++ b/src/GF/evUGF2.f90 @@ -197,7 +197,7 @@ subroutine evUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved if(BSE) then - print*,'!!! BSE2 NYI for evUGF2 !!!' + print*,'!!! BSE@GF2 NYI for evUGF2 !!!' end if diff --git a/src/GF/qsRGF2.f90 b/src/GF/qsRGF2.f90 index 8bd1467..1a0afc0 100644 --- a/src/GF/qsRGF2.f90 +++ b/src/GF/qsRGF2.f90 @@ -291,11 +291,11 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis) -! Perform BSE calculation +! Perform phBSE@GF2 calculation if(dophBSE) then - call RGF2_phBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, nC, nO, & + call RGF2_phBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, nC, nO, & nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE) write(*,*) @@ -310,11 +310,11 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & end if -! Perform ppBSE2 calculation +! Perform ppBSE@GF2 calculation if(doppBSE) then - call RGF2_ppBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & + call RGF2_ppBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE) write(*,*) diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 2aedac1..8ab90a2 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -333,7 +333,7 @@ subroutine qsUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved if(BSE) then - print*,'!!! BSE2 NYI for qsUGF2 !!!' + print*,'!!! BSE@GF(2) NYI for qsUGF2 !!!' end if From fa69888f48a80fd0e73fd7b305c12a1ebe482031 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 6 Sep 2024 13:47:15 +0200 Subject: [PATCH 3/9] ppBSE@GW@GHF --- mol/LiF.xyz | 2 +- src/GW/GG0W0.f90 | 33 ++++++----- src/GW/GGW_ppBSE.f90 | 127 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 144 insertions(+), 18 deletions(-) create mode 100644 src/GW/GGW_ppBSE.f90 diff --git a/mol/LiF.xyz b/mol/LiF.xyz index 1c3c4ff..3dd0df2 100644 --- a/mol/LiF.xyz +++ b/mol/LiF.xyz @@ -1,4 +1,4 @@ 2 Li 0.0000 0.0000 0.0000 -F 0.0000 0.0000 1.5783 +F 0.0000 0.0000 1.58753 diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90 index 5d670d2..c313619 100644 --- a/src/GW/GG0W0.f90 +++ b/src/GW/GG0W0.f90 @@ -1,5 +1,5 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) ! Perform G0W0 calculation implicit none @@ -31,7 +31,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF + double precision,intent(in) :: EGHF double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) @@ -156,7 +156,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Dump results ! !--------------! - call print_GG0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA,EcGM) + call print_GG0W0(nBas,nO,eHF,ENuc,EGHF,SigC,Z,eGW,EcRPA,EcGM) ! Deallocate memory @@ -171,7 +171,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@GHF correlation energy = ',EcBSE,' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@GHF total energy = ',ENuc + ERHF + EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0@GHF total energy = ',ENuc + EGHF + EcBSE,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -198,7 +198,7 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au' ! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (triplet) =',EcBSE(2),' au' ! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy =',ENuc + EGHF + EcBSE(1) + EcBSE(2),' au' ! write(*,*)'-------------------------------------------------------------------------------' ! write(*,*) @@ -206,26 +206,25 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA end if -! if(doppBSE) then + 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) + call GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 correlation energy =',EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0W0 total energy =',ENuc + EGHF + EcBSE,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) -! end if + end if ! Testing zone if(dotest) then - call dump_test_value('G','G0W0 correlation energy',EcRPA) + call dump_test_value('G','RPA@G0W0 correlation energy',EcRPA) + call dump_test_value('G','Tr@ppBSE@G0W0 correlation energy',EcBSE) call dump_test_value('G','G0W0 HOMO energy',eGW(nO)) call dump_test_value('G','G0W0 LUMO energy',eGW(nO+1)) diff --git a/src/GW/GGW_ppBSE.f90 b/src/GW/GGW_ppBSE.f90 new file mode 100644 index 0000000..6decd48 --- /dev/null +++ b/src/GW/GGW_ppBSE.f90 @@ -0,0 +1,127 @@ +subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) + +! Compute the Bethe-Salpeter excitation energies at the pp level based on a GHF reference + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + + double precision,intent(in) :: eta + 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) :: eW(nBas) + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: ispin + integer :: isp_W + + logical :: dRPA = .false. + logical :: dRPA_W = .true. + + integer :: nOO + integer :: nVV + + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + + double precision :: EcRPA + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rho_RPA(:,:,:) + + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) + + double precision,allocatable :: Om1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + + double precision,allocatable :: Om2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) + + double precision,allocatable :: KB_sta(:,:) + double precision,allocatable :: KC_sta(:,:) + double precision,allocatable :: KD_sta(:,:) + +! Output variables + + double precision,intent(out) :: EcBSE + +!-----------------------! +! Compute RPA screening ! +!-----------------------! + + isp_W = 3 + EcRPA = 0d0 + + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + Aph(nS,nS),Bph(nS,nS)) + + 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 phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + + call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + + deallocate(XpY_RPA,XmY_RPA,Aph,Bph) + + ispin = 4 + EcBSE = 0d0 + + nOO = nO*(nO-1)/2 + nVV = nV*(nV-1)/2 + + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), & + KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) + + ! Compute BSE excitation energies + + call RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) + call RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) + if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) + + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + + Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) + Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) + Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) + + call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE) + + call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + + !----------------------------------------------------! + ! Compute the dynamical screening at the ppBSE level ! + !----------------------------------------------------! + +! if(dBSE) & +! call GGW_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) + + + deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) + +end subroutine From c9a269cbb61095df88796ad2b5fa38c51276dc3e Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 6 Sep 2024 15:54:24 +0200 Subject: [PATCH 4/9] ppBSE@GF2@GHF --- src/GF/GG0F2.f90 | 26 ++++---- src/GF/GGF2_ppBSE.f90 | 90 +++++++++++++++++++++++++++ src/GF/RGF2_ppBSE_static_kernel_B.f90 | 40 ++++++++++++ src/GF/RGF2_ppBSE_static_kernel_C.f90 | 42 ++++++++++++- src/GF/RGF2_ppBSE_static_kernel_D.f90 | 40 ++++++++++++ 5 files changed, 224 insertions(+), 14 deletions(-) create mode 100644 src/GF/GGF2_ppBSE.f90 diff --git a/src/GF/GG0F2.f90 b/src/GF/GG0F2.f90 index 1ae1a66..bfd3a91 100644 --- a/src/GF/GG0F2.f90 +++ b/src/GF/GG0F2.f90 @@ -1,5 +1,5 @@ subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) ! Perform a one-shot second-order Green function calculation @@ -25,7 +25,7 @@ subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ENuc - double precision,intent(in) :: ERHF + double precision,intent(in) :: EGHF double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -84,7 +84,7 @@ subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, ! Print results call GMP2(.false.,regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,eGF,Ec) - call print_RG0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec) + call print_RG0F2(nBas,nO,eHF,SigC,eGF,Z,ENuc,EGHF,Ec) ! Perform BSE@GF2 calculation @@ -103,18 +103,18 @@ subroutine GG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,linearize,eta,regularize, ! Perform ppBSE@GF2 calculation -! if(doppBSE) then -! -! call GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) + if(doppBSE) then + + call GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) -! write(*,*) -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 correlation energy =',EcBSE,' au' -! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 total energy =',ENuc + ERHF + EcBSE,' au' -! write(*,*)'-------------------------------------------------------------------------------' -! write(*,*) + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 correlation energy =',EcBSE,' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@GG0F2 total energy =',ENuc + EGHF + EcBSE,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) -! end if + end if ! Testing zone diff --git a/src/GF/GGF2_ppBSE.f90 b/src/GF/GGF2_ppBSE.f90 new file mode 100644 index 0000000..aa89da3 --- /dev/null +++ b/src/GF/GGF2_ppBSE.f90 @@ -0,0 +1,90 @@ +subroutine GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBSE) + +! Compute the Bethe-Salpeter excitation energies at the pp level + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + +! Local variables + + integer :: ispin + + integer :: nOO + integer :: nVV + + double precision,allocatable :: Bpp(:,:) + double precision,allocatable :: Cpp(:,:) + double precision,allocatable :: Dpp(:,:) + + double precision,allocatable :: Om1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + + double precision,allocatable :: Om2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) + + double precision,allocatable :: KB_sta(:,:) + double precision,allocatable :: KC_sta(:,:) + double precision,allocatable :: KD_sta(:,:) + +! Output variables + + double precision,intent(out) :: EcBSE + + ispin = 4 + EcBSE = 0d0 + + nOO = nO*(nO-1)/2 + nVV = nV*(nV-1)/2 + + allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), & + KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) + + ! Compute BSE excitation energies + + if(.not.TDA) call RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta) + call RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta) + call RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta) + + if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp) + call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp) + + Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) + Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) + Dpp(:,:) = Dpp(:,:) + KD_sta(:,:) + + call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE) + + call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + + !----------------------------------------------------! + ! Compute the dynamical screening at the ppBSE level ! + !----------------------------------------------------! + +! if(dBSE) & +! call RGF2_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & +! Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) + + deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) + +end subroutine diff --git a/src/GF/RGF2_ppBSE_static_kernel_B.f90 b/src/GF/RGF2_ppBSE_static_kernel_B.f90 index f8525a5..be2572d 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_B.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_B.f90 @@ -116,4 +116,44 @@ subroutine RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda, end if +! Second-order correlation kernel for the block B of the spinorbital manifold + + if(ispin == 4) 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 = 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) + ERI(a,k,c,i)*ERI(b,c,k,j) + + KB_sta(ab,ij) = KB_sta(ab,ij) + 2d0*num*dem/(dem**2 + eta**2) + + dem = eGF(k) - eGF(c) + num = 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) + ERI(b,k,c,i)*ERI(a,c,k,j) + + KB_sta(ab,ij) = KB_sta(ab,ij) - 2d0*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/RGF2_ppBSE_static_kernel_C.f90 b/src/GF/RGF2_ppBSE_static_kernel_C.f90 index d26dcfb..6c47166 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_C.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_C.f90 @@ -39,7 +39,8 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, KC_sta(:,:) = 0d0 eta2 = eta * eta - allocate(Om_tmp(nO,nV)) + 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) @@ -196,7 +197,46 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, end if + if(ispin == 4) 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 = 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) + + KC_sta(ab,cd) = KC_sta(ab,cd) + 2d0*num*dem/(dem**2 + eta**2) + + dem = eGF(m) - eGF(e) + num = 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*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 spinorbital manifold + ! --- --- --- ! Naive implementation ! --- --- --- diff --git a/src/GF/RGF2_ppBSE_static_kernel_D.f90 b/src/GF/RGF2_ppBSE_static_kernel_D.f90 index ab02ade..4f91bd5 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_D.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_D.f90 @@ -116,4 +116,44 @@ subroutine RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI, end if +! Second-order correlation kernel for the block D of the spinorbital manifold + + if(ispin == 4) 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 = 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) + ERI(i,e,m,k)*ERI(j,m,e,l) + + KD_sta(ij,kl) = KD_sta(ij,kl) + 2d0*num*dem/(dem**2 + eta**2) + + dem = - eGF(e) + eGF(m) + num = 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) + ERI(j,e,m,k)*ERI(i,m,e,l) + + KD_sta(ij,kl) = KD_sta(ij,kl) - 2d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + end subroutine From cd25d07f6f5cf61c9cb1265688f62c14ce002e20 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 6 Sep 2024 16:16:16 +0200 Subject: [PATCH 5/9] ebugging GF2 ppBSE spin adaptation --- src/GF/RGF2_ppBSE_static_kernel_B.f90 | 1 + src/GF/RGF2_ppBSE_static_kernel_C.f90 | 331 +++++++++++++------------- src/GF/RGF2_ppBSE_static_kernel_D.f90 | 5 +- 3 files changed, 171 insertions(+), 166 deletions(-) diff --git a/src/GF/RGF2_ppBSE_static_kernel_B.f90 b/src/GF/RGF2_ppBSE_static_kernel_B.f90 index be2572d..8c1af34 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_B.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_B.f90 @@ -68,6 +68,7 @@ subroutine RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda, end do KB_sta(ab,ij) = 2d0*lambda*KB_sta(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) + end do end do diff --git a/src/GF/RGF2_ppBSE_static_kernel_C.f90 b/src/GF/RGF2_ppBSE_static_kernel_C.f90 index 6c47166..7a5aecd 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_C.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_C.f90 @@ -37,22 +37,22 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, ! Initialization KC_sta(:,:) = 0d0 - eta2 = eta * eta +! eta2 = eta * eta - allocate(Om_tmp(nBas,nBas)) - Om_tmp(:,:) = 0d0 +! 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 +! ! 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 @@ -60,93 +60,94 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, ! 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) & +! !$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 - 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 +! +! 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) - 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) +! +! 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) = 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 +! 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 do +! end do +! !$OMP END DO +! !$OMP END PARALLEL - end if +! end if ! --- --- --- ! Naive implementation ! --- --- --- -! if(ispin == 1) then + + 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 -! 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) - ERI(a,m,e,c)*ERI(b,e,m,d) - -! KC_sta(ab,cd) = KC_sta(ab,cd) + 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) - ERI(b,m,e,c)*ERI(a,e,m,d) - -! KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) - -! 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 - -! end if + 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) - ERI(a,m,e,c)*ERI(b,e,m,d) + + KC_sta(ab,cd) = KC_sta(ab,cd) + 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) - ERI(b,m,e,c)*ERI(a,e,m,d) + + KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) + + 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 + + end if ! Second-order correlation kernel for the block C of the triplet manifold @@ -154,47 +155,89 @@ 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) & +! !$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 - cd = 0 - do c=nO+1,nBas-nR - do d=c+1,nBas-nR - cd = cd + 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 +! +! 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) + +! 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 - - 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) + 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) + ERI(a,m,e,c)*ERI(b,e,m,d) - 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*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) + 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 + KC_sta(ab,cd) = KC_sta(ab,cd) - 2d0*num*dem/(dem**2 + eta**2) + + end do + end do - end do - end do - - end do - end do - !$OMP END DO - !$OMP END PARALLEL + end do + end do + end do + end do + end if if(ispin == 4) then @@ -237,48 +280,8 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, ! Second-order correlation kernel for the block C of the spinorbital manifold -! --- --- --- -! 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(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) - -! KC_sta(ab,cd) = KC_sta(ab,cd) + 2d0*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) + ERI(b,m,e,c)*ERI(a,e,m,d) - -! KC_sta(ab,cd) = KC_sta(ab,cd) - 2d0*num*dem/(dem**2 + eta**2) - -! end do -! end do - -! end do -! end do - -! end do -! end do - - ! end if - deallocate(Om_tmp) +! deallocate(Om_tmp) end subroutine diff --git a/src/GF/RGF2_ppBSE_static_kernel_D.f90 b/src/GF/RGF2_ppBSE_static_kernel_D.f90 index 4f91bd5..0838c71 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_D.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_D.f90 @@ -52,13 +52,13 @@ subroutine RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI, do m=nC+1,nO do e=nO+1,nBas-nR - dem = - eGF(e) + eGF(m) + dem = eGF(m) - eGF(e) 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) - ERI(i,e,m,k)*ERI(j,m,e,l) KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) - dem = - eGF(e) + eGF(m) + dem = eGF(m) - eGF(e) 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) - ERI(j,e,m,k)*ERI(i,m,e,l) @@ -68,6 +68,7 @@ subroutine RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI, end do KD_sta(ij,kl) = 2d0*lambda*KD_sta(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + end do end do From 4134aadac5034550c8b690745129efc85246803c Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Fri, 6 Sep 2024 20:40:54 +0200 Subject: [PATCH 6/9] fixed ppBSE@GF2, still need to add openMP on loops --- src/GF/RGF2_ppBSE_static_kernel_B.f90 | 87 +++++++++++++++------------ src/GF/RGF2_ppBSE_static_kernel_C.f90 | 70 ++++++++++++--------- src/GF/RGF2_ppBSE_static_kernel_D.f90 | 75 +++++++++++++---------- 3 files changed, 136 insertions(+), 96 deletions(-) diff --git a/src/GF/RGF2_ppBSE_static_kernel_B.f90 b/src/GF/RGF2_ppBSE_static_kernel_B.f90 index 8c1af34..47a998e 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_B.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_B.f90 @@ -24,7 +24,7 @@ subroutine RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda, double precision,external :: Kronecker_delta double precision :: dem,num - integer :: i,j,k,a,b,c + integer :: i,j,a,b,m,e integer :: ab,ij ! Output variables @@ -49,25 +49,34 @@ subroutine RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda, do j=i,nO ij = ij + 1 - do k=nC+1,nO - do c=nO+1,nBas-nR + do m=nC+1,nO + do e=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) - ERI(a,k,c,i)*ERI(b,c,k,j) - + dem = eGF(m) - eGF(e) + num = 2d0*ERI(a,m,i,e)*ERI(e,b,m,j) - ERI(a,m,i,e)*ERI(e,b,j,m) & + - ERI(a,m,e,i)*ERI(e,b,m,j) - ERI(a,m,e,i)*ERI(e,b,j,m) + KB_sta(ab,ij) = KB_sta(ab,ij) + 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) - ERI(b,k,c,i)*ERI(a,c,k,j) + num = 2d0*ERI(a,e,i,m)*ERI(m,b,e,j) - ERI(a,e,i,m)*ERI(m,b,j,e) & + - ERI(a,e,m,i)*ERI(m,b,e,j) - ERI(a,e,m,i)*ERI(m,b,j,e) + + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + num = 2d0*ERI(b,m,i,e)*ERI(e,a,m,j) - ERI(b,m,i,e)*ERI(e,a,j,m) & + - ERI(b,m,e,i)*ERI(e,a,m,j) - ERI(b,m,e,i)*ERI(e,a,j,m) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + num = 2d0*ERI(b,e,i,m)*ERI(m,a,e,j) - ERI(b,e,i,m)*ERI(m,a,j,e) & + - ERI(b,e,m,i)*ERI(m,a,e,j) - ERI(b,e,m,i)*ERI(m,a,j,e) + + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) end do end do - KB_sta(ab,ij) = 2d0*lambda*KB_sta(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) + KB_sta(ab,ij) = lambda*KB_sta(ab,ij)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) end do end do @@ -91,20 +100,29 @@ subroutine RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda, do j=i+1,nO ij = ij + 1 - do k=nC+1,nO - do c=nO+1,nBas-nR + do m=nC+1,nO + do e=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) + ERI(a,k,c,i)*ERI(b,c,k,j) + dem = eGF(m) - eGF(e) + num = 2d0*ERI(a,m,i,e)*ERI(e,b,m,j) - ERI(a,m,i,e)*ERI(e,b,j,m) & + - ERI(a,m,e,i)*ERI(e,b,m,j) + ERI(a,m,e,i)*ERI(e,b,j,m) + + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) - KB_sta(ab,ij) = KB_sta(ab,ij) + 2d0*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) + ERI(b,k,c,i)*ERI(a,c,k,j) + num = 2d0*ERI(a,e,i,m)*ERI(m,b,e,j) - ERI(a,e,i,m)*ERI(m,b,j,e) & + - ERI(a,e,m,i)*ERI(m,b,e,j) + ERI(a,e,m,i)*ERI(m,b,j,e) + + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + num = 2d0*ERI(b,m,i,e)*ERI(e,a,m,j) - ERI(b,m,i,e)*ERI(e,a,j,m) & + - ERI(b,m,e,i)*ERI(e,a,m,j) + ERI(b,m,e,i)*ERI(e,a,j,m) + + KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2) - KB_sta(ab,ij) = KB_sta(ab,ij) - 2d0*num*dem/(dem**2 + eta**2) + num = 2d0*ERI(b,e,i,m)*ERI(m,a,e,j) - ERI(b,e,i,m)*ERI(m,a,j,e) & + - ERI(b,e,m,i)*ERI(m,a,e,j) + ERI(b,e,m,i)*ERI(m,a,j,e) + + KB_sta(ab,ij) = KB_sta(ab,ij) - num*dem/(dem**2 + eta**2) end do end do @@ -131,21 +149,16 @@ subroutine RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda, do j=i+1,nO ij = ij + 1 - do k=nC+1,nO - do c=nO+1,nBas-nR + do m=nC+1,nO + do e=nO+1,nBas-nR - dem = eGF(k) - eGF(c) - num = 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) + ERI(a,k,c,i)*ERI(b,c,k,j) - - KB_sta(ab,ij) = KB_sta(ab,ij) + 2d0*num*dem/(dem**2 + eta**2) - - dem = eGF(k) - eGF(c) - num = 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) + ERI(b,k,c,i)*ERI(a,c,k,j) - - KB_sta(ab,ij) = KB_sta(ab,ij) - 2d0*num*dem/(dem**2 + eta**2) - + dem = eGF(m) - eGF(e) + num = (ERI(a,m,i,e) - ERI(a,m,e,i)) * (ERI(e,b,m,j) - ERI(e,b,j,m)) + num = num + (ERI(a,e,i,m) - ERI(a,e,m,i)) * (ERI(m,b,e,j) - ERI(m,b,j,e)) + num = num - (ERI(b,m,i,e) - ERI(b,m,e,i)) * (ERI(e,a,m,j) - ERI(e,a,j,m)) + num = num - (ERI(b,e,i,m) - ERI(b,e,m,i)) * (ERI(m,a,e,j) - ERI(m,a,j,e)) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + end do end do diff --git a/src/GF/RGF2_ppBSE_static_kernel_C.f90 b/src/GF/RGF2_ppBSE_static_kernel_C.f90 index 7a5aecd..340e55a 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_C.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_C.f90 @@ -125,21 +125,30 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, 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) - ERI(a,m,e,c)*ERI(b,e,m,d) - + 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) - - 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) - ERI(b,m,e,c)*ERI(a,e,m,d) - + + 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) = 2d0*lambda*KC_sta(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + KC_sta(ab,cd) = lambda*KC_sta(ab,cd)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -217,17 +226,26 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, 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) + ERI(a,m,e,c)*ERI(b,e,m,d) + 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) - KC_sta(ab,cd) = KC_sta(ab,cd) + 2d0*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) + ERI(b,m,e,c)*ERI(a,e,m,d) + 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) - KC_sta(ab,cd) = KC_sta(ab,cd) - 2d0*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 @@ -256,16 +274,12 @@ subroutine RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI, do e=nO+1,nBas-nR dem = eGF(m) - eGF(e) - num = 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) - - KC_sta(ab,cd) = KC_sta(ab,cd) + 2d0*num*dem/(dem**2 + eta**2) - - dem = eGF(m) - eGF(e) - num = 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*dem/(dem**2 + eta**2) + num = (ERI(a,m,c,e) - ERI(a,m,e,c)) * (ERI(e,b,m,d) - ERI(e,b,d,m)) + num = num + (ERI(a,e,c,m) - ERI(a,e,m,c)) * (ERI(m,b,e,d) - ERI(m,b,d,e)) + num = num - (ERI(b,m,c,e) - ERI(b,m,e,c)) * (ERI(e,a,m,d) - ERI(e,a,d,m)) + num = num - (ERI(b,e,c,m) - ERI(b,e,m,c)) * (ERI(m,a,e,d) - ERI(m,a,d,e)) + + KC_sta(ab,cd) = KC_sta(ab,cd) + num*dem/(dem**2 + eta**2) end do end do diff --git a/src/GF/RGF2_ppBSE_static_kernel_D.f90 b/src/GF/RGF2_ppBSE_static_kernel_D.f90 index 0838c71..0474370 100644 --- a/src/GF/RGF2_ppBSE_static_kernel_D.f90 +++ b/src/GF/RGF2_ppBSE_static_kernel_D.f90 @@ -52,22 +52,31 @@ subroutine RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI, do m=nC+1,nO do e=nO+1,nBas-nR - dem = eGF(m) - eGF(e) - 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) - ERI(i,e,m,k)*ERI(j,m,e,l) - + dem = eGF(m) - eGF(e) + num = 2d0*ERI(i,m,k,e)*ERI(e,j,m,l) - ERI(i,m,k,e)*ERI(e,j,l,m) & + - ERI(i,m,e,k)*ERI(e,j,m,l) - ERI(i,m,e,k)*ERI(e,j,l,m) + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) - - dem = eGF(m) - eGF(e) - 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) - ERI(j,e,m,k)*ERI(i,m,e,l) + num = 2d0*ERI(i,e,k,m)*ERI(m,j,e,l) - ERI(i,e,k,m)*ERI(m,j,l,e) & + - ERI(i,e,m,k)*ERI(m,j,e,l) - ERI(i,e,m,k)*ERI(m,j,l,e) + + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + + num = 2d0*ERI(j,m,k,e)*ERI(e,i,m,l) - ERI(j,m,k,e)*ERI(e,i,l,m) & + - ERI(j,m,e,k)*ERI(e,i,m,l) - ERI(j,m,e,k)*ERI(e,i,l,m) + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + + num = 2d0*ERI(j,e,k,m)*ERI(m,i,e,l) - ERI(j,e,k,m)*ERI(m,i,l,e) & + - ERI(j,e,m,k)*ERI(m,i,e,l) - ERI(j,e,m,k)*ERI(m,i,l,e) + + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) end do end do - KD_sta(ij,kl) = 2d0*lambda*KD_sta(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + KD_sta(ij,kl) = lambda*KD_sta(ij,kl)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) end do end do @@ -94,17 +103,26 @@ subroutine RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI, 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) + ERI(i,e,m,k)*ERI(j,m,e,l) + dem = eGF(m) - eGF(e) + num = 2d0*ERI(i,m,k,e)*ERI(e,j,m,l) - ERI(i,m,k,e)*ERI(e,j,l,m) & + - ERI(i,m,e,k)*ERI(e,j,m,l) + ERI(i,m,e,k)*ERI(e,j,l,m) + + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) - KD_sta(ij,kl) = KD_sta(ij,kl) + 2d0*num*dem/(dem**2 + eta**2) + num = 2d0*ERI(i,e,k,m)*ERI(m,j,e,l) - ERI(i,e,k,m)*ERI(m,j,l,e) & + - ERI(i,e,m,k)*ERI(m,j,e,l) + ERI(i,e,m,k)*ERI(m,j,l,e) + + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + + num = 2d0*ERI(j,m,k,e)*ERI(e,i,m,l) - ERI(j,m,k,e)*ERI(e,i,l,m) & + - ERI(j,m,e,k)*ERI(e,i,m,l) + ERI(j,m,e,k)*ERI(e,i,l,m) + + KD_sta(ij,kl) = KD_sta(ij,kl) - 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) + ERI(j,e,m,k)*ERI(i,m,e,l) - - KD_sta(ij,kl) = KD_sta(ij,kl) - 2d0*num*dem/(dem**2 + eta**2) + num = 2d0*ERI(j,e,k,m)*ERI(m,i,e,l) - ERI(j,e,k,m)*ERI(m,i,l,e) & + - ERI(j,e,m,k)*ERI(m,i,e,l) + ERI(j,e,m,k)*ERI(m,i,l,e) + + KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2) end do end do @@ -134,19 +152,14 @@ subroutine RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI, do m=nC+1,nO do e=nO+1,nBas-nR - dem = - eGF(e) + eGF(m) - num = 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) + ERI(i,e,m,k)*ERI(j,m,e,l) - - KD_sta(ij,kl) = KD_sta(ij,kl) + 2d0*num*dem/(dem**2 + eta**2) - - dem = - eGF(e) + eGF(m) - num = 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) + ERI(j,e,m,k)*ERI(i,m,e,l) - - KD_sta(ij,kl) = KD_sta(ij,kl) - 2d0*num*dem/(dem**2 + eta**2) - - end do + dem = eGF(m) - eGF(e) + num = (ERI(i,m,k,e) - ERI(i,m,e,k)) * (ERI(e,j,m,l) - ERI(e,j,l,m)) + num = num + (ERI(i,e,k,m) - ERI(i,e,m,k)) * (ERI(m,j,e,l) - ERI(m,j,l,e)) + num = num - (ERI(j,m,k,e) - ERI(j,m,e,k)) * (ERI(e,i,m,l) - ERI(e,i,l,m)) + num = num - (ERI(j,e,k,m) - ERI(j,e,m,k)) * (ERI(m,i,e,l) - ERI(m,i,l,e)) + + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + end do end do end do From b9da29a020b78c25cd808a21fb148933b66e4f7d Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Fri, 6 Sep 2024 20:54:12 +0200 Subject: [PATCH 7/9] openmp in ppbse@gf2 --- src/GF/RGF2_ppBSE_static_kernel_C.f90 | 323 +++++++++++++------------- 1 file changed, 163 insertions(+), 160 deletions(-) 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 From 27f14363c71284c1b1c2718070a48e66c991d6e2 Mon Sep 17 00:00:00 2001 From: pfloos Date: Fri, 6 Sep 2024 21:35:42 +0200 Subject: [PATCH 8/9] phGLR subroutines --- src/LR/phGLR.f90 | 82 ++++++++++++++++++++++++++++++++++++++++++++++ src/LR/phGLR_A.f90 | 56 +++++++++++++++++++++++++++++++ src/LR/phGLR_B.f90 | 53 ++++++++++++++++++++++++++++++ 3 files changed, 191 insertions(+) create mode 100644 src/LR/phGLR.f90 create mode 100644 src/LR/phGLR_A.f90 create mode 100644 src/LR/phGLR_B.f90 diff --git a/src/LR/phGLR.f90 b/src/LR/phGLR.f90 new file mode 100644 index 0000000..1492eb3 --- /dev/null +++ b/src/LR/phGLR.f90 @@ -0,0 +1,82 @@ +subroutine phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA + integer,intent(in) :: nS + double precision,intent(in) :: Aph(nS,nS) + double precision,intent(in) :: Bph(nS,nS) + + ! Local variables + + double precision :: trace_matrix + double precision,allocatable :: ApB(:,:) + double precision,allocatable :: AmB(:,:) + double precision,allocatable :: AmBSq(:,:) + double precision,allocatable :: AmBIv(:,:) + double precision,allocatable :: Z(:,:) + double precision,allocatable :: tmp(:,:) + +! Output variables + + double precision,intent(out) :: EcRPA + double precision,intent(out) :: Om(nS) + double precision,intent(out) :: XpY(nS,nS) + double precision,intent(out) :: XmY(nS,nS) + +! Memory allocation + + allocate(ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS),tmp(nS,nS)) + +! Tamm-Dancoff approximation + + if(TDA) then + + XpY(:,:) = Aph(:,:) + call diagonalize_matrix(nS,XpY,Om) + XpY(:,:) = transpose(XpY(:,:)) + XmY(:,:) = XpY(:,:) + + else + + ApB(:,:) = Aph(:,:) + Bph(:,:) + AmB(:,:) = Aph(:,:) - Bph(:,:) + + ! Diagonalize linear response matrix + + call diagonalize_matrix(nS,AmB,Om) + + if(minval(Om) < 0d0) & + call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') + + call ADAt(nS,AmB,1d0*dsqrt(Om),AmBSq) + call ADAt(nS,AmB,1d0/dsqrt(Om),AmBIv) + + call dgemm('N','N',nS,nS,nS,1d0,ApB,size(ApB,1),AmBSq,size(AmBSq,1),0d0,tmp,size(tmp,1)) + call dgemm('N','N',nS,nS,nS,1d0,AmBSq,size(AmBSq,1),tmp,size(tmp,1),0d0,Z,size(Z,1)) + + call diagonalize_matrix(nS,Z,Om) + + if(minval(Om) < 0d0) & + call print_warning('You may have instabilities in linear response: negative excitations!!') + + Om = sqrt(Om) + + call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBSq,size(AmBSq,1),0d0,XpY,size(XpY,1)) + call DA(nS,1d0/dsqrt(Om),XpY) + + call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBIv,size(AmBIv,1),0d0,XmY,size(XmY,1)) + call DA(nS,1d0*dsqrt(Om),XmY) + + end if + + ! Compute the RPA correlation energy + + EcRPA = 0.5d0*(sum(Om) - trace_matrix(nS,Aph)) + +end subroutine diff --git a/src/LR/phGLR_A.f90 b/src/LR/phGLR_A.f90 new file mode 100644 index 0000000..29ae8a8 --- /dev/null +++ b/src/LR/phGLR_A.f90 @@ -0,0 +1,56 @@ +subroutine phRLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) + +! Compute resonant block of the ph channel + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: lambda + double precision,intent(in) :: e(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + +! Local variables + + double precision :: delta_dRPA + double precision,external :: Kronecker_delta + + integer :: i,j,a,b,ia,jb + +! Output variables + + double precision,intent(out) :: Aph(nS,nS) + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +! Build A matrix for spin orbitals + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 + + Aph(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + + lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a) + + end do + end do + end do + end do + +end subroutine diff --git a/src/LR/phGLR_B.f90 b/src/LR/phGLR_B.f90 new file mode 100644 index 0000000..acf94ec --- /dev/null +++ b/src/LR/phGLR_B.f90 @@ -0,0 +1,53 @@ +subroutine phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph) + +! Compute the coupling block of the ph channel + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + +! Local variables + + double precision :: delta_dRPA + + integer :: i,j,a,b,ia,jb + +! Output variables + + double precision,intent(out) :: Bph(nS,nS) + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +! Build B matrix for spin orbitals + + ia = 0 + do i=nC+1,nO + do a=nO+1,nOrb-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nOrb-nR + jb = jb + 1 + + Bph(ia,jb) = lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a) + + end do + end do + end do + end do + +end subroutine From f85bd15e03fa5fb1ab700c0fe7509b5ac8c95828 Mon Sep 17 00:00:00 2001 From: pfloos Date: Fri, 6 Sep 2024 22:14:51 +0200 Subject: [PATCH 9/9] creating new routines gor G branch in RPA, LR, GF and GW --- src/GF/GGF2_phBSE.f90 | 14 ++- src/GF/GGF2_phBSE_static_kernel_A.f90 | 94 +++++++++++++++++++ src/GF/GGF2_phBSE_static_kernel_B.f90 | 94 +++++++++++++++++++ src/GF/GGF2_ppBSE.f90 | 19 ++-- src/GF/GGF2_ppBSE_static_kernel_B.f90 | 68 ++++++++++++++ src/GF/GGF2_ppBSE_static_kernel_C.f90 | 69 ++++++++++++++ src/GF/GGF2_ppBSE_static_kernel_D.f90 | 68 ++++++++++++++ src/GW/GG0W0.f90 | 19 ++-- src/GW/GGW_ppBSE.f90 | 38 ++++---- src/GW/GGW_ppBSE_static_kernel_B.f90 | 66 ++++++++++++++ src/GW/GGW_ppBSE_static_kernel_C.f90 | 67 ++++++++++++++ src/GW/GGW_ppBSE_static_kernel_D.f90 | 65 ++++++++++++++ src/GW/evGGW.f90 | 8 +- src/GW/qsGGW.f90 | 10 +-- src/LR/phGLR_A.f90 | 2 +- src/LR/phLR_oscillator_strength.f90 | 12 +-- src/LR/ppGLR.f90 | 125 ++++++++++++++++++++++++++ src/LR/ppGLR_B.f90 | 47 ++++++++++ src/LR/ppGLR_C.f90 | 61 +++++++++++++ src/LR/ppGLR_D.f90 | 54 +++++++++++ src/LR/ppLR.f90 | 5 +- src/RPA/crGRPA.f90 | 21 ++--- src/RPA/phGRPA.f90 | 21 ++--- src/RPA/phGRPAx.f90 | 21 ++--- src/RPA/ppGRPA.f90 | 11 +-- 25 files changed, 962 insertions(+), 117 deletions(-) create mode 100644 src/GF/GGF2_phBSE_static_kernel_A.f90 create mode 100644 src/GF/GGF2_phBSE_static_kernel_B.f90 create mode 100644 src/GF/GGF2_ppBSE_static_kernel_B.f90 create mode 100644 src/GF/GGF2_ppBSE_static_kernel_C.f90 create mode 100644 src/GF/GGF2_ppBSE_static_kernel_D.f90 create mode 100644 src/GW/GGW_ppBSE_static_kernel_B.f90 create mode 100644 src/GW/GGW_ppBSE_static_kernel_C.f90 create mode 100644 src/GW/GGW_ppBSE_static_kernel_D.f90 create mode 100644 src/LR/ppGLR.f90 create mode 100644 src/LR/ppGLR_B.f90 create mode 100644 src/LR/ppGLR_C.f90 create mode 100644 src/LR/ppGLR_D.f90 diff --git a/src/GF/GGF2_phBSE.f90 b/src/GF/GGF2_phBSE.f90 index f8a75b3..a705ecf 100644 --- a/src/GF/GGF2_phBSE.f90 +++ b/src/GF/GGF2_phBSE.f90 @@ -25,7 +25,6 @@ subroutine GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,E ! Local variables logical :: dRPA = .false. - integer :: ispin double precision,allocatable :: OmBSE(:) double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) @@ -43,16 +42,15 @@ subroutine GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,E allocate(OmBSE(nS),XpY(nS,nS),XmY(nS,nS),A_sta(nS,nS),KA_sta(nS,nS)) allocate(B_sta(nS,nS),KB_sta(nS,nS)) - ispin = 3 EcBSE = 0d0 - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta) - if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) + call phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,A_sta) + if(.not.TDA) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) ! Compute static kernel - call RGF2_phBSE_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) - if(.not.TDA) call RGF2_phBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) + call GGF2_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) + if(.not.TDA) call GGF2_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) @@ -60,12 +58,12 @@ subroutine GGF2_phBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,E ! Compute phBSE@GF2 excitation energies call phLR(TDA,nS,A_sta,B_sta,EcBSE,OmBSE,XpY,XmY) - call print_excitation_energies('phBSE@GGF2','spinorbital',nS,OmBSE) + call print_excitation_energies('phBSE@GGF2','generalized',nS,OmBSE) call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) ! Compute dynamic correction for BSE via perturbation theory ! if(dBSE) & -! call GF2_phBSE_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) +! call GGF2_phBSE_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) end subroutine diff --git a/src/GF/GGF2_phBSE_static_kernel_A.f90 b/src/GF/GGF2_phBSE_static_kernel_A.f90 new file mode 100644 index 0000000..91bf2d5 --- /dev/null +++ b/src/GF/GGF2_phBSE_static_kernel_A.f90 @@ -0,0 +1,94 @@ +subroutine GGF2_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KA_sta) + +! Compute the resonant part of the static BSE@GF2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + 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,l + integer :: a,b,c,d + integer :: ia,jb + +! Output variables + + double precision,intent(out) :: KA_sta(nS,nS) + +! Initialization + + KA_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block A of the spinorbital manifold + + jb = 0 + +!$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) + + 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 = - (eGF(c) - eGF(k)) + num = 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) + ERI(j,k,c,i)*ERI(a,c,k,b) + + KA_sta(ia,jb) = KA_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = + (eGF(c) - eGF(k)) + num = 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) + ERI(j,c,k,i)*ERI(a,k,c,b) + + KA_sta(ia,jb) = KA_sta(ia,jb) + num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - (eGF(c) + eGF(d)) + num = 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) + ERI(a,j,d,c)*ERI(c,d,b,i) + + KA_sta(ia,jb) = KA_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = - (eGF(k) + eGF(l)) + num = 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) + ERI(a,j,l,k)*ERI(k,l,b,i) + + KA_sta(ia,jb) = KA_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do +!$omp end parallel do + +end subroutine diff --git a/src/GF/GGF2_phBSE_static_kernel_B.f90 b/src/GF/GGF2_phBSE_static_kernel_B.f90 new file mode 100644 index 0000000..6980268 --- /dev/null +++ b/src/GF/GGF2_phBSE_static_kernel_B.f90 @@ -0,0 +1,94 @@ +subroutine GGF2_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,KB_sta) + +! Compute the anti-resonant part of the static BSE@GF2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + 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,l + integer :: a,b,c,d + integer :: ia,jb + +! Output variables + + double precision,intent(out) :: KB_sta(nS,nS) + +! Initialization + + KB_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block A of the spinorbital manifold + + jb = 0 + +!$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) + + 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 = + eGF(k) - eGF(c) + num = 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) + ERI(b,k,c,i)*ERI(a,c,k,j) + + KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = - eGF(c) + eGF(k) + num = 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) + ERI(b,c,k,i)*ERI(a,k,c,j) + + KB_sta(ia,jb) = KB_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - eGF(c) - eGF(d) + num = 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) + ERI(a,b,d,c)*ERI(c,d,j,i) + + KB_sta(ia,jb) = KB_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = + eGF(k) + eGF(l) + num = 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) + ERI(a,b,l,k)*ERI(k,l,j,i) + + 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 +!$omp end parallel do + +end subroutine diff --git a/src/GF/GGF2_ppBSE.f90 b/src/GF/GGF2_ppBSE.f90 index aa89da3..d38f45c 100644 --- a/src/GF/GGF2_ppBSE.f90 +++ b/src/GF/GGF2_ppBSE.f90 @@ -23,8 +23,6 @@ subroutine GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBS ! Local variables - integer :: ispin - integer :: nOO integer :: nVV @@ -48,7 +46,8 @@ subroutine GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBS double precision,intent(out) :: EcBSE - ispin = 4 +! Initialization + EcBSE = 0d0 nOO = nO*(nO-1)/2 @@ -61,13 +60,13 @@ subroutine GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBS ! Compute BSE excitation energies - if(.not.TDA) call RGF2_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta) - call RGF2_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta) - call RGF2_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta) + if(.not.TDA) call GGF2_ppBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGF,KB_sta) + call GGF2_ppBSE_static_kernel_C(eta,nBas,nC,nO,nV,nR,nVV,1d0,ERI,eGF,KC_sta) + call GGF2_ppBSE_static_kernel_D(eta,nBas,nC,nO,nV,nR,nOO,1d0,ERI,eGF,KD_sta) - if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp) + if(.not.TDA) call ppGLR_B(nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppGLR_C(nBas,nC,nO,nV,nR,nVV,1d0,eGF,ERI,Cpp) + call ppGLR_D(nBas,nC,nO,nV,nR,nOO,1d0,eGF,ERI,Dpp) Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) @@ -82,7 +81,7 @@ subroutine GGF2_ppBSE(TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,ERI,dipole_int,eGF,EcBS !----------------------------------------------------! ! if(dBSE) & -! call RGF2_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & +! call GGF2_ppBSE_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nOO,nVV,eGF,ERI,dipole_int, & ! Om1,X1,Y1,Om2,X2,Y2,KB_sta,KC_sta,KD_sta) deallocate(Om1,X1,Y1,Om2,X2,Y2,Bpp,Cpp,Dpp,KB_sta,KC_sta,KD_sta) diff --git a/src/GF/GGF2_ppBSE_static_kernel_B.f90 b/src/GF/GGF2_ppBSE_static_kernel_B.f90 new file mode 100644 index 0000000..3c8c10d --- /dev/null +++ b/src/GF/GGF2_ppBSE_static_kernel_B.f90 @@ -0,0 +1,68 @@ +subroutine GGF2_ppBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,eGF,KB_sta) + +! Compute the resonant part of the dynamic BSE@GF2 matrix + + 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) :: nOO + 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,external :: Kronecker_delta + double precision :: dem,num + integer :: i,j,a,b,m,e + 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 spinorbital manifold + + 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 m=nC+1,nO + do e=nO+1,nBas-nR + + dem = eGF(m) - eGF(e) + num = (ERI(a,m,i,e) - ERI(a,m,e,i)) * (ERI(e,b,m,j) - ERI(e,b,j,m)) + num = num + (ERI(a,e,i,m) - ERI(a,e,m,i)) * (ERI(m,b,e,j) - ERI(m,b,j,e)) + num = num - (ERI(b,m,i,e) - ERI(b,m,e,i)) * (ERI(e,a,m,j) - ERI(e,a,j,m)) + num = num - (ERI(b,e,i,m) - ERI(b,e,m,i)) * (ERI(m,a,e,j) - ERI(m,a,j,e)) + KB_sta(ab,ij) = KB_sta(ab,ij) + num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + +end subroutine diff --git a/src/GF/GGF2_ppBSE_static_kernel_C.f90 b/src/GF/GGF2_ppBSE_static_kernel_C.f90 new file mode 100644 index 0000000..25badb9 --- /dev/null +++ b/src/GF/GGF2_ppBSE_static_kernel_C.f90 @@ -0,0 +1,69 @@ +subroutine GGF2_ppBSE_static_kernel_C(eta,nBas,nC,nO,nV,nR,nVV,lambda,ERI,eGF,KC_sta) + +! Compute the resonant part of the static BSE@GF2 matrix + + 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) :: 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,external :: Kronecker_delta + double precision :: dem,num + integer :: m + integer :: a,b,c,d,e + integer :: a0,aa,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 + + 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 = (ERI(a,m,c,e) - ERI(a,m,e,c)) * (ERI(e,b,m,d) - ERI(e,b,d,m)) + num = num + (ERI(a,e,c,m) - ERI(a,e,m,c)) * (ERI(m,b,e,d) - ERI(m,b,d,e)) + num = num - (ERI(b,m,c,e) - ERI(b,m,e,c)) * (ERI(e,a,m,d) - ERI(e,a,d,m)) + num = num - (ERI(b,e,c,m) - ERI(b,e,m,c)) * (ERI(m,a,e,d) - 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 subroutine diff --git a/src/GF/GGF2_ppBSE_static_kernel_D.f90 b/src/GF/GGF2_ppBSE_static_kernel_D.f90 new file mode 100644 index 0000000..cddc221 --- /dev/null +++ b/src/GF/GGF2_ppBSE_static_kernel_D.f90 @@ -0,0 +1,68 @@ +subroutine GGF2_ppBSE_static_kernel_D(eta,nBas,nC,nO,nV,nR,nOO,lambda,ERI,eGF,KD_sta) + +! Compute the resonant part of the static BSE@GF2 matrix + + 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) :: 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) + +! Local variables + + double precision,external :: Kronecker_delta + 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 spinorbital manifold + + 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(m) - eGF(e) + num = (ERI(i,m,k,e) - ERI(i,m,e,k)) * (ERI(e,j,m,l) - ERI(e,j,l,m)) + num = num + (ERI(i,e,k,m) - ERI(i,e,m,k)) * (ERI(m,j,e,l) - ERI(m,j,l,e)) + num = num - (ERI(j,m,k,e) - ERI(j,m,e,k)) * (ERI(e,i,m,l) - ERI(e,i,l,m)) + num = num - (ERI(j,e,k,m) - ERI(j,e,m,k)) * (ERI(m,i,e,l) - ERI(m,i,l,e)) + + KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + end do + end do + + end do + end do + + end do + end do + +end subroutine diff --git a/src/GW/GG0W0.f90 b/src/GW/GG0W0.f90 index c313619..9b1e30f 100644 --- a/src/GW/GG0W0.f90 +++ b/src/GW/GG0W0.f90 @@ -40,7 +40,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA logical :: print_W = .true. logical :: dRPA - integer :: ispin double precision :: EcRPA double precision :: EcBSE double precision :: EcGM @@ -85,10 +84,6 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA write(*,*) end if -! Spin manifold - - ispin = 3 - ! Memory allocation allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & @@ -98,12 +93,12 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute screening ! !-------------------! - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) - if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) + call phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + if(.not.TDA_W) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call phGLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) - if(print_W) call print_excitation_energies('phRPA@GHF','spinorbital',nS,Om) + if(print_W) call print_excitation_energies('phRPA@GHF','generalized',nS,Om) !--------------------------! ! Compute spectral weights ! @@ -147,10 +142,10 @@ subroutine GG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA ! Compute the RPA correlation energy - 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 phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) + if(.not.TDA_W) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call phGLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) !--------------! ! Dump results ! diff --git a/src/GW/GGW_ppBSE.f90 b/src/GW/GGW_ppBSE.f90 index 6decd48..e88cd78 100644 --- a/src/GW/GGW_ppBSE.f90 +++ b/src/GW/GGW_ppBSE.f90 @@ -1,4 +1,4 @@ -subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) +subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies at the pp level based on a GHF reference @@ -13,20 +13,19 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, logical,intent(in) :: dTDA double precision,intent(in) :: eta - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS - double precision,intent(in) :: eW(nBas) - double precision,intent(in) :: eGW(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eW(nOrb) + double precision,intent(in) :: eGW(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables - integer :: ispin integer :: isp_W logical :: dRPA = .false. @@ -71,20 +70,19 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, isp_W = 3 EcRPA = 0d0 - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nOrb,nOrb,nS), & Aph(nS,nS),Bph(nS,nS)) - call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph) + call phLR_A(isp_W,dRPA_W,nOrb,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) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call RGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) deallocate(XpY_RPA,XmY_RPA,Aph,Bph) - ispin = 4 EcBSE = 0d0 nOO = nO*(nO-1)/2 @@ -97,13 +95,13 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, ! Compute BSE excitation energies - call RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) - call RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) - if(.not.TDA) call RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) + call GGW_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nS,nVV,1d0,ERI,OmRPA,rho_RPA,KC_sta) + call GGW_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nS,nOO,1d0,ERI,OmRPA,rho_RPA,KD_sta) + if(.not.TDA) call GGW_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,KB_sta) - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) - if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGW,ERI,Cpp) + call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGW,ERI,Dpp) + if(.not.TDA) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) @@ -111,14 +109,14 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE) - call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) + call ppLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) !----------------------------------------------------! ! Compute the dynamical screening at the ppBSE level ! !----------------------------------------------------! ! if(dBSE) & -! call GGW_ppBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & +! call GGW_ppBSE_dynamic_perturbation(dTDA,eta,nOrb,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) diff --git a/src/GW/GGW_ppBSE_static_kernel_B.f90 b/src/GW/GGW_ppBSE_static_kernel_B.f90 new file mode 100644 index 0000000..ac74f26 --- /dev/null +++ b/src/GW/GGW_ppBSE_static_kernel_B.f90 @@ -0,0 +1,66 @@ +subroutine GGW_ppBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Om,rho,KB) + +! Compute the VVOO block of the static screening W for the pp-BSE + + 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 + integer,intent(in) :: nOO + 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) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision,external :: Kronecker_delta + double precision :: chi + double precision :: eps + integer :: a,b,i,j,ab,ij,m + +! Output variables + + double precision,intent(out) :: KB(nVV,nOO) + +! Initialization + + KB(:,:) = 0d0 + +!---------------! +! SpinOrb block ! +!---------------! + + 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 + + chi = 0d0 + do m=1,nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,a,m)*rho(j,b,m)*Om(m)/eps & + + rho(i,b,m)*rho(a,j,m)*Om(m)/eps + end do + + KB(ab,ij) = 2d0*lambda*chi + + end do + end do + end do + end do + +end subroutine diff --git a/src/GW/GGW_ppBSE_static_kernel_C.f90 b/src/GW/GGW_ppBSE_static_kernel_C.f90 new file mode 100644 index 0000000..0e8cd05 --- /dev/null +++ b/src/GW/GGW_ppBSE_static_kernel_C.f90 @@ -0,0 +1,67 @@ +subroutine GGW_ppBSE_static_kernel_C(eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ERI,Om,rho,KC) + +! Compute the VVVV block of the static screening W for the pp-BSE + + 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 + 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) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + 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(:) + double precision, allocatable :: tmp_m(:,:,:) + double precision, allocatable :: tmp(:,:,:,:) + +! Output variables + + double precision,intent(out) :: KC(nVV,nVV) + +!---------------! +! SpinOrb block ! +!---------------! + + 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 + + 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) = 2d0*lambda*chi + + end do + end do + end do + end do + +end subroutine diff --git a/src/GW/GGW_ppBSE_static_kernel_D.f90 b/src/GW/GGW_ppBSE_static_kernel_D.f90 new file mode 100644 index 0000000..37c0672 --- /dev/null +++ b/src/GW/GGW_ppBSE_static_kernel_D.f90 @@ -0,0 +1,65 @@ +subroutine GGW_ppBSE_static_kernel_D(eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ERI,Om,rho,KD) + +! Compute the OOOO block of the static screening W for the pp-BSE + + 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 + 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) :: Om(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision,external :: Kronecker_delta + double precision :: chi + double precision :: eps + integer :: i,j,k,l,ij,kl,m + +! Output variables + + double precision,intent(out) :: KD(nOO,nOO) + +! Initialization + + KD(:,:) = 0d0 + +!---------------! +! SpinOrb block ! +!---------------! + + 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 + + chi = 0d0 + do m=1,nS + eps = Om(m)**2 + eta**2 + chi = chi - rho(i,k,m)*rho(j,l,m)*Om(m)/eps & + + rho(i,l,m)*rho(j,k,m)*Om(m)/eps + end do + + KD(ij,kl) = 2d0*lambda*chi + + end do + end do + end do + end do + +end subroutine diff --git a/src/GW/evGGW.f90 b/src/GW/evGGW.f90 index b8f18ec..7c862e9 100644 --- a/src/GW/evGGW.f90 +++ b/src/GW/evGGW.f90 @@ -43,7 +43,6 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop logical :: linear_mixing logical :: dRPA = .true. - integer :: ispin integer :: nSCF integer :: n_diis double precision :: rcond @@ -100,7 +99,6 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Initialization nSCF = 0 - ispin = 3 n_diis = 0 Conv = 1d0 e_diis(:,:) = 0d0 @@ -118,10 +116,10 @@ subroutine evGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Compute screening - 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 phGLR_A(dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) + if(.not.TDA_W) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call phGLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) ! Compute spectral weights diff --git a/src/GW/qsGGW.f90 b/src/GW/qsGGW.f90 index c98542a..ce6a64a 100644 --- a/src/GW/qsGGW.f90 +++ b/src/GW/qsGGW.f90 @@ -58,7 +58,6 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop integer :: nSCF integer :: nBasSq integer :: nBas2Sq - integer :: ispin integer :: ixyz integer :: n_diis double precision :: ET,ETaa,ETbb @@ -154,7 +153,6 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop nSCF = -1 n_diis = 0 - ispin = 3 Conv = 1d0 P(:,:) = PHF(:,:) eGW(:) = eHF(:) @@ -253,11 +251,11 @@ subroutine qsGGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! Compute linear response - call phLR_A(ispin,dRPA,nBas2,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) - if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas2,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + call phGLR_A(dRPA,nBas2,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) + if(.not.TDA_W) call phGLR_B(dRPA,nBas2,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) - call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) - if(print_W) call print_excitation_energies('phRPA@GW@GHF','spinorbital',nS,Om) + call phGLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + if(print_W) call print_excitation_energies('phRPA@GW@GHF','generalized',nS,Om) ! Compute correlation part of the self-energy diff --git a/src/LR/phGLR_A.f90 b/src/LR/phGLR_A.f90 index 29ae8a8..d23a6b8 100644 --- a/src/LR/phGLR_A.f90 +++ b/src/LR/phGLR_A.f90 @@ -1,4 +1,4 @@ -subroutine phRLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) +subroutine phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) ! Compute resonant block of the ph channel diff --git a/src/LR/phLR_oscillator_strength.f90 b/src/LR/phLR_oscillator_strength.f90 index f8cb55d..796a6a2 100644 --- a/src/LR/phLR_oscillator_strength.f90 +++ b/src/LR/phLR_oscillator_strength.f90 @@ -1,4 +1,4 @@ -subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,XmY,os) +subroutine phLR_oscillator_strength(nOrb,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,XmY,os) ! Compute linear response @@ -7,14 +7,14 @@ subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,X ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS integer,intent(in) :: maxS - double precision :: dipole_int(nBas,nBas,ncart) + double precision :: dipole_int(nOrb,nOrb,ncart) double precision,intent(in) :: Om(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) @@ -44,7 +44,7 @@ subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,X do ixyz=1,ncart jb = 0 do j=nC+1,nO - do b=nO+1,nBas-nR + do b=nO+1,nOrb-nR jb = jb + 1 f(m,ixyz) = f(m,ixyz) + dipole_int(j,b,ixyz)*XpY(m,jb) end do @@ -68,8 +68,4 @@ subroutine phLR_oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Om,XpY,X write(*,*) '---------------------------------------------------------------' write(*,*) -! do m=1,maxS -! write(*,'(I3,3F12.6)') m,Om(m),os(m) -! end do - end subroutine diff --git a/src/LR/ppGLR.f90 b/src/LR/ppGLR.f90 new file mode 100644 index 0000000..536015b --- /dev/null +++ b/src/LR/ppGLR.f90 @@ -0,0 +1,125 @@ +subroutine ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) + + ! + ! Solve the pp-RPA linear eigenvalue problem + ! + ! right eigen-problem: H R = R w + ! left eigen-problem: H.T L = L w + ! + ! where L.T R = 1 + ! + ! + ! (+C +B) + ! H = ( ) where C = C.T and D = D.T + ! (-B.T -D) + ! + ! (w1 0) (X1 X2) (+X1 +X2) + ! w = ( ), R = ( ) and L = ( ) + ! (0 w2) (Y1 Y2) (-Y1 -Y2) + ! + ! + ! the normalisation condition reduces to + ! + ! X1.T X2 - Y1.T Y2 = 0 + ! X1.T X1 - Y1.T Y1 = 1 + ! X2.T X2 - Y2.T Y2 = 1 + ! + + implicit none + include 'parameters.h' + + logical, intent(in) :: TDA + integer, intent(in) :: nOO, nVV + double precision, intent(in) :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO) + double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV) + double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO) + double precision, intent(out) :: EcRPA + + logical :: imp_bio, verbose + integer :: i, j, N + double precision :: EcRPA1, EcRPA2 + double precision :: thr_d, thr_nd, thr_deg + double precision,allocatable :: M(:,:), Z(:,:), Om(:) + + double precision, external :: trace_matrix + + + + N = nOO + nVV + + allocate(M(N,N), Z(N,N), Om(N)) + + if(TDA) then + + X1(:,:) = +Cpp(:,:) + Y1(:,:) = 0d0 + if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1) + + X2(:,:) = 0d0 + Y2(:,:) = -Dpp(:,:) + if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2) + + else + + ! Diagonal blocks + M( 1:nVV , 1:nVV) = + Cpp(1:nVV,1:nVV) + M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO) + + ! Off-diagonal blocks + M( 1:nVV ,nVV+1:nOO+nVV) = - Bpp(1:nVV,1:nOO) + M(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp(1:nVV,1:nOO)) + + if((nOO .eq. 0) .or. (nVV .eq. 0)) then + + ! Diagonalize the p-p matrix + if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z) + ! Split the various quantities in p-p and h-h parts + call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2) + + else + + thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 + thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 + thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not + imp_bio = .True. ! impose bi-orthogonality + verbose = .False. + call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) + + do i = 1, nOO + Om2(i) = Om(i) + do j = 1, nVV + X2(j,i) = Z(j,i) + enddo + do j = 1, nOO + Y2(j,i) = Z(nVV+j,i) + enddo + enddo + + do i = 1, nVV + Om1(i) = Om(nOO+i) + do j = 1, nVV + X1(j,i) = M(j,nOO+i) + enddo + do j = 1, nOO + Y1(j,i) = M(nVV+j,nOO+i) + enddo + enddo + + endif + + end if + + ! Compute the RPA correlation energy + EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp)) + EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp) + EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp) + + if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then + print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' + endif + + deallocate(M, Z, Om) + +end subroutine + + diff --git a/src/LR/ppGLR_B.f90 b/src/LR/ppGLR_B.f90 new file mode 100644 index 0000000..10b5107 --- /dev/null +++ b/src/LR/ppGLR_B.f90 @@ -0,0 +1,47 @@ +subroutine ppGLR_B(nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) + +! Compute the B matrix of the pp channel + + 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) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision,external :: Kronecker_delta + integer :: a,b,i,j,ab,ij + +! Output variables + + double precision,intent(out) :: Bpp(nVV,nOO) + +! Build the B matrix for the spin-orbital basis + + 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 + + Bpp(ab,ij) = lambda*(ERI(a,b,i,j) - ERI(a,b,j,i)) + + end do + end do + end do + end do + +end subroutine diff --git a/src/LR/ppGLR_C.f90 b/src/LR/ppGLR_C.f90 new file mode 100644 index 0000000..583443a --- /dev/null +++ b/src/LR/ppGLR_C.f90 @@ -0,0 +1,61 @@ +subroutine ppGLR_C(nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) + +! Compute the C matrix of the pp channel + + 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) :: nVV + double precision,intent(in) :: lambda + double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: eF + double precision,external :: Kronecker_delta + + integer :: a,b,c,d,ab,cd + integer :: a0, aa + double precision :: e_ab, tmp_ab, delta_ac, tmp_cd + +! Output variables + + double precision,intent(out) :: Cpp(nVV,nVV) + +! Define the chemical potential + +! eF = e(nO) + e(nO+1) + eF = 0d0 + +! Build C matrix for the singlet manifold + + !$OMP PARALLEL & + !$OMP SHARED(Cpp,lambda,ERI,e,eF,nC,nO,nBas,nR) & + !$OMP PRIVATE(c,d,a,b,ab,cd) & + !$OMP DEFAULT(NONE) + !$OMP DO + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = (c-(nO+1))*(nBas-nR-(nO+1)) - (c-1-(nO+1))*(c-(nO+1))/2 + d - c + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = (a-(nO+1))*(nBas-nR-(nO+1)) - (a-1-(nO+1))*(a-(nO+1))/2 + b - a + + Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + lambda*(ERI(a,b,c,d) - ERI(a,b,d,c)) + + end do + end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +end subroutine diff --git a/src/LR/ppGLR_D.f90 b/src/LR/ppGLR_D.f90 new file mode 100644 index 0000000..b5ecf0b --- /dev/null +++ b/src/LR/ppGLR_D.f90 @@ -0,0 +1,54 @@ +subroutine ppGLR_D(nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp) + +! Compute the D matrix of the pp channel + + 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) :: nOO + double precision,intent(in) :: lambda + double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: eF + double precision,external :: Kronecker_delta + + integer :: i,j,k,l,ij,kl + +! Output variables + + double precision,intent(out) :: Dpp(nOO,nOO) + +! Define the chemical potential + +! eF = e(nO) + e(nO+1) + eF = 0d0 + +! Build the D matrix for the spin-orbital basis + + 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 + + Dpp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + lambda*(ERI(i,j,k,l) - ERI(i,j,l,k)) + + end do + end do + end do + end do + +end subroutine diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 582e28f..d5b4146 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -1,7 +1,4 @@ - -! --- - -subroutine ppLR(TDA, nOO, nVV, Bpp, Cpp, Dpp, Om1, X1, Y1, Om2, X2, Y2, EcRPA) +subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! ! Solve the pp-RPA linear eigenvalue problem diff --git a/src/RPA/crGRPA.f90 b/src/RPA/crGRPA.f90 index 2f0a481..33db3a8 100644 --- a/src/RPA/crGRPA.f90 +++ b/src/RPA/crGRPA.f90 @@ -1,4 +1,4 @@ -subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) +subroutine crGRPA(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) ! Crossed-ring channel of the random phase approximation @@ -11,7 +11,7 @@ subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) logical,intent(in) :: dotest logical,intent(in) :: TDA - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -19,13 +19,12 @@ subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: EGHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables - integer :: ispin logical :: dRPA double precision,allocatable :: Aph(:,:) double precision,allocatable :: Bph(:,:) @@ -59,14 +58,12 @@ subroutine crGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS)) - ispin = 3 + call phLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph) + if(.not.TDA) call phLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,-1d0,ERI,Bph) - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Aph) - if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,-1d0,ERI,Bph) - - call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) call print_excitation_energies('crRPA@GHF','spinorbital',nS,Om) - call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) + call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/RPA/phGRPA.f90 b/src/RPA/phGRPA.f90 index 3907920..5f67e8c 100644 --- a/src/RPA/phGRPA.f90 +++ b/src/RPA/phGRPA.f90 @@ -1,4 +1,4 @@ -subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) +subroutine phGRPA(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) ! Perform a direct random phase approximation calculation @@ -11,7 +11,7 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) logical,intent(in) :: dotest logical,intent(in) :: TDA - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -19,13 +19,12 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: EGHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables - integer :: ispin logical :: dRPA double precision :: lambda @@ -62,14 +61,12 @@ subroutine phGRPA(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS)) - ispin = 3 + call phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph) + if(.not.TDA) call phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph) - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,Aph) - if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) - - call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) call print_excitation_energies('phRPA@GHF','spinorbital',nS,Om) - call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) + call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/RPA/phGRPAx.f90 b/src/RPA/phGRPAx.f90 index 1137322..acf4b0c 100644 --- a/src/RPA/phGRPAx.f90 +++ b/src/RPA/phGRPAx.f90 @@ -1,4 +1,4 @@ -subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) +subroutine phGRPAx(dotest,TDA,nOrb,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) ! Perform random phase approximation calculation with exchange (aka TDHF) @@ -11,7 +11,7 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) logical,intent(in) :: dotest logical,intent(in) :: TDA - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -19,13 +19,12 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: EGHF - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables - integer :: ispin logical :: dRPA double precision,allocatable :: Aph(:,:) double precision,allocatable :: Bph(:,:) @@ -59,14 +58,12 @@ subroutine phGRPAx(dotest,TDA,nBas,nC,nO,nV,nR,nS,ENuc,EGHF,ERI,dipole_int,eHF) allocate(Om(nS),XpY(nS,nS),XmY(nS,nS),Aph(nS,nS),Bph(nS,nS)) - ispin = 3 + call phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) + if(.not.TDA) call phLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Aph) - if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph) - - call phLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + call phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY) call print_excitation_energies('phRPAx@GHF','spinorbital',nS,Om) - call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) + call phLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nS,dipole_int,Om,XpY,XmY) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/RPA/ppGRPA.f90 b/src/RPA/ppGRPA.f90 index 6a9c62a..2747d0c 100644 --- a/src/RPA/ppGRPA.f90 +++ b/src/RPA/ppGRPA.f90 @@ -23,7 +23,6 @@ subroutine ppGRPA(dotest,TDA,nBas,nC,nO,nV,nR,ENuc,EGHF,ERI,dipole_int,eHF) ! Local variables - integer :: ispin integer :: nOO integer :: nVV double precision,allocatable :: Bpp(:,:) @@ -50,19 +49,17 @@ subroutine ppGRPA(dotest,TDA,nBas,nC,nO,nV,nR,ENuc,EGHF,ERI,dipole_int,eHF) EcRPA = 0d0 - ispin = 4 - nOO = nO*(nO-1)/2 nVV = nV*(nV-1)/2 allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) - if(.not.TDA) call ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) - call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp) - call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp) + if(.not.TDA) call ppGLR_B(nBas,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppGLR_C(nBas,nC,nO,nV,nR,nVV,1d0,eHF,ERI,Cpp) + call ppGLR_D(nBas,nC,nO,nV,nR,nOO,1d0,eHF,ERI,Dpp) - call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) + call ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) ! call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2)