diff --git a/input/options b/input/options index 82ee6d1..625c9c8 100644 --- a/input/options +++ b/input/options @@ -4,15 +4,15 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 -# spin: singlet triplet - 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 G0W GW0 lin eta - 256 0.00001 T 5 F F T F F F F T 0.00367493 +# spin: singlet triplet TDA + T T F +# GF: maxSCF thresh DIIS n_diis lin eta renorm + 256 0.00001 T 5 T 0.00367493 3 +# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 + 256 0.00001 T 5 T 0.00367493 F F F F F # ACFDT: AC Kx XBS F F T -# dBSE: dBSE dTDA evDyn - T T T +# BSE: BSE dBSE dTDA evDyn + T T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/BSE2.f90 b/src/QuAcK/BSE2.f90 index 9dc90ad..6e65a2d 100644 --- a/src/QuAcK/BSE2.f90 +++ b/src/QuAcK/BSE2.f90 @@ -1,4 +1,5 @@ -subroutine BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,EcBSE) +subroutine BSE2(TDA,dBSE,dTDA,singlet_manifold,triplet_manifold, & + eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -8,6 +9,8 @@ subroutine BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ER ! Input variables logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold @@ -54,9 +57,9 @@ subroutine BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ER call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) ! Compute dynamic correction for BSE via perturbation theory - - call BSE2_dynamic_perturbation(ispin,eta,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + if(dBSE) & + call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & + ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) end if @@ -77,8 +80,9 @@ subroutine BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ER ! Compute dynamic correction for BSE via perturbation theory - call BSE2_dynamic_perturbation(ispin,eta,nBas,nC,nO,nV,nR,nS, & - ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + if(dBSE) & + call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS, & + ERI(:,:,:,:),eHF(:),eGF(:),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) end if diff --git a/src/QuAcK/BSE2_dynamic_perturbation.f90 b/src/QuAcK/BSE2_dynamic_perturbation.f90 index 37bc8ff..0efb9bc 100644 --- a/src/QuAcK/BSE2_dynamic_perturbation.f90 +++ b/src/QuAcK/BSE2_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine BSE2_dynamic_perturbation(ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,OmBSE,XpY,XmY) +subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -7,6 +7,7 @@ subroutine BSE2_dynamic_perturbation(ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,O ! Input variables + logical,intent(in) :: dTDA integer,intent(in) :: ispin double precision,intent(in) :: eta integer,intent(in) :: nBas @@ -25,7 +26,6 @@ subroutine BSE2_dynamic_perturbation(ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,O ! Local variables - logical :: dTDA = .true. integer :: ia integer,parameter :: maxS = 10 double precision :: gapGF diff --git a/src/QuAcK/G0F2.f90 b/src/QuAcK/G0F2.f90 index 0bd5f46..8213ea2 100644 --- a/src/QuAcK/G0F2.f90 +++ b/src/QuAcK/G0F2.f90 @@ -1,4 +1,4 @@ -subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine G0F2(BSE,TDA,dBSE,dTDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Perform a one-shot second-order Green function calculation @@ -9,6 +9,8 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC, logical,intent(in) :: BSE logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold logical,intent(in) :: linearize @@ -111,7 +113,7 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC, if(BSE) then - call BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE) + call BSE2(TDA,dBSE,dTDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE) end if diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 084319c..bef4d19 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -93,18 +93,19 @@ program QuAcK logical :: singlet_manifold logical :: triplet_manifold + logical :: TDA integer :: maxSCF_GF,n_diis_GF,renormGF double precision :: thresh_GF - logical :: DIIS_GF,linGF,BSE_GF,TDA_GF + logical :: DIIS_GF,linGF double precision :: eta_GF integer :: maxSCF_GW,n_diis_GW double precision :: thresh_GW - logical :: DIIS_GW,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,G0W,GW0,linGW + logical :: DIIS_GW,COHSEX,SOSEX,TDA_W,G0W,GW0,linGW double precision :: eta_GW - logical :: dBSE,dTDA,evDyn + logical :: BSE,dBSE,dTDA,evDyn integer :: nMC,nEq,nWalk,nPrint,iSeed double precision :: dt @@ -145,13 +146,12 @@ program QuAcK call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, & maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - singlet_manifold,triplet_manifold, & - 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,G0W,GW0,linGW,eta_GW, & + singlet_manifold,triplet_manifold,TDA, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & + COHSEX,SOSEX,TDA_W,G0W,GW0, & doACFDT,exchange_kernel,doXBS, & - dBSE,dTDA,evDyn, & + BSE,dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Weird stuff @@ -606,7 +606,7 @@ program QuAcK if(doG0F2) then call cpu_time(start_GF2) - call G0F2(BSE_GF,TDA_GF,singlet_manifold,triplet_manifold,linGF, & + call G0F2(BSE,TDA,dBSE,dTDA,singlet_manifold,triplet_manifold,linGF, & eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) @@ -623,7 +623,8 @@ program QuAcK if(doevGF2) then call cpu_time(start_GF2) - call evGF2(BSE_GF,TDA_GF,maxSCF_GF,thresh_GF,n_diis_GF,singlet_manifold,triplet_manifold,linGF, & + call evGF2(BSE,TDA,dBSE,dTDA,maxSCF_GF,thresh_GF,n_diis_GF, & + singlet_manifold,triplet_manifold,linGF, & eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) @@ -674,7 +675,7 @@ program QuAcK if(doG0W0) then call cpu_time(start_G0W0) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & 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) @@ -692,8 +693,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,G0W,GW0,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, & + call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & + BSE,TDA_W,TDA,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) @@ -710,8 +711,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,G0W,GW0,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, & + call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & + BSE,TDA_W,TDA,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) @@ -730,7 +731,7 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_W,TDA_GW, & + call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA, & 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) @@ -748,8 +749,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,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta_GW, & + call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & + BSE,TDA_W,TDA,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) @@ -865,7 +866,7 @@ 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, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & 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) @@ -879,8 +880,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,dBSE,dTDA,evDyn, & - singlet_manifold,triplet_manifold,linGW,eta_GW, & + call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,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/evGF2.f90 b/src/QuAcK/evGF2.f90 index 9689c40..e5e9a97 100644 --- a/src/QuAcK/evGF2.f90 +++ b/src/QuAcK/evGF2.f90 @@ -1,4 +1,4 @@ -subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold,linearize, & +subroutine evGF2(BSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold,linearize, & eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Perform eigenvalue self-consistent second-order Green function calculation @@ -10,6 +10,8 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol logical,intent(in) :: BSE logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA integer,intent(in) :: maxSCF double precision,intent(in) :: thresh integer,intent(in) :: max_diis @@ -165,7 +167,7 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol if(BSE) then - call BSE2(TDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE) + call BSE2(TDA,dBSE,dTDA,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF2,EcBSE) end if diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 00782bc..005cd22 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -1,12 +1,11 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type, & maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - singlet_manifold,triplet_manifold, & - 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,G0W,GW0,linGW,eta_GW, & + singlet_manifold,triplet_manifold,TDA, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & + COHSEX,SOSEX,TDA_W,G0W,GW0, & doACFDT,exchange_kernel,doXBS, & - dBSE,dTDA,evDyn, & + BSE,dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Read desired methods @@ -29,6 +28,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: singlet_manifold logical,intent(out) :: triplet_manifold + logical,intent(out) :: TDA integer,intent(out) :: maxSCF_GF double precision,intent(out) :: thresh_GF @@ -36,8 +36,6 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t integer,intent(out) :: n_diis_GF logical,intent(out) :: linGF integer,intent(out) :: renormGF - logical,intent(out) :: BSE_GF - logical,intent(out) :: TDA_GF double precision,intent(out) :: eta_GF integer,intent(out) :: maxSCF_GW @@ -46,9 +44,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t integer,intent(out) :: n_diis_GW logical,intent(out) :: COHSEX logical,intent(out) :: SOSEX - logical,intent(out) :: BSE_GW logical,intent(out) :: TDA_W - logical,intent(out) :: TDA_GW logical,intent(out) :: G0W logical,intent(out) :: GW0 logical,intent(out) :: linGW @@ -58,6 +54,7 @@ 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) :: BSE logical,intent(out) :: dBSE logical,intent(out) :: dTDA logical,intent(out) :: evDyn @@ -118,12 +115,14 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t singlet_manifold = .false. triplet_manifold = .false. + TDA = .false. read(1,*) - read(1,*) answer1,answer2 + read(1,*) answer1,answer2,answer3 if(answer1 == 'T') singlet_manifold = .true. if(answer2 == 'T') triplet_manifold = .true. + if(answer3 == 'T') TDA = .true. ! Read Green function options @@ -132,18 +131,14 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t DIIS_GF = .false. n_diis_GF = 5 linGF = .false. - renormGF = 0 - BSE_GF = .false. - TDA_GF = .false. eta_GF = 0d0 + renormGF = 0 read(1,*) - read(1,*) maxSCF_GF,thresh_GW,answer1,n_diis_GF,answer2,renormGF,answer3,answer4,eta_GF + read(1,*) maxSCF_GF,thresh_GF,answer1,n_diis_GF,answer2,eta_GF,renormGF if(answer1 == 'T') DIIS_GF = .true. if(answer2 == 'T') linGF = .true. - if(answer3 == 'T') BSE_GF = .true. - if(answer4 == 'T') TDA_GF = .true. if(.not.DIIS_GF) n_diis_GF = 1 ! Read GW options @@ -154,35 +149,30 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t n_diis_GW = 5 COHSEX = .false. SOSEX = .false. - BSE_GW = .false. TDA_W = .false. - TDA_GW = .false. G0W = .false. GW0 = .false. linGW = .false. eta_GW = 0d0 read(1,*) - read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2, & - answer3,answer4,answer5,answer6,answer7,answer8, & - answer9,eta_GW + read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2,eta_GW, & + answer3,answer4,answer5,answer6,answer7 if(answer1 == 'T') DIIS_GW = .true. - if(answer2 == 'T') COHSEX = .true. - if(answer3 == 'T') SOSEX = .true. - if(answer4 == 'T') BSE_GW = .true. + if(answer2 == 'T') linGW = .true. + if(answer3 == 'T') COHSEX = .true. + if(answer4 == 'T') SOSEX = .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(answer6 == 'T') G0W = .true. + if(answer7 == 'T') GW0 = .true. if(.not.DIIS_GW) n_diis_GW = 1 ! Options for adiabatic connection - doACFDT = .false. + doACFDT = .false. exchange_kernel = .false. - doXBS = .false. + doXBS = .false. read(1,*) read(1,*) answer1,answer2,answer3 @@ -193,16 +183,18 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t ! Options for dynamical BSE calculations + BSE = .false. dBSE = .false. dTDA = .true. evDyn = .false. read(1,*) - read(1,*) answer1,answer2,answer3 + read(1,*) answer1,answer2,answer3,answer4 - if(answer1 == 'T') dBSE = .true. - if(answer2 == 'F') dTDA = .false. - if(answer3 == 'T') evDyn = .true. + if(answer1 == 'T') BSE = .true. + if(answer2 == 'T') dBSE = .true. + if(answer3 == 'F') dTDA = .false. + if(answer4 == '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