From 3ebcd29b197aec3140631be19ffa96add3ee95d0 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 14 Jan 2020 16:36:11 +0100 Subject: [PATCH] GW AC --- input/basis | 21 ++++--- input/methods | 4 +- input/weight | 21 ++++--- src/QuAcK/G0W0.f90 | 114 +++++++----------------------------- src/QuAcK/RPA.f90 | 121 ++++++--------------------------------- src/QuAcK/TDHF.f90 | 119 ++++++-------------------------------- src/QuAcK/evGW.f90 | 48 ++++++++++++---- src/QuAcK/ladder_CCD.f90 | 2 +- src/QuAcK/qsGW.f90 | 76 ++++++++++++++++-------- 9 files changed, 176 insertions(+), 350 deletions(-) diff --git a/input/basis b/input/basis index b9ca7b5..4c61da7 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,16 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 6 +S 4 1.00 + 234.0000000 0.0025870 + 35.1600000 0.0195330 + 7.9890000 0.0909980 + 2.2120000 0.2720500 S 1 1.00 - 0.2976000 1.0000000 + 0.6669000 1.0000000 +S 1 1.00 + 0.2089000 1.0000000 P 1 1.00 - 1.2750000 1.0000000 + 3.0440000 1.0000000 +P 1 1.00 + 0.7580000 1.0000000 +D 1 1.00 + 1.9650000 1.0000000 diff --git a/input/methods b/input/methods index 9e9a79e..02bb6a3 100644 --- a/input/methods +++ b/input/methods @@ -3,9 +3,9 @@ # MP2 MP3 MP2-F12 F F F # CCD CCSD CCSD(T) ringCCD ladderCCD - F F F T T + F F F F F # CIS RPA TDHF ppRPA ADC - F T T T F + F T T F F # GF2 GF3 F F # G0W0 evGW qsGW diff --git a/input/weight b/input/weight index b9ca7b5..4c61da7 100644 --- a/input/weight +++ b/input/weight @@ -1,9 +1,16 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 6 +S 4 1.00 + 234.0000000 0.0025870 + 35.1600000 0.0195330 + 7.9890000 0.0909980 + 2.2120000 0.2720500 S 1 1.00 - 0.2976000 1.0000000 + 0.6669000 1.0000000 +S 1 1.00 + 0.2089000 1.0000000 P 1 1.00 - 1.2750000 1.0000000 + 3.0440000 1.0000000 +P 1 1.00 + 0.7580000 1.0000000 +D 1 1.00 + 1.9650000 1.0000000 diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 0343877..cdf110f 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -29,8 +29,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & ! Local variables - logical :: dRPA - integer :: ispin,jspin + integer :: ispin double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcGM @@ -42,11 +41,8 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rhox(:,:,:,:) - logical :: AC - integer :: iAC - double precision :: lambda - double precision,allocatable :: EcACBSE(:,:) - double precision,allocatable :: EcAC(:,:) + logical :: adiabatic_connection + logical :: scaled_screening ! Output variables @@ -70,10 +66,6 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & if(COHSEX) write(*,*) 'COHSEX approximation activated!' write(*,*) -! Switch off exchange for G0W0 - - dRPA = .true. - ! Spin manifold ispin = 1 @@ -83,12 +75,9 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & allocate(SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), & rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin)) - AC = .true. - allocate(EcACBSE(nAC,nspin),EcAC(nAC,nspin)) - ! Compute linear response - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & + 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)) ! Compute correlation part of the self-energy @@ -129,11 +118,11 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & ispin = 1 EcBSE(ispin) = 0d0 - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & + 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,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, & + 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)) @@ -146,11 +135,11 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & ispin = 2 EcBSE(ispin) = 0d0 - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & + 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,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI, & + 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)) @@ -158,98 +147,35 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy (singlet) =',EcBSE(1) - write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy (triplet) =',EcBSE(2) - write(*,'(2X,A40,F15.6)') 'BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2) - write(*,'(2X,A40,F15.6)') 'BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) ! Compute the BSE correlation energy via the adiabatic connection - if(AC) then + adiabatic_connection = .true. + scaled_screening = .true. + + if(adiabatic_connection) then write(*,*) '------------------------------------------------------' write(*,*) 'Adiabatic connection version of BSE correlation energy' write(*,*) '------------------------------------------------------' write(*,*) - - if(singlet_manifold) then - ispin = 1 - EcACBSE(:,ispin) = 0d0 + if(scaled_screening) then - write(*,*) '--------------' - write(*,*) 'Singlet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcBSE(lambda)','Tr(V x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,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,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, & - rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin) - - end do - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(BSE) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin)) - write(*,*) '-----------------------------------------------------------------------------------' + write(*,*) '*** scaled screening version (extended BSE) ***' write(*,*) end if - - if(triplet_manifold) then - ispin = 2 - EcACBSE(:,ispin) = 0d0 + call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI,eG0W0,Omega,XpY,XmY,rho) - write(*,*) '--------------' - write(*,*) 'Triplet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcBSE(lambda)','Tr(V x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,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,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,eG0W0,ERI, & - rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin) - - end do - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(BSE) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin)) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - end if end if diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/RPA.f90 index 2a5765b..ff0a135 100644 --- a/src/QuAcK/RPA.f90 +++ b/src/QuAcK/RPA.f90 @@ -23,9 +23,6 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E ! Local variables - logical :: dRPA - logical :: TDA - logical :: BSE integer :: ispin double precision,allocatable :: Omega(:,:) double precision,allocatable :: XpY(:,:,:) @@ -34,11 +31,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E double precision :: rho double precision :: EcRPA(nspin) - logical :: AC - integer :: iAC - double precision :: lambda - double precision,allocatable :: EcACRPA(:,:) - double precision,allocatable :: EcAC(:,:) + logical :: adiabatic_connection ! Hello world @@ -52,32 +45,17 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E EcRPA(:) = 0d0 -! Switch off exchange for RPA - - dRPA = .true. - -! Switch off Tamm-Dancoff approximation for RPA - - TDA = .false. - -! Switch off Bethe-Salpeter equation for RPA - - BSE = .false. - ! Memory allocation allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) - AC = .true. - allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin)) - ! Singlet manifold if(singlet_manifold) then ispin = 1 - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + 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)) @@ -89,7 +67,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E ispin = 2 - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + 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)) @@ -97,90 +75,27 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A40,F15.6)') 'RPA@RPA correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A40,F15.6)') 'RPA@RPA correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A40,F15.6)') 'RPA@RPA correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A40,F15.6)') 'RPA@RPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) + write(*,'(2X,A40,F15.6)') 'Tr@RPA correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A40,F15.6)') 'Tr@RPA correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A40,F15.6)') 'Tr@RPA correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A40,F15.6)') 'Tr@RPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! Compute the correlation energy via the adiabatic connection +! Compute the correlation energy via the adiabatic connection - if(AC) then + adiabatic_connection = .true. - write(*,*) '------------------------------------------------------' - write(*,*) 'Adiabatic connection version of RPA correlation energy' - write(*,*) '------------------------------------------------------' - write(*,*) + if(adiabatic_connection) then + + write(*,*) '------------------------------------------------------' + write(*,*) 'Adiabatic connection version of RPA correlation energy' + write(*,*) '------------------------------------------------------' + write(*,*) - if(singlet_manifold) then + call ACDFT(.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho) - ispin = 1 - EcACRPA(:,ispin) = 0d0 - - write(*,*) '--------------' - write(*,*) 'Singlet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcRPA(lambda)','Tr(V x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, & - EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin) - - end do - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(RPA) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin)) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - - if(triplet_manifold) then - - ispin = 2 - EcACRPA(:,ispin) = 0d0 - - write(*,*) '--------------' - write(*,*) 'Triplet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcRPA(lambda)','Tr(V x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho, & - EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin) - - end do - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(RPA) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin)) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - - end if + end if end subroutine RPA diff --git a/src/QuAcK/TDHF.f90 b/src/QuAcK/TDHF.f90 index c52ecee..545f7f8 100644 --- a/src/QuAcK/TDHF.f90 +++ b/src/QuAcK/TDHF.f90 @@ -23,9 +23,6 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, ! Local variables - logical :: dRPA - logical :: TDA - logical :: BSE integer :: ispin double precision,allocatable :: Omega(:,:) double precision,allocatable :: XpY(:,:,:) @@ -34,11 +31,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, double precision :: rho double precision :: EcRPA(nspin) - logical :: AC - integer :: iAC - double precision :: lambda - double precision,allocatable :: EcACRPA(:,:) - double precision,allocatable :: EcAC(:,:) + logical :: adiabatic_connection ! Hello world @@ -52,32 +45,17 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, EcRPA(:) = 0d0 -! Switch on exchange for TDHF - - dRPA = .false. - -! Switch off Tamm-Dancoff approximation for TDHF - - TDA = .false. - -! Switch off Bethe-Salpeter equation for TDHF - - BSE = .false. - ! Memory allocation allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) - AC = .true. - allocate(EcACRPA(nAC,nspin),EcAC(nAC,nspin)) - ! Singlet manifold if(singlet_manifold) then ispin = 1 - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('TDHF ',ispin,nS,Omega(:,ispin)) @@ -89,7 +67,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, ispin = 2 - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + call linear_response(ispin,.false.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('TDHF ',ispin,nS,Omega(:,ispin)) @@ -97,90 +75,27 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A40,F15.6)') 'RPA@TDHF correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A40,F15.6)') 'RPA@TDHF correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A40,F15.6)') 'RPA@TDHF correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A40,F15.6)') 'RPA@TDHF total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) + write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy (singlet) =',EcRPA(1) + write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy (triplet) =',EcRPA(2) + write(*,'(2X,A40,F15.6)') 'Tr@TDHF correlation energy =',EcRPA(1) + EcRPA(2) + write(*,'(2X,A40,F15.6)') 'Tr@TDHF total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! Compute the correlation energy via the adiabatic connection +! Compute the correlation energy via the adiabatic connection - if(AC) then + adiabatic_connection = .true. - write(*,*) '------------------------------------------------------' - write(*,*) 'Adiabatic connection version of RPA correlation energy' - write(*,*) '------------------------------------------------------' - write(*,*) - - if(singlet_manifold) then + if(adiabatic_connection) then - ispin = 1 - EcACRPA(:,ispin) = 0d0 + write(*,*) '-------------------------------------------------------' + write(*,*) 'Adiabatic connection version of TDHF correlation energy' + write(*,*) '-------------------------------------------------------' + write(*,*) - write(*,*) '--------------' - write(*,*) 'Singlet states' - write(*,*) '--------------' - write(*,*) + call ACDFT(.false.,.false.,.false.,.false.,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcRPA(lambda)','Tr(V x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin) - - end do - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(RPA) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin)) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - - if(triplet_manifold) then - - ispin = 2 - EcACRPA(:,ispin) = 0d0 - - write(*,*) '--------------' - write(*,*) 'Triplet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','EcRPA(lambda)','Tr(V x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),EcAC(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACRPA(iAC,ispin),EcAC(iAC,ispin) - - end do - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(RPA) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin)) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - - end if + end if end subroutine TDHF diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 7dbb03f..dd307eb 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -53,9 +53,13 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani double precision,allocatable :: SigC(:) double precision,allocatable :: Omega(:,:) double precision,allocatable :: XpY(:,:,:) + double precision,allocatable :: XmY(:,:,:) double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rhox(:,:,:,:) + logical :: adiabatic_connection + logical :: scaled_screening + ! Hello world write(*,*) @@ -86,7 +90,7 @@ 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),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,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 @@ -112,7 +116,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, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) endif @@ -221,11 +225,11 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani 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)) + 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)) + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif @@ -238,24 +242,48 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani 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)) + 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)) + rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) endif write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy (singlet) =',EcBSE(1) - write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy (triplet) =',EcBSE(2) - write(*,'(2X,A40,F15.6)') 'BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) - write(*,'(2X,A40,F15.6)') 'BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) +! Compute the BSE correlation energy via the adiabatic connection + + adiabatic_connection = .true. + scaled_screening = .true. + + if(adiabatic_connection) then + + write(*,*) '------------------------------------------------------' + write(*,*) 'Adiabatic connection version of BSE correlation energy' + write(*,*) '------------------------------------------------------' + write(*,*) + + if(scaled_screening) then + + write(*,*) '*** scaled screening version (extended BSE) ***' + write(*,*) + + end if + + call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho) + + end if + endif end subroutine evGW diff --git a/src/QuAcK/ladder_CCD.f90 b/src/QuAcK/ladder_CCD.f90 index fd3adc6..43c3586 100644 --- a/src/QuAcK/ladder_CCD.f90 +++ b/src/QuAcK/ladder_CCD.f90 @@ -187,7 +187,7 @@ subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) write(*,*)' ladder-CCD energy ' write(*,*)'----------------------------------------------------' write(*,'(1X,A30,1X,F15.10)')' E(ladder-CCD) = ',ECCD - write(*,'(1X,A30,1X,F10.6)')' Ec(ladder-CCSD) = ',EcCCD + write(*,'(1X,A30,1X,F15.10)')' Ec(ladder-CCSD) = ',EcCCD write(*,*)'----------------------------------------------------' write(*,*) diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 7c34a49..eb4a3f2 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -52,11 +52,12 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani double precision,allocatable :: F_diis(:,:) double precision,allocatable :: Omega(:,:) double precision,allocatable :: XpY(:,:,:) + double precision,allocatable :: XmY(:,:,:) double precision,allocatable :: rho(:,:,:,:) double precision,allocatable :: rhox(:,:,:,:) double precision,allocatable :: c(:,:) double precision,allocatable :: cp(:,:) - double precision,allocatable :: e(:) + double precision,allocatable :: eGW(:) double precision,allocatable :: P(:,:) double precision,allocatable :: F(:,:) double precision,allocatable :: Fp(:,:) @@ -68,6 +69,9 @@ 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(*,*) @@ -101,9 +105,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Memory allocation - allocate(e(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & + allocate(eGW(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & - Omega(nS,nspin),XpY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & + Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -113,7 +117,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ispin = 1 Conv = 1d0 P(:,:) = PHF(:,:) - e(:) = eHF(:) + eGW(:) = eHF(:) c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 @@ -140,8 +144,8 @@ 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,e,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + 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)) endif @@ -160,9 +164,9 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani else - call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e, & + call self_energy_correlation(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC) - call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,e, & + call renormalization_factor(COHSEX,SOSEX,eta,nBas,nC,nO,nV,nR,nS,eGW, & Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z) endif @@ -196,7 +200,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,e) + call diagonalize_matrix(nBas,cp,eGW) c = matmul(X,cp) ! Compute new density matrix in the AO basis @@ -206,7 +210,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_qsGW(nBas,nO,nSCF,Conv,thresh,eHF,e,c,ENuc,P,T,V,Hc,J,K,F,SigCp,Z,EcRPA(ispin),EcGM,EqsGW) + 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 @@ -219,7 +223,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ! Compute second-order correction of the Hermitization error - call qsGW_PT(nBas,nC,nO,nV,nR,nS,e,SigCm) + call qsGW_PT(nBas,nC,nO,nV,nR,nS,eGW,SigCm) ! Compute the overlap between HF and GW orbitals @@ -253,12 +257,12 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ispin = 1 EcBSE(ispin) = 0d0 - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + 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,e,ERI_MO_basis, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,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 @@ -269,25 +273,49 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani ispin = 2 EcBSE(ispin) = 0d0 - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO_basis, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + 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,e,ERI_MO_basis, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,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 write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy (singlet) =',EcBSE(1) - write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy (triplet) =',EcBSE(2) - write(*,'(2X,A40,F15.6)') 'BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2) - write(*,'(2X,A40,F15.6)') 'BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A40,F15.6)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) - endif +! Compute the BSE correlation energy via the adiabatic connection + + adiabatic_connection = .true. + scaled_screening = .true. + + if(adiabatic_connection) then + + write(*,*) '------------------------------------------------------' + write(*,*) 'Adiabatic connection version of BSE correlation energy' + write(*,*) '------------------------------------------------------' + write(*,*) + + if(scaled_screening) then + + write(*,*) '*** scaled screening version (extended BSE) ***' + write(*,*) + + end if + + call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & + nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho) + + end if + + end if end subroutine qsGW