LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-10-31 19:23:52 +01:00
Code Releases Activity

GIC functional

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-04-06 19:40:42 +02:00
parent 5a80b38685
commit 192a6345de
26 changed files with 602 additions and 283 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
include/parameters.h
View File

@ -18,6 +18,5 @@
double precision,parameter :: CxLDA = - (3d0/4d0)*(3d0/pi)**(1d0/3d0)
double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0)
! double precision,parameter :: Cx1 = - 0.913d0*(4d0/3d0)*(1d0/pi)**(1d0/3d0)
double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)

48
input/basis
View File

@ -1,27 +1,43 @@
1 5
1 9
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 33.8700000 0.0060680
2 5.0950000 0.0453080
3 1.1590000 0.2028220
S 1
1 0.1220000 1.0000000
1 0.3258000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.1027000 1.0000000
S 1
1 0.0252600 1.0000000
P 1
1 0.7270000 1.0000000
1 1.4070000 1.0000000
P 1
1 0.1410000 1.0000000
2 5
1 0.3880000 1.0000000
P 1
1 0.1020000 1.0000000
D 1
1 1.0570000 1.0000000
D 1
1 0.2470000 1.0000000
2 9
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 33.8700000 0.0060680
2 5.0950000 0.0453080
3 1.1590000 0.2028220
S 1
1 0.1220000 1.0000000
1 0.3258000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.1027000 1.0000000
S 1
1 0.0252600 1.0000000
P 1
1 0.7270000 1.0000000
1 1.4070000 1.0000000
P 1
1 0.1410000 1.0000000
1 0.3880000 1.0000000
P 1
1 0.1020000 1.0000000
D 1
1 1.0570000 1.0000000
D 1
1 0.2470000 1.0000000

4
input/dft
View File

@ -2,11 +2,11 @@
GOK-RKS
# exchange rung:
# Hartree = 0
# LDA = 1: RS51,RMFL20
# LDA = 1: RS51,RMFL20,RGIC
# GGA = 2: RB88
# Hybrid = 4
# Hartree-Fock = 666
666 HF
1 RGIC
# correlation rung:
# Hartree = 0
# LDA = 1: RVWN5,RMFL20

48
input/weight
View File

@ -1,27 +1,43 @@
1 5
1 9
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 33.8700000 0.0060680
2 5.0950000 0.0453080
3 1.1590000 0.2028220
S 1
1 0.1220000 1.0000000
1 0.3258000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.1027000 1.0000000
S 1
1 0.0252600 1.0000000
P 1
1 0.7270000 1.0000000
1 1.4070000 1.0000000
P 1
1 0.1410000 1.0000000
2 5
1 0.3880000 1.0000000
P 1
1 0.1020000 1.0000000
D 1
1 1.0570000 1.0000000
D 1
1 0.2470000 1.0000000
2 9
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 33.8700000 0.0060680
2 5.0950000 0.0453080
3 1.1590000 0.2028220
S 1
1 0.1220000 1.0000000
1 0.3258000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.1027000 1.0000000
S 1
1 0.0252600 1.0000000
P 1
1 0.7270000 1.0000000
1 1.4070000 1.0000000
P 1
1 0.1410000 1.0000000
1 0.3880000 1.0000000
P 1
1 0.1020000 1.0000000
D 1
1 1.0570000 1.0000000
D 1
1 0.2470000 1.0000000

69
src/QuAcK/G0T0.f90
View File

