diff --git a/src/QuAcK/ACDFT.f90 b/src/QuAcK/ACDFT.f90 deleted file mode 100644 index 624d1fe..0000000 --- a/src/QuAcK/ACDFT.f90 +++ /dev/null @@ -1,127 +0,0 @@ -subroutine ACDFT(scaled_screening,dRPA,TDA,BSE,singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho) - -! Compute the correlation energy via the adiabatic connection dissipation fluctuation theorem - - implicit none - include 'parameters.h' - include 'quadrature.h' - -! Input variables - - logical,intent(in) :: scaled_screening - logical,intent(in) :: dRPA - logical,intent(in) :: TDA - logical,intent(in) :: BSE - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold - - integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: e(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 - logical :: adiabatic_connection - integer :: iAC - double precision :: lambda - double precision,allocatable :: Ec(:,:) - double precision,allocatable :: EcAC(:,:) - -! Memory allocation - - allocate(Ec(nAC,nspin),EcAC(nAC,nspin)) - - if(singlet_manifold) then - - ispin = 1 - Ec(:,ispin) = 0d0 - EcAC(:,ispin) = 0d0 - - write(*,*) '--------------' - write(*,*) 'Singlet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(V x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - if(scaled_screening) then - - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho(:,:,:,ispin),Ec(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - - end if - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho(:,:,:,ispin),Ec(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,Ec(iAC,ispin),EcAC(iAC,ispin) - - end do - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin)) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - - if(triplet_manifold) then - - ispin = 2 - Ec(:,ispin) = 0d0 - EcAC(:,ispin) = 0d0 - - write(*,*) '--------------' - write(*,*) 'Triplet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(V x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - if(scaled_screening) then - - call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho(:,:,:,ispin),Ec(iAC,ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) - - end if - - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, & - rho(:,:,:,ispin),Ec(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,Ec(iAC,ispin),EcAC(iAC,ispin) - - end do - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',0.5d0*dot_product(wAC,EcAC(:,ispin)) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - -end subroutine ACDFT diff --git a/src/QuAcK/Ec_AC.f90 b/src/QuAcK/Ec_AC.f90 deleted file mode 100644 index b7f2b7c..0000000 --- a/src/QuAcK/Ec_AC.f90 +++ /dev/null @@ -1,92 +0,0 @@ -subroutine Ec_AC(ispin,dRPA,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,EcAC) - -! Compute the correlation energy via the adiabatic connection formula - - implicit none - include 'parameters.h' - -! 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) - double precision,intent(in) :: XmY(nS,nS) - -! Local variables - - integer :: i,j,a,b - integer :: ia,jb,kc - double precision :: delta_spin - double precision :: delta_dRPA - double precision,allocatable :: P(:,:) - double precision,allocatable :: Ap(:,:) - double precision,allocatable :: Bp(:,:) - double precision,allocatable :: X(:,:) - double precision,allocatable :: Y(:,:) - double precision,external :: trace_matrix - -! Output variables - - double precision,intent(out) :: EcAC - -! Singlet or triplet manifold? - - delta_spin = 0d0 - 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),Ap(nS,nS),Bp(nS,nS),X(nS,nS),Y(nS,nS)) - -! Compute P = (X+Y)(X+Y) - 1 - - P(:,:) = matmul(transpose(XpY),XpY) - - do ia=1,nS - P(ia,ia) = P(ia,ia) - 1d0 - enddo - -! Compute Aiajb = (ia|bj) and Biajb = (ia|jb) - - ia = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR - ia = ia + 1 - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - - Ap(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) - Bp(ia,jb) = (1d0 + delta_spin)*ERI(i,j,b,a) - - enddo - enddo - enddo - enddo - -! Compute Tr(A x P) - -! EcAC = trace_matrix(nS,matmul(Ap,P)) - -! print*,'EcAC =',EcAC - - X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:)) - Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:)) - - EcAC = trace_matrix(nS,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) & - + trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) & - - trace_matrix(nS,Ap) - -! print*,'EcAC =',EcAC - -end subroutine Ec_AC - diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index cdf110f..774dcb4 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -173,7 +173,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold,eta, & end if - call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & + call ACFDT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & nBas,nC,nO,nV,nR,nS,ERI,eG0W0,Omega,XpY,XmY,rho) end if diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index ef2e2c4..44e063a 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -8,7 +8,7 @@ program QuAcK logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doCCSD,doCCSDT logical :: do_ring_CCD,do_ladder_CCD - logical :: doCIS,doRPA,doTDHF + logical :: doCIS,doRPA,doRPAx logical :: doppRPA,doADC logical :: doGF2,doGF3 logical :: doG0W0,doevGW,doqsGW @@ -47,7 +47,7 @@ program QuAcK 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_RPAx ,end_RPAx ,t_RPAx double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA double precision :: start_ADC ,end_ADC ,t_ADC double precision :: start_GF2 ,end_GF2 ,t_GF2 @@ -113,7 +113,7 @@ program QuAcK doMP2,doMP3,doMP2F12, & doCCD,doCCSD,doCCSDT, & do_ring_CCD,do_ladder_CCD, & - doCIS,doRPA,doTDHF, & + doCIS,doRPA,doRPAx, & doppRPA,doADC, & doGF2,doGF3, & doG0W0,doevGW,doqsGW, & @@ -444,17 +444,17 @@ program QuAcK end if !------------------------------------------------------------------------ -! Compute TDHF excitations +! Compute RPAx excitations !------------------------------------------------------------------------ - if(doTDHF) then + if(doRPAx) then - call cpu_time(start_TDHF) - call TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) - call cpu_time(end_TDHF) + call cpu_time(start_RPAx) + call RPAx(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) + call cpu_time(end_RPAx) - t_TDHF = end_TDHF - start_TDHF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for TDHF = ',t_TDHF,' seconds' + t_RPAx = end_RPAx - start_RPAx + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RPAx = ',t_RPAx,' seconds' write(*,*) end if diff --git a/src/QuAcK/RPA.f90 b/src/QuAcK/RPA.f90 index ff0a135..451c166 100644 --- a/src/QuAcK/RPA.f90 +++ b/src/QuAcK/RPA.f90 @@ -93,7 +93,7 @@ subroutine RPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,E write(*,*) '------------------------------------------------------' write(*,*) - call ACDFT(.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold, & + call ACFDT(.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold, & nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho) end if diff --git a/src/QuAcK/TDHF.f90 b/src/QuAcK/TDHF.f90 deleted file mode 100644 index 545f7f8..0000000 --- a/src/QuAcK/TDHF.f90 +++ /dev/null @@ -1,101 +0,0 @@ -subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) - -! Perform 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 - - integer :: ispin - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - - double precision :: rho - double precision :: EcRPA(nspin) - - logical :: adiabatic_connection - -! Hello world - - write(*,*) - write(*,*)'************************************************' - write(*,*)'| Time-dependent Hartree-Fock calculation |' - write(*,*)'************************************************' - write(*,*) - -! Initialization - - EcRPA(:) = 0d0 - -! Memory allocation - - allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) - -! Singlet manifold - - if(singlet_manifold) then - - ispin = 1 - - 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)) - - endif - -! Triplet manifold - - if(triplet_manifold) then - - ispin = 2 - - 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)) - - endif - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - 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 - - adiabatic_connection = .true. - - if(adiabatic_connection) then - - write(*,*) '-------------------------------------------------------' - write(*,*) 'Adiabatic connection version of TDHF correlation energy' - write(*,*) '-------------------------------------------------------' - write(*,*) - - call ACDFT(.false.,.false.,.false.,.false.,singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho) - - end if - -end subroutine TDHF diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index dd307eb..bcad887 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -279,7 +279,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani end if - call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & + call ACFDT(scaled_screening,.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 7fd4536..89baebf 100644 --- a/src/QuAcK/ppRPA.f90 +++ b/src/QuAcK/ppRPA.f90 @@ -21,7 +21,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER ! Local variables - logical :: BSE integer :: ispin integer :: nOO integer :: nVV @@ -46,10 +45,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER Ec_ppRPA(:) = 0d0 -! Switch off Bethe-Salpeter equation for TDHF - - BSE = .false. - ! Singlet manifold if(singlet_manifold) then @@ -66,9 +61,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) - call linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & + Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & + Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & Ec_ppRPA(ispin)) call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) @@ -95,9 +90,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) - call linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & + Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & + Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & Ec_ppRPA(ispin)) call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) @@ -109,10 +104,10 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy (singlet) =',Ec_ppRPA(1) - write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy (triplet) =',3d0*Ec_ppRPA(2) - write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) - write(*,'(2X,A40,F15.6)') 'pp-RPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) + write(*,'(2X,A40,F15.6)') 'Tr@ppRPA correlation energy (singlet) =',Ec_ppRPA(1) + write(*,'(2X,A40,F15.6)') 'Tr@ppRPA correlation energy (triplet) =',3d0*Ec_ppRPA(2) + write(*,'(2X,A40,F15.6)') 'Tr@ppRPA correlation energy =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) + write(*,'(2X,A40,F15.6)') 'Tr@ppRPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index eb4a3f2..e539c43 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -311,7 +311,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani end if - call ACDFT(scaled_screening,.true.,TDA,BSE,singlet_manifold,triplet_manifold, & + call ACFDT(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 diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 638dbf6..df6991f 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -2,7 +2,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & doMP2,doMP3,doMP2F12, & doCCD,doCCSD,doCCSDT, & do_ring_CCD,do_ladder_CCD, & - doCIS,doRPA,doTDHF, & + doCIS,doRPA,doRPAx, & doppRPA,doADC, & doGF2,doGF3, & doG0W0,doevGW,doqsGW, & @@ -19,7 +19,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & logical,intent(out) :: doMP2,doMP3,doMP2F12 logical,intent(out) :: doCCD,doCCSD,doCCSDT logical,intent(out) :: do_ring_CCD,do_ladder_CCD - logical,intent(out) :: doCIS,doRPA,doTDHF,doppRPA,doADC + logical,intent(out) :: doCIS,doRPA,doRPAx,doppRPA,doADC logical,intent(out) :: doGF2,doGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW logical,intent(out) :: doG0T0,doevGT,doqsGT @@ -52,7 +52,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & doCIS = .false. doRPA = .false. - doTDHF = .false. + doRPAx = .false. doppRPA = .false. doADC = .false. @@ -101,7 +101,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & read(1,*) answer1,answer2,answer3,answer4,answer5 if(answer1 == 'T') doCIS = .true. if(answer2 == 'T') doRPA = .true. - if(answer3 == 'T') doTDHF = .true. + if(answer3 == 'T') doRPAx = .true. if(answer4 == 'T') doppRPA = .true. if(answer5 == 'T') doADC = .true.