pp-RPA

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

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
 !------------------------------------------------------------------------

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 

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(    1:nVV    ,nVV+1:nOO+nVV) = -           B(1:nVV,1:nOO)
-  M(nOO+1:nOO+nVV,1:nOO) = -          B(1:nVV,1:nOO)
+  M(nVV+1:nOO+nVV,    1:nVV)     = + transpose(B(1:nVV,1:nOO))
 
-  print*, 'pp-RPA matrix'
+! print*, 'pp-RPA matrix'
-  call matout(nOO+nVV,nOO+nVV,M(:,:))
+! 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
 

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