@ -45,6 +45,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
double precision,allocatable :: X2s(:,:),X2t(:,:)
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
double precision,allocatable :: rho1st(:,:,:),rho2st(:,:,:)
double precision,allocatable :: SigT(:)
double precision,allocatable :: Z(:)
@ -81,6 +82,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nO,nVVt),rho2t(nBas,nV,nOOt), &
rho1st(nBas,nO,nVVt),rho2st(nBas,nV,nOOt), &
SigT(nBas),Z(nBas),eG0T0(nBas))
!----------------------------------------------
@ -100,14 +102,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s(:))
call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s(:))
! Compute excitation densities for the T-matrix
rho1s(:,:,:) = 0d0
rho2s(:,:,:) = 0d0
call excitation_density_Tmatrix(ispin,1d0,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), &
X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:))
!----------------------------------------------
! Triplet manifold
!----------------------------------------------
@ -125,72 +119,33 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t(:))
call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t(:))
! Compute excitation densities for the T-matrix
!-----------------------------------------------------------
! Compute excitation densities for the T-matrix self-energy
!-----------------------------------------------------------
rho1t(:,:,:) = 0d0
rho2t(:,:,:) = 0d0
call excitation_density_Tmatrix(ispin,1d0,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), &
call excitation_density_Tmatrix(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,ERI(:,:,:,:),rho1st(:,:,:),rho2st(:,:,:), &
X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:), &
X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:))
!----------------------------------------------
! Compute T-matrix version of the self-energy
!----------------------------------------------
SigT(:) = 0d0
rho2s(:,:,:) = 0d0
rho2t(:,:,:) = 0d0
! rho2s(:,:,:) = 0d0
! rho2t(:,:,:) = 0d0
! rho2st(:,:,:) = 0d0
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), &
Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), &
Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), &
SigT(:))
!----------------------------------------------
! Singlet manifold
!----------------------------------------------
ispin = 1
! Compute excitation densities for the T-matrix
rho1s(:,:,:) = 0d0
rho2s(:,:,:) = 0d0
call excitation_density_Tmatrix(ispin,0d0,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), &
X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:))
!----------------------------------------------
! Triplet manifold
!----------------------------------------------
ispin = 2
! Compute excitation densities for the T-matrix
rho1t(:,:,:) = 0d0
rho2t(:,:,:) = 0d0
call excitation_density_Tmatrix(ispin,0d0,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), &
X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:))
!----------------------------------------------
! Compute T-matrix version of the self-energy
!----------------------------------------------
rho2s(:,:,:) = 0d0
rho2t(:,:,:) = 0d0
call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), &
Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), &
Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), &
SigT(:))
rho1st(:,:,:),rho2st(:,:,:),SigT(:))
! Compute renormalization factor for T-matrix self-energy
call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), &
Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), &
Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), &
Z(:))
rho1st(:,:,:),rho2st(:,:,:),Z(:))
!----------------------------------------------
! Solve the quasi-particle equation

302
src/QuAcK/excitation_density_Tmatrix.f90
View File

