From 3d3806da641e6a9ae9f092c5b471692f6dd5ef98 Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Mon, 13 Dec 2021 14:28:05 +0100
Subject: [PATCH] ppURPA

---
 input/methods      |   4 +-
 src/RPA/ppURPA.f90 | 133 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 135 insertions(+), 2 deletions(-)
 create mode 100644 src/RPA/ppURPA.f90

diff --git a/input/methods b/input/methods
index 77c21ce..0444112 100644
--- a/input/methods
+++ b/input/methods
@@ -13,9 +13,9 @@
 # G0F2* evGF2* qsGF2* G0F3 evGF3
   F     F     F      F    F
 # G0W0* evGW* qsGW* ufG0W0 ufGW
-  F    F    T       F      F
+  F    F    F       F      F
 # G0T0 evGT qsGT 
-  F    F    F
+  T    F    F
 # MCMP2
   F
 # * unrestricted version available
diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90
new file mode 100644
index 0000000..ac09ec1
--- /dev/null
+++ b/src/RPA/ppURPA.f90
@@ -0,0 +1,133 @@
+subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb,e)
+
+! Perform unrestricted pp-RPA calculations
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: TDA
+  logical,intent(in)            :: doACFDT
+  logical,intent(in)            :: spin_conserved
+  logical,intent(in)            :: spin_flip
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nR(nspin)
+  double precision,intent(in)   :: ENuc
+  double precision,intent(in)   :: EUHF
+  double precision,intent(in)   :: e(nBas,nspin)
+  double precision,intent(in)   :: ERI_aaaa(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_aabb(nBas,nBas,nBas,nBas)
+  double precision,intent(in)   :: ERI_bbbb(nBas,nBas,nBas,nBas)
+
+! Local variables
+
+  integer                       :: ispin
+  integer                       :: nS
+  integer                       :: nOOs,nOOt
+  integer                       :: nVVs,nVVt
+  double precision,allocatable  :: Omega1s(:),Omega1t(:)
+  double precision,allocatable  :: X1s(:,:),X1t(:,:)
+  double precision,allocatable  :: Y1s(:,:),Y1t(:,:)
+  double precision,allocatable  :: Omega2s(:),Omega2t(:)
+  double precision,allocatable  :: X2s(:,:),X2t(:,:)
+  double precision,allocatable  :: Y2s(:,:),Y2t(:,:)
+
+  double precision              :: Ec_ppRPA(nspin)
+  double precision              :: EcAC(nspin)
+
+! Hello world
+
+  write(*,*)
+  write(*,*)'****************************************'
+  write(*,*)'|  particle-particle URPA calculation  |'
+  write(*,*)'****************************************'
+  write(*,*)
+
+! Initialization
+
+  Ec_ppRPA(:) = 0d0
+  EcAC(:)   = 0d0
+
+! Useful quantities
+
+! nS   = nO*nV
+
+! nOOs = nO*(nO+1)/2
+! nVVs = nV*(nV+1)/2
+
+! nOOt = nO*(nO-1)/2
+! nVVt = nV*(nV-1)/2
+
+ ! Memory allocation
+
+! allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & 
+!          Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs))
+
+! allocate(Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & 
+!          Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt))
+
+! Spin-conserved manifold
+
+  if(spin_conserved) then 
+
+    ispin = 1
+
+!   call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,e,ERI, & 
+!                           Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,Ec_ppRPA(ispin))
+
+!   call print_excitation('pp-RPA (N+2)',5,nVVs,Omega1s)
+!   call print_excitation('pp-RPA (N-2)',5,nOOs,Omega2s)
+
+  endif
+
+! Spin-flip manifold 
+
+  if(spin_flip) then 
+
+    ispin = 2
+
+!   call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,e,ERI, &
+!                           Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,Ec_ppRPA(ispin))
+
+!   call print_excitation('pp-RPA (N+2)',6,nVVt,Omega1t)
+!   call print_excitation('pp-RPA (N-2)',6,nOOt,Omega2t)
+
+  endif
+
+  write(*,*)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-conserved) =',Ec_ppRPA(1)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-flip)      =',3d0*Ec_ppRPA(2)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy                  =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
+  write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy                        =',ENuc + EUHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2)
+  write(*,*)'-------------------------------------------------------------------------------'
+  write(*,*)
+
+! Compute the correlation energy via the adiabatic connection 
+
+! if(doACFDT) then
+
+!   write(*,*) '---------------------------------------------------------'
+!   write(*,*) 'Adiabatic connection version of pp-RPA correlation energy'
+!   write(*,*) '---------------------------------------------------------'
+!   write(*,*)
+
+!   call ACFDT_pp(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC)
+
+!   write(*,*)
+!   write(*,*)'-------------------------------------------------------------------------------'
+!   write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcAC(1),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcAC(2),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy           =',EcAC(1) + EcAC(2),' au'
+!   write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy                 =',ENuc + ERHF + EcAC(1) + EcAC(2),' au'
+!   write(*,*)'-------------------------------------------------------------------------------'
+!   write(*,*)
+
+! end if
+
+
+end subroutine ppURPA