LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-18 11:15:30 +02:00
Code Releases Activity

GF clean up

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-19 10:21:18 +01:00
parent d0020765fc
commit bcbeba2a7c
29 changed files with 904 additions and 356 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.H2
View File

@ -2,4 +2,4 @@
2 1 1 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 1.399
H 0. 0. 1.402

2
examples/molecule.LiH
View File

@ -2,4 +2,4 @@
2 2 2 0 0
# Znuc x y z
Li 0. 0. 0.
H 0. 0. 3.5
H 0. 0. 3.017

55
input/basis
View File

@ -1,30 +1,39 @@
1 6
1 10
S 8
1 2940.0000000 0.0006800
2 441.2000000 0.0052360
3 100.5000000 0.0266060
4 28.4300000 0.0999930
5 9.1690000 0.2697020
6 3.1960000 0.4514690
7 1.1590000 0.2950740
8 0.1811000 0.0125870
1 24350.0000000 0.0005020
2 3650.0000000 0.0038810
3 829.6000000 0.0199970
4 234.0000000 0.0784180
5 75.6100000 0.2296760
6 26.7300000 0.4327220
7 9.9270000 0.3506420
8 1.1020000 -0.0076450
S 8
1 2940.0000000 -0.0001230
2 441.2000000 -0.0009660
3 100.5000000 -0.0048310
4 28.4300000 -0.0193140
5 9.1690000 -0.0532800
6 3.1960000 -0.1207230
7 1.1590000 -0.1334350
8 0.1811000 0.5307670
1 24350.0000000 -0.0001180
2 3650.0000000 -0.0009150
3 829.6000000 -0.0047370
4 234.0000000 -0.0192330
5 75.6100000 -0.0603690
6 26.7300000 -0.1425080
7 9.9270000 -0.1777100
8 1.1020000 0.6058360
S 1
1 0.0589000 1.0000000
1 2.8360000 1.0000000
S 1
1 0.3782000 1.0000000
P 3
1 3.6190000 0.0291110
2 0.7110000 0.1693650
3 0.1951000 0.5134580
1 54.7000000 0.0171510
2 12.4300000 0.1076560
3 3.6790000 0.3216810
P 1
1 0.0601800 1.0000000
1 1.1430000 1.0000000
P 1
1 0.3300000 1.0000000
D 1
1 0.2380000 1.0000000
1 4.0140000 1.0000000
D 1
1 1.0960000 1.0000000
F 1
1 2.5440000 1.0000000

4
input/methods
View File

@ -7,9 +7,9 @@
# CIS RPA RPAx ppRPA ADC
F F F F F
# GF2 GF3
F F
T F
# G0W0 evGW qsGW
T F F
F F F
# G0T0 evGT qsGT
F F F
# MCMP2

4
input/molecule
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd
1 2 2 0 0
1 5 5 0 0
# Znuc x y z
Be 0.0 0.0 0.0
Ne 0.0 0.0 0.0

2
input/molecule.xyz
View File

@ -1,3 +1,3 @@
1
Be 0.0000000000 0.0000000000 0.0000000000
Ne 0.0000000000 0.0000000000 0.0000000000

2
input/options
View File

@ -11,6 +11,6 @@
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 5 F F T F F F T 0.000
# ACFDT: AC Kx XBS
T F T
T F F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

55
input/weight
View File