@ -1,4 +1,5 @@
subroutine excitation_density_Tmatrix(ispin,db,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
subroutine excitation_density_Tmatrix(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,ERI,rho1st,rho2st, &
X1s,Y1s,rho1s,X2s,Y2s,rho2s,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
! Compute excitation densities for T-matrix self-energy
@ -6,14 +7,24 @@ subroutine excitation_density_Tmatrix(ispin,db,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y
! Input variables
integer,intent(in) :: ispin
double precision,intent(in) :: db
integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: X1(nVV,nVV)
double precision,intent(in) :: Y1(nOO,nVV)
double precision,intent(in) :: X2(nVV,nOO)
double precision,intent(in) :: Y2(nOO,nOO)
integer,intent(in) :: nOOs
integer,intent(in) :: nOOt
double precision,intent(in) :: X1s(nVVs,nVVs)
double precision,intent(in) :: Y1s(nOOs,nVVs)
double precision,intent(in) :: X2s(nVVs,nOOs)
double precision,intent(in) :: Y2s(nOOs,nOOs)
integer,intent(in) :: nVVs
integer,intent(in) :: nVVt
double precision,intent(in) :: X1t(nVVt,nVVt)
double precision,intent(in) :: Y1t(nOOt,nVVt)
double precision,intent(in) :: X2t(nVVt,nOOt)
double precision,intent(in) :: Y2t(nOOt,nOOt)
! Local variables
@ -25,187 +36,176 @@ subroutine excitation_density_Tmatrix(ispin,db,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y
! Output variables
double precision,intent(out) :: rho1(nBas,nO,nVV)
double precision,intent(out) :: rho2(nBas,nV,nOO)
double precision,intent(out) :: rho1s(nBas,nO,nVVs)
double precision,intent(out) :: rho2s(nBas,nV,nOOs)
double precision,intent(out) :: rho1t(nBas,nO,nVVt)
double precision,intent(out) :: rho2t(nBas,nV,nOOt)
double precision,intent(out) :: rho1st(nBas,nO,nVVt)
double precision,intent(out) :: rho2st(nBas,nV,nOOt)
! Initialization
rho1s(:,:,:) = 0d0
rho2s(:,:,:) = 0d0
rho1t(:,:,:) = 0d0
rho2t(:,:,:) = 0d0
rho1st(:,:,:) = 0d0
rho2st(:,:,:) = 0d0
!----------------------------------------------
! Singlet manifold
!----------------------------------------------
if(ispin == 1) then
do p=nC+1,nBas-nR
do p=nC+1,nBas-nR
do i=nC+1,nO
do ab=1,nVVs
do i=nC+1,nO
do ab=1,nVV
cd = 0
do c=nO+1,nBas-nR
do d=c,nBas-nR
cd = cd + 1
rho1(p,i,ab) = rho1(p,i,ab) &
! + db*(ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) &
! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) &
+ 0d0*db*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) &
+ (1d0 - db)*ERI(p,i,c,d)*X1(cd,ab)
end do
cd = 0
do c=nO+1,nBas-nR
do d=c,nBas-nR
cd = cd + 1
rho1s(p,i,ab) = rho1s(p,i,ab) + 2.0d0*ERI(p,i,c,d)*X1s(cd,ab)
end do
kl = 0
do k=nC+1,nO
do l=k,nO
kl = kl + 1
rho1(p,i,ab) = rho1(p,i,ab) &
! + db*(ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) &
! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) &
+ 0d0*db*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) &
+ (1d0 - db)*ERI(p,i,k,l)*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,nBas-nR
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) + db*(ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
end do
kl = 0
do k=nC+1,nO
do l=k,nO
kl = kl + 1
rho1s(p,i,ab) = rho1s(p,i,ab) + 2.0d0*ERI(p,i,k,l)*Y1s(kl,ab)
end do
kl = 0
do k=nC+1,nO
do l=k,nO
kl = kl + 1
rho2(p,a,ij) = rho2(p,a,ij) + db*(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 do
end if
do a=1,nV-nR
do ij=1,nOOs
cd = 0
do c=nO+1,nBas-nR
do d=c,nBas-nR
cd = cd + 1
rho2s(p,a,ij) = rho2s(p,a,ij) + 2.0d0*ERI(p,nO+a,c,d)*X2s(cd,ij)
end do
end do
kl = 0
do k=nC+1,nO
do l=k,nO
kl = kl + 1
rho2s(p,a,ij) = rho2s(p,a,ij) + 2.0d0*ERI(p,nO+a,k,l)*Y2s(kl,ij)
end do
end do
end do
end do
end do
!----------------------------------------------
! Triplet manifold
!----------------------------------------------
if(ispin == 2) then
do p=nC+1,nBas-nR
do p=nC+1,nBas-nR
do i=nC+1,nO
do ab=1,nVVt
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) &
+ 1.0d0*db*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) &
+ (1d0-db)*0d0*(ERI(p,i,c,d))*X1(cd,ab)
end do
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho1t(p,i,ab) = rho1t(p,i,ab) + 1d0*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1t(cd,ab)
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) &
+ 1.0d0*db*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) &
+ (1d0-db)*0d0*(ERI(p,i,k,l))*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) + 1d0*(ERI(p,nO+a,c,d) - db*ERI(p,nO+a,d,c))*X2(cd,ij)
end do
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho1t(p,i,ab) = rho1t(p,i,ab) + 1d0*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1t(kl,ab)
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) + 1d0*(ERI(p,nO+a,k,l) - db*ERI(p,nO+a,l,k))*Y2(kl,ij)
end do
end do
end do
end do
end do
end do
end if
do a=1,nV-nR
do ij=1,nOOt
!----------------------------------------------
! 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
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho2t(p,a,ij) = rho2t(p,a,ij) + 1d0*(ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2t(cd,ij)
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
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho2t(p,a,ij) = rho2t(p,a,ij) + 1d0*(ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2t(kl,ij)
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
end do
end if
end do
!----------------------------------------------
! Singlet-triplet crossed term
!----------------------------------------------
do p=nC+1,nBas-nR
do i=nC+1,nO
do ab=1,nVVt
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho1st(p,i,ab) = rho1st(p,i,ab) + 2d0*ERI(p,i,c,d)*X1t(cd,ab)
end do
end do
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho1st(p,i,ab) = rho1st(p,i,ab) + 2d0*ERI(p,i,k,l)*Y1t(kl,ab)
end do
end do
end do
end do
do a=1,nV-nR
do ij=1,nOOt
cd = 0
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho2st(p,a,ij) = rho2st(p,a,ij) + 2d0*ERI(p,nO+a,c,d)*X2t(cd,ij)
end do
end do
kl = 0
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho2st(p,a,ij) = rho2st(p,a,ij) + 2d0*ERI(p,nO+a,k,l)*Y2t(kl,ij)
end do
end do
end do
end do
end do
end subroutine excitation_density_Tmatrix

82
src/QuAcK/excitation_density_Tmatrix_so.f90 Normal file
View File

@ -0,0 +1,82 @@
subroutine excitation_density_Tmatrix_so(nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
! Compute excitation densities for T-matrix self-energy
implicit none
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: X1(nVV,nVV)
double precision,intent(in) :: Y1(nOO,nVV)
double precision,intent(in) :: X2(nVV,nOO)
double precision,intent(in) :: Y2(nOO,nOO)
! Local variables
integer :: i,j,k,l
integer :: a,b,c,d
integer :: p
integer :: ab,cd,ij,kl
double precision,external :: Kronecker_delta
! Output variables
double precision,intent(out) :: rho1(nBas,nO,nVV)
double precision,intent(out) :: rho2(nBas,nV,nOO)
! Initialization
rho1(:,:,:) = 0d0
rho2(:,:,:) = 0d0
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
end subroutine excitation_density_Tmatrix_so

2
src/QuAcK/linear_response_A_matrix.f90
View File

@ -49,8 +49,6 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ (1d0 + delta_spin)*lambda*ERI(i,b,a,j) &
- (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
! + (1d0 + delta_spin)*lambda*ERI(i,j,a,b) &
! - (1d0 - delta_dRPA)*lambda*ERI(i,a,j,b)
enddo
enddo

2
src/QuAcK/linear_response_B_matrix.f90
View File

@ -46,8 +46,6 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,j,a,b) &
- (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
! B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,b,a,j) &
! - (1d0 - delta_dRPA)*lambda*ERI(i,a,b,j)
enddo
enddo

8
src/QuAcK/renormalization_factor_Tmatrix.f90
View File

@ -1,6 +1,6 @@
subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, &
Omega1s,rho1s,Omega2s,rho2s,Omega1t,rho1t,Omega2t,rho2t, &
Z)
rho1st,rho2st,Z)
! Compute renormalization factor of the T-matrix self-energy
@ -15,9 +15,9 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
integer,intent(in) :: nVVs,nVVt
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1s(nVVs),Omega1t(nVVt)
double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt)
double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt),rho1st(nBas,nO,nVVt)
double precision,intent(in) :: Omega2s(nOOs),Omega2t(nOOt)
double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt)
double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt),rho2st(nBas,nV,nOOt)
! Local variables
@ -69,6 +69,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do cd=1,nVVt
eps = e(p) + e(i) - Omega1t(cd)
Z(p) = Z(p) + (rho1t(p,i,cd)/eps)**2
Z(p) = Z(p) + (rho1st(p,i,cd)/eps)**2
enddo
enddo
enddo
@ -80,6 +81,7 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV
do kl=1,nOOt
eps = e(p) + e(nO+a) - Omega2t(kl)
Z(p) = Z(p) + (rho2t(p,a,kl)/eps)**2
Z(p) = Z(p) + (rho2st(p,a,kl)/eps)**2
enddo
enddo
enddo

