diff --git a/input/options b/input/options index 8414408..82ee6d1 100644 --- a/input/options +++ b/input/options @@ -8,9 +8,11 @@ T T # GF: maxSCF thresh DIIS n_diis lin renorm BSE TDA eta 256 0.00001 T 5 T 3 T T 0.00367493 -# GW/GT: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA_W TDA dTDA G0W GW0 lin eta - 256 0.00001 T 5 F F T F F F F F T 0.00367493 +# GW/GT: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA_W TDA G0W GW0 lin eta + 256 0.00001 T 5 F F T F F F F T 0.00367493 # ACFDT: AC Kx XBS - T F T + F F T +# dBSE: dBSE dTDA evDyn + T T T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index ca5eac8..9e4b26a 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & +subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -10,7 +10,9 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dBSE logical,intent(in) :: dTDA + logical,intent(in) :: evDyn logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold @@ -27,7 +29,6 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, ! Local variables - logical :: evDyn = .true. logical :: W_BSE = .false. integer :: ispin integer :: isp_W @@ -70,24 +71,30 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) + !------------------------------------------------- ! Compute the dynamical screening at the BSE level + !------------------------------------------------- - if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + if(dBSE) then - ! Compute dynamic correction for BSE via perturbation theory - - if(evDyn) then - - call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) - else - - if(W_BSE) then - call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + + ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) + + if(evDyn) then + + call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + + if(W_BSE) then + call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + else + call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + end if + end if end if @@ -118,24 +125,30 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, rho_RPA(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) + !------------------------------------------------- ! Compute the dynamical screening at the BSE level + !------------------------------------------------- - if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + if(dBSE) then - ! Compute dynamic correction for BSE via perturbation theory + if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) - if(evDyn) then + ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) - call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) - else - - if(W_BSE) then - call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + if(evDyn) then + + call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + + if(W_BSE) then + call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) + else + call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) + end if + end if end if diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index c36d140..231c842 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -1,4 +1,4 @@ -subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold, & +subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -14,7 +14,9 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dTDA,singlet_manifol logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dBSE logical,intent(in) :: dTDA + logical,intent(in) :: evDyn logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold logical,intent(in) :: linearize @@ -212,7 +214,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dTDA,singlet_manifol allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin)) - call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eHF,eG0T0,Omega,XpY,XmY,rho,EcRPA,EcBSE) if(exchange_kernel) then diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 4c3e87d..a330cdf 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -1,5 +1,5 @@ -subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dTDA, & - singlet_manifold,triplet_manifold,linearize,eta, & +subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & + dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linearize,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW) ! Perform G0W0 calculation @@ -18,7 +18,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dTDA, & logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dBSE logical,intent(in) :: dTDA + logical,intent(in) :: evDyn logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold logical,intent(in) :: linearize @@ -164,7 +166,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dTDA, & if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) if(exchange_kernel) then diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index d1c3472..084319c 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -101,9 +101,11 @@ program QuAcK integer :: maxSCF_GW,n_diis_GW double precision :: thresh_GW - logical :: DIIS_GW,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,linGW + logical :: DIIS_GW,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,G0W,GW0,linGW double precision :: eta_GW + logical :: dBSE,dTDA,evDyn + integer :: nMC,nEq,nWalk,nPrint,iSeed double precision :: dt logical :: doDrift @@ -147,8 +149,9 @@ program QuAcK maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,renormGF, & BSE_GF,TDA_GF,eta_GF, & maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, & - COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,linGW,eta_GW, & + COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,G0W,GW0,linGW,eta_GW, & doACFDT,exchange_kernel,doXBS, & + dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Weird stuff @@ -671,8 +674,8 @@ program QuAcK if(doG0W0) then call cpu_time(start_G0W0) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA, & - singlet_manifold,triplet_manifold,linGW,eta_GW, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW, & + dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) @@ -689,8 +692,8 @@ program QuAcK if(doevGW) then call cpu_time(start_evGW) - call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & - BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & + call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & + BSE_GW,TDA_W,TDA_GW,G0W,GW0,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_evGW) @@ -707,8 +710,8 @@ program QuAcK if(doqsGW) then call cpu_time(start_qsGW) - call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & - BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & + call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & + BSE_GW,TDA_W,TDA_GW,G0W,GW0,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF) call cpu_time(end_qsGW) @@ -727,8 +730,8 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_W,TDA_GW,dTDA, & - singlet_manifold,triplet_manifold,linGW,eta_GW, & + call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_W,TDA_GW, & + dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_G0T0) @@ -745,8 +748,8 @@ program QuAcK if(doevGT) then call cpu_time(start_evGT) - call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & - BSE_GW,TDA_W,TDA_GW,dTDA,singlet_manifold,triplet_manifold,eta_GW, & + call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & + BSE_GW,TDA_W,TDA_GW,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_evGT) @@ -862,8 +865,8 @@ program QuAcK ! Long-range G0W0 calculation call cpu_time(start_G0W0) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA, & - singlet_manifold,triplet_manifold,linGW,eta_GW, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW, & + dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) @@ -876,8 +879,8 @@ program QuAcK ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:) call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_W,TDA_GW,dTDA, & - singlet_manifold,triplet_manifold,linGW,eta_GW, & + call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_W,TDA_GW,dBSE,dTDA,evDyn, & + singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0) call cpu_time(end_G0T0) diff --git a/src/QuAcK/evGT.f90 b/src/QuAcK/evGT.f90 index 436379a..1e8c030 100644 --- a/src/QuAcK/evGT.f90 +++ b/src/QuAcK/evGT.f90 @@ -1,5 +1,5 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & - BSE,TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold, & + BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) @@ -18,7 +18,9 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dBSE logical,intent(in) :: dTDA + logical,intent(in) :: evDyn logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold double precision,intent(in) :: eta @@ -260,7 +262,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin)) - call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcRPA,EcBSE) if(exchange_kernel) then diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 2d25f4e..625c1fa 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,4 +1,4 @@ -subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dTDA,G0W,GW0, & +subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0, & singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -21,7 +21,9 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dBSE logical,intent(in) :: dTDA + logical,intent(in) :: evDyn logical,intent(in) :: G0W logical,intent(in) :: GW0 logical,intent(in) :: singlet_manifold @@ -228,7 +230,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) write(*,*) diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 0e61ea7..8191c4b 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -1,5 +1,5 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & - COHSEX,SOSEX,BSE,TDA_W,TDA,dTDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, & + COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,G0W,GW0,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -20,7 +20,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dBSE logical,intent(in) :: dTDA + logical,intent(in) :: evDyn logical,intent(in) :: G0W logical,intent(in) :: GW0 logical,intent(in) :: singlet_manifold @@ -271,7 +273,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) write(*,*) diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 9f4bc05..00782bc 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -4,8 +4,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,renormGF, & BSE_GF,TDA_GF,eta_GF, & maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, & - COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,linGW,eta_GW, & + COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,G0W,GW0,linGW,eta_GW, & doACFDT,exchange_kernel,doXBS, & + dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Read desired methods @@ -48,7 +49,6 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: BSE_GW logical,intent(out) :: TDA_W logical,intent(out) :: TDA_GW - logical,intent(out) :: dTDA logical,intent(out) :: G0W logical,intent(out) :: GW0 logical,intent(out) :: linGW @@ -58,6 +58,10 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: exchange_kernel logical,intent(out) :: doXBS + logical,intent(out) :: dBSE + logical,intent(out) :: dTDA + logical,intent(out) :: evDyn + integer,intent(out) :: nMC integer,intent(out) :: nEq integer,intent(out) :: nWalk @@ -153,7 +157,6 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t BSE_GW = .false. TDA_W = .false. TDA_GW = .false. - dTDA = .false. G0W = .false. GW0 = .false. linGW = .false. @@ -162,18 +165,17 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t read(1,*) read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2, & answer3,answer4,answer5,answer6,answer7,answer8, & - answer9,answer10,eta_GW + answer9,eta_GW - if(answer1 == 'T') DIIS_GW = .true. - if(answer2 == 'T') COHSEX = .true. - if(answer3 == 'T') SOSEX = .true. - if(answer4 == 'T') BSE_GW = .true. - if(answer5 == 'T') TDA_W = .true. - if(answer6 == 'T') TDA_GW = .true. - if(answer7 == 'T') dTDA = .true. - if(answer8 == 'T') G0W = .true. - if(answer9 == 'T') GW0 = .true. - if(answer10 == 'T') linGW = .true. + if(answer1 == 'T') DIIS_GW = .true. + if(answer2 == 'T') COHSEX = .true. + if(answer3 == 'T') SOSEX = .true. + if(answer4 == 'T') BSE_GW = .true. + if(answer5 == 'T') TDA_W = .true. + if(answer6 == 'T') TDA_GW = .true. + if(answer7 == 'T') G0W = .true. + if(answer8 == 'T') GW0 = .true. + if(answer9 == 'T') linGW = .true. if(.not.DIIS_GW) n_diis_GW = 1 ! Options for adiabatic connection @@ -189,6 +191,19 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t if(answer2 == 'T') exchange_kernel = .true. if(answer3 == 'T') doXBS = .true. +! Options for dynamical BSE calculations + + dBSE = .false. + dTDA = .true. + evDyn = .false. + + read(1,*) + read(1,*) answer1,answer2,answer3 + + if(answer1 == 'T') dBSE = .true. + if(answer2 == 'F') dTDA = .false. + if(answer3 == 'T') evDyn = .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