@ -1,30 +1,39 @@
1 6
1 10
S 8
1 2940.0000000 0.0006800
2 441.2000000 0.0052360
3 100.5000000 0.0266060
4 28.4300000 0.0999930
5 9.1690000 0.2697020
6 3.1960000 0.4514690
7 1.1590000 0.2950740
8 0.1811000 0.0125870
1 24350.0000000 0.0005020
2 3650.0000000 0.0038810
3 829.6000000 0.0199970
4 234.0000000 0.0784180
5 75.6100000 0.2296760
6 26.7300000 0.4327220
7 9.9270000 0.3506420
8 1.1020000 -0.0076450
S 8
1 2940.0000000 -0.0001230
2 441.2000000 -0.0009660
3 100.5000000 -0.0048310
4 28.4300000 -0.0193140
5 9.1690000 -0.0532800
6 3.1960000 -0.1207230
7 1.1590000 -0.1334350
8 0.1811000 0.5307670
1 24350.0000000 -0.0001180
2 3650.0000000 -0.0009150
3 829.6000000 -0.0047370
4 234.0000000 -0.0192330
5 75.6100000 -0.0603690
6 26.7300000 -0.1425080
7 9.9270000 -0.1777100
8 1.1020000 0.6058360
S 1
1 0.0589000 1.0000000
1 2.8360000 1.0000000
S 1
1 0.3782000 1.0000000
P 3
1 3.6190000 0.0291110
2 0.7110000 0.1693650
3 0.1951000 0.5134580
1 54.7000000 0.0171510
2 12.4300000 0.1076560
3 3.6790000 0.3216810
P 1
1 0.0601800 1.0000000
1 1.1430000 1.0000000
P 1
1 0.3300000 1.0000000
D 1
1 0.2380000 1.0000000
1 4.0140000 1.0000000
D 1
1 1.0960000 1.0000000
F 1
1 2.5440000 1.0000000

122
src/QuAcK/G0F2.f90 Normal file
View File

@ -0,0 +1,122 @@
subroutine G0F2(linearize,nBas,nC,nO,nV,nR,V,e0)
! Perform a one-shot second-order Green function calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: linearize
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: nC
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: e0(nBas)
double precision,intent(in) :: V(nBas,nBas,nBas,nBas)
! Local variables
double precision :: eps
double precision,allocatable :: eGF2(:)
double precision,allocatable :: Bpp(:,:)
double precision,allocatable :: Z(:)
integer :: i,j,a,b,p
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| One-shot second-order Green function |'
write(*,*)'************************************************'
write(*,*)
! Memory allocation
allocate(Bpp(nBas,2),Z(nBas),eGF2(nBas))
! 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 = eGF2(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
end do
end do
end do
end do
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(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
end do
end do
end do
end do
! Compute the renormalization factor
Z(:) = 0d0
do p=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + e0(a) - e0(i) - e0(j)
Z(p) = Z(p) - (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps**2
end do
end do
end do
end do
do p=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + e0(i) - e0(a) - e0(b)
Z(p) = Z(p) - (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps**2
end do
end do
end do
end do
Z(:) = 1d0/(1d0 - Z(:))
if(linearize) then
eGF2(:) = e0(:) + Z(:)*(Bpp(:,1) + Bpp(:,2))
else
eGF2(:) = e0(:) + Bpp(:,1) + Bpp(:,2)
end if
! Print results
call print_G0F2(nBas,nO,e0,eGF2)
end subroutine G0F2

433
src/QuAcK/G0F3.f90 Normal file
View File

@ -0,0 +1,433 @@
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

4
src/QuAcK/G0T0.f90
View File

@ -77,6 +77,8 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
Omega2s(:),X2s(:,:),Y2s(:,:), &
EcRPA(ispin))
EcRPA(ispin) = 1d0*EcRPA(ispin)
call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s(:))
call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s(:))
@ -100,6 +102,8 @@ subroutine G0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
Omega2t(:),X2t(:,:),Y2t(:,:), &
EcRPA(ispin))
EcRPA(ispin) = 3d0*EcRPA(ispin)
call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t(:))
call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t(:))

137
src/QuAcK/GF2.f90
View File

