10
1
mirror of https://github.com/pfloos/quack synced 2024-10-20 06:48:15 +02:00

remove redundant code

This commit is contained in:
Pierre-Francois Loos 2024-09-08 16:19:03 +02:00
parent 8c01e47080
commit d5a396200e
17 changed files with 22 additions and 285 deletions

View File

@ -135,39 +135,4 @@ 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 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 if
end subroutine

View File

@ -260,45 +260,5 @@ end if
! end do
! 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(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 if
! Second-order correlation kernel for the block C of the spinorbital manifold
! deallocate(Om_tmp)
end subroutine

View File

@ -135,39 +135,4 @@ 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(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 if
end subroutine

View File

@ -30,9 +30,6 @@ subroutine GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,di
logical :: dRPA = .false.
logical :: dRPA_W = .true.
integer :: ispin
integer :: isp_W
double precision :: EcRPA
double precision,allocatable :: OmRPA(:)
double precision,allocatable :: XpY_RPA(:,:)
@ -65,25 +62,23 @@ subroutine GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,di
! Compute (singlet) RPA screening
!---------------------------------
isp_W = 3
EcRPA = 0d0
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 phGLR_A(dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
if(.not.TDA_W) call phGLR_B(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 phGLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)
call GGW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA)
call GGW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta)
call GGW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta)
ispin = 3
EcBSE = 0d0
! Compute BSE excitation energies
call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph)
if(.not.TDA) 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) call phGLR_B(dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,Bph)
! Second-order BSE static kernel
@ -107,7 +102,7 @@ subroutine GGW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,eta,nBas,nC,nO,nV,nR,nS,ERI,di
Aph(:,:) = Aph(:,:) + KA_sta(:,:)
if(.not.TDA) Bph(:,:) = Bph(:,:) + KB_sta(:,:)
call phLR(TDA,nS,Aph,Bph,EcBSE,OmBSE,XpY_BSE,XmY_BSE)
call phGLR(TDA,nS,Aph,Bph,EcBSE,OmBSE,XpY_BSE,XmY_BSE)
call print_excitation_energies('phBSE@GW@GHF','spinorbital',nS,OmBSE)
call phLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE)

View File

@ -26,8 +26,6 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,
! Local variables
integer :: isp_W
logical :: dRPA = .false.
logical :: dRPA_W = .true.
@ -67,15 +65,13 @@ subroutine GGW_ppBSE(TDA_W,TDA,dBSE,dTDA,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,
! Compute RPA screening !
!-----------------------!
isp_W = 3
EcRPA = 0d0
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,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phGLR_A(dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,eW,ERI,Aph)
if(.not.TDA_W) call phGLR_B(dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph)
call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA)

View File

@ -99,35 +99,4 @@ subroutine RGW_ppBSE_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambd
end if
!---------------!
! SpinOrb block !
!---------------!
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
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 if
end subroutine

View File

@ -348,35 +348,4 @@ subroutine RGW_ppBSE_static_kernel_C(ispin,eta,nBas,nC,nO,nV,nR,nS,nVV,lambda,ER
end if
!---------------!
! SpinOrb block !
!---------------!
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
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 if
end subroutine

View File

@ -98,35 +98,4 @@ subroutine RGW_ppBSE_static_kernel_D(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,lambda,ER
end if
!---------------!
! SpinOrb block !
!---------------!
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
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 if
end subroutine

View File

@ -48,7 +48,6 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
integer,parameter :: maxS = 20
integer :: ia,i,a,mu
integer :: ispin
double precision,allocatable :: Aph(:,:)
double precision,allocatable :: Bph(:,:)
@ -152,10 +151,8 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Stability analysis: Real GHF -> Real GHF
!-------------------------------------------------------------!
ispin = 3
call phLR_A(ispin,.false.,nBas2,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
call phLR_B(ispin,.false.,nBas2,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
call phGLR_A(.false.,nBas2,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph)
call phGLR_B(.false.,nBas2,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph)
AB(:,:) = Aph(:,:) + Bph(:,:)

View File

@ -21,7 +21,6 @@ subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
integer,parameter :: maxS = 20
integer :: ia
integer :: ispin
double precision,allocatable :: A(:,:)
double precision,allocatable :: B(:,:)
@ -36,10 +35,8 @@ subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI)
! Stability analysis: Real GHF -> Real GHF
!-------------------------------------------------------------!
ispin = 3
call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
call phGLR_A(.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,A)
call phGLR_B(.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI,B)
AB(:,:) = A(:,:) + B(:,:)

View File

@ -81,27 +81,4 @@ subroutine phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
end if
! Build A matrix for spin orbitals
if(ispin == 3) then
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-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 if
end subroutine

View File

@ -71,26 +71,4 @@ subroutine phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph)
end if
! Build B matrix for spin orbitals
if(ispin == 3) then
ia = 0
do i=nC+1,nO
do a=nO+1,nBas-nR
ia = ia + 1
jb = 0
do j=nC+1,nO
do b=nO+1,nBas-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 if
end subroutine

View File

@ -49,9 +49,9 @@ subroutine ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp)
end if
! Build the B matrix for the triplet manifold, or alpha-alpha, or in the spin-orbital basis
! Build the B matrix for the triplet or alpha-alpha manifold
if(ispin == 2 .or. ispin == 4) then
if(ispin == 2) then
ab = 0
do a=nO+1,nOrb-nR

View File

@ -106,9 +106,9 @@ subroutine ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp)
end if
! Build C matrix for the triplet manifold, or alpha-alpha block, or in the spin-orbital basis
! Build C matrix for the triplet or alpha-alpha manifold
if(ispin == 2 .or. ispin == 4) then
if(ispin == 2) then
!$OMP PARALLEL &
!$OMP SHARED(Cpp,lambda,ERI,e,eF,nC,nO,nOrb,nR) &
!$OMP PRIVATE(c,d,a,b,ab,cd) &

View File

@ -44,9 +44,9 @@ subroutine ppLR_C_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Cpp)
end if
! Build C matrix for the triplet manifold, or alpha-alpha block, or in the spin-orbital basis
! Build C matrix for the triplet or alpha-alpha manifold
if(ispin == 2 .or. ispin == 4) then
if(ispin == 2) then
ab = 0
do a=nO+1,nBas-nR

View File

@ -56,9 +56,9 @@ subroutine ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp)
end if
! Build the D matrix for the triplet manifold, the alpha-alpha block, or in the spin-orbital basis
! Build the D matrix for the triplet or alpha-alpha manifold
if(ispin == 2 .or. ispin == 4) then
if(ispin == 2) then
ij = 0
do i=nC+1,nO

View File

@ -44,9 +44,9 @@ subroutine ppLR_D_od(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Dpp)
end if
! Build the D matrix for the triplet manifold, the alpha-alpha block, or in the spin-orbital basis
! Build the D matrix for the triplet or alpha-alpha manifold
if(ispin == 2 .or. ispin == 4) then
if(ispin == 2) then
ij = 0
do i=nC+1,nO