LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:50 +01:00
Code Releases Activity

more cleaning

Browse Source
This commit is contained in:
Pierre-Francois Loos 2024-10-30 10:44:39 +01:00
parent 1a4ad3e225
commit 127aeabf72
5 changed files with 74 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
src/GT/GGTpp_ppBSE.f90
View File

@ -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)

1
src/GT/GGTpp_ppBSE_static_kernel_D.f90
View File

@ -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

63
src/GT/RGT_Tmatrix.f90
View File

@ -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

45
src/GT/RGTpp_ppBSE.f90
View File

@ -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(:,:)

7
src/GT/RGTpp_ppBSE_static_kernel_D.f90
View File

@ -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