@ -1,137 +0,0 @@
subroutine GF2(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,V,e0)
! Perform second-order Green function calculation in diagonal approximation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: maxSCF
double precision,intent(in) :: thresh
integer,intent(in) :: max_diis
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e0(nBas),V(nBas,nBas,nBas,nBas)
! Local variables
integer :: nSCF,n_diis
double precision :: eps
double precision :: Conv
double precision :: rcond
double precision,allocatable :: eGF2(:),eOld(:),Bpp(:,:,:),error_diis(:,:),e_diis(:,:)
integer :: i,j,a,b,p,q
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Second-order Green function calculation |'
write(*,*)'************************************************'
write(*,*)
! Memory allocation
allocate(Bpp(nBas,nBas,2),eGF2(nBas),eOld(nBas), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization
Conv = 1d0
nSCF = 0
n_diis = 0
e_diis(:,:) = 0d0
error_diis(:,:) = 0d0
eGF2(:) = e0(:)
eOld(:) = e0(:)
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
do while(Conv > thresh .and. nSCF < maxSCF)
! Frequency-dependent second-order contribution
Bpp(:,:,:) = 0d0
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do j=nC+1,nO
do a=nO+1,nBas-nR
eps = eGF2(p) + e0(a) - e0(i) - e0(j)
Bpp(p,q,1) = Bpp(p,q,1) &
+ (2d0*V(p,a,i,j) - V(p,a,j,i))*V(q,a,i,j)/eps
enddo
enddo
enddo
enddo
enddo
do p=nC+1,nBas-nR
do q=nC+1,nBas-nR
do i=nC+1,nO
do a=nO+1,nBas-nR
do b=nO+1,nBas-nR
eps = eGF2(p) + e0(i) - e0(a) - e0(b)
Bpp(p,q,2) = Bpp(p,q,2) &
+ (2d0*V(p,i,a,b) - V(p,i,b,a))*V(q,i,a,b)/eps
enddo
enddo
enddo
enddo
enddo
print*,'Sig2 in GF2'
call matout(nBas,nBas,Bpp(:,:,1) + Bpp(:,:,2))
! eGF2(:) = e0(:) &
! + Bpp(:,1) + Bpp(:,2)
Conv = maxval(abs(eGF2 - eOld))
! DIIS extrapolation
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGF2-eOld,eGF2)
! Reset DIIS if required
if(abs(rcond) < 1d-15) n_diis = 0
eOld = eGF2
! Print results
call print_GF2(nBas,nO,nSCF,Conv,e0,eGF2)
! Increment
nSCF = nSCF + 1
enddo
!------------------------------------------------------------------------
! End main SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF+1) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
endif
end subroutine GF2

64
src/QuAcK/QuAcK.f90
View File

