From c6ec8f91e964a641a3295d53733e8141bafdca03 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 13 Jan 2020 23:08:03 +0100 Subject: [PATCH] RPA and ACRPA --- input/basis | 21 ++--- input/methods | 8 +- input/weight | 21 ++--- src/QuAcK/Ec_AC.f90 | 34 +++---- src/QuAcK/G0W0.f90 | 16 +++- src/QuAcK/QuAcK.f90 | 37 ++++++-- src/QuAcK/RPA.f90 | 185 +++++++++++++++++++++++++++++++++++++ src/QuAcK/TDHF.f90 | 106 +++++++++++++++++++-- src/QuAcK/read_methods.f90 | 29 +++--- 9 files changed, 377 insertions(+), 80 deletions(-) create mode 100644 src/QuAcK/RPA.f90 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 6379c0b..7530736 100644 --- a/input/methods +++ b/input/methods @@ -3,13 +3,13 @@ # MP2 MP3 MP2-F12 F F F # CCD CCSD CCSD(T) - F F F -# CIS TDHF ppRPA ADC - F F F F + F F T +# CIS RPA TDHF ppRPA ADC + F T T F F # GF2 GF3 F F # G0W0 evGW qsGW - T F F + F F F # G0T0 evGT qsGT F F F # MCMP2 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/Ec_AC.f90 b/src/QuAcK/Ec_AC.f90 index 21c94ce..a4aa102 100644 --- a/src/QuAcK/Ec_AC.f90 +++ b/src/QuAcK/Ec_AC.f90 @@ -1,4 +1,4 @@ -subroutine Ec_AC(ispin,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC) +subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC) ! Compute the correlation energy via the adiabatic connection formula @@ -8,6 +8,7 @@ subroutine Ec_AC(ispin,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC) ! Input variables integer,intent(in) :: ispin + logical,intent(in) :: dRPA integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: XpY(nS,nS) @@ -17,7 +18,10 @@ subroutine Ec_AC(ispin,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC) integer :: i,j,a,b integer :: ia,jb,kc double precision :: delta_spin + double precision :: delta_dRPA double precision,allocatable :: P(:,:) + double precision,allocatable :: V(:,:) + double precision,external :: trace_matrix ! Output variables @@ -29,30 +33,24 @@ subroutine Ec_AC(ispin,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC) if(ispin == 1) delta_spin = +1d0 if(ispin == 2) delta_spin = -1d0 +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + ! Memory allocation - allocate(P(nS,nS)) + allocate(P(nS,nS),V(nS,nS)) ! Compute P = (X+Y)(X+Y) - 1 - P(:,:) = 0d0 + P(:,:) = matmul(transpose(XpY),XpY) do ia=1,nS - do jb=1,nS - do kc=1,nS - - P(ia,jb) = P(ia,jb) + XpY(ia,kc)*XpY(kc,jb) - - enddo - enddo - P(ia,ia) = P(ia,ia) - 1d0 - enddo -! Compute Tr[VP] - - EcAC = 0d0 +! Compute Viajb = (ia|bj) ia = 0 do i=nC+1,nO @@ -63,12 +61,16 @@ subroutine Ec_AC(ispin,nBas,nC,nO,nV,nR,nS,ERI,XpY,EcAC) do b=nO+1,nBas-nR jb = jb + 1 - EcAC = EcAC + (1d0 + delta_spin)*ERI(i,b,a,j)*P(jb,ia) + V(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) enddo enddo enddo enddo +! Compute Tr(VP) + + EcAC = trace_matrix(nS,matmul(V,P)) + end subroutine Ec_AC diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 152fed6..2d6df3a 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -191,10 +191,14 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & 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)) + 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)) + rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) - call Ec_AC(ispin,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) + call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin) @@ -225,10 +229,14 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & 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)) + 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)) + rho(:,:,:,ispin),EcACBSE(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) - call Ec_AC(ispin,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) + call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),EcAC(iAC,ispin)) write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcACBSE(iAC,ispin),EcAC(iAC,ispin) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index c4b6976..2df031c 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -7,7 +7,8 @@ program QuAcK logical :: doRHF,doUHF,doMOM logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doCCSD,doCCSDT - logical :: doCIS,doTDHF,doppRPA,doADC + logical :: doCIS,doRPA,doTDHF + logical :: doppRPA,doADC logical :: doGF2,doGF3 logical :: doG0W0,doevGW,doqsGW logical :: doG0T0,doevGT,doqsGT @@ -44,6 +45,7 @@ program QuAcK double precision :: start_CCD ,end_CCD ,t_CCD double precision :: start_CCSD ,end_CCSD ,t_CCSD double precision :: start_CIS ,end_CIS ,t_CIS + double precision :: start_RPA ,end_RPA ,t_RPA double precision :: start_TDHF ,end_TDHF ,t_TDHF double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA double precision :: start_ADC ,end_ADC ,t_ADC @@ -106,13 +108,14 @@ program QuAcK ! Which calculations do you want to do? - call read_methods(doRHF,doUHF,doMOM, & - doMP2,doMP3,doMP2F12, & - doCCD,doCCSD,doCCSDT, & - doCIS,doTDHF,doppRPA,doADC, & - doGF2,doGF3, & - doG0W0,doevGW,doqsGW, & - doG0T0,doevGT,doqsGT, & + call read_methods(doRHF,doUHF,doMOM, & + doMP2,doMP3,doMP2F12, & + doCCD,doCCSD,doCCSDT, & + doCIS,doRPA,doTDHF, & + doppRPA,doADC, & + doGF2,doGF3, & + doG0W0,doevGW,doqsGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read options for methods @@ -363,6 +366,8 @@ program QuAcK ! Perform CCSD or CCSD(T) calculation !------------------------------------------------------------------------ + if(doCCSDT) doCCSD = .true. + if(doCCSD) then call cpu_time(start_CCSD) @@ -391,6 +396,22 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Compute (direct) RPA excitations +!------------------------------------------------------------------------ + + 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 cpu_time(end_RPA) + + t_RPA = end_RPA - start_RPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPA = ',t_RPA,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Compute TDHF excitations !------------------------------------------------------------------------ diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/RPA.f90 new file mode 100644 index 0000000..4b1c561 --- /dev/null +++ b/src/QuAcK/RPA.f90 @@ -0,0 +1,185 @@ +subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) + +! Perform a direct random phase approximation calculation + + implicit none + include 'parameters.h' + include 'quadrature.h' + +! Input variables + + logical,intent(in) :: singlet_manifold + logical,intent(in) :: triplet_manifold + 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) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + logical :: dRPA + logical :: TDA + logical :: BSE + integer :: ispin + double precision,allocatable :: Omega(:,:) + double precision,allocatable :: XpY(:,:,:) + + double precision :: rho + double precision :: EcRPA(nspin) + + logical :: AC + integer :: iAC + double precision :: lambda + double precision,allocatable :: EcACRPA(:,:) + double precision,allocatable :: EcAC(:,:) + +! Hello world + + write(*,*) + write(*,*)'***********************************************' + write(*,*)'| random-phase approximation calculation |' + write(*,*)'***********************************************' + write(*,*) + +! Initialization + + 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)) + + 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, & + EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + + endif + +! Triplet manifold + + if(triplet_manifold) then + + ispin = 2 + + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & + EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + call print_excitation('RPA ',ispin,nS,Omega(:,ispin)) + + endif + + 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(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the correlation energy via the adiabatic connection + + if(AC) then + + write(*,*) '------------------------------------------------------' + write(*,*) 'Adiabatic connection version of RPA correlation energy' + write(*,*) '------------------------------------------------------' + write(*,*) + + if(singlet_manifold) then + + 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)) + + call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,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)) + + call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,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 subroutine RPA diff --git a/src/QuAcK/TDHF.f90 b/src/QuAcK/TDHF.f90 index 5a2b085..50f22e2 100644 --- a/src/QuAcK/TDHF.f90 +++ b/src/QuAcK/TDHF.f90 @@ -4,6 +4,7 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, implicit none include 'parameters.h' + include 'quadrature.h' ! Input variables @@ -26,13 +27,18 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, logical :: TDA logical :: BSE integer :: ispin - double precision :: lambda double precision,allocatable :: Omega(:,:) double precision,allocatable :: XpY(:,:,:) double precision :: rho double precision :: EcRPA(nspin) + logical :: AC + integer :: iAC + double precision :: lambda + double precision,allocatable :: EcACRPA(:,:) + double precision,allocatable :: EcAC(:,:) + ! Hello world write(*,*) @@ -45,10 +51,6 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, EcRPA(:) = 0d0 -! Adiabatic connection scaling - - lambda = 1d0 - ! Switch on exchange for TDHF dRPA = .false. @@ -65,13 +67,16 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, allocate(Omega(nS,nspin),XpY(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,lambda,e,ERI,rho, & + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('TDHF ',ispin,nS,Omega(:,ispin)) @@ -83,7 +88,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,lambda,e,ERI,rho, & + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('TDHF ',ispin,nS,Omega(:,ispin)) @@ -98,4 +103,91 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, write(*,*)'-------------------------------------------------------------------------------' write(*,*) +! Compute the correlation energy via the adiabatic connection + + if(AC) then + + write(*,*) '------------------------------------------------------' + write(*,*) 'Adiabatic connection version of RPA correlation energy' + write(*,*) '------------------------------------------------------' + write(*,*) + + if(singlet_manifold) then + + 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,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, & +! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,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,e,ERI, & + rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) + + call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,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,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI, & +! rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,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,e,ERI, & + rho,EcACRPA(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin)) + + call Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,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 subroutine TDHF diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index a8ea50e..602f0b1 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -1,10 +1,11 @@ -subroutine read_methods(doRHF,doUHF,doMOM, & - doMP2,doMP3,doMP2F12, & - doCCD,doCCSD,doCCSDT, & - doCIS,doTDHF,doppRPA,doADC, & - doGF2,doGF3, & - doG0W0,doevGW,doqsGW, & - doG0T0,doevGT,doqsGT, & +subroutine read_methods(doRHF,doUHF,doMOM, & + doMP2,doMP3,doMP2F12, & + doCCD,doCCSD,doCCSDT, & + doCIS,doRPA,doTDHF, & + doppRPA,doADC, & + doGF2,doGF3, & + doG0W0,doevGW,doqsGW, & + doG0T0,doevGT,doqsGT, & doMCMP2) ! Read desired methods @@ -16,7 +17,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & logical,intent(out) :: doRHF,doUHF,doMOM logical,intent(out) :: doMP2,doMP3,doMP2F12 logical,intent(out) :: doCCD,doCCSD,doCCSDT - logical,intent(out) :: doCIS,doTDHF,doppRPA,doADC + logical,intent(out) :: doCIS,doRPA,doTDHF,doppRPA,doADC logical,intent(out) :: doGF2,doGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW logical,intent(out) :: doG0T0,doevGT,doqsGT @@ -24,7 +25,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & ! Local variables - character(len=1) :: answer1,answer2,answer3,answer4 + character(len=1) :: answer1,answer2,answer3,answer4,answer5 ! Open file with method specification @@ -45,6 +46,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & doCCSDT = .false. doCIS = .false. + doRPA = .false. doTDHF = .false. doppRPA = .false. doADC = .false. @@ -89,11 +91,12 @@ subroutine read_methods(doRHF,doUHF,doMOM, & ! Read excited state methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4 + read(1,*) answer1,answer2,answer3,answer4,answer5 if(answer1 == 'T') doCIS = .true. - if(answer2 == 'T') doTDHF = .true. - if(answer3 == 'T') doppRPA = .true. - if(answer4 == 'T') doADC = .true. + if(answer2 == 'T') doRPA = .true. + if(answer3 == 'T') doTDHF = .true. + if(answer4 == 'T') doppRPA = .true. + if(answer5 == 'T') doADC = .true. ! Read Green function methods