From ac0b9336dfb495129ad714053ab5f769a9109115 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 19 Jul 2023 11:03:55 +0200 Subject: [PATCH] clean up LR in GW mess --- src/CC/CCD.f90 | 13 +- src/CC/CCSD.f90 | 13 +- src/CC/rCCD.f90 | 17 +-- src/GF/GF2_phBSE2.f90 | 48 ++++--- src/GF/GF2_phBSE2_dynamic_perturbation.f90 | 33 +++-- src/GW/Bethe_Salpeter_A_matrix.f90 | 53 -------- src/GW/Bethe_Salpeter_B_matrix.f90 | 50 -------- src/GW/G0W0.f90 | 2 +- src/GW/GW_phACFDT.f90 | 47 ++++--- src/GW/GW_phBSE.f90 | 81 ++++++++---- ...mic.f90 => GW_phBSE2_dynamic_kernel_A.f90} | 4 +- ...l_KA.f90 => GW_phBSE2_static_kernel_A.f90} | 4 +- ...l_KB.f90 => GW_phBSE2_static_kernel_B.f90} | 4 +- ...ynamic.f90 => GW_phBSE_dynamic_kernel.f90} | 4 +- ...amic.f90 => GW_phBSE_dynamic_kernel_A.f90} | 4 +- ...amic.f90 => GW_phBSE_dynamic_kernel_Z.f90} | 5 +- ...mic.f90 => GW_phBSE_dynamic_kernel_ZA.f90} | 4 +- ....f90 => GW_phBSE_dynamic_perturbation.f90} | 19 ++- ..._phBSE_dynamic_perturbation_iterative.f90} | 14 +- ...el_KA.f90 => GW_phBSE_static_kernel_A.f90} | 18 ++- ...el_KB.f90 => GW_phBSE_static_kernel_B.f90} | 16 +-- src/GW/soBSE.f90 | 121 ------------------ src/GW/static_screening_WA_so.f90 | 53 -------- src/GW/static_screening_WB_so.f90 | 58 --------- src/QuAcK/QuAcK.f90 | 6 +- 25 files changed, 186 insertions(+), 505 deletions(-) delete mode 100644 src/GW/Bethe_Salpeter_A_matrix.f90 delete mode 100644 src/GW/Bethe_Salpeter_B_matrix.f90 rename src/GW/{BSE2_GW_A_matrix_dynamic.f90 => GW_phBSE2_dynamic_kernel_A.f90} (95%) rename src/GW/{BSE2_static_kernel_KA.f90 => GW_phBSE2_static_kernel_A.f90} (95%) rename src/GW/{BSE2_static_kernel_KB.f90 => GW_phBSE2_static_kernel_B.f90} (95%) rename src/GW/{Bethe_Salpeter_AB_matrix_dynamic.f90 => GW_phBSE_dynamic_kernel.f90} (95%) rename src/GW/{Bethe_Salpeter_A_matrix_dynamic.f90 => GW_phBSE_dynamic_kernel_A.f90} (93%) rename src/GW/{Bethe_Salpeter_ZAB_matrix_dynamic.f90 => GW_phBSE_dynamic_kernel_Z.f90} (93%) rename src/GW/{Bethe_Salpeter_ZA_matrix_dynamic.f90 => GW_phBSE_dynamic_kernel_ZA.f90} (90%) rename src/GW/{Bethe_Salpeter_dynamic_perturbation.f90 => GW_phBSE_dynamic_perturbation.f90} (82%) rename src/GW/{Bethe_Salpeter_dynamic_perturbation_iterative.f90 => GW_phBSE_dynamic_perturbation_iterative.f90} (87%) rename src/GW/{BSE_static_kernel_KA.f90 => GW_phBSE_static_kernel_A.f90} (65%) rename src/GW/{BSE_static_kernel_KB.f90 => GW_phBSE_static_kernel_B.f90} (73%) delete mode 100644 src/GW/soBSE.f90 delete mode 100644 src/GW/static_screening_WA_so.f90 delete mode 100644 src/GW/static_screening_WB_so.f90 diff --git a/src/CC/CCD.f90 b/src/CC/CCD.f90 index e418dac..8200dc7 100644 --- a/src/CC/CCD.f90 +++ b/src/CC/CCD.f90 @@ -1,4 +1,4 @@ -subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF) +subroutine CCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF) ! CCD module @@ -6,7 +6,6 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER ! Input variables - logical,intent(in) :: BSE integer,intent(in) :: maxSCF integer,intent(in) :: max_diis double precision,intent(in) :: thresh @@ -93,15 +92,7 @@ subroutine CCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ER allocate(dbERI(nBas,nBas,nBas,nBas)) - if(BSE) then - - call static_screening(nBas,nC,nO,nV,nR,seHF,sERI,dbERI) - - else - - call antisymmetrize_ERI(2,nBas,sERI,dbERI) - - end if + call antisymmetrize_ERI(2,nBas,sERI,dbERI) deallocate(sERI) diff --git a/src/CC/CCSD.f90 b/src/CC/CCSD.f90 index 0324992..7707659 100644 --- a/src/CC/CCSD.f90 +++ b/src/CC/CCSD.f90 @@ -1,4 +1,4 @@ -subroutine CCSD(BSE,maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF) +subroutine CCSD(maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF) ! CCSD module @@ -6,7 +6,6 @@ subroutine CCSD(BSE,maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ER ! Input variables - logical,intent(in) :: BSE integer,intent(in) :: maxSCF integer,intent(in) :: max_diis double precision,intent(in) :: thresh @@ -105,15 +104,7 @@ subroutine CCSD(BSE,maxSCF,thresh,max_diis,doCCSDT,nBasin,nCin,nOin,nVin,nRin,ER allocate(dbERI(nBas,nBas,nBas,nBas)) - if(BSE) then - - call static_screening(nBas,nC,nO,nV,nR,seHF,sERI,dbERI) - - else - - call antisymmetrize_ERI(2,nBas,sERI,dbERI) - - end if + call antisymmetrize_ERI(2,nBas,sERI,dbERI) deallocate(sERI) diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90 index f30ebd8..766c4bf 100644 --- a/src/CC/rCCD.f90 +++ b/src/CC/rCCD.f90 @@ -1,4 +1,4 @@ -subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF,eGW) +subroutine rCCD(maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,ERHF,eHF,eGW) ! Ring CCD module @@ -6,7 +6,6 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E ! Input variables - logical,intent(in) :: BSE integer,intent(in) :: maxSCF integer,intent(in) :: max_diis double precision,intent(in) :: thresh @@ -77,23 +76,11 @@ subroutine rCCD(BSE,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,E call spatial_to_spin_MO_energy(nBasin,eGW,nBas,seGW) call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI) -! call soRPAx(.false.,nBas,nC,nO,nV,nR,nO*nV,ENuc,ERHF,sERI,seHF) - - call soBSE(.false.,.false.,0.0,nBas,nC,nO,nV,nR,nO*nV,sERI,seHF,seGW) - ! Antysymmetrize ERIs allocate(dbERI(nBas,nBas,nBas,nBas)) - if(BSE) then - - call static_screening(nBas,nC,nO,nV,nR,seHF,sERI,dbERI) - - else - - call antisymmetrize_ERI(2,nBas,sERI,dbERI) - - end if + call antisymmetrize_ERI(2,nBas,sERI,dbERI) deallocate(sERI) diff --git a/src/GF/GF2_phBSE2.f90 b/src/GF/GF2_phBSE2.f90 index 5b6dee7..3cb031a 100644 --- a/src/GF/GF2_phBSE2.f90 +++ b/src/GF/GF2_phBSE2.f90 @@ -28,13 +28,15 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n ! Local variables + logical :: dRPA = .false. integer :: ispin double precision,allocatable :: OmBSE(:) double precision,allocatable :: XpY(:,:) double precision,allocatable :: XmY(:,:) - double precision :: rho double precision,allocatable :: A_sta(:,:) double precision,allocatable :: B_sta(:,:) + double precision,allocatable :: KA_sta(:,:) + double precision,allocatable :: KB_sta(:,:) ! Output variables @@ -42,7 +44,8 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n ! Memory allocation - allocate(OmBSE(nS),XpY(nS,nS),XmY(nS,nS),A_sta(nS,nS),B_sta(nS,nS)) + 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)) !------------------- ! Singlet manifold @@ -53,18 +56,23 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n ispin = 1 EcBSE(ispin) = 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) + ! Compute static kernel - call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta) - if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta) + call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) + if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) - ! Compute BSE2 excitation energies + A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) + if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) - call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta,-B_sta,EcBSE(ispin),OmBSE,XpY,XmY) - call print_excitation('BSE2 ',ispin,nS,OmBSE) + ! Compute phBSE2@GF2 excitation energies + + call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY) + call print_excitation('phBSE2@GF2 ',ispin,nS,OmBSE) call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) - ! Compute dynamic correction for BSE via perturbation theory if(dBSE) then @@ -72,10 +80,10 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n if(evDyn) then call GF2_phBSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - A_sta,B_sta,OmBSE,XpY,XmY) + KA_sta,KB_sta,OmBSE,XpY,XmY) else - call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY) + call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) end if @@ -92,15 +100,21 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n ispin = 2 EcBSE(ispin) = 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) + ! Compute static kernel - call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta) - if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta) + call GF2_phBSE2_static_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KA_sta) + if(.not.TDA) call GF2_phBSE2_static_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,KB_sta) - ! Compute BSE2 excitation energies + A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) + if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) - call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta,-B_sta,EcBSE(ispin),OmBSE,XpY,XmY) - call print_excitation('BSE2 ',ispin,nS,OmBSE) + ! Compute phBSE2@GF2 excitation energies + + call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY,XmY) + call print_excitation('phBSE2@GF2 ',ispin,nS,OmBSE) call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY,XmY) ! Compute dynamic correction for BSE via perturbation theory @@ -110,10 +124,10 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n if(evDyn) then call GF2_phBSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - A_sta,B_sta,OmBSE,XpY,XmY) + KA_sta,KB_sta,OmBSE,XpY,XmY) else - call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY) + call GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) end if diff --git a/src/GF/GF2_phBSE2_dynamic_perturbation.f90 b/src/GF/GF2_phBSE2_dynamic_perturbation.f90 index 4b3ba2b..9486868 100644 --- a/src/GF/GF2_phBSE2_dynamic_perturbation.f90 +++ b/src/GF/GF2_phBSE2_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY) +subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,KA_sta,KB_sta,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -21,8 +21,8 @@ subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ER double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF(nBas) - double precision,intent(in) :: A_sta(nS,nS) - double precision,intent(in) :: B_sta(nS,nS) + double precision,intent(in) :: KA_sta(nS,nS) + double precision,intent(in) :: KB_sta(nS,nS) double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) @@ -38,18 +38,17 @@ subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ER double precision,allocatable :: X(:) double precision,allocatable :: Y(:) - double precision,allocatable :: Ap_dyn(:,:) - double precision,allocatable :: Am_dyn(:,:) + double precision,allocatable :: KAp_dyn(:,:) + double precision,allocatable :: KAm_dyn(:,:) double precision,allocatable :: ZAp_dyn(:,:) double precision,allocatable :: ZAm_dyn(:,:) - double precision,allocatable :: B_dyn(:,:) + double precision,allocatable :: KB_dyn(:,:) ! Memory allocation - allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) - - if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),B_dyn(nS,nS)) + allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),KAp_dyn(nS,nS),ZAp_dyn(nS,nS)) + allocate(KAm_dyn(nS,nS),ZAm_dyn(nS,nS),KB_dyn(nS,nS)) if(dTDA) then write(*,*) @@ -78,28 +77,28 @@ subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ER ! Resonant part of the BSE correction for dynamical TDA - call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),Ap_dyn,ZAp_dyn) + call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmBSE(ia),KAp_dyn,ZAp_dyn) if(dTDA) then ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) - OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) + OmDyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) else ! Second part of the resonant and anti-resonant part of the BSE correction (frequency independent) - call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),Am_dyn,ZAm_dyn) + call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),KAm_dyn,ZAm_dyn) - call GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn) + call GF2_phBSE2_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)) - OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) & - - dot_product(Y,matmul(Am_dyn - A_sta,Y)) & - + dot_product(X,matmul(B_dyn - B_sta,Y)) & - - dot_product(Y,matmul(B_dyn - B_sta,X)) + OmDyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) & + - dot_product(Y,matmul(KAm_dyn - KA_sta,Y)) & + + dot_product(X,matmul(KB_dyn - KB_sta,Y)) & + - dot_product(Y,matmul(KB_dyn - KB_sta,X)) end if diff --git a/src/GW/Bethe_Salpeter_A_matrix.f90 b/src/GW/Bethe_Salpeter_A_matrix.f90 deleted file mode 100644 index 38de0f2..0000000 --- a/src/GW/Bethe_Salpeter_A_matrix.f90 +++ /dev/null @@ -1,53 +0,0 @@ -subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A_lr) - -! Compute the extra term for Bethe-Salpeter equation for linear response - - 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) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ia,jb,kc - -! Output variables - - double precision,intent(out) :: A_lr(nS,nS) - - jb = 0 -!$omp parallel do default(private) shared(A_lr,ERI,Omega,rho,nO,nBas,nS,chi,eps,eta,nC,nR,lambda) - 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) - - chi = 0d0 - do kc=1,nS - eps = Omega(kc)**2 + eta**2 -! chi = chi + lambda*rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps - chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps - enddo - - A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI(i,b,j,a) + 4d0*lambda*chi - - enddo - enddo - enddo - enddo -!$omp end parallel do - -end subroutine Bethe_Salpeter_A_matrix diff --git a/src/GW/Bethe_Salpeter_B_matrix.f90 b/src/GW/Bethe_Salpeter_B_matrix.f90 deleted file mode 100644 index 97f60d5..0000000 --- a/src/GW/Bethe_Salpeter_B_matrix.f90 +++ /dev/null @@ -1,50 +0,0 @@ -subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,B_lr) - -! Compute the extra term for Bethe-Salpeter equation for linear response - - 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) :: OmRPA(nS) - double precision,intent(in) :: rho_RPA(nBas,nBas,nS) - -! Local variables - - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ia,jb,kc - -! Output variables - - double precision,intent(out) :: B_lr(nS,nS) - - 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 - - chi = 0d0 - do kc=1,nS - eps = OmRPA(kc)**2 + eta**2 -! chi = chi + lambda*rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/eps - chi = chi + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/eps - enddo - - B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,j,b,a) + 4d0*lambda*chi - - enddo - enddo - enddo - enddo - -end subroutine Bethe_Salpeter_B_matrix diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 9e00a42..ff40f9c 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -107,7 +107,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT ! Compute screening ! !-------------------! - call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO,Aph) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO,Aph) if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) diff --git a/src/GW/GW_phACFDT.f90 b/src/GW/GW_phACFDT.f90 index dc43e32..697e8a8 100644 --- a/src/GW/GW_phACFDT.f90 +++ b/src/GW/GW_phACFDT.f90 @@ -30,6 +30,7 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e ! Local variables + logical :: dRPA_W = .true. integer :: ispin integer :: isp_W integer :: iAC @@ -57,8 +58,9 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e ! Memory allocation allocate(Ec(nAC,nspin)) - allocate(KA(nS,nS),KB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) - allocate(Om(nS),XpY(nS,nS),XmY(nS,nS)) + allocate(Aph(nS,nS),KA(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & + rho_RPA(nBas,nBas,nS),Om(nS),XpY(nS,nS),XmY(nS,nS)) + if(.not.TDA) allocate(Aph(nS,nS),KB(nS,nS)) ! Antisymmetrized kernel version @@ -78,11 +80,14 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e isp_W = 1 EcRPA = 0d0 - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + 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 GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB) + call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA) + call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB) ! Singlet manifold @@ -105,15 +110,21 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e if(doXBS) then - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) + + call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) + + call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) + call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) end if - call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY) + Aph(:,:) = Aph(:,:) + KA(:,:) + if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:) + + call phLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY) call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) @@ -153,15 +164,21 @@ subroutine GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,e if(doXBS) then - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Aph) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,lambda,ERI,Bph) + + call phLR(TDA_W,nS,Aph,Bph,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) + + call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KA) + call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,KB) end if - call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,KA,KB,EcAC(ispin),Om,XpY,XmY) + Aph(:,:) = Aph(:,:) + KA(:,:) + if(.not.TDA) Bph(:,:) = Bph(:,:) + KB(:,:) + + call phLR(TDA,nS,Aph,Bph,EcAC(ispin),Om,XpY,XmY) call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,Ec(iAC,ispin)) diff --git a/src/GW/GW_phBSE.f90 b/src/GW/GW_phBSE.f90 index 7c96331..fd767da 100644 --- a/src/GW/GW_phBSE.f90 +++ b/src/GW/GW_phBSE.f90 @@ -30,6 +30,9 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n ! Local variables + logical :: dRPA = .false. + logical :: dRPA_W = .true. + integer :: ispin integer :: isp_W @@ -43,6 +46,9 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n double precision,allocatable :: XpY_BSE(:,:) double precision,allocatable :: XmY_BSE(:,:) + double precision,allocatable :: A_sta(:,:) + double precision,allocatable :: B_sta(:,:) + double precision,allocatable :: KA_sta(:,:) double precision,allocatable :: KB_sta(:,:) @@ -57,7 +63,8 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n ! Memory allocation allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & - KA_sta(nS,nS),KB_sta(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) + A_sta(nS,nS),KA_sta(nS,nS),KA2_sta(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) + allocate(B_sta(nS,nS),KB_sta(nS,nS),KB2_sta(nS,nS)) !--------------------------------- ! Compute (singlet) RPA screening @@ -66,11 +73,14 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n isp_W = 1 EcRPA = 0d0 - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phLR_A(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,A_sta) + if(.not.TDA_W) call phLR_B(isp_W,dRPA_W,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) + + call phLR(TDA_W,nS,A_sta,B_sta,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) - call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) + call GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KA_sta) + call GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,KB_sta) !------------------- ! Singlet manifold @@ -81,27 +91,39 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n ispin = 1 EcBSE(ispin) = 0d0 + ! Compute BSE excitation energies + + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,A_sta) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) + + A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) + if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) + ! Second-order BSE static kernel - allocate(W(nBas,nBas,nBas,nBas),KA2_sta(nS,nS),KB2_sta(nS,nS)) - KA2_sta(:,:) = 0d0 - KB2_sta(:,:) = 0d0 - if(BSE2) then - write(*,*) '*** Second-order BSE static kernel activated! ***' - call static_kernel_W(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W) - call BSE2_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KA2_sta) + allocate(W(nBas,nBas,nBas,nBas)) - if(.not.TDA) call BSE2_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KB2_sta) + write(*,*) + write(*,*) '*** Second-order BSE static kernel activated! ***' + write(*,*) + + call static_kernel_W(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W) + call GW_phBSE2_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KA2_sta) + + if(.not.TDA) call GW_phBSE2_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KB2_sta) + + deallocate(W) + + A_sta(:,:) = A_sta(:,:) + KA2_sta(:,:) + if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB2_sta(:,:) end if - ! Compute BSE excitation energies + call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) - call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,KA_sta-KA2_sta,KB_sta-KB2_sta, & - EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) - call print_excitation('BSE@GW ',ispin,nS,OmBSE) + call print_excitation('phBSE@GW ',ispin,nS,OmBSE) call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) !------------------------------------------------- @@ -114,12 +136,12 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & - OmBSE,XpY_BSE,XmY_BSE) + call GW_phBSE_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & + OmBSE,XpY_BSE,XmY_BSE) else - call Bethe_Salpeter_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, & - OmBSE,XpY_BSE,XmY_BSE,W,KA2_sta) + call GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, & + OmBSE,XpY_BSE,XmY_BSE,W,KA2_sta) end if end if @@ -137,8 +159,15 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n ! Compute BSE excitation energies - call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,KA_sta,KB_sta,EcBSE,OmBSE,XpY_BSE,XmY_BSE) - call print_excitation('BSE@GW ',ispin,nS,OmBSE) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,A_sta) + if(.not.TDA) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI,B_sta) + + A_sta(:,:) = A_sta(:,:) + KA_sta(:,:) + if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB_sta(:,:) + + call phLR(TDA,nS,A_sta,B_sta,EcBSE(ispin),OmBSE,XpY_BSE,XmY_BSE) + + call print_excitation('phBSE@GW ',ispin,nS,OmBSE) call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) !------------------------------------------------- @@ -151,12 +180,12 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & - OmBSE,XpY_BSE,XmY_BSE) + call GW_phBSE_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & + OmBSE,XpY_BSE,XmY_BSE) else - call Bethe_Salpeter_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, & - OmBSE,XpY_BSE,XmY_BSE,W,KA2_sta) + call GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int,OmRPA,rho_RPA, & + OmBSE,XpY_BSE,XmY_BSE,W,KA2_sta) end if end if diff --git a/src/GW/BSE2_GW_A_matrix_dynamic.f90 b/src/GW/GW_phBSE2_dynamic_kernel_A.f90 similarity index 95% rename from src/GW/BSE2_GW_A_matrix_dynamic.f90 rename to src/GW/GW_phBSE2_dynamic_kernel_A.f90 index b4e1051..609396e 100644 --- a/src/GW/BSE2_GW_A_matrix_dynamic.f90 +++ b/src/GW/GW_phBSE2_dynamic_kernel_A.f90 @@ -1,4 +1,4 @@ -subroutine BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,ZA_dyn) +subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,ZA_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -92,4 +92,4 @@ subroutine BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn,ZA enddo !$omp end parallel do -end subroutine BSE2_GW_A_matrix_dynamic +end subroutine diff --git a/src/GW/BSE2_static_kernel_KA.f90 b/src/GW/GW_phBSE2_static_kernel_A.f90 similarity index 95% rename from src/GW/BSE2_static_kernel_KA.f90 rename to src/GW/GW_phBSE2_static_kernel_A.f90 index 5780c68..709c542 100644 --- a/src/GW/BSE2_static_kernel_KA.f90 +++ b/src/GW/GW_phBSE2_static_kernel_A.f90 @@ -1,4 +1,4 @@ -subroutine BSE2_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KA2_sta) +subroutine GW_phBSE2_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KA2_sta) ! Compute the second-order static BSE kernel for the resonant block (only for singlets!) @@ -91,4 +91,4 @@ subroutine BSE2_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KA2_sta) !$omp end parallel do -end subroutine BSE2_static_kernel_KA +end subroutine diff --git a/src/GW/BSE2_static_kernel_KB.f90 b/src/GW/GW_phBSE2_static_kernel_B.f90 similarity index 95% rename from src/GW/BSE2_static_kernel_KB.f90 rename to src/GW/GW_phBSE2_static_kernel_B.f90 index 04a4e54..dfe2618 100644 --- a/src/GW/BSE2_static_kernel_KB.f90 +++ b/src/GW/GW_phBSE2_static_kernel_B.f90 @@ -1,4 +1,4 @@ -subroutine BSE2_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KB2_sta) +subroutine GW_phBSE2_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KB2_sta) ! Compute the second-order static BSE kernel for the coupling block (only for singlets!) @@ -91,4 +91,4 @@ subroutine BSE2_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,W,KB2_sta) !$omp end parallel do -end subroutine BSE2_static_kernel_KB +end subroutine diff --git a/src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/GW/GW_phBSE_dynamic_kernel.f90 similarity index 95% rename from src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90 rename to src/GW/GW_phBSE_dynamic_kernel.f90 index 69b026b..516049b 100644 --- a/src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rhO_RPA,OmBSE,Ap,Am,Bp,Bm) +subroutine GW_phBSE_dynamic_kernel(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rhO_RPA,OmBSE,Ap,Am,Bp,Bm) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -115,4 +115,4 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O enddo enddo -end subroutine Bethe_Salpeter_AB_matrix_dynamic +end subroutine diff --git a/src/GW/Bethe_Salpeter_A_matrix_dynamic.f90 b/src/GW/GW_phBSE_dynamic_kernel_A.f90 similarity index 93% rename from src/GW/Bethe_Salpeter_A_matrix_dynamic.f90 rename to src/GW/GW_phBSE_dynamic_kernel_A.f90 index a804a13..81ee6ba 100644 --- a/src/GW/Bethe_Salpeter_A_matrix_dynamic.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel_A.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn,ZA_dyn) +subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,A_dyn,ZA_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -92,4 +92,4 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om !$omp end parallel do -end subroutine Bethe_Salpeter_A_matrix_dynamic +end subroutine diff --git a/src/GW/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/GW/GW_phBSE_dynamic_kernel_Z.f90 similarity index 93% rename from src/GW/Bethe_Salpeter_ZAB_matrix_dynamic.f90 rename to src/GW/GW_phBSE_dynamic_kernel_Z.f90 index c7550da..5d7a7e3 100644 --- a/src/GW/Bethe_Salpeter_ZAB_matrix_dynamic.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel_Z.f90 @@ -1,5 +1,4 @@ -subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE, & - ZAp,ZAm,ZBp,ZBm) +subroutine GW_phBSE_dynamic_kernel_Z(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZAp,ZAm,ZBp,ZBm) ! Compute the dynamic part of the renormalization for the Bethe-Salpeter equation matrices @@ -99,4 +98,4 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, enddo enddo -end subroutine Bethe_Salpeter_ZAB_matrix_dynamic +end subroutine diff --git a/src/GW/Bethe_Salpeter_ZA_matrix_dynamic.f90 b/src/GW/GW_phBSE_dynamic_kernel_ZA.f90 similarity index 90% rename from src/GW/Bethe_Salpeter_ZA_matrix_dynamic.f90 rename to src/GW/GW_phBSE_dynamic_kernel_ZA.f90 index e15677d..1ef836a 100644 --- a/src/GW/Bethe_Salpeter_ZA_matrix_dynamic.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel_ZA.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZA_dyn) +subroutine GW_phBSE_dynamic_kernel_ZA(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rho_RPA,OmBSE,ZA_dyn) ! Compute the dynamic part of the Bethe-Salpeter equation matrices @@ -63,4 +63,4 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O enddo enddo -end subroutine Bethe_Salpeter_ZA_matrix_dynamic +end subroutine diff --git a/src/GW/Bethe_Salpeter_dynamic_perturbation.f90 b/src/GW/GW_phBSE_dynamic_perturbation.f90 similarity index 82% rename from src/GW/Bethe_Salpeter_dynamic_perturbation.f90 rename to src/GW/GW_phBSE_dynamic_perturbation.f90 index e12df44..1285c78 100644 --- a/src/GW/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/GW/GW_phBSE_dynamic_perturbation.f90 @@ -1,5 +1,5 @@ -subroutine Bethe_Salpeter_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW, & -dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY,W,A_stat) +subroutine GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,dipole_int, & + OmRPA,rho_RPA,OmBSE,XpY,XmY,W,A_stat) ! Compute dynamical effects via perturbation theory for BSE @@ -57,8 +57,7 @@ dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY,W,A_stat) maxS = min(nS,maxS) allocate(OmDyn(maxS),ZDyn(maxS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) - - if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS)) + allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS)) if(dTDA) then write(*,*) @@ -87,9 +86,9 @@ dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY,W,A_stat) ! Resonant part of the BSE correction for dynamical TDA - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn,ZAp_dyn) + call GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn,ZAp_dyn) -if(BSE2) call BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),Ap_dyn,ZAp_dyn,W) + if(BSE2) call GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),Ap_dyn,ZAp_dyn,W) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) OmDyn(ia) = dot_product(X,matmul( Ap_dyn - A_stat,X)) @@ -98,13 +97,11 @@ if(BSE2) call BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),A ! Resonant and anti-resonant part of the BSE correction - call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), & - Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) + call GW_phBSE_dynamic_kernel(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) ! Renormalization factor of the resonant and anti-resonant parts - call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia), & - ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn) + call GW_phBSE_dynamic_kernel_Z(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmBSE(ia),ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & - dot_product(Y,matmul(ZAm_dyn,Y)) & @@ -128,4 +125,4 @@ if(BSE2) call BSE2_GW_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),A write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) -end subroutine Bethe_Salpeter_dynamic_perturbation +end subroutine diff --git a/src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/GW/GW_phBSE_dynamic_perturbation_iterative.f90 similarity index 87% rename from src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90 rename to src/GW/GW_phBSE_dynamic_perturbation_iterative.f90 index 26ba94d..4d39011 100644 --- a/src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90 +++ b/src/GW/GW_phBSE_dynamic_perturbation_iterative.f90 @@ -1,5 +1,4 @@ -subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int, & - OmRPA,rho_RPA,OmBSE,XpY,XmY) +subroutine GW_phBSE_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA,OmBSE,XpY,XmY) ! Compute self-consistently the dynamical effects via perturbation theory for BSE @@ -57,9 +56,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, maxS = min(nS,maxS) allocate(OmDyn(maxS),OmOld(maxS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) - - if(.not.dTDA) & - allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS)) + allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS)) if(dTDA) then write(*,*) @@ -101,7 +98,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, ! Resonant part of the BSE correction - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn,ZAp_dyn) + call GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn,ZAp_dyn) OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) @@ -109,8 +106,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, ! Anti-resonant part of the BSE correction - call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia), & - Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) + call GW_phBSE_dynamic_kernel(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,OmOld(ia),Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & @@ -146,4 +142,4 @@ subroutine Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV, endif -end subroutine Bethe_Salpeter_dynamic_perturbation_iterative +end subroutine diff --git a/src/GW/BSE_static_kernel_KA.f90 b/src/GW/GW_phBSE_static_kernel_A.f90 similarity index 65% rename from src/GW/BSE_static_kernel_KA.f90 rename to src/GW/GW_phBSE_static_kernel_A.f90 index 725c73e..2747173 100644 --- a/src/GW/BSE_static_kernel_KA.f90 +++ b/src/GW/GW_phBSE_static_kernel_A.f90 @@ -1,4 +1,4 @@ -subroutine BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA) +subroutine GW_phBSE_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,KA) ! Compute the BSE static kernel for the resonant block @@ -11,7 +11,7 @@ subroutine BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA) double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -22,11 +22,9 @@ subroutine BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA) ! Output variables - double precision,intent(out) :: WA(nS,nS) + double precision,intent(out) :: KA(nS,nS) -! Initialize - - WA(:,:) = 0d0 +! Compute static kernel ia = 0 do i=nC+1,nO @@ -39,15 +37,15 @@ subroutine BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA) chi = 0d0 do kc=1,nS - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,j,kc)*rho(a,b,kc)*Om(kc)/eps enddo - WA(ia,jb) = WA(ia,jb) + lambda*ERI(i,b,j,a) - 4d0*lambda*chi + KA(ia,jb) = 4d0*lambda*chi enddo enddo enddo enddo -end subroutine BSE_static_kernel_KA +end subroutine diff --git a/src/GW/BSE_static_kernel_KB.f90 b/src/GW/GW_phBSE_static_kernel_B.f90 similarity index 73% rename from src/GW/BSE_static_kernel_KB.f90 rename to src/GW/GW_phBSE_static_kernel_B.f90 index deb7336..3554de1 100644 --- a/src/GW/BSE_static_kernel_KB.f90 +++ b/src/GW/GW_phBSE_static_kernel_B.f90 @@ -1,4 +1,4 @@ -subroutine BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB) +subroutine GW_phBSE_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Om,rho,KB) ! Compute the BSE static kernel for the coupling block @@ -16,7 +16,7 @@ subroutine BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB) double precision,intent(in) :: eta double precision,intent(in) :: lambda double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: Om(nS) double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables @@ -29,9 +29,7 @@ subroutine BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB) double precision,intent(out) :: KB(nS,nS) -! Initialize - - KB(:,:) = 0d0 +! Compute static kernel ia = 0 do i=nC+1,nO @@ -44,15 +42,15 @@ subroutine BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB) chi = 0d0 do kc=1,nS - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps + eps = Om(kc)**2 + eta**2 + chi = chi + rho(i,b,kc)*rho(a,j,kc)*Om(kc)/eps enddo - KB(ia,jb) = KB(ia,jb) + lambda*ERI(i,j,b,a) - 4d0*lambda*chi + KB(ia,jb) = 4d0*lambda*chi enddo enddo enddo enddo -end subroutine BSE_static_kernel_KB +end subroutine diff --git a/src/GW/soBSE.f90 b/src/GW/soBSE.f90 deleted file mode 100644 index fd095f5..0000000 --- a/src/GW/soBSE.f90 +++ /dev/null @@ -1,121 +0,0 @@ -subroutine soBSE(TDA_W,TDA,eta,nBas,nC,nO,nV,nR,nS,ERI,eW,eGW) - -! Compute the Bethe-Salpeter excitation energies - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: TDA_W - logical,intent(in) :: TDA - - 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) - -! Local variables - - integer :: ispin - integer :: isp_W - - 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 :: OmBSE(:) - double precision,allocatable :: XpY_BSE(:,:) - double precision,allocatable :: XmY_BSE(:,:) - - double precision,allocatable :: WA_sta(:,:) - double precision,allocatable :: WB_sta(:,:) - - double precision,allocatable :: X(:,:) - double precision,allocatable :: Y(:,:) - double precision,allocatable :: Xinv(:,:) - double precision,allocatable :: t(:,:,:,:) - - integer ::i,a,j,b,k,c,ia,jb,kc - - double precision :: EcBSE - -! Memory allocation - - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & - WA_sta(nS,nS),WB_sta(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) - -!--------------------------------- -! Compute (singlet) RPA screening -!--------------------------------- - - isp_W = 3 - EcRPA = 0d0 - - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - - call static_screening_WA_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA_sta) - call static_screening_WB_so(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB_sta) - - - ispin = 3 - EcBSE = 0d0 - - ! Compute BSE excitation energies - - call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, & - EcBSE,OmBSE,XpY_BSE,XmY_BSE) - call print_excitation('soBSE@GW ',ispin,nS,OmBSE) - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@G0W0 correlation energy =',0.5d0*EcBSE,' au' - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - -! allocate(X(nS,nS),Y(nS,nS),Xinv(nS,nS),t(nO,nO,nV,nV)) - -! X(:,:) = transpose(0.5d0*(XpY_BSE(:,:) + XmY_BSE(:,:))) -! Y(:,:) = transpose(0.5d0*(XpY_BSE(:,:) - XmY_BSE(:,:))) - -! call matout(nS,nS,matmul(X,transpose(X))-matmul(Y,transpose(Y))) - -! call inverse_matrix(nS,X,Xinv) - -! t = 0d0 -! ia = 0 -! do i=1,nO -! do a=1,nV -! ia = ia + 1 - -! jb = 0 -! do j=1,nO -! do b=1,nV -! jb = jb + 1 - -! kc = 0 -! do k=1,nO -! do c=1,nV -! kc = kc + 1 - -! t(i,j,a,b) = t(i,j,a,b) + Y(ia,kc)*Xinv(kc,jb) - -! end do -! end do -! end do -! end do -! end do -! end do - -! call matout(nO*nO,nV*nV,t) - -end subroutine diff --git a/src/GW/static_screening_WA_so.f90 b/src/GW/static_screening_WA_so.f90 deleted file mode 100644 index 951a013..0000000 --- a/src/GW/static_screening_WA_so.f90 +++ /dev/null @@ -1,53 +0,0 @@ -subroutine static_screening_WA_so(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA) - -! Compute the OOVV block of the static screening W for the resonant block - - 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) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ia,jb,kc - -! Output variables - - double precision,intent(out) :: WA(nS,nS) - -! Initialize - - WA(:,:) = 0d0 - - 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 - - chi = 0d0 - do kc=1,nS - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps - enddo - - WA(ia,jb) = WA(ia,jb) + lambda*ERI(i,b,j,a) - 2d0*lambda*chi - - enddo - enddo - enddo - enddo - -end subroutine static_screening_WA_so diff --git a/src/GW/static_screening_WB_so.f90 b/src/GW/static_screening_WB_so.f90 deleted file mode 100644 index 849f489..0000000 --- a/src/GW/static_screening_WB_so.f90 +++ /dev/null @@ -1,58 +0,0 @@ -subroutine static_screening_WB_so(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WB) - -! Compute the static screening W for the coupling block - - 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 - double precision,intent(in) :: eta - double precision,intent(in) :: lambda - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ia,jb,kc - -! Output variables - - double precision,intent(out) :: WB(nS,nS) - -! Initialize - - WB(:,:) = 0d0 - - 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 - - chi = 0d0 - do kc=1,nS - eps = Omega(kc)**2 + eta**2 - chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps - enddo - - WB(ia,jb) = WB(ia,jb) + lambda*ERI(i,j,b,a) - 2d0*lambda*chi - - enddo - enddo - enddo - enddo - -end subroutine static_screening_WB_so diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 64a8033..ef566ed 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -384,7 +384,7 @@ program QuAcK if(doCCD) then call cpu_time(start_CC) - call CCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) + call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) call cpu_time(end_CC) t_CC = end_CC - start_CC @@ -419,7 +419,7 @@ program QuAcK if(doCCSD) then call cpu_time(start_CC) - call CCSD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) + call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF) call cpu_time(end_CC) t_CC = end_CC - start_CC @@ -451,7 +451,7 @@ program QuAcK if(do_rCCD) then call cpu_time(start_CC) - call rCCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF,epsHF) + call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF,epsHF) call cpu_time(end_CC) t_CC = end_CC - start_CC