diff --git a/input/basis b/input/basis index 4c61da7..b9ca7b5 100644 --- a/input/basis +++ b/input/basis @@ -1,16 +1,9 @@ -1 6 -S 4 1.00 - 234.0000000 0.0025870 - 35.1600000 0.0195330 - 7.9890000 0.0909980 - 2.2120000 0.2720500 +1 3 +S 3 1.00 + 38.3600000 0.0238090 + 5.7700000 0.1548910 + 1.2400000 0.4699870 S 1 1.00 - 0.6669000 1.0000000 -S 1 1.00 - 0.2089000 1.0000000 + 0.2976000 1.0000000 P 1 1.00 - 3.0440000 1.0000000 -P 1 1.00 - 0.7580000 1.0000000 -D 1 1.00 - 1.9650000 1.0000000 + 1.2750000 1.0000000 diff --git a/input/methods b/input/methods index 02bb6a3..4d309d6 100644 --- a/input/methods +++ b/input/methods @@ -4,12 +4,12 @@ F F F # CCD CCSD CCSD(T) ringCCD ladderCCD F F F F F -# CIS RPA TDHF ppRPA ADC +# CIS RPA RPAx ppRPA ADC F T T F F # GF2 GF3 F F # G0W0 evGW qsGW - T F F + T T T # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index d04d036..07e4f0b 100644 --- a/input/options +++ b/input/options @@ -8,7 +8,9 @@ T T # GF: maxSCF thresh DIIS n_diis renormalization 1 0.00001 T 10 3 -# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta - 256 0.00001 T 5 F F T F F F F 0.000 +# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta + 256 0.00001 T 5 F F T F F F F 0.000 +# ACFDT: AC XBS + T T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/input/weight b/input/weight index 4c61da7..b9ca7b5 100644 --- a/input/weight +++ b/input/weight @@ -1,16 +1,9 @@ -1 6 -S 4 1.00 - 234.0000000 0.0025870 - 35.1600000 0.0195330 - 7.9890000 0.0909980 - 2.2120000 0.2720500 +1 3 +S 3 1.00 + 38.3600000 0.0238090 + 5.7700000 0.1548910 + 1.2400000 0.4699870 S 1 1.00 - 0.6669000 1.0000000 -S 1 1.00 - 0.2089000 1.0000000 + 0.2976000 1.0000000 P 1 1.00 - 3.0440000 1.0000000 -P 1 1.00 - 0.7580000 1.0000000 -D 1 1.00 - 1.9650000 1.0000000 + 1.2750000 1.0000000 diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 new file mode 100644 index 0000000..80cf51e --- /dev/null +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -0,0 +1,68 @@ +subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + +! Compute the Bethe-Salpeter excitation energies + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: singlet_manifold + logical,intent(in) :: triplet_manifold + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eW(nBas) + double precision,intent(in) :: eGW(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + + double precision :: Omega(nS,nspin) + double precision :: XpY(nS,nS,nspin) + double precision :: XmY(nS,nS,nspin) + double precision :: rho(nBas,nBas,nS,nspin) + +! Local variables + + integer :: ispin + +! Output variables + + double precision,intent(out) :: EcRPA(nspin) + double precision,intent(out) :: EcBSE(nspin) + +! Singlet manifold + + if(singlet_manifold) then + + ispin = 1 + EcBSE(ispin) = 0d0 + + call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + + call linear_response(ispin,.true.,TDA,.true.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) + + end if + +! Triplet manifold + + if(triplet_manifold) then + + ispin = 2 + EcBSE(ispin) = 0d0 + + call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + + call linear_response(ispin,.true.,TDA,.true.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) + + end if + +end subroutine Bethe_Salpeter diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 774dcb4..b9d4a4e 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -1,5 +1,5 @@ -subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) +subroutine G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW) ! Perform G0W0 calculation @@ -9,6 +9,8 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & ! Input variables + logical,intent(in) :: doACFDT + logical,intent(in) :: doXBS logical,intent(in) :: COHSEX logical,intent(in) :: SOSEX logical,intent(in) :: BSE @@ -41,12 +43,9 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rhox(:,:,:,:) - logical :: adiabatic_connection - logical :: scaled_screening - ! Output variables - double precision :: eG0W0(nBas) + double precision :: eGW(nBas) ! Hello world @@ -96,54 +95,23 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & ! Solve the quasi-particle equation - eG0W0(:) = eHF(:) + Z(:)*SigC(:) + eGW(:) = eHF(:) + Z(:)*SigC(:) ! Dump results - call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) - call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eG0W0,EcRPA(ispin),EcGM) + call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM) ! Plot stuff - call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eG0W0,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin)) + call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin)) ! Perform BSE calculation if(BSE) then - ! Singlet manifold - - if(singlet_manifold) then - - ispin = 1 - EcBSE(ispin) = 0d0 - - call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - - call linear_response(ispin,.true.,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) - - end if - - ! Triplet manifold - - if(triplet_manifold) then - - ispin = 2 - EcBSE(ispin) = 0d0 - - call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - - call linear_response(ispin,.true.,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) - - end if + call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -156,25 +124,22 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & ! Compute the BSE correlation energy via the adiabatic connection - adiabatic_connection = .true. - scaled_screening = .true. - - if(adiabatic_connection) then + if(doACFDT) then write(*,*) '------------------------------------------------------' write(*,*) 'Adiabatic connection version of BSE correlation energy' write(*,*) '------------------------------------------------------' write(*,*) - if(scaled_screening) then + if(doXBS) then - write(*,*) '*** scaled screening version (extended BSE) ***' + write(*,*) '*** scaled screening version (XBS) ***' write(*,*) end if - call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,nS,ERI,eG0W0,Omega,XpY,XmY,rho) + call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho) end if diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 44e063a..00cc4cd 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -24,9 +24,13 @@ program QuAcK double precision,allocatable :: ZNuc(:),rNuc(:,:) double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:) + double precision,allocatable :: eG0W0(:) double precision,allocatable :: eG0T0(:) + logical :: doACFDT + logical :: doXBS + integer :: nShell integer,allocatable :: TotAngMomShell(:),KShell(:) double precision,allocatable :: CenterShell(:,:),DShell(:,:),ExpShell(:,:) @@ -122,11 +126,13 @@ program QuAcK ! Read options for methods - 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,renormalization, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, & + 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,renormalization, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, & + COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, & + doACFDT,doXBS, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Weird stuff @@ -434,7 +440,7 @@ program QuAcK if(doRPA) then call cpu_time(start_RPA) - call RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) + call RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) call cpu_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -450,7 +456,7 @@ program QuAcK if(doRPAx) then call cpu_time(start_RPAx) - call RPAx(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) + call RPAx(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) call cpu_time(end_RPAx) t_RPAx = end_RPAx - start_RPAx @@ -466,7 +472,7 @@ program QuAcK if(doppRPA) then call cpu_time(start_ppRPA) - call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF) + call ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO_basis,eHF) call cpu_time(end_ppRPA) t_ppRPA = end_ppRPA - start_ppRPA @@ -532,7 +538,7 @@ program QuAcK if(doG0W0) then call cpu_time(start_G0W0) - call G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & + call G0W0(doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) @@ -549,7 +555,8 @@ program QuAcK if(doevGW) then call cpu_time(start_evGW) - call evGW(maxSCF_GW,thresh_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize,eta, & + call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, & + singlet_manifold,triplet_manifold,linearize,eta, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0) call cpu_time(end_evGW) @@ -566,8 +573,8 @@ program QuAcK if(doqsGW) then call cpu_time(start_qsGW) - call qsGW(maxSCF_GW,thresh_GW,n_diis_GW, & - COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, & + call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, & + singlet_manifold,triplet_manifold,eta, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF) call cpu_time(end_qsGW) diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/RPA.f90 index 451c166..bf9bac9 100644 --- a/src/QuAcK/RPA.f90 +++ b/src/QuAcK/RPA.f90 @@ -1,4 +1,4 @@ -subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) +subroutine RPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) ! Perform a direct random phase approximation calculation @@ -8,6 +8,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E ! Input variables + logical,intent(in) :: doACFDT logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold integer,intent(in) :: nBas @@ -31,8 +32,6 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E double precision :: rho double precision :: EcRPA(nspin) - logical :: adiabatic_connection - ! Hello world write(*,*) @@ -57,7 +56,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E call linear_response(ispin,.true.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) endif @@ -69,7 +68,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E call linear_response(ispin,.true.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) endif @@ -84,9 +83,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E ! Compute the correlation energy via the adiabatic connection - adiabatic_connection = .true. - - if(adiabatic_connection) then + if(doACFDT) then write(*,*) '------------------------------------------------------' write(*,*) 'Adiabatic connection version of RPA correlation energy' diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index bcad887..91c39fa 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,4 +1,5 @@ -subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,linearize,eta, & +subroutine evGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, & + singlet_manifold,triplet_manifold,linearize,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -13,6 +14,8 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani double precision,intent(in) :: thresh double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF + logical,intent(in) :: doACFDT + logical,intent(in) :: doXBS logical,intent(in) :: COHSEX logical,intent(in) :: SOSEX logical,intent(in) :: BSE @@ -34,7 +37,6 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Local variables - logical :: dRPA logical :: linear_mixing integer :: ispin integer :: nSCF @@ -57,9 +59,6 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rhox(:,:,:,:) - logical :: adiabatic_connection - logical :: scaled_screening - ! Hello world write(*,*) @@ -80,8 +79,6 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Switch off exchange for G0W0 - dRPA = .true. - ! Linear mixing linear_mixing = .false. @@ -89,8 +86,9 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Memory allocation - allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), & - XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & + allocate(eGW(nBas),eOld(nBas),Z(nBas),SigC(nBas),Omega(nS,nspin), & + XpY(nS,nS,nspin),XmY(nS,nS,nspin), & + rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & error_diis(nBas,max_diis),e_diis(nBas,max_diis)) ! Initialization @@ -115,7 +113,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(.not. GW0 .or. nSCF == 0) then - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI, & + call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif @@ -162,7 +160,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Print results - call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM) ! Linear mixing or DIIS extrapolation @@ -217,39 +215,8 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(BSE) then - ! Singlet manifold - - if(singlet_manifold) then - - ispin = 1 - EcBSE(ispin) = 0d0 - - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) - - endif - - ! Triplet manifold - - if(triplet_manifold) then - - ispin = 2 - EcBSE(ispin) = 0d0 - - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) - - endif + call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -262,24 +229,21 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Compute the BSE correlation energy via the adiabatic connection - adiabatic_connection = .true. - scaled_screening = .true. - - if(adiabatic_connection) then + if(doACFDT) then write(*,*) '------------------------------------------------------' write(*,*) 'Adiabatic connection version of BSE correlation energy' write(*,*) '------------------------------------------------------' write(*,*) - if(scaled_screening) then + if(doXBS) then - write(*,*) '*** scaled screening version (extended BSE) ***' + write(*,*) '*** scaled screening version (XBS) ***' write(*,*) end if - call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & + call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho) end if diff --git a/src/QuAcK/ppRPA.f90 b/src/QuAcK/ppRPA.f90 index 89baebf..0d45d81 100644 --- a/src/QuAcK/ppRPA.f90 +++ b/src/QuAcK/ppRPA.f90 @@ -1,4 +1,4 @@ -subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) +subroutine ppRPA(doACFDT,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,e) ! Perform pp-RPA calculation @@ -7,6 +7,7 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER ! Input variables + logical,intent(in) :: doACFDT logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold integer,intent(in) :: nBas diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index e539c43..375039d 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -1,7 +1,7 @@ -subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, & +subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,doXBS,COHSEX,SOSEX,BSE,TDA,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) -! Compute linear response +! Perform a quasiparticle self-consistent GW calculation implicit none include 'parameters.h' @@ -11,6 +11,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani integer,intent(in) :: maxSCF integer,intent(in) :: max_diis double precision,intent(in) :: thresh + logical,intent(in) :: doACFDT + logical,intent(in) :: doXBS logical,intent(in) :: COHSEX logical,intent(in) :: SOSEX logical,intent(in) :: BSE @@ -36,7 +38,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Local variables - logical :: dRPA integer :: nSCF integer :: nBasSq integer :: ispin @@ -69,9 +70,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani double precision,allocatable :: Z(:) double precision,allocatable :: error(:,:) - logical :: adiabatic_connection - logical :: scaled_screening - ! Hello world write(*,*) @@ -99,10 +97,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(COHSEX) write(*,*) 'COHSEX approximation activated!' write(*,*) -! Switch off exchange for qsGW - - dRPA = .true. - ! Memory allocation allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & @@ -144,7 +138,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(.not. GW0 .or. nSCF == 0) then - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, & + call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, & rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif @@ -209,7 +203,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Print results - call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) call print_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW) ! Increment @@ -251,37 +245,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani if(BSE) then - ! Singlet manifold - if(singlet_manifold) then - - ispin = 1 - EcBSE(ispin) = 0d0 - - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) - - endif - - ! Triplet manifold - if(triplet_manifold) then - - ispin = 2 - EcBSE(ispin) = 0d0 - - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin)) - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) - - endif + call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -294,24 +259,21 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Compute the BSE correlation energy via the adiabatic connection - adiabatic_connection = .true. - scaled_screening = .true. - - if(adiabatic_connection) then + if(doACFDT) then write(*,*) '------------------------------------------------------' write(*,*) 'Adiabatic connection version of BSE correlation energy' write(*,*) '------------------------------------------------------' write(*,*) - if(scaled_screening) then + if(doXBS) then - write(*,*) '*** scaled screening version (extended BSE) ***' + write(*,*) '*** scaled screening version (XBS) ***' write(*,*) end if - call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & + call ACFDT(doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho) end if diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 43d1455..0f67374 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -1,8 +1,10 @@ -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,renormalization, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, & +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,renormalization, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, & + COHSEX,SOSEX,BSE,TDA,G0W,GW0,linearize,eta, & + doACFDT,doXBS, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Read desired methods @@ -45,6 +47,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: linearize double precision,intent(out) :: eta + logical,intent(out) :: doACFDT + logical,intent(out) :: doXBS + integer,intent(out) :: nMC integer,intent(out) :: nEq integer,intent(out) :: nWalk @@ -55,7 +60,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t ! Local variables - character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7,answer8 + character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7,answer8,answer9 ! Open file with method specification @@ -150,6 +155,17 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t if(answer8 == 'T') linearize = .true. if(.not.DIIS_GW) n_diis_GW = 1 +! Options for adiabatic connection + + doACFDT = .false. + doXBS = .false. + + read(1,*) + read(1,*) answer1,answer2 + + if(answer1 == 'T') doACFDT = .true. + if(answer2 == 'T') doXBS = .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