From 127aeabf727b99dd151abb75204549107683211c Mon Sep 17 00:00:00 2001 From: pfloos Date: Wed, 30 Oct 2024 10:44:39 +0100 Subject: [PATCH] more cleaning --- src/GT/GGTpp_ppBSE.f90 | 15 +++--- src/GT/GGTpp_ppBSE_static_kernel_D.f90 | 1 + src/GT/RGT_Tmatrix.f90 | 63 +++++++++++++++----------- src/GT/RGTpp_ppBSE.f90 | 45 ++++++++++-------- src/GT/RGTpp_ppBSE_static_kernel_D.f90 | 7 +-- 5 files changed, 74 insertions(+), 57 deletions(-) diff --git a/src/GT/GGTpp_ppBSE.f90 b/src/GT/GGTpp_ppBSE.f90 index dbc7e15..3f0be29 100644 --- a/src/GT/GGTpp_ppBSE.f90 +++ b/src/GT/GGTpp_ppBSE.f90 @@ -73,9 +73,9 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV,ERI,dipo allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO)) - if(.not.TDA_T) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) - call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGT,ERI,Cpp) - call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGT,ERI,Dpp) + if(.not.TDA) call ppGLR_B(nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,Bpp) + call ppGLR_C(nOrb,nC,nO,nV,nR,nVV,1d0,eGT,ERI,Cpp) + call ppGLR_D(nOrb,nC,nO,nV,nR,nOO,1d0,eGT,ERI,Dpp) !---------------------------------------------- ! Compute T matrix tensor @@ -91,10 +91,9 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV,ERI,dipo allocate(KB_sta(nVV,nOO),KC_sta(nVV,nVV),KD_sta(nOO,nOO)) - call GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KC_sta) - call GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KD_sta) - if(.not.TDA_T) & - call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KB_sta) + call GGTpp_ppBSE_static_kernel_C(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KC_sta) + call GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KD_sta) + if(.not.TDA) call GGTpp_ppBSE_static_kernel_B(eta,nOrb,nC,nO,nV,nR,nOO,nVV,1d0,ERI,eGT,Om1,rho1,Om2,rho2,T,KB_sta) deallocate(Om1,Om2,rho1,rho2) ! Deallocate the 4-tensor T @@ -110,7 +109,7 @@ subroutine GGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nOO,nVV,ERI,dipo allocate(Om1(nVV),X1(nVV,nVV),Y1(nOO,nVV),Om2(nOO),X2(nVV,nOO),Y2(nOO,nOO)) - call ppGLR(TDA_T,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE) + call ppGLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE) call ppLR_transition_vectors(.true.,nOrb,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) diff --git a/src/GT/GGTpp_ppBSE_static_kernel_D.f90 b/src/GT/GGTpp_ppBSE_static_kernel_D.f90 index 15e2e4f..5052c6b 100644 --- a/src/GT/GGTpp_ppBSE_static_kernel_D.f90 +++ b/src/GT/GGTpp_ppBSE_static_kernel_D.f90 @@ -58,6 +58,7 @@ subroutine GGTpp_ppBSE_static_kernel_D(eta,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,e num = num - T(j,m,k,e) * T(e,i,m,l) - T(j,e,k,m) * T(m,i,e,l) KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) + end do end do diff --git a/src/GT/RGT_Tmatrix.f90 b/src/GT/RGT_Tmatrix.f90 index b3eb153..6a8b3d5 100644 --- a/src/GT/RGT_Tmatrix.f90 +++ b/src/GT/RGT_Tmatrix.f90 @@ -38,10 +38,12 @@ subroutine RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,lambda,ERI,Om1 double precision,intent(out) :: T(nBas,nBas,nBas,nBas) ! Initialization + T(:,:,:,:) = 0d0 - ! This corresponds to the alpha alpha alpha alpha elements of T - if (isp_T == 1) then + ! Elements aaaa + + if(isp_T == 1) then !$OMP PARALLEL & !$OMP SHARED(nC,nO,nBas,nR,T,ERI,rho1t,rho2t,Om1t,Om2t) & @@ -52,18 +54,20 @@ subroutine RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,lambda,ERI,Om1 do r=nC+1,nBas-nR do q=nC+1,nBas-nR do p=nC+1,nBas-nR + T(p,q,r,s) = ERI(p,q,r,s) - ERI(p,q,s,r) cd = 0 - do c = nO+1, nBas-nR - do d = c+1, nBas-nR + do c=nO+1,nBas-nR + do d=c+1,nBas-nR cd = cd + 1 T(p,q,r,s) = T(p,q,r,s) - rho1t(p,q,cd) * rho1t(r,s,cd) / Om1t(cd) end do ! d end do ! c + kl = 0 - do k = nC+1, nO - do l = k+1, nO + do k=nC+1,nO + do l=k+1,nO kl = kl + 1 T(p,q,r,s) = T(p,q,r,s) + rho2t(p,q,kl) * rho2t(r,s,kl) / Om2t(kl) enddo ! l @@ -78,8 +82,9 @@ subroutine RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,lambda,ERI,Om1 endif - ! This corresponds to the alpha beta alpha beta elements of T - if (isp_T == 2) then + ! Elements abab + + if(isp_T == 2) then !$OMP PARALLEL & !$OMP SHARED(nC,nO,nBas,nR,T,ERI,rho1s,rho2s,Om1s,Om2s,rho1t,rho2t,Om1t,Om2t) & !$OMP PRIVATE(p,q,r,s,c,d,cd,k,l,kl) & @@ -89,35 +94,36 @@ subroutine RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,lambda,ERI,Om1 do r=nC+1,nBas-nR do q=nC+1,nBas-nR do p=nC+1,nBas-nR + T(p,q,r,s) = ERI(p,q,r,s) cd = 0 - do c = nO+1, nBas-nR - do d = c, nBas-nR + do c=nO+1,nBas-nR + do d=c,nBas-nR cd = cd + 1 T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho1s(p,q,cd) * rho1s(r,s,cd) / Om1s(cd) end do ! d end do ! c cd = 0 - do c = nO+1, nBas-nR - do d = c+1, nBas-nR + do c=nO+1,nBas-nR + do d=c+1,nBas-nR cd = cd + 1 T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho1t(p,q,cd) * rho1t(r,s,cd) / Om1t(cd) end do ! d end do ! c kl = 0 - do k = nC+1, nO - do l = k, nO + do k=nC+1,nO + do l=k,nO kl = kl + 1 T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho2s(p,q,kl) * rho2s(r,s,kl) / Om2s(kl) enddo ! l enddo ! k kl = 0 - do k = nC+1, nO - do l = k+1, nO + do k=nC+1,nO + do l=k+1,nO kl = kl + 1 T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho2t(p,q,kl) * rho2t(r,s,kl) / Om2t(kl) enddo ! l @@ -132,8 +138,10 @@ subroutine RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,lambda,ERI,Om1 endif - ! This corresponds to the beta alpha alpha beta elements of T - if (isp_T == 3) then + ! Elements baab + + if(isp_T == 3) then + !$OMP PARALLEL & !$OMP SHARED(nC,nO,nBas,nR,T,ERI,rho1s,rho2s,Om1s,Om2s,rho1t,rho2t,Om1t,Om2t) & !$OMP PRIVATE(p,q,r,s,c,d,cd,k,l,kl) & @@ -143,33 +151,36 @@ subroutine RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,lambda,ERI,Om1 do r=nC+1,nBas-nR do q=nC+1,nBas-nR do p=nC+1,nBas-nR + T(p,q,r,s) = - ERI(p,q,s,r) cd = 0 - do c = nO+1, nBas-nR - do d = c+1, nBas-nR + do c=nO+1,nBas-nR + do d=c+1,nBas-nR cd = cd + 1 T(p,q,r,s) = T(p,q,r,s) + 0.5d0 * rho1t(p,q,cd) * rho1s(r,s,cd) / Om1t(cd) end do ! d end do ! c + cd = 0 - do c = nO+1, nBas-nR - do d = c, nBas-nR + do c=nO+1,nBas-nR + do d=c,nBas-nR cd = cd + 1 T(p,q,r,s) = T(p,q,r,s) - (1d0 - Kronecker_delta(c,d)) * 0.5d0 * rho1s(p,q,cd) * rho1t(r,s,cd) / Om1s(cd) end do ! d end do ! c kl = 0 - do k = nC+1, nO - do l = k+1, nO + do k=nC+1,nO + do l=k+1,nO kl = kl + 1 T(p,q,r,s) = T(p,q,r,s) - 0.5d0 * rho2t(p,q,kl) * rho2s(r,s,kl) / Om2t(kl) enddo ! l enddo ! k + kl = 0 - do k = nC+1, nO - do l = k, nO + do k=nC+1,nO + do l=k,nO kl = kl + 1 T(p,q,r,s) = T(p,q,r,s) + (1d0 - Kronecker_delta(k,l)) * 0.5d0 * rho2s(p,q,kl) * rho2t(r,s,kl) / Om2s(kl) enddo ! l diff --git a/src/GT/RGTpp_ppBSE.f90 b/src/GT/RGTpp_ppBSE.f90 index 1751075..2da003d 100644 --- a/src/GT/RGTpp_ppBSE.f90 +++ b/src/GT/RGTpp_ppBSE.f90 @@ -59,7 +59,10 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, ! Compute ppRPA excitation density !--------------------------------- + ! Singlet contribution + isp_T = 1 + allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) @@ -68,28 +71,34 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs)) allocate(Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs)) + call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) ! call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOOs,nVVs,dipole_int,Om1s,X1s,Y1s,Om2s,X2s,Y2s) - allocate(rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs)) + call RGTpp_excitation_density(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) deallocate(X1s,Y1s,X2s,Y2s,Bpp,Cpp,Dpp) + ! Triplet contribution + isp_T = 2 - allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt)) - allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt)) + allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) if(.not.TDA_T) call ppLR_B(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) call ppLR_C(isp_T,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) call ppLR_D(isp_T,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) + allocate(Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt)) + allocate(Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt)) + call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) ! call ppLR_transition_vectors(.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,dipole_int,Om1t,X1t,Y1t,Om2t,X2t,Y2t) allocate(rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) + call RGTpp_excitation_density(isp_T,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) deallocate(X1t,Y1t,X2t,Y2t,Bpp,Cpp,Dpp) @@ -98,19 +107,25 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, ! Compute T matrix elements !--------------------------------- - ! This correspond to the alpha-alpha-alpha-alpha elements + ! Elements aaaa + isp_T = 1 allocate(Taaaa(nBas,nBas,nBas,nBas)) + call RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,1d0,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,Taaaa) - ! This correspond to the alpha-alpha-alpha-alpha elements + ! Elements abab + isp_T = 2 allocate(Tabab(nBas,nBas,nBas,nBas)) + call RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,1d0,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,Tabab) - ! This correspond to the alpha-alpha-alpha-alpha elements + ! Elements baab + isp_T = 3 allocate(Tbaab(nBas,nBas,nBas,nBas)) + call RGT_Tmatrix(isp_T,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,1d0,ERI,Om1s,rho1s,Om2s,rho2s,Om1t,rho1t,Om2t,rho2t,Tbaab) deallocate(Om1s,Om2s,Om1t,Om2t,rho1s,rho2s,rho1t,rho2t) @@ -135,12 +150,9 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) - if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT, & - Taaaa,Tabab,Tbaab,KB_sta) - call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT, & - Taaaa,Tabab,Tbaab,KC_sta) - call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT, & - Taaaa,Tabab,Tbaab,KD_sta) + if(.not.TDA) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,Taaaa,Tabab,Tbaab,KB_sta) + call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,Taaaa,Tabab,Tbaab,KC_sta) + call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,Taaaa,Tabab,Tbaab,KD_sta) Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) @@ -176,12 +188,9 @@ subroutine RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR, call ppLR_C(ispin,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) call ppLR_D(ispin,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) - if(.not.TDA_T) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT, & - Taaaa,Tabab,Tbaab,KB_sta) - call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT, & - Taaaa,Tabab,Tbaab,KC_sta) - call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT, & - Taaaa,Tabab,Tbaab,KD_sta) + if(.not.TDA) call RGTpp_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,Taaaa,Tabab,Tbaab,KB_sta) + call RGTpp_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,Taaaa,Tabab,Tbaab,KC_sta) + call RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,Taaaa,Tabab,Tbaab,KD_sta) Bpp(:,:) = Bpp(:,:) + KB_sta(:,:) Cpp(:,:) = Cpp(:,:) + KC_sta(:,:) diff --git a/src/GT/RGTpp_ppBSE_static_kernel_D.f90 b/src/GT/RGTpp_ppBSE_static_kernel_D.f90 index d3c91f4..5e5c41f 100644 --- a/src/GT/RGTpp_ppBSE_static_kernel_D.f90 +++ b/src/GT/RGTpp_ppBSE_static_kernel_D.f90 @@ -104,22 +104,19 @@ subroutine RGTpp_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nOO,nVV,lambda do m=nC+1,nO do e=nO+1,nBas-nR dem = eGF(m) - eGF(e) + num = Taaaa(i,m,k,e)*Taaaa(e,j,m,l) + Tabab(i,m,k,e)*Tabab(e,j,m,l) - KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) num = Taaaa(i,e,k,m)*Taaaa(m,j,e,l) + Tabab(i,e,k,m)*Tabab(m,j,e,l) - KD_sta(ij,kl) = KD_sta(ij,kl) + num*dem/(dem**2 + eta**2) num = Taaaa(j,m,k,e)*Taaaa(e,i,m,l) + Tabab(j,m,k,e)*Tabab(e,i,m,l) - KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2) num = Taaaa(j,e,k,m)*Taaaa(m,i,e,l) + Tabab(j,e,k,m)*Tabab(m,i,e,l) - KD_sta(ij,kl) = KD_sta(ij,kl) - num*dem/(dem**2 + eta**2) - + end do end do