phGLR subroutines

src/LR/phGLR.f90

@@ -0,0 +1,82 @@
+subroutine phGLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
+
+! Compute linear response
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: TDA
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: Aph(nS,nS)
+  double precision,intent(in)   :: Bph(nS,nS)
+
+  ! Local variables
+
+  double precision              :: trace_matrix
+  double precision,allocatable  :: ApB(:,:)
+  double precision,allocatable  :: AmB(:,:)
+  double precision,allocatable  :: AmBSq(:,:)
+  double precision,allocatable  :: AmBIv(:,:)
+  double precision,allocatable  :: Z(:,:)
+  double precision,allocatable  :: tmp(:,:)
+
+! Output variables
+
+  double precision,intent(out)  :: EcRPA
+  double precision,intent(out)  :: Om(nS)
+  double precision,intent(out)  :: XpY(nS,nS)
+  double precision,intent(out)  :: XmY(nS,nS)
+
+! Memory allocation
+
+  allocate(ApB(nS,nS),AmB(nS,nS),AmBSq(nS,nS),AmBIv(nS,nS),Z(nS,nS),tmp(nS,nS))
+
+! Tamm-Dancoff approximation
+
+  if(TDA) then
+ 
+    XpY(:,:) = Aph(:,:)
+    call diagonalize_matrix(nS,XpY,Om)
+    XpY(:,:) = transpose(XpY(:,:))
+    XmY(:,:) = XpY(:,:)
+
+  else
+
+    ApB(:,:) = Aph(:,:) + Bph(:,:)
+    AmB(:,:) = Aph(:,:) - Bph(:,:)
+
+  ! Diagonalize linear response matrix
+
+    call diagonalize_matrix(nS,AmB,Om)
+
+    if(minval(Om) < 0d0) &
+      call print_warning('You may have instabilities in linear response: A-B is not positive definite!!')
+
+    call ADAt(nS,AmB,1d0*dsqrt(Om),AmBSq)
+    call ADAt(nS,AmB,1d0/dsqrt(Om),AmBIv)
+
+    call dgemm('N','N',nS,nS,nS,1d0,ApB,size(ApB,1),AmBSq,size(AmBSq,1),0d0,tmp,size(tmp,1))
+    call dgemm('N','N',nS,nS,nS,1d0,AmBSq,size(AmBSq,1),tmp,size(tmp,1),0d0,Z,size(Z,1))
+
+    call diagonalize_matrix(nS,Z,Om)
+
+    if(minval(Om) < 0d0) & 
+      call print_warning('You may have instabilities in linear response: negative excitations!!')
+ 
+    Om = sqrt(Om)
+
+    call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBSq,size(AmBSq,1),0d0,XpY,size(XpY,1))
+    call DA(nS,1d0/dsqrt(Om),XpY)
+
+    call dgemm('T','N',nS,nS,nS,1d0,Z,size(Z,1),AmBIv,size(AmBIv,1),0d0,XmY,size(XmY,1))
+    call DA(nS,1d0*dsqrt(Om),XmY)
+ 
+  end if
+
+  ! Compute the RPA correlation energy
+
+  EcRPA = 0.5d0*(sum(Om) - trace_matrix(nS,Aph))
+
+end subroutine 

src/LR/phGLR_A.f90

@@ -0,0 +1,56 @@
+subroutine phRLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,e,ERI,Aph)
+
+! Compute resonant block of the ph channel
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: dRPA
+  integer,intent(in)            :: nOrb
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: e(nOrb)
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb) 
+  
+! Local variables
+
+  double precision              :: delta_dRPA
+  double precision,external     :: Kronecker_delta
+
+  integer                       :: i,j,a,b,ia,jb
+
+! Output variables
+
+  double precision,intent(out)  :: Aph(nS,nS)
+
+! Direct RPA
+
+  delta_dRPA = 0d0
+  if(dRPA) delta_dRPA = 1d0
+
+! Build A matrix for spin orbitals
+
+  ia = 0
+  do i=nC+1,nO
+    do a=nO+1,nOrb-nR
+      ia = ia + 1
+      jb = 0
+      do j=nC+1,nO
+        do b=nO+1,nOrb-nR
+          jb = jb + 1
+
+          Aph(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+                     + lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
+
+        end  do
+      end  do
+    end  do
+  end  do
+
+end subroutine 

src/LR/phGLR_B.f90

@@ -0,0 +1,53 @@
+subroutine phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph)
+
+! Compute the coupling block of the ph channel
+
+  implicit none
+  include 'parameters.h'
+
+! Input variables
+
+  logical,intent(in)            :: dRPA
+  integer,intent(in)            :: nOrb
+  integer,intent(in)            :: nC
+  integer,intent(in)            :: nO
+  integer,intent(in)            :: nV
+  integer,intent(in)            :: nR
+  integer,intent(in)            :: nS
+  double precision,intent(in)   :: lambda
+  double precision,intent(in)   :: ERI(nOrb,nOrb,nOrb,nOrb)
+  
+! Local variables
+
+  double precision              :: delta_dRPA
+
+  integer                       :: i,j,a,b,ia,jb
+
+! Output variables
+
+  double precision,intent(out)  :: Bph(nS,nS)
+
+! Direct RPA
+
+  delta_dRPA = 0d0
+  if(dRPA) delta_dRPA = 1d0
+
+! Build B matrix for spin orbitals 
+
+  ia = 0
+  do i=nC+1,nO
+    do a=nO+1,nOrb-nR
+      ia = ia + 1
+      jb = 0
+      do j=nC+1,nO
+        do b=nO+1,nOrb-nR
+          jb = jb + 1
+
+          Bph(ia,jb) = lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
+
+        end do
+      end do
+    end do
+  end do
+
+end subroutine