diff --git a/input/options b/input/options index a36a2b8..dfa37b1 100644 --- a/input/options +++ b/input/options @@ -13,6 +13,6 @@ # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 256 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS - T T F -# BSE: phBSE phBSE2 ppBSE dBSE dTDA evDyn - T T T T F F + F F T +# BSE: phBSE phBSE2 ppBSE dBSE dTDA + T T T T F diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90 index 9cb7aba..db1d272 100644 --- a/src/GF/G0F2.f90 +++ b/src/GF/G0F2.f90 @@ -1,4 +1,4 @@ -subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize, & +subroutine G0F2(BSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform a one-shot second-order Green function calculation @@ -12,7 +12,6 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet logical,intent(in) :: linearize @@ -95,7 +94,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize if(BSE) then - call GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) + call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) end if diff --git a/src/GF/GF2_phBSE2.f90 b/src/GF/GF2_phBSE2.f90 index 3cb031a..0ebe1d3 100644 --- a/src/GF/GF2_phBSE2.f90 +++ b/src/GF/GF2_phBSE2.f90 @@ -1,4 +1,4 @@ -subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,EcBSE) +subroutine GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,EcBSE) ! Compute the second-order Bethe-Salpeter excitation energies @@ -10,7 +10,6 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet @@ -75,19 +74,8 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n ! Compute dynamic correction for BSE via perturbation theory - if(dBSE) then - - if(evDyn) then - - call GF2_phBSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - KA_sta,KB_sta,OmBSE,XpY,XmY) - else - + if(dBSE) & 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 - - end if end if @@ -119,19 +107,8 @@ subroutine GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,n ! Compute dynamic correction for BSE via perturbation theory - if(dBSE) then - - if(evDyn) then - - call GF2_phBSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - KA_sta,KB_sta,OmBSE,XpY,XmY) - else - + if(dBSE) & 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 - - end if end if diff --git a/src/GF/GF2_phBSE2_dynamic_perturbation.f90 b/src/GF/GF2_phBSE2_dynamic_perturbation.f90 index 9486868..dd91664 100644 --- a/src/GF/GF2_phBSE2_dynamic_perturbation.f90 +++ b/src/GF/GF2_phBSE2_dynamic_perturbation.f90 @@ -43,12 +43,12 @@ subroutine GF2_phBSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ER double precision,allocatable :: ZAp_dyn(:,:) double precision,allocatable :: ZAm_dyn(:,:) - double precision,allocatable :: KB_dyn(:,:) + double precision,allocatable :: KB_dyn(:,:) ! Memory allocation - 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)) + allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),KAp_dyn(nS,nS),ZAp_dyn(nS,nS), & + KAm_dyn(nS,nS),ZAm_dyn(nS,nS),KB_dyn(nS,nS)) if(dTDA) then write(*,*) diff --git a/src/GF/GF2_phBSE2_dynamic_perturbation_iterative.f90 b/src/GF/GF2_phBSE2_dynamic_perturbation_iterative.f90 deleted file mode 100644 index a030d54..0000000 --- a/src/GF/GF2_phBSE2_dynamic_perturbation_iterative.f90 +++ /dev/null @@ -1,155 +0,0 @@ -subroutine 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) - -! Compute self-consistently the dynamical effects via perturbation theory for BSE2 - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: dTDA - integer,intent(in) :: ispin - 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) :: ERI(nBas,nBas,nBas,nBas) - 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) :: OmBSE(nS) - double precision,intent(in) :: XpY(nS,nS) - double precision,intent(in) :: XmY(nS,nS) - -! Local variables - - integer :: ia - - integer,parameter :: maxS = 10 - double precision :: gapGF - - integer :: nSCF - integer :: maxSCF = 10 - double precision :: Conv - double precision :: thresh = 1d-3 - - - double precision,allocatable :: OmDyn(:) - double precision,allocatable :: ZDyn(:) - double precision,allocatable :: OmOld(:) - double precision,allocatable :: X(:) - double precision,allocatable :: Y(:) - - double precision,allocatable :: Ap_dyn(:,:) - double precision,allocatable :: Am_dyn(:,:) - double precision,allocatable :: ZAp_dyn(:,:) - double precision,allocatable :: ZAm_dyn(:,:) - - double precision,allocatable :: B_dyn(:,:) - -! Memory allocation - - allocate(OmDyn(nS),OmOld(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)) - -! Print main components of transition vectors - - call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) - - if(dTDA) then - write(*,*) - write(*,*) '*** dynamical TDA activated ***' - write(*,*) - end if - - gapGF = eGF(nO+1) - eGF(nO) - - Conv = 1d0 - nSCF = 0 - OmOld(:) = OmBSE(:) - - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies ' - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(A57,F10.6,A3)') ' BSE2 neutral excitation must be lower than the GF gap = ',gapGF*HaToeV,' eV' - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) - - do while(Conv > thresh .and. nSCF < maxSCF) - - nSCF = nSCF + 1 - - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(2X,A15,I3)') 'Iteration n.',nSCF - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Convergence (eV)' - write(*,*) '---------------------------------------------------------------------------------------------------' - - do ia=1,min(nS,maxS) - - - X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) - Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - - ! Resonant part of the BSE correction - call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,+OmOld(ia),Ap_dyn,ZAp_dyn) - - if(dTDA) then - - OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) - ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) - - else - - ! Anti-resonant part of the BSE correction - - call GF2_phBSE2_dynamic_kernel_A(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmOld(ia),Am_dyn,ZAm_dyn) - - call GF2_phBSE2_dynamic_kernel_B(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_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)) - - end if - - write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & - ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,(OmBSE(ia) + OmDyn(ia) - OmOld(ia))*HaToeV - - end do - - Conv = maxval(abs(OmBSE(:) + OmDyn(:) - OmOld(:)))*HaToeV - OmOld(:) = OmBSE(:) + OmDyn(:) - - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(2X,A20,1X,F10.6)') ' Convergence = ',Conv - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) - - end do - -! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) - - endif - -end subroutine diff --git a/src/GF/UG0F2.f90 b/src/GF/UG0F2.f90 index 327d8e4..3019e55 100644 --- a/src/GF/UG0F2.f90 +++ b/src/GF/UG0F2.f90 @@ -1,4 +1,4 @@ -subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & +subroutine UG0F2(BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF) ! Perform unrestricted G0W0 calculation @@ -13,7 +13,6 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta, logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip logical,intent(in) :: linearize diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90 index f9131bc..d72e2e0 100644 --- a/src/GF/evGF2.f90 +++ b/src/GF/evGF2.f90 @@ -1,4 +1,4 @@ -subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, & +subroutine evGF2(BSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform eigenvalue self-consistent second-order Green function calculation @@ -12,7 +12,6 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn integer,intent(in) :: maxSCF double precision,intent(in) :: thresh integer,intent(in) :: max_diis @@ -143,7 +142,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, if(BSE) then - call GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) + call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) end if diff --git a/src/GF/evUGF2.f90 b/src/GF/evUGF2.f90 index 6849320..68ad156 100644 --- a/src/GF/evUGF2.f90 +++ b/src/GF/evUGF2.f90 @@ -1,4 +1,4 @@ -subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & +subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip, & eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,cHF,eHF) @@ -18,7 +18,6 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index 8f80d52..8d27d52 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -1,4 +1,4 @@ -subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, & +subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,singlet,triplet, & eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) @@ -16,7 +16,6 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta @@ -262,7 +261,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, if(BSE) then - call GF2_phBSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGF2,EcBSE) + call GF2_phBSE2(TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGF2,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 11e8e8c..327495a 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -1,4 +1,4 @@ -subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize, & +subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,regularize, & nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) @@ -16,7 +16,6 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta diff --git a/src/GT/G0T0eh.f90 b/src/GT/G0T0eh.f90 index 44df49f..2810c3d 100644 --- a/src/GT/G0T0eh.f90 +++ b/src/GT/G0T0eh.f90 @@ -1,4 +1,4 @@ -subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,evDyn,doppBSE, & +subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) @@ -20,7 +20,6 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet logical,intent(in) :: linearize @@ -186,7 +185,7 @@ subroutine G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE, ! if(BSE) then -! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) +! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) ! if(exchange_kernel) then ! diff --git a/src/GT/G0T0pp.f90 b/src/GT/G0T0pp.f90 index 385014a..75a4676 100644 --- a/src/GT/G0T0pp.f90 +++ b/src/GT/G0T0pp.f90 @@ -1,4 +1,4 @@ -subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & +subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,ppBSE,singlet,triplet, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -17,7 +17,6 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet logical,intent(in) :: linearize @@ -249,7 +248,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp if(BSE) then - call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & + call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & ERI_MO,dipole_int,eHF,eGT,EcBSE) @@ -311,7 +310,7 @@ subroutine G0T0pp(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,pp if(ppBSE) then - call GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOaa,nVVaa, & + call GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & ERI_MO,dipole_int,eHF,eGT,EcppBSE) diff --git a/src/GT/GTpp_phBSE.f90 b/src/GT/GTpp_phBSE.f90 index edf9462..e315d56 100644 --- a/src/GT/GTpp_phBSE.f90 +++ b/src/GT/GTpp_phBSE.f90 @@ -1,4 +1,4 @@ -subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & +subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & ERI,dipole_int,eT,eGT,EcBSE) @@ -13,7 +13,6 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet @@ -143,22 +142,12 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n call print_excitation('phBSE@GTpp ',ispin,nS,OmBSE) call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) - if(dBSE) then - - ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) - - if(evDyn) then - - print*, ' Iterative dynamical correction for BSE@GT NYI' + ! Compute dynamic correction for BSE via renormalized perturbation theory - else - + if(dBSE) & call GTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) - end if - - end if end if @@ -183,23 +172,12 @@ subroutine GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n call print_excitation('phBSE@GTpp ',ispin,nS,OmBSE) call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE,XpY_BSE,XmY_BSE) - if(dBSE) then - - ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) - - if(evDyn) then - - print*, ' Iterative dynamical correction for BSE@GT NYI' + ! Compute dynamic correction for BSE via renormalized perturbation theory - else - + if(dBSE) & call GTpp_phBSE_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,Om2ab,Om1aa,Om2aa,rho1ab,rho2ab,rho1aa,rho2aa,eT,eGT, & dipole_int,OmBSE,XpY_BSE,XmY_BSE,TAab,TAaa) - end if - - end if - end if end subroutine diff --git a/src/GT/GTpp_ppBSE.f90 b/src/GT/GTpp_ppBSE.f90 index ad1e552..758770c 100644 --- a/src/GT/GTpp_ppBSE.f90 +++ b/src/GT/GTpp_ppBSE.f90 @@ -1,4 +1,4 @@ -subroutine GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOaa,nVVaa, & +subroutine GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOab,nVVab,nOOaa,nVVaa, & Om1ab,X1ab,Y1ab,Om2ab,X2ab,Y2ab,rho1ab,rho2ab,Om1aa,X1aa,Y1aa,Om2aa,X2aa,Y2aa,rho1aa,rho2aa, & ERI,dipole_int,eT,eGT,EcBSE) @@ -13,7 +13,6 @@ subroutine GTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,n logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet diff --git a/src/GT/UG0T0.f90 b/src/GT/UG0T0.f90 index 56bfad6..a79f411 100644 --- a/src/GT/UG0T0.f90 +++ b/src/GT/UG0T0.f90 @@ -1,4 +1,4 @@ -subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & +subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV, & nR,nS,ENuc,EUHF,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0T0) @@ -18,7 +18,6 @@ subroutine UG0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip logical,intent(in) :: linearize diff --git a/src/GT/evGTeh.f90 b/src/GT/evGTeh.f90 index c9ebdea..7d2ef01 100644 --- a/src/GT/evGTeh.f90 +++ b/src/GT/evGTeh.f90 @@ -1,4 +1,4 @@ -subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,evDyn,doppBSE, & +subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, & cHF,eHF,Vxc) @@ -23,7 +23,6 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: doppBSE logical,intent(in) :: singlet logical,intent(in) :: triplet @@ -237,7 +236,7 @@ subroutine evGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,d ! if(BSE) then -! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE) +! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE) ! if(exchange_kernel) then diff --git a/src/GT/evGTpp.f90 b/src/GT/evGTpp.f90 index a58a4ef..c564808 100644 --- a/src/GT/evGTpp.f90 +++ b/src/GT/evGTpp.f90 @@ -1,4 +1,4 @@ -subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn, & +subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA, & singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) @@ -19,7 +19,6 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta @@ -233,7 +232,7 @@ subroutine evGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T if(BSE) then - call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI_MO,dipole_int,eGT,eGT,EcBSE) diff --git a/src/GT/evUGT.f90 b/src/GT/evUGT.f90 index 212bf2d..78efda5 100644 --- a/src/GT/evUGT.f90 +++ b/src/GT/evUGT.f90 @@ -1,5 +1,5 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & - TDA_T,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,& + TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,ERI_AO,ERI_aaaa, & ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF, & Vxc,eG0T0) @@ -21,7 +21,6 @@ subroutine evUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta diff --git a/src/GT/qsGTeh.f90 b/src/GT/qsGTeh.f90 index c206f87..dc6ebd4 100644 --- a/src/GT/qsGTeh.f90 +++ b/src/GT/qsGTeh.f90 @@ -1,5 +1,5 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA, & - dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & + dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GTeh calculation @@ -21,7 +21,6 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta @@ -292,7 +291,7 @@ subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2, ! if(BSE) then -! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & +! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & ! eGT,eGT,EcBSE) ! if(exchange_kernel) then diff --git a/src/GT/qsGTpp.f90 b/src/GT/qsGTpp.f90 index 0ed0487..0c0d349 100644 --- a/src/GT/qsGTpp.f90 +++ b/src/GT/qsGTpp.f90 @@ -1,5 +1,5 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA, & - dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & + dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GT calculation @@ -20,7 +20,6 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta @@ -355,7 +354,7 @@ subroutine qsGTpp(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T if(BSE) then - call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI_MO,dipole_int_MO,eGT,eGT,EcBSE) diff --git a/src/GT/qsUGT.f90 b/src/GT/qsUGT.f90 index 05a3ce5..137325e 100644 --- a/src/GT/qsUGT.f90 +++ b/src/GT/qsUGT.f90 @@ -1,5 +1,5 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & - TDA_T,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,& + TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& eta,regularize,nBas,nC,nO,nV,nR,nS,nNuc,ZNuc,rNuc,ENuc,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,& ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) @@ -20,7 +20,6 @@ subroutine qsUGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE, & logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index fd1a3e7..8a12163 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -1,4 +1,4 @@ -subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE, & +subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO, & dipole_int,PHF,cHF,eHF,Vxc) @@ -20,7 +20,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet logical,intent(in) :: linearize @@ -182,7 +181,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dT if(dophBSE) then - call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) + call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/GW_phBSE.f90 b/src/GW/GW_phBSE.f90 index fd767da..afc48c6 100644 --- a/src/GW/GW_phBSE.f90 +++ b/src/GW/GW_phBSE.f90 @@ -1,4 +1,4 @@ -subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) +subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -12,7 +12,6 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet @@ -53,8 +52,6 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n double precision,allocatable :: KB_sta(:,:) double precision,allocatable :: W(:,:,:,:) - double precision,allocatable :: KA2_sta(:,:) - double precision,allocatable :: KB2_sta(:,:) ! Output variables @@ -63,8 +60,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), & - 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)) + A_sta(nS,nS),B_sta(nS,nS),KA_sta(nS,nS),KB_sta(nS,nS), & + OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) !--------------------------------- ! Compute (singlet) RPA screening @@ -96,9 +93,6 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n 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 if(BSE2) then @@ -110,17 +104,17 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n 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) + call GW_phBSE2_static_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KA_sta) - if(.not.TDA) call GW_phBSE2_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KB2_sta) + if(.not.TDA) call GW_phBSE2_static_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,W,KB_sta) deallocate(W) - A_sta(:,:) = A_sta(:,:) + KA2_sta(:,:) - if(.not.TDA) B_sta(:,:) = B_sta(:,:) + KB2_sta(:,:) - end if + 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) @@ -130,22 +124,10 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n ! Compute the dynamical screening at the BSE level !------------------------------------------------- - if(dBSE) then + if(dBSE) & + call GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & + OmBSE,XpY_BSE,XmY_BSE,KA_sta,KB_sta) - ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) - - if(evDyn) then - - 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 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 - end if !------------------- @@ -174,21 +156,9 @@ subroutine GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,n ! Compute the dynamical screening at the BSE level !------------------------------------------------- - if(dBSE) then - - ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) - - if(evDyn) then - - 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 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 + if(dBSE) & + call GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int,OmRPA,rho_RPA, & + OmBSE,XpY_BSE,XmY_BSE,KA_sta,KB_sta) end if diff --git a/src/GW/GW_phBSE2_dynamic_kernel_A.f90 b/src/GW/GW_phBSE2_dynamic_kernel_A.f90 index 609396e..1d2e4d8 100644 --- a/src/GW/GW_phBSE2_dynamic_kernel_A.f90 +++ b/src/GW/GW_phBSE2_dynamic_kernel_A.f90 @@ -24,16 +24,10 @@ subroutine GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE,A_dyn, double precision,intent(out) :: A_dyn(nS,nS) double precision,intent(out) :: ZA_dyn(nS,nS) -! Initialization - -! A_dyn(:,:) = 0d0 -! ZA_dyn(:,:) = 0d0 - -! Number of poles taken into account - ! Build dynamic A matrix jb = 0 + !$omp parallel do default(private) shared(A_dyn,ZA_dyn,OmBSE,W,num,dem,eGW,nO,nBas,eta,nC,nR) do j=nC+1,nO do b=nO+1,nBas-nR diff --git a/src/GW/GW_phBSE_dynamic_kernel.f90 b/src/GW/GW_phBSE_dynamic_kernel.f90 deleted file mode 100644 index 516049b..0000000 --- a/src/GW/GW_phBSE_dynamic_kernel.f90 +++ /dev/null @@ -1,118 +0,0 @@ -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 - - 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) :: eGW(nBas) - double precision,intent(in) :: OmRPA(nS) - double precision,intent(in) :: rho_RPA(nBas,nBas,nS) - double precision,intent(in) :: OmBSE - -! Local variables - - integer :: maxS - double precision :: chi_A,chi_B - double precision :: chi_Ap,chi_Bp,eps_Ap,eps_Bp - double precision :: chi_Am,chi_Bm,eps_Am,eps_Bm - integer :: i,j,a,b,ia,jb,kc - -! Output variables - - double precision,intent(out) :: Ap(nS,nS) - double precision,intent(out) :: Am(nS,nS) - - double precision,intent(out) :: Bp(nS,nS) - double precision,intent(out) :: Bm(nS,nS) - -! Initialization - - Ap(:,:) = 0d0 - Am(:,:) = 0d0 - - Bp(:,:) = 0d0 - Bm(:,:) = 0d0 - -! Number of poles taken into account - - maxS = nS - -! Build dynamic A matrix - - 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_A = 0d0 - chi_B = 0d0 - - do kc=1,maxS - - chi_A = chi_A + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) - chi_B = chi_B + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) - - enddo - - Ap(ia,jb) = Ap(ia,jb) - 4d0*lambda*chi_A - Am(ia,jb) = Am(ia,jb) - 4d0*lambda*chi_A - - Bp(ia,jb) = Bp(ia,jb) - 4d0*lambda*chi_B - Bm(ia,jb) = Bm(ia,jb) - 4d0*lambda*chi_B - - chi_Ap = 0d0 - chi_Am = 0d0 - - chi_Bp = 0d0 - chi_Bm = 0d0 - - do kc=1,maxS - - eps_Ap = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) - - eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Ap = chi_Ap + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Ap/(eps_Ap**2 + eta**2) - - eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) - chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) - - eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) - chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) - - eps_Bp = - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) - - eps_Bp = - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) - - eps_Bm = - OmRPA(kc) - (eGW(a) - eGW(b)) - chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) - - eps_Bm = - OmRPA(kc) - (eGW(j) - eGW(i)) - chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) - - enddo - - Ap(ia,jb) = Ap(ia,jb) - 2d0*lambda*chi_Ap - Am(ia,jb) = Am(ia,jb) - 2d0*lambda*chi_Am - - Bp(ia,jb) = Bp(ia,jb) - 2d0*lambda*chi_Bp - Bm(ia,jb) = Bm(ia,jb) - 2d0*lambda*chi_Bm - - enddo - enddo - enddo - enddo - -end subroutine diff --git a/src/GW/GW_phBSE_dynamic_kernel_A.f90 b/src/GW/GW_phBSE_dynamic_kernel_A.f90 index 81ee6ba..ea5a244 100644 --- a/src/GW/GW_phBSE_dynamic_kernel_A.f90 +++ b/src/GW/GW_phBSE_dynamic_kernel_A.f90 @@ -17,7 +17,6 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh ! Local variables - integer :: maxS double precision :: chi double precision :: eps integer :: i,j,a,b,ia,jb,kc @@ -32,14 +31,10 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh A_dyn(:,:) = 0d0 ZA_dyn(:,:) = 0d0 -! Number of poles taken into account - - maxS = nS - ! Build dynamic A matrix jb = 0 -!$omp parallel do default(private) shared(A_dyn,ZA_dyn,OmRPA,OmBSE,eGW,rho_RPA,nO,nBas,maxS,chi,eps,eta,nC,nR,lambda) +!$omp parallel do default(private) shared(A_dyn,ZA_dyn,OmRPA,OmBSE,eGW,rho_RPA,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) @@ -51,7 +46,7 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh chi = 0d0 - do kc=1,maxS + do kc=1,nS eps = OmRPA(kc) chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2) enddo @@ -60,7 +55,7 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh chi = 0d0 - do kc=1,maxS + do kc=1,nS eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps/(eps**2 + eta**2) @@ -73,7 +68,7 @@ subroutine GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rh A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi chi = 0d0 - do kc=1,maxS + do kc=1,nS eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) chi = chi + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 diff --git a/src/GW/GW_phBSE_dynamic_kernel_B.f90 b/src/GW/GW_phBSE_dynamic_kernel_B.f90 new file mode 100644 index 0000000..29365e3 --- /dev/null +++ b/src/GW/GW_phBSE_dynamic_kernel_B.f90 @@ -0,0 +1,60 @@ +subroutine GW_phBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,rhO,OmBSE,KB) + +! Compute the dynamic part of the Bethe-Salpeter equation matrices + + 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) :: eGW(nBas) + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + double precision,intent(in) :: OmBSE + +! Local variables + + double precision :: chi,eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: KB(nS,nS) + +! Initialization + + KB(:,:) = 0d0 + +! Build dynamic A matrix + + 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) - (eGW(a) - eGW(b)) + chi = chi + rho(i,b,kc)*rho(a,j,kc)*eps/(eps**2 + eta**2) + + eps = - OmRPA(kc) - (eGW(j) - eGW(i)) + chi = chi + rho(i,b,kc)*rho(a,j,kc)*eps/(eps**2 + eta**2) + + enddo + + KB(ia,jb) = KB(ia,jb) - 2d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine diff --git a/src/GW/GW_phBSE_dynamic_perturbation.f90 b/src/GW/GW_phBSE_dynamic_perturbation.f90 index 1285c78..fb2fc9b 100644 --- a/src/GW/GW_phBSE_dynamic_perturbation.f90 +++ b/src/GW/GW_phBSE_dynamic_perturbation.f90 @@ -1,5 +1,5 @@ -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) +subroutine GW_phBSE_dynamic_perturbation(dophBSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eGW,ERI,dipole_int, & + OmRPA,rho_RPA,OmBSE,XpY,XmY,KA_sta,KB_sta) ! Compute dynamical effects via perturbation theory for BSE @@ -8,7 +8,7 @@ subroutine GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eG ! Input variables - logical,intent(in) :: BSE2 + logical,intent(in) :: dophBSE2 logical,intent(in) :: dTDA double precision,intent(in) :: eta integer,intent(in) :: nBas @@ -18,6 +18,7 @@ subroutine GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eG integer,intent(in) :: nR integer,intent(in) :: nS + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: eW(nBas) double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: dipole_int(nBas,nBas,ncart) @@ -26,8 +27,8 @@ subroutine GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eG double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) - double precision,intent(in) :: W(nBas,nBas,nBas,nBas) - double precision,intent(in) :: A_stat(nS,nS) + double precision,intent(in) :: KA_sta(nS,nS) + double precision,intent(in) :: KB_sta(nS,nS) ! Local variables @@ -41,23 +42,19 @@ subroutine GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eG double precision,allocatable :: X(:) double precision,allocatable :: Y(:) - double precision,allocatable :: Ap_dyn(:,:) + double precision,allocatable :: KAp_dyn(:,:) + double precision,allocatable :: KAm_dyn(:,:) double precision,allocatable :: ZAp_dyn(:,:) - - double precision,allocatable :: Bp_dyn(:,:) - double precision,allocatable :: ZBp_dyn(:,:) - - double precision,allocatable :: Am_dyn(:,:) double precision,allocatable :: ZAm_dyn(:,:) - double precision,allocatable :: Bm_dyn(:,:) - double precision,allocatable :: ZBm_dyn(:,:) + double precision,allocatable :: KB_dyn(:,:) + + double precision,allocatable :: W(:,:,:,:) ! Memory allocation - maxS = min(nS,maxS) - allocate(OmDyn(maxS),ZDyn(maxS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_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)) + allocate(OmDyn(maxS),ZDyn(maxS),X(nS),Y(nS),KAp_dyn(nS,nS),ZAp_dyn(nS,nS), & + KAm_dyn(nS,nS),ZAm_dyn(nS,nS),KB_dyn(nS,nS)) if(dTDA) then write(*,*) @@ -65,6 +62,17 @@ subroutine GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eG write(*,*) end if + if(dophBSE2) then + + write(*,*) + write(*,*) '*** Second-order BSE static kernel activated! ***' + write(*,*) + + allocate(W(nBas,nBas,nBas,nBas)) + call static_kernel_W(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,W) + + end if + gapGW = eGW(nO+1) - eGW(nO) write(*,*) '---------------------------------------------------------------------------------------------------' @@ -75,7 +83,7 @@ subroutine GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eG write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' - do ia=1,maxS + do ia=1,min(nS,maxS) X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) @@ -86,32 +94,29 @@ subroutine GW_phBSE_dynamic_perturbation(BSE2,dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,eG ! Resonant part of the BSE correction for dynamical TDA - call GW_phBSE_dynamic_kernel_A(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),KAp_dyn,ZAp_dyn) - if(BSE2) call GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),Ap_dyn,ZAp_dyn,W) + if(dophBSE2) call GW_phBSE2_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,eGW,W,OmBSE(ia),KAp_dyn,ZAp_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) - OmDyn(ia) = dot_product(X,matmul( Ap_dyn - A_stat,X)) + OmDyn(ia) = dot_product(X,matmul(KAp_dyn - KA_sta,X)) else ! Resonant and anti-resonant part of the BSE correction - 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) + call GW_phBSE_dynamic_kernel_A(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,-OmBSE(ia),KAm_dyn,ZAm_dyn) + call GW_phBSE_dynamic_kernel_B(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,rho_RPA,KB_dyn) ! Renormalization factor of the resonant and anti-resonant parts - 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)) & - + dot_product(X,matmul(ZBp_dyn,Y)) & - - dot_product(Y,matmul(ZBm_dyn,X)) + + dot_product(Y,matmul(ZAm_dyn,Y)) - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & - - dot_product(Y,matmul(Am_dyn,Y)) & - + dot_product(X,matmul(Bp_dyn,Y)) & - - dot_product(Y,matmul(Bm_dyn,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/GW_phBSE_dynamic_perturbation_iterative.f90 b/src/GW/GW_phBSE_dynamic_perturbation_iterative.f90 deleted file mode 100644 index 4d39011..0000000 --- a/src/GW/GW_phBSE_dynamic_perturbation_iterative.f90 +++ /dev/null @@ -1,145 +0,0 @@ -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 - - implicit none - include 'parameters.h' - -! Input variables - - 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) :: eGW(nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: OmRPA(nS) - double precision,intent(in) :: rho_RPA(nBas,nBas,nS) - double precision,intent(in) :: OmBSE(nS) - double precision,intent(in) :: XpY(nS,nS) - double precision,intent(in) :: XmY(nS,nS) - -! Local variables - - integer :: ia - - integer :: maxS = 10 - double precision :: gapGW - - integer :: nSCF - integer :: maxSCF = 10 - double precision :: Conv - double precision :: thresh = 1d-3 - - - double precision,allocatable :: OmDyn(:) - double precision,allocatable :: OmOld(:) - double precision,allocatable :: X(:) - double precision,allocatable :: Y(:) - - double precision,allocatable :: Ap_dyn(:,:) - double precision,allocatable :: ZAp_dyn(:,:) - - double precision,allocatable :: Bp_dyn(:,:) - - double precision,allocatable :: Am_dyn(:,:) - double precision,allocatable :: ZAm_dyn(:,:) - - double precision,allocatable :: Bm_dyn(:,:) - -! Memory allocation - - maxS = min(nS,maxS) - allocate(OmDyn(maxS),OmOld(maxS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) - allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),Bm_dyn(nS,nS)) - - if(dTDA) then - write(*,*) - write(*,*) '*** dynamical TDA activated ***' - write(*,*) - end if - - gapGW = eGW(nO+1) - eGW(nO) - - Conv = 1d0 - nSCF = 0 - OmOld(1:maxS) = OmBSE(1:maxS) - - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies ' - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV' - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) - - do while(Conv > thresh .and. nSCF < maxSCF) - - nSCF = nSCF + 1 - - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(2X,A15,I3)') 'Iteration n.',nSCF - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Convergence (eV)' - write(*,*) '---------------------------------------------------------------------------------------------------' - - do ia=1,maxS - - X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) - Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - - ! First-order correction - - if(dTDA) then - - ! Resonant part of the BSE correction - - 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(:))) - - else - - ! Anti-resonant part of the BSE correction - - 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(:))) & - + dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) - - end if - - write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & - ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,(OmBSE(ia) + OmDyn(ia) - OmOld(ia))*HaToeV - - end do - - Conv = maxval(abs(OmBSE(1:maxS) + OmDyn(:) - OmOld(:)))*HaToeV - OmOld(:) = OmBSE(1:maxS) + OmDyn(:) - - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(2X,A20,1X,F10.6)') ' Convergence = ',Conv - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) - - end do - -! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) - - endif - -end subroutine diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index 62d26af..df88aef 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -1,5 +1,5 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & + dBSE,dTDA,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -21,7 +21,6 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta @@ -308,7 +307,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE if(BSE) then - call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) + call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 975e754..3862679 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -1,4 +1,4 @@ -subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & +subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eGW) @@ -19,7 +19,6 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip logical,intent(in) :: linearize @@ -193,7 +192,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev if(BSE) then - call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, & + call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,S, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eHF,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 9d07238..1c875c9 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -1,4 +1,4 @@ -subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE, & +subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, & cHF,eHF,Vxc) @@ -23,7 +23,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: doppBSE logical,intent(in) :: singlet logical,intent(in) :: triplet @@ -234,7 +233,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop if(dophBSE) then - call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE) + call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index 29f74f8..43fe783 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -1,5 +1,5 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, & + dBSE,dTDA,spin_conserved,spin_flip,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, & EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -23,7 +23,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta @@ -241,7 +240,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE if(BSE) then - call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & + call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,cHF,eGW,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index 9062b32..e486e58 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -1,4 +1,4 @@ -subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE, & +subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO, & ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) @@ -21,7 +21,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: doppBSE logical,intent(in) :: singlet logical,intent(in) :: triplet @@ -292,7 +291,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dop if(dophBSE) then - call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) + call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index 189f4d5..1dfce97 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -1,5 +1,5 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & + dBSE,dTDA,spin_conserved,spin_flip,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa, & dipole_int_bb,PHF,cHF,eHF) @@ -22,7 +22,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta @@ -371,7 +370,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE if(BSE) then - call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & + call unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS, & S,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,c,eGW,eGW,EcBSE) if(exchange_kernel) then diff --git a/src/GW/unrestricted_Bethe_Salpeter.f90 b/src/GW/unrestricted_Bethe_Salpeter.f90 index a483c0e..cd50b1e 100644 --- a/src/GW/unrestricted_Bethe_Salpeter.f90 +++ b/src/GW/unrestricted_Bethe_Salpeter.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & +subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta, & nBas,nC,nO,nV,nR,nS,S,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,cW,eW,eGW,EcBSE) @@ -13,7 +13,6 @@ subroutine unrestricted_Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved, logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA - logical,intent(in) :: evDyn logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index a0ee3dd..025b590 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -101,7 +101,7 @@ program QuAcK logical :: DIIS_GT,TDA_T,linGT,regGT double precision :: eta_GT - logical :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA,evDyn + logical :: dophBSE,dophBSE2,doppBSE,dBSE,dTDA ! Hello World @@ -145,7 +145,7 @@ program QuAcK maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW,COHSEX,TDA_W, & maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & doACFDT,exchange_kernel,doXBS, & - dophBSE,dBSE,dTDA,evDyn,doppBSE,dophBSE2) + dophBSE,dBSE,dTDA,doppBSE,dophBSE2) !------------------------------------------------------------------------ ! Read input information @@ -667,13 +667,13 @@ program QuAcK if(unrestricted) then - call UG0F2(dophBSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linGF,eta_GF,regGF, & + call UG0F2(dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,linGF,eta_GF,regGF, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & dipole_int_aa,dipole_int_bb,epsHF) else - call G0F2(dophBSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF,eta_GF,regGF, & + call G0F2(dophBSE,TDA,dBSE,dTDA,singlet,triplet,linGF,eta_GF,regGF, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF) end if @@ -696,13 +696,13 @@ program QuAcK if(unrestricted) then - call evUGF2(maxSCF_GF,thresh_GF,n_diis_GF,dophBSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & + call evUGF2(maxSCF_GF,thresh_GF,n_diis_GF,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip, & eta_GF,regGF,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & dipole_int_aa,dipole_int_bb,cHF,epsHF) else - call evGF2(dophBSE,TDA,dBSE,dTDA,evDyn,maxSCF_GF,thresh_GF,n_diis_GF, & + call evGF2(dophBSE,TDA,dBSE,dTDA,maxSCF_GF,thresh_GF,n_diis_GF, & singlet,triplet,linGF,eta_GF,regGF,nBas,nC,nO,nV,nR,nS,ENuc,EHF, & ERI_MO,dipole_int_MO,epsHF) @@ -726,13 +726,13 @@ program QuAcK if(unrestricted) then - call qsUGF2(maxSCF_GF,thresh_GF,n_diis_GF,dophBSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GF,regGF, & + call qsUGF2(maxSCF_GF,thresh_GF,n_diis_GF,dophBSE,TDA,dBSE,dTDA,spin_conserved,spin_flip,eta_GF,regGF, & nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO, & ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) else - call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,dophBSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GF,regGF,nNuc,ZNuc,rNuc,ENuc, & + call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,dophBSE,TDA,dBSE,dTDA,singlet,triplet,eta_GF,regGF,nNuc,ZNuc,rNuc,ENuc, & nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) end if @@ -806,12 +806,12 @@ program QuAcK call cpu_time(start_GW) if(unrestricted) then - call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & + call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,TDA_W,TDA,dBSE,dTDA,spin_conserved,spin_flip, & linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,EHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF,Vxc) else - call G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, & + call G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) end if @@ -833,14 +833,14 @@ program QuAcK if(unrestricted) then call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc, & + dBSE,dTDA,spin_conserved,spin_flip,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc, & EHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & PHF,cHF,epsHF,Vxc) else call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet,linGW,eta_GW,regGW, & + dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet,linGW,eta_GW,regGW, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) end if call cpu_time(end_GW) @@ -862,14 +862,14 @@ program QuAcK if(unrestricted) then call qsUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,dophBSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,regGW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & + dBSE,dTDA,spin_conserved,spin_flip,eta_GW,regGW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO, & dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) else call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,regGW,nNuc,ZNuc,rNuc,ENuc, & + dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta_GW,regGW,nNuc,ZNuc,rNuc,ENuc, & nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) end if @@ -896,7 +896,7 @@ program QuAcK else - call SRG_qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn, & + call SRG_qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA, & singlet,triplet,eta_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO, & dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) @@ -952,14 +952,14 @@ program QuAcK if(unrestricted) then - call UG0T0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn, & + call UG0T0(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA, & spin_conserved,spin_flip,linGT,eta_GT,regGT,nBas,nC,nO,nV, & nR,nS,ENuc,EHF,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF,Vxc) else - call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, & + call G0T0pp(doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) end if @@ -983,7 +983,7 @@ program QuAcK if(unrestricted) then call evUGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & - dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,& + dophBSE,TDA_T,TDA,dBSE,dTDA,spin_conserved,spin_flip,& eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO, & ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa, & dipole_int_bb,PHF,cHF,epsHF,Vxc) @@ -991,7 +991,7 @@ program QuAcK else call evGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & - dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT, & + dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta_GT,regGT, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO, & PHF,cHF,epsHF,Vxc) @@ -1016,13 +1016,13 @@ program QuAcK if(unrestricted) then call qsUGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T, & - TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GT,regGT,nBas,nC,nO,nV,& + TDA,dBSE,dTDA,spin_conserved,spin_flip,eta_GT,regGT,nBas,nC,nO,nV,& nR,nS,nNuc,ZNuc,rNuc,ENuc,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,& ERI_MO_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,epsHF) else call qsGTpp(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & - dophBSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT, & + dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta_GT,regGT, & nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc, & ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) @@ -1050,7 +1050,7 @@ program QuAcK else - call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet, & + call G0T0eh(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) end if @@ -1075,7 +1075,7 @@ program QuAcK else call evGTeh(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,evDyn,doppBSE,singlet,triplet,linGT,eta_GT,regGT, & + dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet,linGT,eta_GT,regGT, & nBas,nC,nO,nV,nR,nS,ENuc,EHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,epsHF,Vxc) end if call cpu_time(end_GT) @@ -1099,7 +1099,7 @@ program QuAcK else call qsGTeh(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT,nNuc,ZNuc,rNuc,ENuc, & + dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta_GT,regGT,nNuc,ZNuc,rNuc,ENuc, & nBas,nC,nO,nV,nR,nS,EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF) end if diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index fbcde29..f3fe5df 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -6,7 +6,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW,COHSEX,TDA_W, & maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & doACFDT,exchange_kernel,doXBS, & - dophBSE,dophBSE2,doppBSE,dBSE,dTDA,evDyn) + dophBSE,dophBSE2,doppBSE,dBSE,dTDA) ! Read desired methods @@ -74,7 +74,6 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: doppBSE logical,intent(out) :: dBSE logical,intent(out) :: dTDA - logical,intent(out) :: evDyn ! Local variables @@ -227,17 +226,15 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t doppBSE = .false. dBSE = .false. dTDA = .true. - evDyn = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 + read(1,*) answer1,answer2,answer3,answer4,answer5 if(answer1 == 'T') dophBSE = .true. if(answer2 == 'T') dophBSE2 = .true. if(answer3 == 'T') doppBSE = .true. if(answer4 == 'T') dBSE = .true. if(answer5 == 'F') dTDA = .false. - if(answer6 == 'T') evDyn = .true. ! Close file with options