12
src/QuAcK/self_energy_Tmatrix_diag.f90
View File

@ -1,6 +1,6 @@
subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, &
Omega1s,rho1s,Omega2s,rho2s,Omega1t,rho1t,Omega2t,rho2t, &
SigT)
rho1st,rho2st,SigT)
! Compute diagonal of the correlation part of the T-matrix self-energy
@ -19,9 +19,9 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
integer,intent(in) :: nVVs,nVVt
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: Omega1s(nVVs),Omega1t(nVVt)
double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt)
double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt),rho1st(nBas,nO,nVVt)
double precision,intent(in) :: Omega2s(nOOs),Omega2t(nOOt)
double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt)
double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt),rho2st(nBas,nV,nOOt)
! Local variables
@ -32,6 +32,10 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
double precision,intent(out) :: SigT(nBas)
! Initialization
SigT(:) = 0d0
!----------------------------------------------
! Singlet part of the T-matrix self-energy
!----------------------------------------------
@ -69,6 +73,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do cd=1,nVVt
eps = e(p) + e(i) - Omega1t(cd)
SigT(p) = SigT(p) + rho1t(p,i,cd)**2/eps
SigT(p) = SigT(p) + rho1st(p,i,cd)**2/eps
enddo
enddo
enddo
@ -80,6 +85,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
do kl=1,nOOt
eps = e(p) + e(nO+a) - Omega2t(kl)
SigT(p) = SigT(p) + rho2t(p,a,kl)**2/eps
SigT(p) = SigT(p) + rho2st(p,a,kl)**2/eps
enddo
enddo
enddo

