From 494204113b784bf8c87e00bc2eff358fd337543c Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 22 Oct 2021 21:10:31 +0200 Subject: [PATCH] new version of SOSEX (need serious debugging) --- include/parameters.h | 2 +- input/dft | 15 +- input/methods | 4 +- input/options | 4 +- src/MBPT/G0W0.f90 | 10 +- src/MBPT/G0W0_SOSEX.f90 | 196 ++++++++++++++++++ src/MBPT/evGW.f90 | 10 +- src/MBPT/excitation_density_SOSEX.f90 | 52 ++++- src/MBPT/print_SOSEX.f90 | 55 +++++ src/MBPT/qsGW.f90 | 10 +- src/MBPT/qsUGW.f90 | 10 +- src/MBPT/renormalization_factor_SOSEX.f90 | 73 +++++++ .../self_energy_correlation_SOSEX_diag.f90 | 80 +++++++ src/QuAcK/QuAcK.f90 | 12 +- src/eDFT/read_options_dft.f90 | 23 +- src/eDFT/unrestricted_auxiliary_energy.f90 | 4 - ...stricted_correlation_individual_energy.f90 | 1 - 17 files changed, 483 insertions(+), 78 deletions(-) create mode 100644 src/MBPT/G0W0_SOSEX.f90 create mode 100644 src/MBPT/print_SOSEX.f90 create mode 100644 src/MBPT/renormalization_factor_SOSEX.f90 create mode 100644 src/MBPT/self_energy_correlation_SOSEX_diag.f90 diff --git a/include/parameters.h b/include/parameters.h index d0e1202..89c2c93 100644 --- a/include/parameters.h +++ b/include/parameters.h @@ -1,7 +1,7 @@ integer,parameter :: ncart = 3 integer,parameter :: nspin = 2 integer,parameter :: nsp = 3 - integer,parameter :: maxEns = 3 + integer,parameter :: maxEns = 4 integer,parameter :: maxShell = 512 integer,parameter :: maxL = 7 integer,parameter :: n1eInt = 3 diff --git a/input/dft b/input/dft index c95d3f4..796cb10 100644 --- a/input/dft +++ b/input/dft @@ -6,31 +6,34 @@ # GGA = 2: B88,G96,PBE # MGGA = 3: # Hybrid = 4: HF,B3LYP,PBE - 4 HF + 1 S51 # correlation rung: # Hartree = 0: H # LDA = 1: PW92,VWN3,VWN5,eVWN5 # GGA = 2: LYP,PBE # MGGA = 3: # Hybrid = 4: HF,B3LYP,PBE - 0 H + 1 VWN5 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) - 1 + 4 # occupation numbers 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + + 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.0 0.0 + 0.25 0.25 0.25 # N-centered? - F + T # Parameters for CC weight-dependent exchange functional 0.0 0.0 0.0 0.0 0.0 0.0 diff --git a/input/methods b/input/methods index ddcffcb..a05c2f1 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 F + T F F F F # G0T0 evGT qsGT - T F F + F F F # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index 061ca17..31e9ef5 100644 --- a/input/options +++ b/input/options @@ -9,10 +9,10 @@ # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 - 256 0.00001 T 5 T 0.0 F F F F F + 256 0.00001 T 5 T 0.0 F T F F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - T T T F + F T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/MBPT/G0W0.f90 b/src/MBPT/G0W0.f90 index 2f3f2fd..9065661 100644 --- a/src/MBPT/G0W0.f90 +++ b/src/MBPT/G0W0.f90 @@ -1,4 +1,4 @@ -subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & +subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0) @@ -14,7 +14,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA @@ -76,13 +75,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & EcRPA = 0d0 -! SOSEX correction - - if(SOSEX) then - write(*,*) 'SOSEX correction activated but BUG!' - stop - end if - ! COHSEX approximation if(COHSEX) then diff --git a/src/MBPT/G0W0_SOSEX.f90 b/src/MBPT/G0W0_SOSEX.f90 new file mode 100644 index 0000000..fa1eeb5 --- /dev/null +++ b/src/MBPT/G0W0_SOSEX.f90 @@ -0,0 +1,196 @@ +subroutine G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eSOSEX) + +! Perform the SOSEX extension of G0W0 + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: 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 + logical,intent(in) :: triplet + 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) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: Vxc(nBas) + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: PHF(nBas,nBas) + +! Local variables + + logical :: print_W = .true. + integer :: ispin + double precision :: EcRPA(nspin) + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: EcGM + double precision,allocatable :: SigX(:) + double precision,allocatable :: SigC(:) + double precision,allocatable :: Z(:) + double precision,allocatable :: OmRPA(:,:) + double precision,allocatable :: XpY_RPA(:,:,:) + double precision,allocatable :: XmY_RPA(:,:,:) + double precision,allocatable :: rho_RPA(:,:,:,:) + +! Output variables + + double precision :: eSOSEX(nBas) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| One-shot SOSEX calculation |' + write(*,*)'************************************************' + write(*,*) + +! Initialization + + EcRPA = 0d0 + +! TDA for W + + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Memory allocation + + allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS,nspin),XpY_RPA(nS,nS,nspin),XmY_RPA(nS,nS,nspin),rho_RPA(nBas,nBas,nS,nspin)) + +!-------------------! +! Compute screening ! +!-------------------! + + do ispin=1,nspin + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO, & + OmRPA(:,ispin),rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) + if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) + end do + +!--------------------------! +! Compute spectral weights ! +!--------------------------! + + call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) + +!------------------------! +! Compute GW self-energy ! +!------------------------! + + call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) + call self_energy_correlation_SOSEX_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + +!--------------------------------! +! Compute renormalization factor ! +!--------------------------------! + + call renormalization_factor_SOSEX(eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + +!-----------------------------------! +! Solve the quasi-particle equation ! +!-----------------------------------! + + eSOSEX(:) = eHF(:) + Z(:)*(SigX(:) + SigC(:) - Vxc(:)) + +! Compute the RPA correlation energy + + do ispin=1,nspin + call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eSOSEX,ERI_MO,OmRPA(:,ispin), & + rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) + end do + +!--------------! +! Dump results ! +!--------------! + + call print_SOSEX(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eSOSEX,EcRPA,EcGM) + +! Deallocate memory + + deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA) + +! Perform BSE calculation + + 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,eSOSEX,EcBSE) + + if(exchange_kernel) then + + EcBSE(1) = 0.5d0*EcBSE(1) + EcBSE(2) = 1.5d0*EcBSE(2) + + end if + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@BSE@SOSEX correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@SOSEX correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@SOSEX correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@SOSEX total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + + if(doACFDT) then + + write(*,*) '--------------------------------------------------------------' + write(*,*) ' Adiabatic connection version of BSE@SOSEX correlation energy ' + write(*,*) '--------------------------------------------------------------' + write(*,*) + + if(doXBS) then + + write(*,*) '*** scaled screening version (XBS) ***' + write(*,*) + + end if + + call ACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eSOSEX,EcAC) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'AC@BSE@SOSEX correlation energy (singlet) =',EcAC(1) + write(*,'(2X,A50,F20.10)') 'AC@BSE@SOSEX correlation energy (triplet) =',EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@SOSEX correlation energy =',EcAC(1) + EcAC(2) + write(*,'(2X,A50,F20.10)') 'AC@BSE@SOSEX total energy =',ENuc + ERHF + EcAC(1) + EcAC(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + + end if + +end subroutine G0W0_SOSEX diff --git a/src/MBPT/evGW.f90 b/src/MBPT/evGW.f90 index 56aeef6..9a74734 100644 --- a/src/MBPT/evGW.f90 +++ b/src/MBPT/evGW.f90 @@ -1,4 +1,4 @@ -subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & +subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0) @@ -18,7 +18,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA @@ -78,13 +77,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE write(*,*)'************************************************' write(*,*) -! SOSEX correction - - if(SOSEX) then - write(*,*) 'SOSEX correction activated but BUG!' - stop - end if - ! COHSEX approximation if(COHSEX) then diff --git a/src/MBPT/excitation_density_SOSEX.f90 b/src/MBPT/excitation_density_SOSEX.f90 index 6d78f48..651ef58 100644 --- a/src/MBPT/excitation_density_SOSEX.f90 +++ b/src/MBPT/excitation_density_SOSEX.f90 @@ -1,31 +1,63 @@ -subroutine excitation_density_SOSEX(nBas,nC,nO,nR,nS,G,XpY,rho) +subroutine excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI,XpY,rho) -! Compute excitation densities +! Compute excitation densities for SOSEX implicit none + include 'parameters.h' ! Input variables - integer,intent(in) :: nBas,nC,nO,nR,nS - double precision,intent(in) :: G(nBas,nBas,nBas,nBas),XpY(nS,nS) + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: XpY(nS,nS,nspin) ! Local variables - integer :: ia,jb,x,y,j,b + integer :: ispin + integer :: p,q + integer :: i,a,j,b + integer :: ia,jb ! Output variables - double precision,intent(out) :: rho(nBas,nBas,nS) + double precision,intent(out) :: rho(nBas,nBas,nS,nspin) + + rho(:,:,:,:) = 0d0 + +! Singlet part + + ispin = 1 - rho(:,:,:) = 0d0 do ia=1,nS - do x=nC+1,nBas-nR - do y=nC+1,nBas-nR + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR jb = 0 do j=nC+1,nO do b=nO+1,nBas-nR jb = jb + 1 - rho(x,y,ia) = rho(x,y,ia) + G(x,y,b,j)*XpY(ia,jb) + rho(p,q,ia,ispin) = rho(p,q,ia,ispin) + ERI(p,j,q,b)*XpY(ia,jb,ispin) + enddo + enddo + enddo + enddo + enddo + +! Triplet part + + ispin = 2 + + do ia=1,nS + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + rho(p,q,ia,ispin) = rho(p,q,ia,ispin) + (ERI(p,j,q,b) - ERI(p,j,b,q))*XpY(ia,jb,ispin) enddo enddo enddo diff --git a/src/MBPT/print_SOSEX.f90 b/src/MBPT/print_SOSEX.f90 new file mode 100644 index 0000000..ffcbd93 --- /dev/null +++ b/src/MBPT/print_SOSEX.f90 @@ -0,0 +1,55 @@ +subroutine print_SOSEX(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eSOSEX,EcRPA,EcGM) + +! Print one-electron energies and other stuff for SOSEX + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas,nO + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: EcRPA + double precision,intent(in) :: EcGM + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: SigC(nBas) + double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: eSOSEX(nBas) + + integer :: p,HOMO,LUMO + double precision :: Gap + +! HOMO and LUMO + + HOMO = nO + LUMO = HOMO + 1 + Gap = eSOSEX(LUMO)-eSOSEX(HOMO) + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + write(*,*)' One-shot SOSEX calculation ' + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & + '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do p=1,nBas + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',p,'|',eHF(p)*HaToeV,'|',SigC(p)*HaToeV,'|',Z(p),'|',eSOSEX(p)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') 'SOSEX HOMO energy:',eSOSEX(HOMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'SOSEX LUMO energy:',eSOSEX(LUMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'SOSEX HOMO-LUMO gap :',Gap*HaToeV,' eV' + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') 'RPA@SOSEX total energy :',ENuc + ERHF + EcRPA,' au' + write(*,'(2X,A30,F15.6,A3)') 'RPA@SOSEX correlation energy:',EcRPA,' au' + write(*,'(2X,A30,F15.6,A3)') 'GM@SOSEX total energy :',ENuc + ERHF + EcGM,' au' + write(*,'(2X,A30,F15.6,A3)') 'GM@SOSEX correlation energy:',EcGM,' au' + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +end subroutine print_SOSEX + + diff --git a/src/MBPT/qsGW.f90 b/src/MBPT/qsGW.f90 index 8281d25..d6bc2bb 100644 --- a/src/MBPT/qsGW.f90 +++ b/src/MBPT/qsGW.f90 @@ -1,4 +1,4 @@ -subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & +subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,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) @@ -16,7 +16,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA @@ -113,13 +112,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE nBasSq = nBas*nBas -! SOSEX correction - - if(SOSEX) then - write(*,*) 'SOSEX correction activated but BUG!' - stop - end if - ! COHSEX approximation if(COHSEX) then diff --git a/src/MBPT/qsUGW.f90 b/src/MBPT/qsUGW.f90 index 0718811..b9d055d 100644 --- a/src/MBPT/qsUGW.f90 +++ b/src/MBPT/qsUGW.f90 @@ -1,4 +1,4 @@ -subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & +subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,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) @@ -17,7 +17,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS logical,intent(in) :: COHSEX - logical,intent(in) :: SOSEX logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA @@ -121,13 +120,6 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOS nBasSq = nBas*nBas -! SOSEX correction - - if(SOSEX) then - write(*,*) 'SOSEX correction activated but BUG!' - stop - end if - ! COHSEX approximation if(COHSEX) then diff --git a/src/MBPT/renormalization_factor_SOSEX.f90 b/src/MBPT/renormalization_factor_SOSEX.f90 new file mode 100644 index 0000000..dea2d08 --- /dev/null +++ b/src/MBPT/renormalization_factor_SOSEX.f90 @@ -0,0 +1,73 @@ +subroutine renormalization_factor_SOSEX(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) + +! Compute renormalization factor for the SOSEX version of GW + + 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) :: e(nBas) + double precision,intent(in) :: Omega(nS,nspin) + double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + +! Local variables + + integer :: ispin + integer :: i,j,a,b,p,jb + double precision :: eps + +! Output variables + + double precision,intent(out) :: Z(nBas) + +! Initialize + + Z(:) = 0d0 + + ! Occupied part of the correlation self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + do ispin=1,nspin + eps = e(p) - e(i) + Omega(jb,ispin) + Z(p) = Z(p) - 2d0*rho(p,i,jb,ispin)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + do ispin=1,nspin + eps = e(p) - e(a) - Omega(jb,ispin) + Z(p) = Z(p) - 2d0*rho(p,a,jb,ispin)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + end do + end do + +! Compute renormalization factor from derivative of SigC + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine renormalization_factor_SOSEX diff --git a/src/MBPT/self_energy_correlation_SOSEX_diag.f90 b/src/MBPT/self_energy_correlation_SOSEX_diag.f90 new file mode 100644 index 0000000..a5ea6ff --- /dev/null +++ b/src/MBPT/self_energy_correlation_SOSEX_diag.f90 @@ -0,0 +1,80 @@ +subroutine self_energy_correlation_SOSEX_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) + +! Compute diagonal of the correlation part of the self-energy + + 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) :: e(nBas) + double precision,intent(in) :: Omega(nS,nspin) + double precision,intent(in) :: rho(nBas,nBas,nS,nspin) + +! Local variables + + integer :: ispin + integer :: i,a,p,q,jb + double precision :: eps + +! Output variables + + double precision,intent(out) :: SigC(nBas) + double precision,intent(out) :: EcGM + +! Initialize + + SigC(:) = 0d0 + +!----------------------------- +! SOSEX self-energy +!----------------------------- + + ! Occupied part of the correlation self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do jb=1,nS + do ispin=1,nspin + eps = e(p) - e(i) + Omega(jb,ispin) + SigC(p) = SigC(p) + 2d0*rho(p,i,jb,ispin)**2*eps/(eps**2 + eta**2) + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do jb=1,nS + do ispin=1,nspin + eps = e(p) - e(a) - Omega(jb,ispin) + SigC(p) = SigC(p) + 2d0*rho(p,a,jb,ispin)**2*eps/(eps**2 + eta**2) + end do + end do + end do + end do + + ! GM correlation energy + + EcGM = 0d0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do jb=1,nS + do ispin=1,nspin + eps = e(a) - e(i) + Omega(jb,ispin) + EcGM = EcGM - 4d0*rho(a,i,jb,ispin)*rho(a,i,jb,ispin)*eps/(eps**2 + eta**2) + end do + end do + end do + end do + +end subroutine self_energy_correlation_SOSEX_diag diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 8758535..aad3ba8 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -953,8 +953,10 @@ program QuAcK dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0) else - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) +! call G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,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) ! call ufBSE(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) @@ -984,7 +986,7 @@ program QuAcK else - call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & + call evGW(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, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) end if @@ -1006,14 +1008,14 @@ program QuAcK if(unrestricted) then - call qsUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & + call qsUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & nV,nR,nS,EUHF,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,eHF) else - call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & + 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,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) @@ -1220,7 +1222,7 @@ program QuAcK ! ! Long-range G0W0 calculation ! call cpu_time(start_G0W0) -! call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & +! call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & ! dBSE,dTDA,evDyn,singlet,triplet,linGW,eta_GW, & ! nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,dipole_int_MO,eHF,eG0W0) ! call cpu_time(end_G0W0) diff --git a/src/eDFT/read_options_dft.f90 b/src/eDFT/read_options_dft.f90 index 1d052ea..4352a21 100644 --- a/src/eDFT/read_options_dft.f90 +++ b/src/eDFT/read_options_dft.f90 @@ -22,7 +22,7 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, character(len=8),intent(out) :: method integer,intent(out) :: x_rung,c_rung - character(len=12),intent(out) :: x_DFA, c_DFA + character(len=12),intent(out) :: x_DFA,c_DFA integer,intent(out) :: SGn integer,intent(out) :: nEns logical,intent(out) :: doNcentered @@ -39,11 +39,11 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, ! Default values - method = 'GOK-RKS' + method = 'eDFT-UKS' x_rung = 1 c_rung = 1 - x_DFA = 'RS51' - c_DFA = 'RVWN5' + x_DFA = 'S51' + c_DFA = 'VWN5' SGn = 0 wEns(:) = 0d0 @@ -134,15 +134,16 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, if(answer == 'T') doNcentered = .true. + wEns(1) = 1d0 + do iEns=2,nEns + wEns(1) = wEns(1) - wEns(iEns) + end do + if (doNcentered) then - wEns(1) = 1d0 - wEns(2) - wEns(3) - wEns(2) = (nEl(1)/nEl(2))*wEns(2) - wEns(3) = (nEl(1)/nEl(3))*wEns(3) - - else - - wEns(1) = 1d0 - wEns(2) - wEns(3) + do iEns=2,nEns + wEns(iEns) = (nEl(1)/nEl(iEns))*wEns(iEns) + end do end if diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/unrestricted_auxiliary_energy.f90 index a5f2b9a..cb7c987 100644 --- a/src/eDFT/unrestricted_auxiliary_energy.f90 +++ b/src/eDFT/unrestricted_auxiliary_energy.f90 @@ -46,10 +46,6 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux) end do end do -! if (doNcentered .NE. 0) then -! Eaux(:,iEns) = Eaux(:,iEns)*(nEl(iEns)/nEl(1)) -! end if - end do end subroutine unrestricted_auxiliary_energy diff --git a/src/eDFT/unrestricted_correlation_individual_energy.f90 b/src/eDFT/unrestricted_correlation_individual_energy.f90 index 4e63d1b..53814a6 100644 --- a/src/eDFT/unrestricted_correlation_individual_energy.f90 +++ b/src/eDFT/unrestricted_correlation_individual_energy.f90 @@ -26,7 +26,6 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns double precision :: EcLDA(nsp) double precision :: EcGGA(nsp) - double precision :: aC ! Output variables