@ -10,7 +10,7 @@ program QuAcK
logical :: do_ring_CCD,do_ladder_CCD
logical :: doCIS,doRPA,doRPAx
logical :: doppRPA,doADC
logical :: doGF2,doGF3
logical :: doG0F2,doevGF2,doG0F3,doevGF3
logical :: doG0W0,doevGW,doqsGW
logical :: doG0T0,doevGT,doqsGT
logical :: doMCMP2,doMinMCMP2
@ -114,15 +114,15 @@ program QuAcK
! Which calculations do you want to do?
call read_methods(doRHF,doUHF,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, &
do_ring_CCD,do_ladder_CCD, &
doCIS,doRPA,doRPAx, &
doppRPA,doADC, &
doGF2,doGF3, &
doG0W0,doevGW,doqsGW, &
doG0T0,doevGT,doqsGT, &
call read_methods(doRHF,doUHF,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, &
do_ring_CCD,do_ladder_CCD, &
doCIS,doRPA,doRPAx, &
doppRPA,doADC, &
doG0F2,doevGF2,doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doG0T0,doevGT,doqsGT, &
doMCMP2)
! Read options for methods
@ -507,13 +507,13 @@ program QuAcK
end if
!------------------------------------------------------------------------
! Compute GF2 electronic binding energies
! Compute G0F2 electronic binding energies
!------------------------------------------------------------------------
if(doGF2) then
if(doG0F2) then
call cpu_time(start_GF2)
call GF2_diag(maxSCF_GF,thresh_GF,n_diis_GF,linearize,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call G0F2(maxSCF_GF,thresh_GF,n_diis_GF,linearize,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2
@ -523,13 +523,45 @@ program QuAcK
end if
!------------------------------------------------------------------------
! Compute GF3 electronic binding energies
! Compute evGF2 electronic binding energies
!------------------------------------------------------------------------
if(doGF3) then
if(doevGF2) then
call cpu_time(start_GF2)
call evGF2(maxSCF_GF,thresh_GF,n_diis_GF,linearize,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF2,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute G0F3 electronic binding energies
!------------------------------------------------------------------------
if(doG0F3) then
call cpu_time(start_GF3)
call GF3_diag(maxSCF_GF,thresh_GF,n_diis_GF,renormalization,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call G0F3(renormalization,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call cpu_time(end_GF3)
t_GF3 = end_GF3 - start_GF3
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF3,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------
! Compute evGF3 electronic binding energies
!------------------------------------------------------------------------
if(doevGF3) then
call cpu_time(start_GF3)
call evGF3(maxSCF_GF,thresh_GF,n_diis_GF,renormalization,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call cpu_time(end_GF3)
t_GF3 = end_GF3 - start_GF3

8
src/QuAcK/GF2_diag.f90 → src/QuAcK/evGF2.f90
View File

@ -1,6 +1,6 @@
subroutine GF2_diag(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
subroutine evGF2(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
! Perform second-order Green function calculation in diagonal approximation
! Perform eigenvalue self-consistent second-order Green function calculation
implicit none
include 'parameters.h'
@ -145,7 +145,7 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
! Print results
call print_GF2(nBas,nO,nSCF,Conv,e0,eGF2)
call print_evGF2(nBas,nO,nSCF,Conv,e0,eGF2)
! DIIS extrapolation
@ -177,4 +177,4 @@ subroutine GF2_diag(maxSCF,thresh,max_diis,linearize,nBas,nC,nO,nV,nR,V,e0)
end if
end subroutine GF2_diag
end subroutine evGF2

6
src/QuAcK/GF3_diag.f90 → src/QuAcK/evGF3.f90
View File

@ -1,4 +1,4 @@
subroutine GF3_diag(maxSCF,thresh,max_diis,renormalization,nBas,nC,nO,nV,nR,V,e0)
subroutine evGF3(maxSCF,thresh,max_diis,renormalization,nBas,nC,nO,nV,nR,V,e0)
! Perform third-order Green function calculation in diagonal approximation
@ -455,7 +455,7 @@
! Print results
call print_GF3(nBas,nO,nSCF,Conv,e0,Z,eGF3)
call print_evGF3(nBas,nO,nSCF,Conv,e0,Z,eGF3)
! DIIS extrapolation
@ -489,4 +489,4 @@
endif
end subroutine GF3_diag
end subroutine evGF3

128
src/QuAcK/excitation_density_Tmatrix.f90
View File

@ -43,7 +43,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c
do d=c,nBas-nR
cd = cd + 1
rho1(p,i,ab) = rho1(p,i,ab) &
! + ERI(p,i,c,d)*X1(cd,ab)
@ -54,7 +54,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
kl = 0
do k=nC+1,nO
do l=nC+1,k
do l=k,nO
kl = kl + 1
rho1(p,i,ab) = rho1(p,i,ab) &
! + ERI(p,i,k,l)*Y1(kl,ab)
@ -71,7 +71,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c
do d=c,nBas-nR
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) &
! + ERI(p,nO+a,c,d)*X2(cd,ij)
@ -83,7 +83,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
kl = 0
do k=nC+1,nO
do l=nC+1,k
do l=k,nO
kl = kl + 1
rho2(p,a,ij) = rho2(p,a,ij) &
! + ERI(p,nO+a,k,l)*Y2(kl,ij)
@ -105,65 +105,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
if(ispin == 2) then
do p=nC+1,nBas-nR
do i=nC+1,nO
do ab=1,nVV
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c-1
cd = cd + 1
rho1(p,i,ab) = rho1(p,i,ab) &
+ (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
end do
end do
kl = 0
do k=nC+1,nO
do l=nC+1,k-1
kl = kl + 1
rho1(p,i,ab) = rho1(p,i,ab) &
+ (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
end do
end do
end do
end do
do a=1,nV-nR
do ij=1,nOO
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c-1
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
end do
end do
kl = 0
do k=nC+1,nO
do l=nC+1,k-1
kl = kl + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
end do
end do
end do
end do
end do
end if
!----------------------------------------------
! Spinorbital basis
!----------------------------------------------
if(ispin == 3) then
do p=nC+1,nBas-nR
do i=nC+1,nO
@ -195,7 +136,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
@ -213,6 +154,65 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
end do
end do
end do
end if
!----------------------------------------------
! Spinorbital basis
!----------------------------------------------
if(ispin == 3) then
do p=nC+1,nBas-nR
do i=nC+1,nO
do ab=1,nVV
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho1(p,i,ab) = rho1(p,i,ab) &
+ (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
end do
end do
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho1(p,i,ab) = rho1(p,i,ab) &
+ (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
end do
end do
end do
end do
do a=1,nV-nR
do ij=1,nOO
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
end do
end do
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
end do
end do
end do
end do
end do

8
src/QuAcK/linear_response_B_pp.f90
View File

@ -26,11 +26,11 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
ab = 0
do a=nO+1,nBas-nR
do b=nO+1,a
do b=a,nBas-nR
ab = ab + 1
ij = 0
do i=nC+1,nO
do j=nC+1,i
do j=i,nO
ij = ij + 1
B_pp(ab,ij) = (ERI(a,b,i,j) + ERI(a,b,j,i))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j)))
@ -48,11 +48,11 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp)
ab = 0
do a=nO+1,nBas-nR
do b=nO+1,a-1
do b=a+1,nBas-nR
ab = ab + 1
ij = 0
do i=nC+1,nO
do j=nC+1,i-1
do j=i+1,nO
ij = ij + 1
B_pp(ab,ij) = ERI(a,b,i,j) - ERI(a,b,j,i)

8
src/QuAcK/linear_response_C_pp.f90
View File

@ -32,11 +32,11 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
ab = 0
do a=nO+1,nBas-nR
do b=nO+1,a
do b=a,nBas-nR
ab = ab + 1
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c
do d=c,nBas-nR
cd = cd + 1
C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &
@ -55,11 +55,11 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
ab = 0
do a=nO+1,nBas-nR
do b=nO+1,a-1
do b=a+1,nBas-nR
ab = ab + 1
cd = 0
do c=nO+1,nBas-nR
do d=nO+1,c-1
do d=c+1,nBas-nR
cd = cd + 1
C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) &

8
src/QuAcK/linear_response_D_pp.f90
View File

@ -32,11 +32,11 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
ij = 0
do i=nC+1,nO
do j=nC+1,i
do j=i,nO
ij = ij + 1
kl = 0
do k=nC+1,nO
do l=nC+1,k
do l=k,nO
kl = kl + 1
D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &
@ -55,11 +55,11 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
ij = 0
do i=nC+1,nO
do j=nC+1,i-1
do j=i+1,nO
ij = ij + 1
kl = 0
do k=nC+1,nO
do l=nC+1,k-1
do l=k+1,nO
kl = kl + 1
D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) &

45
src/QuAcK/print_G0F2.f90 Normal file
View File

@ -0,0 +1,45 @@
subroutine print_G0F2(nBas,nO,eHF,eGF2)
! Print one-electron energies and other stuff for G0F2
implicit none
include 'parameters.h'
integer,intent(in) :: nBas
integer,intent(in) :: nO
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: eGF2(nBas)
integer :: x
integer :: HOMO
integer :: LUMO
double precision :: Gap
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eGF2(LUMO)-eGF2(HOMO)
! Dump results
write(*,*)'-------------------------------------------'
write(*,*)' Frequency-dependent G0F2 calculation'
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','e_G0F2 (eV)','|'
write(*,*)'-------------------------------------------'
do x=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',eHF(x)*HaToeV,'|',eGF2(x)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------'
write(*,'(2X,A27,F15.6)') 'G0F2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0F2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0F2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------'
write(*,*)
end subroutine print_G0F2

41
src/QuAcK/print_G0F3.f90 Normal file
View File

@ -0,0 +1,41 @@
subroutine print_G0F3(nBas,nO,eHF,Z,eGF3)
! Print one-electron energies and other stuff for GF3
implicit none
include 'parameters.h'
integer,intent(in) :: nBas,nO
double precision,intent(in) :: eHF(nBas),eGF3(nBas),Z(nBas)
integer :: x,HOMO,LUMO
double precision :: Gap
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eGF3(LUMO)-eGF3(HOMO)
! Dump results
write(*,*)'-------------------------------------------------------------'
write(*,*)' Frequency-dependent G0F3 calculation'
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','Z','|','e_G0F3 (eV)','|'
write(*,*)'-------------------------------------------------------------'
do x=1,nBas
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',x,'|',eHF(x)*HaToeV,'|',Z(x),'|',eGF3(x)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'G0F3 HOMO energy (eV):',eGF3(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0F3 LUMO energy (eV):',eGF3(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'G0F3 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------'
write(*,*)
end subroutine print_G0F3

14
src/QuAcK/print_GF2.f90 → src/QuAcK/print_evGF2.f90
View File

@ -1,4 +1,4 @@
subroutine print_GF2(nBas,nO,nSCF,Conv,eHF,eGF2)
subroutine print_evGF2(nBas,nO,nSCF,Conv,eHF,eGF2)
! Print one-electron energies and other stuff for GF2
@ -20,10 +20,10 @@ subroutine print_GF2(nBas,nO,nSCF,Conv,eHF,eGF2)
! Dump results
write(*,*)'-------------------------------------------'
write(*,*)' Frequency-dependent diagonal GF2 calculation'
write(*,*)' Frequency-dependent evGF2 calculation'
write(*,*)'-------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','e_GF2 (eV)','|'
'|','#','|','e_HF (eV)','|','e_evGF2 (eV)','|'
write(*,*)'-------------------------------------------'
do x=1,nBas
@ -35,10 +35,10 @@ subroutine print_GF2(nBas,nO,nSCF,Conv,eHF,eGF2)
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'-------------------------------------------'
write(*,'(2X,A27,F15.6)') 'GF2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'GF2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'GF2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF2 HOMO energy (eV):',eGF2(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF2 LUMO energy (eV):',eGF2(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF2 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------'
write(*,*)
end subroutine print_GF2
end subroutine print_evGF2

14
src/QuAcK/print_GF3.f90 → src/QuAcK/print_evGF3.f90
View File

@ -1,4 +1,4 @@
subroutine print_GF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
subroutine print_evGF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
! Print one-electron energies and other stuff for GF3
@ -20,10 +20,10 @@ subroutine print_GF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
! Dump results
write(*,*)'-------------------------------------------------------------'
write(*,*)' Frequency-dependent diagonal GF3 calculation'
write(*,*)' Frequency-dependent diagonal evGF3 calculation'
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,1X,A15,1X,A1,1X,A15,1X,A1,1X)') &
'|','#','|','e_HF (eV)','|','Z','|','e_GF3 (eV)','|'
'|','#','|','e_HF (eV)','|','Z','|','e_evGF3 (eV)','|'
write(*,*)'-------------------------------------------------------------'
do x=1,nBas
@ -35,10 +35,10 @@ subroutine print_GF3(nBas,nO,nSCF,Conv,eHF,Z,eGF3)
write(*,'(2X,A10,I3)') 'Iteration ',nSCF
write(*,'(2X,A14,F15.5)')'Convergence = ',Conv
write(*,*)'-------------------------------------------------------------'
write(*,'(2X,A27,F15.6)') 'GF3 HOMO energy (eV):',eGF3(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'GF3 LUMO energy (eV):',eGF3(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'GF3 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF3 HOMO energy (eV):',eGF3(HOMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF3 LUMO energy (eV):',eGF3(LUMO)*HaToeV
write(*,'(2X,A27,F15.6)') 'evGF3 HOMO-LUMO gap (eV):',Gap*HaToeV
write(*,*)'-------------------------------------------------------------'
write(*,*)
end subroutine print_GF3
end subroutine print_evGF3

34
src/QuAcK/read_methods.f90
View File

@ -1,12 +1,12 @@
subroutine read_methods(doRHF,doUHF,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, &
do_ring_CCD,do_ladder_CCD, &
doCIS,doRPA,doRPAx, &
doppRPA,doADC, &
doGF2,doGF3, &
doG0W0,doevGW,doqsGW, &
doG0T0,doevGT,doqsGT, &
subroutine read_methods(doRHF,doUHF,doMOM, &
doMP2,doMP3,doMP2F12, &
doCCD,doCCSD,doCCSDT, &
do_ring_CCD,do_ladder_CCD, &
doCIS,doRPA,doRPAx, &
doppRPA,doADC, &
doG0F2,doevGF2,doG0F3,doevGF3, &
doG0W0,doevGW,doqsGW, &
doG0T0,doevGT,doqsGT, &
doMCMP2)
! Read desired methods
@ -20,7 +20,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
logical,intent(out) :: doCCD,doCCSD,doCCSDT
logical,intent(out) :: do_ring_CCD,do_ladder_CCD
logical,intent(out) :: doCIS,doRPA,doRPAx,doppRPA,doADC
logical,intent(out) :: doGF2,doGF3
logical,intent(out) :: doG0F2,doevGF2,doG0F3,doevGF3
logical,intent(out) :: doG0W0,doevGW,doqsGW
logical,intent(out) :: doG0T0,doevGT,doqsGT
logical,intent(out) :: doMCMP2
@ -56,8 +56,10 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
doppRPA = .false.
doADC = .false.
doGF2 = .false.
doGF3 = .false.
doG0F2 = .false.
doevGF2 = .false.
doG0F3 = .false.
doevGF3 = .false.
doG0W0 = .false.
doevGT = .false.
@ -108,9 +110,11 @@ subroutine read_methods(doRHF,doUHF,doMOM, &
! Read Green function methods
read(1,*)
read(1,*) answer1,answer2
if(answer1 == 'T') doGF2 = .true.
if(answer2 == 'T') doGF3 = .true.
read(1,*) answer1,answer2,answer3,answer4
if(answer1 == 'T') doG0F2 = .true.
if(answer2 == 'T') doevGF2 = .true.
if(answer3 == 'T') doG0F3 = .true.
if(answer4 == 'T') doevGF3 = .true.
! Read GW methods

10
src/QuAcK/renormalization_factor_Tmatrix.f90
View File

@ -42,7 +42,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do i=nC+1,nO
do cd=1,nVVs
eps = e(p) + e(i) - Omega1s(cd)
Z(p) = Z(p) - 2d0*rho1s(p,i,cd)**2/eps**2
Z(p) = Z(p) + 1d0*(rho1s(p,i,cd)/eps)**2
enddo
enddo
enddo
@ -53,7 +53,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do a=1,nV-nR
do kl=1,nOOs
eps = e(p) + e(nO+a) - Omega2s(kl)
Z(p) = Z(p) - 2d0*rho2s(p,a,kl)**2/eps**2
Z(p) = Z(p) + 1d0*(rho2s(p,a,kl)/eps)**2
enddo
enddo
enddo
@ -68,7 +68,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do i=nC+1,nO
do cd=1,nVVt
eps = e(p) + e(i) - Omega1t(cd)
Z(p) = Z(p) - 2d0*rho1t(p,i,cd)**2/eps**2
Z(p) = Z(p) + 1d0*(rho1t(p,i,cd)/eps)**2
enddo
enddo
enddo
@ -79,13 +79,13 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do a=1,nV-nR
do kl=1,nOOt
eps = e(p) + e(nO+a) - Omega2t(kl)
Z(p) = Z(p) - 2d0*rho2t(p,a,kl)**2/eps**2
Z(p) = Z(p) + 1d0*(rho2t(p,a,kl)/eps)**2
enddo
enddo
enddo
! Compute renormalization factor from derivative of SigT
Z(:) = 1d0/(1d0 - Z(:))
Z(:) = 1d0/(1d0 + Z(:))
end subroutine renormalization_factor_Tmatrix

18
src/QuAcK/renormalization_factor_Tmatrix_so.f90
View File

@ -39,12 +39,8 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
do p=nC+1,nBas-nR
do i=nC+1,nO
do cd=1,nVV
! do c=nO+1,nBas-nR
! do d=c+1,nBas-nR
! cd = cd + 1
eps = e(p) + e(i) - Omega1(cd)
Z(p) = Z(p) - rho1(p,i,cd)**2/eps**2
! enddo
eps = e(p) + e(i) - Omega1(cd)
Z(p) = Z(p) + (rho1(p,i,cd)/eps)**2
enddo
enddo
enddo
@ -54,18 +50,14 @@ subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omeg
do p=nC+1,nBas-nR
do a=1,nV-nR
do kl=1,nOO
! do k=nC+1,nO
! do l=k+1,nO
! kl = kl + 1
eps = e(p) + e(nO+a) - Omega2(kl)
Z(p) = Z(p) - rho2(p,a,kl)**2/eps**2
! enddo
eps = e(p) + e(nO+a) - Omega2(kl)
Z(p) = Z(p) + (rho2(p,a,kl)/eps)**2
enddo
enddo
enddo
! Compute renormalization factor from derivative of SigT
Z(:) = 1d0/(1d0 - Z(:))
Z(:) = 1d0/(1d0 + Z(:))
end subroutine renormalization_factor_Tmatrix_so

8
src/QuAcK/self_energy_Tmatrix_diag.f90
View File

@ -46,7 +46,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do i=nC+1,nO
do cd=1,nVVs
eps = e(p) + e(i) - Omega1s(cd)
SigT(p) = SigT(p) + 2d0*rho1s(p,i,cd)**2/eps
SigT(p) = SigT(p) + 1d0*rho1s(p,i,cd)**2/eps
enddo
enddo
enddo
@ -57,7 +57,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do a=1,nV-nR
do kl=1,nOOs
eps = e(p) + e(nO+a) - Omega2s(kl)
SigT(p) = SigT(p) + 2d0*rho2s(p,a,kl)**2/eps
SigT(p) = SigT(p) + 1d0*rho2s(p,a,kl)**2/eps
enddo
enddo
enddo
@ -72,7 +72,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do i=nC+1,nO
do cd=1,nVVt
eps = e(p) + e(i) - Omega1t(cd)
SigT(p) = SigT(p) + 2d0*rho1t(p,i,cd)**2/eps
SigT(p) = SigT(p) + 1d0*rho1t(p,i,cd)**2/eps
enddo
enddo
enddo
@ -83,7 +83,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do a=1,nV-nR
do kl=1,nOOt
eps = e(p) + e(nO+a) - Omega2t(kl)
SigT(p) = SigT(p) + 2d0*rho2t(p,a,kl)**2/eps
SigT(p) = SigT(p) + 1d0*rho2t(p,a,kl)**2/eps
enddo
enddo
enddo

16
src/QuAcK/self_energy_Tmatrix_diag_so.f90
View File

@ -43,12 +43,8 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
do p=nC+1,nBas-nR
do i=nC+1,nO
do cd=1,nVV
! do c=nO+1,nBas-nR
! do d=c+1,nBas-nR
! cd = cd + 1
eps = e(p) + e(i) - Omega1(cd)
SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps
! enddo
eps = e(p) + e(i) - Omega1(cd)
SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps
enddo
enddo
enddo
@ -58,12 +54,8 @@ subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho
do p=nC+1,nBas-nR
do a=1,nV-nR
do kl=1,nOO
! do k=nC+1,nO
! do l=k+1,nO
! kl = kl + 1
eps = e(p) + e(nO+a) - Omega2(kl)
SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps
! enddo
eps = e(p) + e(nO+a) - Omega2(kl)
SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps
enddo
enddo
enddo

8
src/QuAcK/soG0T0.f90
View File

@ -89,7 +89,7 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
! Compute excitation densities for the T-matrix
call excitation_density_Tmatrix(ispin,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), &
X1(:,:),Y1(:,:),rho1(:,:,:), &
X1(:,:),Y1(:,:),rho1(:,:,:), &
X2(:,:),Y2(:,:),rho2(:,:,:))
!----------------------------------------------
@ -97,13 +97,15 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
!----------------------------------------------
call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), &
Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), &
Omega1(:),rho1(:,:,:), &
Omega2(:),rho2(:,:,:), &
SigT(:))
! Compute renormalization factor for T-matrix self-energy
call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), &
Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), &
Omega1(:),rho1(:,:,:), &
Omega2(:),rho2(:,:,:), &
Z(:))
!----------------------------------------------