4
src/QuAcK/soG0T0.f90
View File

@ -91,14 +91,14 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF)
rho1(:,:,:) = 0d0
rho2(:,:,:) = 0d0
call excitation_density_Tmatrix(ispin,1d0,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), &
call excitation_density_Tmatrix_so(nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), &
X1(:,:),Y1(:,:),rho1(:,:,:),X2(:,:),Y2(:,:),rho2(:,:,:))
!----------------------------------------------
! Compute T-matrix version of the self-energy
!----------------------------------------------
rho2(:,:,:) = 0d0
! rho2(:,:,:) = 0d0
call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), &
Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), &

2
src/eDFT/Makefile
View File

@ -3,7 +3,7 @@ BDIR =../../bin
ODIR = obj
OODIR = ../IntPak/obj
SDIR =.
FC = gfortran -I$(IDIR) -Wall
FC = gfortran -I$(IDIR)
ifeq ($(DEBUG),1)
FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
else

65
src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 Normal file
View File

@ -0,0 +1,65 @@
subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,ExDD)
! Compute the restricted version of the GIC exchange individual energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
! Local variables
integer :: iEns,jEns
integer :: iG
double precision :: r
double precision :: dExdw(nEns)
double precision,external :: Kronecker_delta
double precision :: a,b,c,w
double precision :: dCxGICdw
! Output variables
double precision,intent(out) :: ExDD(nEns)
! Compute correlation energy for ground- and doubly-excited states
a = + 0.5751782560799208d0
b = - 0.021108186591137282d0
c = - 0.36718902716347124d0
w = wEns(2)
dCxGICdw = CxLDA*(0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0)))
dExdw(:) = 0d0
do iG=1,nGrid
r = max(0d0,rhow(iG))
if(r > threshold) then
dExdw(1) = 0d0
dExdw(2) = dExdw(2) + weight(iG)*dCxGICdw*r**(4d0/3d0)
end if
end do
ExDD(:) = 0d0
do iEns=1,nEns
do jEns=2,nEns
ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns)
end do
end do
end subroutine RGIC_lda_exchange_derivative_discontinuity

