From 27f14363c71284c1b1c2718070a48e66c991d6e2 Mon Sep 17 00:00:00 2001 From: pfloos Date: Fri, 6 Sep 2024 21:35:42 +0200 Subject: [PATCH] phGLR subroutines --- src/LR/phGLR.f90 | 82 ++++++++++++++++++++++++++++++++++++++++++++++ src/LR/phGLR_A.f90 | 56 +++++++++++++++++++++++++++++++ src/LR/phGLR_B.f90 | 53 ++++++++++++++++++++++++++++++ 3 files changed, 191 insertions(+) create mode 100644 src/LR/phGLR.f90 create mode 100644 src/LR/phGLR_A.f90 create mode 100644 src/LR/phGLR_B.f90 diff --git a/src/LR/phGLR.f90 b/src/LR/phGLR.f90 new file mode 100644 index 0000000..1492eb3 --- /dev/null +++ b/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 diff --git a/src/LR/phGLR_A.f90 b/src/LR/phGLR_A.f90 new file mode 100644 index 0000000..29ae8a8 --- /dev/null +++ b/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 diff --git a/src/LR/phGLR_B.f90 b/src/LR/phGLR_B.f90 new file mode 100644 index 0000000..acf94ec --- /dev/null +++ b/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