From 923d1a8bd6d7ba44eb3a484fea099e11bfaa0787 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sat, 5 Oct 2019 23:09:20 +0200 Subject: [PATCH] pp-RPA --- examples/molecule.LiH | 2 +- src/QuAcK/QuAcK.f90 | 21 ++++++- src/QuAcK/TDHF.f90 | 6 -- src/QuAcK/linear_response_pp.f90 | 17 +++--- src/QuAcK/ppRPA.f90 | 96 ++++++++++++++++++++++++++++++++ 5 files changed, 125 insertions(+), 17 deletions(-) create mode 100644 src/QuAcK/ppRPA.f90 diff --git a/examples/molecule.LiH b/examples/molecule.LiH index 95a3003..6878c9e 100644 --- a/examples/molecule.LiH +++ b/examples/molecule.LiH @@ -2,4 +2,4 @@ 2 2 2 0 0 # Znuc x y z Li 0. 0. 0. - H 0. 0. 3.9 + H 0. 0. 3.015 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 8bbe609..d81486f 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -7,7 +7,7 @@ program QuAcK logical :: doRHF,doUHF,doMOM logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doCCSD,doCCSDT - logical :: doCIS,doTDHF,doADC + logical :: doCIS,doTDHF,doppRPA,doADC logical :: doGF2,doGF3 logical :: doG0W0,doevGW,doqsGW logical :: doMCMP2,doMinMCMP2 @@ -43,6 +43,7 @@ program QuAcK double precision :: start_CCSD ,end_CCSD ,t_CCSD double precision :: start_CIS ,end_CIS ,t_CIS double precision :: start_TDHF ,end_TDHF ,t_TDHF + double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA double precision :: start_ADC ,end_ADC ,t_ADC double precision :: start_GF2 ,end_GF2 ,t_GF2 double precision :: start_GF3 ,end_GF3 ,t_GF3 @@ -104,7 +105,7 @@ program QuAcK call read_methods(doRHF,doUHF,doMOM, & doMP2,doMP3,doMP2F12, & doCCD,doCCSD,doCCSDT, & - doCIS,doTDHF,doADC, & + doCIS,doTDHF,doppRPA,doADC, & doGF2,doGF3, & doG0W0,doevGW,doqsGW, & doMCMP2) @@ -400,6 +401,22 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Compute pp-RPA excitations +!------------------------------------------------------------------------ + + if(doppRPA) then + + call cpu_time(start_ppRPA) + call ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF) + call cpu_time(end_ppRPA) + + t_ppRPA = end_ppRPA - start_ppRPA + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pp-RPA = ',t_ppRPA,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Compute ADC excitations !------------------------------------------------------------------------ diff --git a/src/QuAcK/TDHF.f90 b/src/QuAcK/TDHF.f90 index adf4d3a..e40e9b8 100644 --- a/src/QuAcK/TDHF.f90 +++ b/src/QuAcK/TDHF.f90 @@ -70,12 +70,6 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('TDHF ',ispin,nS,Omega(:,ispin)) - call linear_response_ph(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, & - EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) - call print_excitation('p-h TDHF ',ispin,nS,Omega(:,ispin)) - - call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho) - endif ! Triplet manifold diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index 437a3a3..3b9874a 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho) +subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,Ec_ppRPA) ! Compute the p-p channel of the linear response @@ -25,10 +25,11 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho) double precision,allocatable :: D(:,:) double precision,allocatable :: M(:,:) double precision,allocatable :: w(:) - double precision :: Ec_ppRPA ! Output variables + double precision,intent(out) :: Ec_ppRPA + ! Useful quantities nOO = nO*nO @@ -60,11 +61,11 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho) ! Off-diagonal blocks - M(1:nOO,nOO+1:nOO+nVV) = +transpose(B(1:nVV,1:nOO)) - M(nOO+1:nOO+nVV,1:nOO) = - B(1:nVV,1:nOO) + M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO) + M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO)) - print*, 'pp-RPA matrix' - call matout(nOO+nVV,nOO+nVV,M(:,:)) +! print*, 'pp-RPA matrix' +! call matout(nOO+nVV,nOO+nVV,M(:,:)) ! Diagonalize the p-h matrix @@ -78,11 +79,11 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho) ! call matout(nS,nS,XpY) print*,'pp-RPA excitation energies' - call matout(2*nS,1,w) + call matout(nOO+nVV,1,w) ! Compute the RPA correlation energy - Ec_ppRPA = 0.5d0*(sum(abs(w)) - trace_matrix(nVV,C) - trace_matrix(nOO,D)) + Ec_ppRPA = 0.5d0*(sum(abs(w(:))) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:))) print*,'Ec(pp-RPA) = ',Ec_ppRPA diff --git a/src/QuAcK/ppRPA.f90 b/src/QuAcK/ppRPA.f90 new file mode 100644 index 0000000..594a556 --- /dev/null +++ b/src/QuAcK/ppRPA.f90 @@ -0,0 +1,96 @@ +subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e) + +! Perform pp-RPA calculation + + implicit none + include 'parameters.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 :: Ec_ppRPA(nspin) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Time-dependent Hartree-Fock calculation |' + write(*,*)'************************************************' + write(*,*) + +! Initialization + + Ec_ppRPA(:) = 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)) + +! Singlet manifold + + if(singlet_manifold) then + + ispin = 1 + + call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, & + Ec_ppRPA) +! call print_excitation('pp-RPA ',ispin,nS,Omega(:,ispin)) + + endif + +! Triplet manifold + + if(triplet_manifold) then + + ispin = 2 + + call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, & + Ec_ppRPA) +! call print_excitation('pp-RPA ',ispin,nS,Omega(:,ispin)) + + endif + + 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) =',Ec_ppRPA(2) + write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy =',Ec_ppRPA(1) + Ec_ppRPA(2) + write(*,'(2X,A40,F15.6)') 'pp-RPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + Ec_ppRPA(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +end subroutine ppRPA