mirror of
https://github.com/pfloos/quack
synced 2024-11-07 14:43:58 +01:00
434 lines
11 KiB
Fortran
434 lines
11 KiB
Fortran
|
subroutine G0F3(renormalization,nBas,nC,nO,nV,nR,V,e0)
|
||
|
|
|
||
|
|
! Perform third-order Green function calculation in diagonal approximation
|
||
|
|
|
||
|
|
implicit none
|
||
|
|
include 'parameters.h'
|
||
|
|
|
||
|
|
! Input variables
|
||
|
|
|
||
|
|
integer,intent(in) :: renormalization
|
||
|
|
integer,intent(in) :: nBas,nC,nO,nV,nR
|
||
|
|
double precision,intent(in) :: e0(nBas),V(nBas,nBas,nBas,nBas)
|
||
|
|
|
||
|
|
! Local variables
|
||
|
|
|
||
|
|
double precision :: eps,eps1,eps2
|
||
|
|
double precision,allocatable :: Sig2(:),SigInf(:),Sig3(:),eGF3(:),eOld(:)
|
||
|
|
double precision,allocatable :: App(:,:),Bpp(:,:),Cpp(:,:),Dpp(:,:)
|
||
|
|
double precision,allocatable :: Z(:),X2h1p(:),X1h2p(:),Sig2h1p(:),Sig1h2p(:)
|
||
|
|
|
||
|
|
integer :: i,j,k,l,a,b,c,d,p
|
||
|
|
|
||
|
|
! Hello world
|
||
|
|
|
||
|
|
write(*,*)
|
||
|
|
write(*,*)'************************************************'
|
||
|
|
write(*,*)'| Third-order Green function calculation |'
|
||
|
|
write(*,*)'************************************************'
|
||
|
|
write(*,*)
|
||
|
|
|
||
|
|
! Memory allocation
|
||
|
|
|
||
|
|
allocate(eGF3(nBas),Sig2(nBas),SigInf(nBas),Sig3(nBas), &
|
||
|
|
App(nBas,6),Bpp(nBas,2),Cpp(nBas,6),Dpp(nBas,6), &
|
||
|
|
Z(nBas),X2h1p(nBas),X1h2p(nBas),Sig2h1p(nBas),Sig1h2p(nBas))
|
||
|
|
|
||
|
|
!------------------------------------------------------------------------
|
||
|
|
! Compute third-order frequency-independent contribution
|
||
|
|
!------------------------------------------------------------------------
|
||
|
|
|
||
|
|
App(:,:) = 0d0
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do k=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
|
||
|
|
eps2 = e0(k) + e0(i) - e0(a) - e0(b)
|
||
|
|
|
||
|
|
App(p,1) = App(p,1) &
|
||
|
|
- (2d0*V(p,k,p,j) - V(p,k,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,k,i)/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
do c=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
|
||
|
|
eps2 = e0(j) + e0(i) - e0(a) - e0(c)
|
||
|
|
|
||
|
|
App(p,2) = App(p,2) &
|
||
|
|
+ (2d0*V(p,c,p,b) - V(p,c,b,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(j,i,c,a)/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
do c=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
|
||
|
|
eps2 = e0(j) - e0(c)
|
||
|
|
|
||
|
|
App(p,3) = App(p,3) &
|
||
|
|
+ (2d0*V(p,c,p,j) - V(p,c,j,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(a,b,c,i)/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
App(:,4) = App(:,3)
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do k=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = e0(j) + e0(i) - e0(a) - e0(b)
|
||
|
|
eps2 = e0(k) - e0(b)
|
||
|
|
|
||
|
|
App(p,5) = App(p,5) &
|
||
|
|
- (2d0*V(p,b,p,k) - V(p,b,k,p))*(2d0*V(j,i,a,b) - V(j,i,b,a))*V(i,j,k,a)/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
App(:,6) = App(:,5)
|
||
|
|
|
||
|
|
! Frequency-independent part of the third-order self-energy
|
||
|
|
|
||
|
|
SigInf(:) = App(:,1) + App(:,2) + App(:,3) + App(:,4) + App(:,5) + App(:,6)
|
||
|
|
|
||
|
|
! Frequency-dependent second-order contribution
|
||
|
|
|
||
|
|
Bpp(:,:) = 0d0
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps = eGF3(p) + e0(a) - e0(i) - e0(j)
|
||
|
|
|
||
|
|
Bpp(p,1) = Bpp(p,1) &
|
||
|
|
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps = eGF3(p) + e0(i) - e0(a) - e0(b)
|
||
|
|
|
||
|
|
Bpp(p,2) = Bpp(p,2) &
|
||
|
|
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
! Total second-order Green function
|
||
|
|
|
||
|
|
Sig2(:) = Bpp(:,1) + Bpp(:,2)
|
||
|
|
|
||
|
|
! Frequency-dependent third-order contribution: "C" terms
|
||
|
|
|
||
|
|
Cpp(:,:) = 0d0
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
do c=nO+1,nBas-nR
|
||
|
|
do d=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
|
||
|
|
eps2 = eGF3(p) + e0(i) - e0(c) - e0(d)
|
||
|
|
|
||
|
|
Cpp(p,1) = Cpp(p,1) &
|
||
|
|
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,c,d)*V(p,i,c,d)/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do k=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
|
||
|
|
eps2 = e0(j) + e0(k) - e0(a) - e0(b)
|
||
|
|
|
||
|
|
Cpp(p,2) = Cpp(p,2) &
|
||
|
|
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(a,b,j,k)*V(p,i,j,k)/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
Cpp(:,3) = Cpp(:,2)
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
do c=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
|
||
|
|
eps2 = e0(i) + e0(j) - e0(b) - e0(c)
|
||
|
|
|
||
|
|
Cpp(p,4) = Cpp(p,4) &
|
||
|
|
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(i,j,b,c)*V(p,a,b,c)/(eps1*eps2)
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
Cpp(:,5) = Cpp(:,4)
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do k=nC+1,nO
|
||
|
|
do l=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = eGF3(p) + e0(a) - e0(i) - e0(j)
|
||
|
|
eps2 = eGF3(p) + e0(a) - e0(k) - e0(l)
|
||
|
|
|
||
|
|
Cpp(p,6) = Cpp(p,6) &
|
||
|
|
- (2d0*V(p,a,k,l) - V(p,a,l,k))*V(k,l,i,j)*V(p,a,i,j)/(eps1*eps2)
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
! Frequency-dependent third-order contribution: "D" terms
|
||
|
|
|
||
|
|
Dpp(:,:) = 0d0
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
do c=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = eGF3(p) + e0(i) - e0(a) - e0(b)
|
||
|
|
eps2 = eGF3(p) + e0(j) - e0(b) - e0(c)
|
||
|
|
|
||
|
|
Dpp(p,1) = Dpp(p,1) &
|
||
|
|
+ V(p,i,a,b)*(V(a,j,i,c)*( V(p,j,c,b) - 2d0*V(p,j,b,c)) &
|
||
|
|
+ V(a,j,c,i)*( V(p,j,b,c) - 2d0*V(p,j,c,b)))/(eps1*eps2)
|
||
|
|
|
||
|
|
Dpp(p,1) = Dpp(p,1) &
|
||
|
|
+ V(p,i,b,a)*(V(a,j,i,c)*(4d0*V(p,j,b,c) - 2d0*V(p,j,c,b)) &
|
||
|
|
+ V(a,j,c,i)*( V(p,j,c,b) - 2d0*V(p,j,b,c)))/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
do c=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = eGF3(p) + e0(i) - e0(a) - e0(c)
|
||
|
|
eps2 = e0(i) + e0(j) - e0(a) - e0(b)
|
||
|
|
|
||
|
|
Dpp(p,2) = Dpp(p,2) &
|
||
|
|
+ V(p,i,c,a)*(V(a,b,i,j)*(4d0*V(p,b,c,j) - 2d0*V(p,b,j,c)) &
|
||
|
|
+ V(a,b,j,i)*( V(p,b,j,c) - 2d0*V(p,b,c,j)))/(eps1*eps2)
|
||
|
|
|
||
|
|
Dpp(p,2) = Dpp(p,2) &
|
||
|
|
+ V(p,i,a,c)*(V(a,b,i,j)*( V(p,b,j,c) - 2d0*V(p,b,c,j)) &
|
||
|
|
+ V(a,b,j,i)*( V(p,b,c,j) - 2d0*V(p,b,j,c)))/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
Dpp(:,3) = Dpp(:,2)
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do k=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = eGF3(p) + e0(a) - e0(j) - e0(k)
|
||
|
|
eps2 = e0(i) + e0(j) - e0(a) - e0(b)
|
||
|
|
|
||
|
|
Dpp(p,4) = Dpp(p,4) &
|
||
|
|
+ V(p,a,k,j)*(V(j,i,a,b)*(4d0*V(p,i,k,b) - 2d0*V(p,i,b,k)) &
|
||
|
|
+ V(j,i,b,a)*( V(p,i,b,k) - 2d0*V(p,i,k,b)))/(eps1*eps2)
|
||
|
|
|
||
|
|
Dpp(p,4) = Dpp(p,4) &
|
||
|
|
+ V(p,a,j,k)*(V(j,i,a,b)*( V(p,i,b,k) - 2d0*V(p,i,k,b)) &
|
||
|
|
+ V(j,i,b,a)*( V(p,i,k,b) - 2d0*V(p,i,b,k)))/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
Dpp(:,5) = Dpp(:,4)
|
||
|
|
|
||
|
|
do p=nC+1,nBas-nR
|
||
|
|
do i=nC+1,nO
|
||
|
|
do j=nC+1,nO
|
||
|
|
do k=nC+1,nO
|
||
|
|
do a=nO+1,nBas-nR
|
||
|
|
do b=nO+1,nBas-nR
|
||
|
|
|
||
|
|
eps1 = eGF3(p) + e0(a) - e0(i) - e0(k)
|
||
|
|
eps2 = eGF3(p) + e0(b) - e0(j) - e0(k)
|
||
|
|
|
||
|
|
Dpp(p,6) = Dpp(p,6) &
|
||
|
|
- V(p,a,k,i)*(V(i,b,a,j)*(4d0*V(p,b,k,j) - 2d0*V(p,b,j,k)) &
|
||
|
|
+ V(i,b,j,a)*( V(p,b,j,k) - 2d0*V(p,b,k,j)))/(eps1*eps2)
|
||
|
|
|
||
|
|
Dpp(p,6) = Dpp(p,6) &
|
||
|
|
- V(p,a,i,k)*(V(i,b,a,j)*( V(p,b,j,k) - 2d0*V(p,b,k,j)) &
|
||
|
|
+ V(i,b,j,a)*( V(p,b,k,j) - 2d0*V(p,b,j,k)))/(eps1*eps2)
|
||
|
|
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
|
||
|
|
! Compute renormalization factor (if required)
|
||
|
|
|
||
|
|
Z(:) = 1d0
|
||
|
|
|
||
|
|
if(renormalization == 0) then
|
||
|
|
|
||
|
|
Sig3(:) = SigInf(:) &
|
||
|
|
+ Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) + Cpp(:,6) &
|
||
|
|
+ Dpp(:,1) + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
|
||
|
|
|
||
|
|
elseif(renormalization == 1) then
|
||
|
|
|
||
|
|
Sig3(:) = SigInf(:) &
|
||
|
|
+ Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) + Cpp(:,6) &
|
||
|
|
+ Dpp(:,1) + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
|
||
|
|
|
||
|
|
Z(:) = Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) &
|
||
|
|
+ Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5)
|
||
|
|
|
||
|
|
Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/Sig2(nC+1:nBas-nR)
|
||
|
|
Z(:) = 1d0/(1d0 - Z(:))
|
||
|
|
|
||
|
|
Sig3(:) = Z(:)*Sig3(:)
|
||
|
|
|
||
|
|
elseif(renormalization == 2) then
|
||
|
|
|
||
|
|
Sig2h1p(:) = Cpp(:,4) + Cpp(:,5) + Cpp(:,6) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
|
||
|
|
Sig1h2p(:) = Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Dpp(:,1) + Dpp(:,2) + Dpp(:,3)
|
||
|
|
|
||
|
|
X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
|
||
|
|
X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
|
||
|
|
|
||
|
|
X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
|
||
|
|
X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
|
||
|
|
|
||
|
|
Sig3(:) = SigInf(:) + &
|
||
|
|
+ 1d0/(1d0 - X2h1p(:))*Sig2h1p(:) &
|
||
|
|
+ 1d0/(1d0 - X1h2p(:))*Sig1h2p(:)
|
||
|
|
|
||
|
|
elseif(renormalization == 3) then
|
||
|
|
|
||
|
|
Sig3(:) = SigInf(:) &
|
||
|
|
+ Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Cpp(:,4) + Cpp(:,5) + Cpp(:,6) &
|
||
|
|
+ Dpp(:,1) + Dpp(:,2) + Dpp(:,3) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
|
||
|
|
|
||
|
|
Sig2h1p(:) = Cpp(:,4) + Cpp(:,5) + Cpp(:,6) + Dpp(:,4) + Dpp(:,5) + Dpp(:,6)
|
||
|
|
Sig1h2p(:) = Cpp(:,1) + Cpp(:,2) + Cpp(:,3) + Dpp(:,1) + Dpp(:,2) + Dpp(:,3)
|
||
|
|
|
||
|
|
X2h1p(:) = Cpp(:,4) + Cpp(:,5) + Dpp(:,4) + Dpp(:,5)
|
||
|
|
X1h2p(:) = Cpp(:,2) + Cpp(:,3) + Dpp(:,2) + Dpp(:,3)
|
||
|
|
|
||
|
|
X2h1p(nC+1:nBas-nR) = X2h1p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,1)
|
||
|
|
X1h2p(nC+1:nBas-nR) = X1h2p(nC+1:nBas-nR)/Bpp(nC+1:nBas-nR,2)
|
||
|
|
|
||
|
|
Z(:) = X2h1p(:)*Sig2h1p(:) + X1h2p(:)*Sig1h2p(:)
|
||
|
|
Z(nC+1:nBas-nR) = Z(nC+1:nBas-nR)/(Sig3(nC+1:nBas-nR) - SigInf(nC+1:nBas-nR))
|
||
|
|
Z(:) = 1d0/(1d0 - Z(:))
|
||
|
|
|
||
|
|
Sig3(:) = Z(:)*Sig3(:)
|
||
|
|
|
||
|
|
endif
|
||
|
|
|
||
|
|
! Total third-order Green function
|
||
|
|
|
||
|
|
eGF3(:) = e0(:) + Sig2(:) + Sig3(:)
|
||
|
|
|
||
|
|
! Print results
|
||
|
|
|
||
|
|
call print_G0F3(nBas,nO,e0,Z,eGF3)
|
||
|
|
|
||
|
|
end subroutine G0F3
|