diff --git a/input/methods b/input/methods index 422d8d3..688d56a 100644 --- a/input/methods +++ b/input/methods @@ -13,9 +13,9 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - F F F F T + T F F F F # G0T0 evGT qsGT - f F F + T F F # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index f8084f4..6bc78e4 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability - 512 0.0000001 T 5 2 1 F 0.0 F + 512 0.0000001 T 0 1 1 F 0.0 F # MP: # CC: maxSCF thresh DIIS n_diis @@ -11,7 +11,7 @@ # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg 256 0.00001 T 5 T 0.0 F F F F F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg - 10 0.00001 T 5 T 0.0 T F + 10 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS F T T # BSE: BSE dBSE dTDA evDyn ppBSE diff --git a/int/Nuc.Hu.dat b/int/Nuc.Hu.dat index 9772f82..48f4be0 100644 --- a/int/Nuc.Hu.dat +++ b/int/Nuc.Hu.dat @@ -1,2 +1,2 @@ 1 1 0.0 -2 2 0.1 +2 2 0.0 diff --git a/src/GT/excitation_density_Tmatrix.f90 b/src/GT/excitation_density_Tmatrix.f90 index 6ee28f6..c125ade 100644 --- a/src/GT/excitation_density_Tmatrix.f90 +++ b/src/GT/excitation_density_Tmatrix.f90 @@ -228,4 +228,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r end if +! call matout(nBas**2,nVV,rho1) +! call matout(nBas**2,nOO,rho2) + end subroutine excitation_density_Tmatrix diff --git a/src/GW/BSE2_static_kernel.f90 b/src/GW/BSE2_static_kernel.f90 new file mode 100644 index 0000000..32d7be3 --- /dev/null +++ b/src/GW/BSE2_static_kernel.f90 @@ -0,0 +1,123 @@ +subroutine BSE2_static_kernel(eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,Om,rho,A_sta) + +! Compute the second-order static BSE kernel (only for singlets!) + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: lambda + + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eW(nBas) + double precision,intent(in) :: Om(nS) + + double precision,intent(in) :: rho(nBas,nBas,nS) + + +! Local variables + + double precision :: chi + integer :: p,q,r,s + integer :: m + + double precision :: dem,num + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb + + double precision,allocatable :: W(:,:,:,:) + +! Output variables + + double precision,intent(inout) :: A_sta(nS,nS) + +! memory allocation + + allocate(W(nBas,nBas,nBas,nBas)) + +!------------------------------------------------ +! Compute static screening (physicist's notation) +!------------------------------------------------ + + do p=1,nBas + do q=1,nBas + do r=1,nBas + do s=1,nBas + + chi = 0d0 + do m=1,nS + dem = Om(m)**2 + eta**2 + chi = chi + rho(p,q,m)*rho(r,s,m)*Om(m)/dem + enddo + + W(p,s,q,r) = - ERI(p,s,q,r) + 4d0*chi + + enddo + enddo + enddo + enddo + + 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 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = - (eW(c) - eW(k)) + num = 2d0*W(j,k,i,c)*W(a,c,b,k) + + A_sta(ia,jb) = A_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = + (eW(c) - eW(k)) + num = 2d0*W(j,c,i,k)*W(a,k,b,c) + + A_sta(ia,jb) = A_sta(ia,jb) + num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - (eW(c) + eW(d)) + num = 2d0*W(a,j,c,d)*W(c,d,i,b) + + A_sta(ia,jb) = A_sta(ia,jb) + num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = - (eW(k) + eW(l)) + num = 2d0*W(a,j,k,l)*W(k,l,i,b) + + A_sta(ia,jb) = A_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + +end subroutine BSE2_static_kernel diff --git a/src/GW/static_screening_WA.f90 b/src/GW/BSE_static_kernel_KA.f90 similarity index 84% rename from src/GW/static_screening_WA.f90 rename to src/GW/BSE_static_kernel_KA.f90 index dff07d7..725c73e 100644 --- a/src/GW/static_screening_WA.f90 +++ b/src/GW/BSE_static_kernel_KA.f90 @@ -1,6 +1,6 @@ -subroutine static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA) +subroutine BSE_static_kernel_KA(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 +! Compute the BSE static kernel for the resonant block implicit none include 'parameters.h' @@ -50,4 +50,4 @@ subroutine static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WA) enddo enddo -end subroutine static_screening_WA +end subroutine BSE_static_kernel_KA diff --git a/src/GW/static_screening_WB.f90 b/src/GW/BSE_static_kernel_KB.f90 similarity index 78% rename from src/GW/static_screening_WB.f90 rename to src/GW/BSE_static_kernel_KB.f90 index ee7ac48..deb7336 100644 --- a/src/GW/static_screening_WB.f90 +++ b/src/GW/BSE_static_kernel_KB.f90 @@ -1,6 +1,6 @@ -subroutine static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WB) +subroutine BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,KB) -! Compute the static screening W for the coupling block +! Compute the BSE static kernel for the coupling block implicit none include 'parameters.h' @@ -27,11 +27,11 @@ subroutine static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WB) ! Output variables - double precision,intent(out) :: WB(nS,nS) + double precision,intent(out) :: KB(nS,nS) ! Initialize - WB(:,:) = 0d0 + KB(:,:) = 0d0 ia = 0 do i=nC+1,nO @@ -48,11 +48,11 @@ subroutine static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,WB) 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) - 4d0*lambda*chi + KB(ia,jb) = KB(ia,jb) + lambda*ERI(i,j,b,a) - 4d0*lambda*chi enddo enddo enddo enddo -end subroutine static_screening_WB +end subroutine BSE_static_kernel_KB diff --git a/src/GW/Bethe_Salpeter.f90 b/src/GW/Bethe_Salpeter.f90 index 2440585..0fc3f57 100644 --- a/src/GW/Bethe_Salpeter.f90 +++ b/src/GW/Bethe_Salpeter.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) +subroutine Bethe_Salpeter(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -7,6 +7,7 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, ! Input variables + logical,intent(in) :: BSE2 logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -42,8 +43,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, double precision,allocatable :: XpY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:) - double precision,allocatable :: WA_sta(:,:) - double precision,allocatable :: WB_sta(:,:) + double precision,allocatable :: KA_sta(:,:) + double precision,allocatable :: KB_sta(:,:) ! Output variables @@ -52,7 +53,7 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, ! 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,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) + KA_sta(nS,nS),KB_sta(nS,nS),OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) !--------------------------------- ! Compute (singlet) RPA screening @@ -65,8 +66,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA_sta) - call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB_sta) + 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) !------------------- ! Singlet manifold @@ -77,9 +78,13 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, ispin = 1 EcBSE(ispin) = 0d0 + ! Second-order BSE staic kernel + + if(BSE2) call BSE2_static_kernel(eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA,rho_RPA,KA_sta) + ! 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, & + call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,KA_sta,KB_sta, & EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & @@ -118,7 +123,7 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, ! 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, & + call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,KA_sta,KB_sta, & EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 099e923..0b82089 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -1,4 +1,4 @@ -subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, & +subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGW) @@ -15,6 +15,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evD logical,intent(in) :: doXBS logical,intent(in) :: COHSEX logical,intent(in) :: BSE + logical,intent(in) :: BSE2 logical,intent(in) :: ppBSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA @@ -200,7 +201,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evD if(BSE) then - call Bethe_Salpeter(TDA_W,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_W,TDA,dBSE,dTDA,evDyn,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/Sangalli_dynamic_perturbation.f90 b/src/GW/Sangalli_dynamic_perturbation.f90 deleted file mode 100644 index 9718415..0000000 --- a/src/GW/Sangalli_dynamic_perturbation.f90 +++ /dev/null @@ -1,124 +0,0 @@ -subroutine Sangalli_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) - -! Compute dynamical effects via perturbation theory for Sangalli's kernel - - 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) :: OmRPA(nS) - double precision,intent(in) :: OmBSE(nS) - double precision,intent(in) :: XpY(nS,nS) - double precision,intent(in) :: XmY(nS,nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: ia - - integer,parameter :: maxS = 10 - double precision :: gapGW - - double precision,allocatable :: OmDyn(:) - double precision,allocatable :: ZDyn(:) - 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 :: ZBp_dyn(:,:) - - double precision,allocatable :: Am_dyn(:,:) - double precision,allocatable :: ZAm_dyn(:,:) - - double precision,allocatable :: Bm_dyn(:,:) - double precision,allocatable :: ZBm_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),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS)) - -! Print main components of transition vectors - - call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) - - if(dTDA) then - write(*,*) - write(*,*) '*** dynamical TDA activated ***' - write(*,*) - end if - - gapGW = eGW(nO+1) - eGW(nO) - - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) ' First-order dynamical correction to static BSE excitation energies with Sangalli kernel ' - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV' - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' - write(*,*) '---------------------------------------------------------------------------------------------------' - - do ia=1,min(nS,maxS) - - X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) - Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - - if(dTDA) then - -! call Sangalli_A_matrix_dynamic (eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),Ap_dyn(:,:)) -! call Sangalli_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),ZAp_dyn(:,:)) - - ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) - OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) - - else - -! call Sangalli_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),Ap_dyn(:,:)) -! call Sangalli_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),Am_dyn(:,:)) - -! call Sangalli_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),Bp_dyn(:,:)) -! call Sangalli_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),Bm_dyn(:,:)) - -! call Sangalli_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),ZAp_dyn(:,:)) -! call Sangalli_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),ZAm_dyn(:,:)) -! call Sangalli_ZB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),+OmBSE(ia),rho(:,:,:),ZBp_dyn(:,)) -! call Sangalli_ZB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),-OmBSE(ia),rho(:,:,:),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(:))) - - 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 - - ZDyn(ia) = 1d0/(1d0 - ZDyn(ia)) - OmDyn(ia) = ZDyn(ia)*OmDyn(ia) - - 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,ZDyn(ia) - - end do - write(*,*) '---------------------------------------------------------------------------------------------------' - write(*,*) - -end subroutine Sangalli_dynamic_perturbation diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 64586df..5218a3c 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,COHSEX,BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,ppBSE, & +subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,ppBSE, & singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -18,6 +18,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, logical,intent(in) :: doXBS logical,intent(in) :: COHSEX logical,intent(in) :: BSE + logical,intent(in) :: BSE2 logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -270,7 +271,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(BSE) then - call Bethe_Salpeter(TDA_W,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_W,TDA,dBSE,dTDA,evDyn,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/qsGW.f90 b/src/GW/qsGW.f90 index 9f1d33b..ec35dc0 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,COHSEX,BSE,TDA_W,TDA, & +subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,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) @@ -17,6 +17,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, logical,intent(in) :: doXBS logical,intent(in) :: COHSEX logical,intent(in) :: BSE + logical,intent(in) :: BSE2 logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -330,7 +331,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & + call Bethe_Salpeter(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,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/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 3a22bd1..d7ce9aa 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -139,7 +139,7 @@ program QuAcK logical :: DIIS_GT,TDA_T,linGT,regGT double precision :: eta_GT - logical :: BSE,dBSE,dTDA,evDyn,ppBSE + logical :: BSE,dBSE,dTDA,evDyn,ppBSE,BSE2 integer :: nMC,nEq,nWalk,nPrint,iSeed double precision :: dt @@ -188,7 +188,7 @@ program QuAcK COHSEX,SOSEX,TDA_W,G0W,GW0, & maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & doACFDT,exchange_kernel,doXBS, & - BSE,dBSE,dTDA,evDyn,ppBSE, & + BSE,dBSE,dTDA,evDyn,ppBSE,BSE2, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Weird stuff @@ -1028,10 +1028,10 @@ program QuAcK ! SOSEX extrension of GW if(SOSEX) then - call G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & + call G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) else - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) ! call soG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet, & ! linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) @@ -1064,7 +1064,7 @@ program QuAcK else call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX, & - BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,eta_GW,regGW, & + BSE,BSE2,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,eta_GW,regGW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) end if call cpu_time(end_evGW) @@ -1093,7 +1093,7 @@ program QuAcK else call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX, & - BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,regGW,nNuc,ZNuc,rNuc,ENuc, & + BSE,BSE2,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,regGW,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) end if diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index e3f1e99..2478088 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 COHSEX,SOSEX,TDA_W,G0W,GW0, & maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & doACFDT,exchange_kernel,doXBS, & - BSE,dBSE,dTDA,evDyn,ppBSE, & + BSE,dBSE,dTDA,evDyn,ppBSE,BSE2, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Read desired methods @@ -76,6 +76,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: dTDA logical,intent(out) :: evDyn logical,intent(out) :: ppBSE + logical,intent(out) :: BSE2 integer,intent(out) :: nMC integer,intent(out) :: nEq @@ -239,15 +240,17 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t dTDA = .true. evDyn = .false. ppBSE = .false. + BSE2 = .false. read(1,*) - read(1,*) answer1,answer2,answer3,answer4,answer5 + read(1,*) answer1,answer2,answer3,answer4,answer5,answer6 if(answer1 == 'T') BSE = .true. if(answer2 == 'T') dBSE = .true. if(answer3 == 'F') dTDA = .false. if(answer4 == 'T') evDyn = .true. if(answer5 == 'T') ppBSE = .true. + if(answer6 == 'T') BSE2 = .true. ! Read options for MC-MP2: Monte Carlo steps, number of equilibration steps, number of walkers, ! Monte Carlo time step, frequency of output results, and seed for random number generator diff --git a/src/RPA/ACFDT.f90 b/src/RPA/ACFDT.f90 index 4a7f33a..3a55c37 100644 --- a/src/RPA/ACFDT.f90 +++ b/src/RPA/ACFDT.f90 @@ -75,8 +75,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA) - call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) + call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA) + call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) ! Singlet manifold @@ -104,8 +104,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) ! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) - call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) - call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) end if @@ -155,8 +155,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) - call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) end if diff --git a/src/RPA/ACFDT_cr.f90 b/src/RPA/ACFDT_cr.f90 index 7d128a6..4df5b02 100644 --- a/src/RPA/ACFDT_cr.f90 +++ b/src/RPA/ACFDT_cr.f90 @@ -76,8 +76,8 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA) - call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) + call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA) + call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) ! Singlet manifold @@ -105,8 +105,8 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) ! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) - call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) - call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) end if @@ -156,8 +156,8 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) - call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + call BSE_static_kernel_KA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call BSE_static_kernel_KB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) end if diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index 3d2cace..4028c7f 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -216,4 +216,10 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) X2(1:nVV,1:nOO) = + Z2( 1: nVV,1:nOO) Y2(1:nOO,1:nOO) = - Z2(nVV+1:nOO+nVV,1:nOO) +! call matout(nVV,nVV,X1) +! call matout(nOO,nVV,Y1) + +! call matout(nVV,nOO,X2) +! call matout(nOO,nOO,Y2) + end subroutine sort_ppRPA