52
src/eDFT/RGIC_lda_exchange_energy.f90 Normal file
View File

@ -0,0 +1,52 @@
subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
! Compute the restricted version of the GIC exchange functional
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: iG
double precision :: r
double precision :: a,b,c,w
double precision :: CxGIC
! Output variables
double precision :: Ex
! Weight-denepdent Cx coefficient
a = + 0.5751782560799208d0
b = - 0.021108186591137282d0
c = - 0.36718902716347124d0
w = wEns(2)
CxGIC = CxLDA*w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)
! Compute GIC-LDA exchange energy
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
Ex = Ex + weight(iG)*CxLDA*r**(4d0/3d0)
Ex = Ex + weight(iG)*CxGIC*r**(4d0/3d0)
endif
enddo
end subroutine RGIC_lda_exchange_energy

61
src/eDFT/RGIC_lda_exchange_individual_energy.f90 Normal file
View File

@ -0,0 +1,61 @@
subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex)
! Compute the restricted version of the GIC exchange functional
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: iG
double precision :: CxGIC
double precision :: r,rI
double precision :: e_p,dedr
double precision :: a,b,c,w
! Output variables
double precision,intent(out) :: Ex
! Weight-dependent Cx coefficient for RMFL20 exchange functional
a = + 0.5751782560799208d0
b = - 0.021108186591137282d0
c = - 0.36718902716347124d0
w = wEns(2)
CxGIC = CxLDA*w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)
! Compute LDA exchange matrix in the AO basis
Ex = 0d0
do iG=1,nGrid
r = max(0d0,rhow(iG))
rI = max(0d0,rho(iG))
if(r > threshold .and. rI > threshold) then
e_p = CxLDA*r**(1d0/3d0)
dedr = 1d0/3d0*CxLDA*r**(-2d0/3d0)
Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r)
e_p = CxGIC*r**(1d0/3d0)
dedr = 1d0/3d0*CxGIC*r**(-2d0/3d0)
Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r)
endif
enddo
end subroutine RGIC_lda_exchange_individual_energy

61
src/eDFT/RGIC_lda_exchange_potential.f90 Normal file
View File

@ -0,0 +1,61 @@
subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
! Compute the restricted version of the GIC exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: mu,nu,iG
double precision :: r,vAO
double precision :: a,b,c,w
double precision :: CxGIC
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Weight-dependent Cx coefficient for RMFL20 exchange functional
a = + 0.5751782560799208d0
b = - 0.021108186591137282d0
c = - 0.36718902716347124d0
w = wEns(2)
CxGIC = CxLDA*w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)
! Compute LDA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxLDA*r**(1d0/3d0)
Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*CxGIC*r**(1d0/3d0)
endif
enddo
enddo
enddo
end subroutine RGIC_lda_exchange_potential

10
src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90
View File

@ -34,8 +34,10 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r
r = max(0d0,rhow(iG))
if(r > threshold) then
dExdw(1) = dExdw(1) + weight(iG)*Cx0*r**(4d0/3d0)
dExdw(2) = dExdw(2) + weight(iG)*Cx1*r**(4d0/3d0)
dExdw(1) = 0d0
dExdw(2) = dExdw(2) + weight(iG)*(Cx1 - Cx0)*r**(4d0/3d0)
end if
end do
@ -43,9 +45,9 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r
ExDD(:) = 0d0
do iEns=1,nEns
do jEns=1,nEns
do jEns=2,nEns
ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*(dExdw(jEns) - dExdw(1))
ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns)
end do
end do

2
src/eDFT/RMFL20_lda_exchange_energy.f90
View File

@ -28,7 +28,7 @@ subroutine RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex
! Weight-denepdent Cx coefficient
if(LDA_centered) then
Cxw = CxLDA + (Cx1 - Cx0)*wEns(2)
Cxw = CxLDA + (Cx1 - Cx0)*wEns(2)
else
Cxw = wEns(1)*Cx0 + wEns(2)*Cx1
end if

6
src/eDFT/RMFL20_lda_exchange_individual_energy.f90
View File

@ -20,7 +20,7 @@ subroutine RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,we
integer :: iG
double precision :: Cxw
double precision :: r,rI
double precision :: e,dedr
double precision :: e_p,dedr
! Output variables
@ -44,9 +44,9 @@ subroutine RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,we
if(r > threshold .and. rI > threshold) then
e = Cxw*r**(1d0/3d0)
e_p = Cxw*r**(1d0/3d0)
dedr = 1d0/3d0*Cxw*r**(-2d0/3d0)
Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r)
Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r)
endif

6
src/eDFT/fock_exchange_individual_energy.f90
View File

@ -24,8 +24,8 @@ subroutine fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex)
allocate(Fx(nBas,nBas))
call fock_exchange_potential(nBas,P(:,:),ERI(:,:,:,:),Fx(:,:))
Ex = 0.5d0*trace_matrix(nBas,matmul(P(:,:),Fx(:,:))) !&
call fock_exchange_potential(nBas,Pw(:,:),ERI(:,:,:,:),Fx(:,:))
Ex = trace_matrix(nBas,matmul(P(:,:),Fx(:,:))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),Fx(:,:)))
end subroutine fock_exchange_individual_energy

10
src/eDFT/hartree_individual_energy.f90
View File

@ -28,8 +28,9 @@ subroutine hartree_individual_energy(rung,nBas,ERI,J,Pw,P,EJ)
case(0)
call hartree_coulomb(nBas,P(:,:),ERI(:,:,:,:),J(:,:))
EJ = 0.5d0*trace_matrix(nBas,matmul(P(:,:),J(:,:)))
call hartree_coulomb(nBas,Pw(:,:),ERI(:,:,:,:),J(:,:))
EJ = trace_matrix(nBas,matmul(P(:,:),J(:,:))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),J(:,:)))
! LDA functionals
@ -57,8 +58,9 @@ subroutine hartree_individual_energy(rung,nBas,ERI,J,Pw,P,EJ)
case(666)
call hartree_coulomb(nBas,P(:,:),ERI(:,:,:,:),J(:,:))
EJ = 0.5d0*trace_matrix(nBas,matmul(P(:,:),J(:,:)))
call hartree_coulomb(nBas,Pw(:,:),ERI(:,:,:,:),J(:,:))
EJ = trace_matrix(nBas,matmul(P(:,:),J(:,:))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),J(:,:)))
end select

4
src/eDFT/lda_exchange_derivative_discontinuity.f90
View File

@ -37,6 +37,10 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow
call RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:))
case ('RGIC')
call RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:))
case default
call print_warning('!!! LDA exchange functional not available !!!')

10
src/eDFT/lda_exchange_energy.f90
View File

@ -23,24 +23,22 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex)
select case (DFA)
! Slater's LDA correlation functional
case ('S51')
call S51_lda_exchange_energy(nGrid,weight,rho,Ex)
! Restricted version of Slater's LDA correlation functional
case ('RS51')
call RS51_lda_exchange_energy(nGrid,weight,rho,Ex)
! Restricted version of the weight-dependent Marut-Fromager-Loos 2020 exchange functional
case ('RMFL20')
call RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex)
case ('RGIC')
call RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
case default
call print_warning('!!! LDA exchange functional not available !!!')

4
src/eDFT/lda_exchange_individual_energy.f90
View File

@ -32,6 +32,10 @@ subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weigh
call RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex)
case ('RGIC')
call RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex)
case default
call print_warning('!!! LDA correlation functional not available !!!')

10
src/eDFT/lda_exchange_potential.f90
View File

@ -26,24 +26,22 @@ subroutine lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,A
select case (DFA)
! Restricted version of Slater's LDA correlation functional
case ('RS51')
call RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
! Slater's LDA correlation functional
case ('S51')
call S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
! Restricted version of the weight-dependent Marut-Fromager-Loos 2020 functional
case ('RMFL20')
call RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
case ('RGIC')
call RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
case default
call print_warning('!!! LDA exchange functional not available !!!')