LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

starting implementing B88

This commit is contained in:
Pierre-Francois Loos 2020-03-28 22:52:45 +01:00
parent 59d798afc9
commit 3e72d87adf
17 changed files with 332 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
include/parameters.h
View File

@ -18,5 +18,6 @@
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.904d0*(4d0/3d0)*(1d0/pi)**(1d0/3d0)
double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0)

69
input/basis
View File

@ -1,27 +1,62 @@
1 5
1 14
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
S 1
1 0.0236300 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
2 5
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
P 1
1 0.0848000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
D 1
1 0.1900000 1.0000000
F 1
1 1.3970000 1.0000000
F 1
1 0.3600000 1.0000000
2 14
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
S 1
1 0.0236300 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
P 1
1 0.0848000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
D 1
1 0.1900000 1.0000000
F 1
1 1.3970000 1.0000000
F 1
1 0.3600000 1.0000000

8
input/dft
View File

@ -3,10 +3,10 @@
# exchange rung:
# Hartree = 0
# LDA = 1: RS51,RMFL20
# GGA = 2:
# GGA = 2: RB88
# Hybrid = 4
# Hartree-Fock = 666
666 HF
1 RMFL20
# correlation rung:
# Hartree = 0
# LDA = 1: RVWN5,RMFL20
@ -19,6 +19,6 @@
# Number of states in ensemble (nEns)
2
# Ensemble weights: wEns(1),...,wEns(nEns-1)
0.00000 0.00000
1.00000 0.00000
# GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type
32 0.00001 T 5 1 1
32 0.00001 T 5 1 1

69
input/weight
View File

@ -1,27 +1,62 @@
1 5
1 14
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
S 1
1 0.0236300 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
2 5
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
P 1
1 0.0848000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
D 1
1 0.1900000 1.0000000
F 1
1 1.3970000 1.0000000
F 1
1 0.3600000 1.0000000
2 14
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 82.6400000 0.0020060
2 12.4100000 0.0153430
3 2.8240000 0.0755790
S 1
1 0.1220000 1.0000000
1 0.7977000 1.0000000
S 1
1 0.0297400 1.0000000
1 0.2581000 1.0000000
S 1
1 0.0898900 1.0000000
S 1
1 0.0236300 1.0000000
P 1
1 0.7270000 1.0000000
1 2.2920000 1.0000000
P 1
1 0.1410000 1.0000000
1 0.8380000 1.0000000
P 1
1 0.2920000 1.0000000
P 1
1 0.0848000 1.0000000
D 1
1 2.0620000 1.0000000
D 1
1 0.6620000 1.0000000
D 1
1 0.1900000 1.0000000
F 1
1 1.3970000 1.0000000
F 1
1 0.3600000 1.0000000

2
src/IntPak/IntPak.f90
View File

@ -311,7 +311,7 @@ program IntPak
allocate(DG(1:KG),ExpG(1:KG))
DG = (/ 1d0 /)
ExpG = (/ 1d0 /)
ExpS = ExpS*ExpS
! ExpS = ExpS*ExpS
call cpu_time(start_2eInt(iType))
call Compute2eInt(debug,chemist_notation,iType,nShell, &

50
src/QuAcK/G0T0.f90
View File

@ -102,7 +102,10 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
! Compute excitation densities for the T-matrix
call excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), &
rho1s(:,:,:) = 0d0
rho2s(:,:,:) = 0d0
call excitation_density_Tmatrix(ispin,1d0,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), &
X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:))
!----------------------------------------------
@ -124,7 +127,50 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
! Compute excitation densities for the T-matrix
call excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), &
rho1t(:,:,:) = 0d0
rho2t(:,:,:) = 0d0
call excitation_density_Tmatrix(ispin,1d0,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), &
X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:))
!----------------------------------------------
! Compute T-matrix version of the self-energy
!----------------------------------------------
SigT(:) = 0d0
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(:))
!----------------------------------------------
! 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(:,:,:))
!----------------------------------------------

42
src/QuAcK/excitation_density_Tmatrix.f90
View File

@ -1,4 +1,4 @@
subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
subroutine excitation_density_Tmatrix(ispin,db,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2)
! Compute excitation densities for T-matrix self-energy
@ -7,6 +7,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
! Input variables
integer,intent(in) :: ispin
double precision,intent(in) :: db
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)
@ -27,9 +28,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
double precision,intent(out) :: rho1(nBas,nO,nVV)
double precision,intent(out) :: rho2(nBas,nV,nOO)
rho1(:,:,:) = 0d0
rho2(:,:,:) = 0d0
!----------------------------------------------
! Singlet manifold
!----------------------------------------------
@ -45,10 +43,11 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do c=nO+1,nBas-nR
do d=c,nBas-nR
cd = cd + 1
! rho1(p,i,ab) = rho1(p,i,ab) + ERI(p,i,c,d)*X1(cd,ab)
rho1(p,i,ab) = rho1(p,i,ab) &
+ (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) &
/sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d)))
! + 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
end do
@ -56,10 +55,11 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do k=nC+1,nO
do l=k,nO
kl = kl + 1
! rho1(p,i,ab) = rho1(p,i,ab) + ERI(p,i,k,l)*Y1(kl,ab)
rho1(p,i,ab) = rho1(p,i,ab) &
+ (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) &
/sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l)))
! + 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
@ -73,10 +73,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do c=nO+1,nBas-nR
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)
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) &
/sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d)))
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
end do
@ -84,10 +81,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do k=nC+1,nO
do l=k,nO
kl = kl + 1
! rho2(p,a,ij) = rho2(p,a,ij) + ERI(p,nO+a,k,l)*Y2(kl,ij)
rho2(p,a,ij) = rho2(p,a,ij) &
+ (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) &
/sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l)))
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
@ -113,7 +107,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho1(p,i,ab) = rho1(p,i,ab) + 1.5d0*(ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab)
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
end do
@ -121,7 +117,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho1(p,i,ab) = rho1(p,i,ab) + 1.5d0*(ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab)
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
@ -135,7 +133,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do c=nO+1,nBas-nR
do d=c+1,nBas-nR
cd = cd + 1
rho2(p,a,ij) = rho2(p,a,ij) + 0.5d0*(ERI(p,nO+a,c,d) - ERI(p,nO+a,d,c))*X2(cd,ij)
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
end do
@ -143,7 +141,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r
do k=nC+1,nO
do l=k+1,nO
kl = kl + 1
rho2(p,a,ij) = rho2(p,a,ij) + 0.5d0*(ERI(p,nO+a,k,l) - ERI(p,nO+a,l,k))*Y2(kl,ij)
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

4
src/QuAcK/self_energy_Tmatrix_diag.f90
View File

@ -32,10 +32,6 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e,
double precision,intent(out) :: SigT(nBas)
! Initialize
SigT(:) = 0d0
!----------------------------------------------
! Singlet part of the T-matrix self-energy
!----------------------------------------------

5
src/QuAcK/soG0T0.f90
View File

@ -88,7 +88,10 @@ 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(:,:,:,:), &
rho1(:,:,:) = 0d0
rho2(:,:,:) = 0d0
call excitation_density_Tmatrix(ispin,1d0,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), &
X1(:,:),Y1(:,:),rho1(:,:,:),X2(:,:),Y2(:,:),rho2(:,:,:))
!----------------------------------------------

9
src/eDFT/B88_gga_exchange_energy.f90
View File

@ -1,6 +1,6 @@
subroutine B88_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
subroutine B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Compute Becke's 96 GGA exchange energy
! Compute Becke's 88 GGA exchange energy
implicit none
@ -8,9 +8,6 @@ subroutine B88_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
! Input variables
character(len=12),intent(in) :: DFA
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)
@ -19,7 +16,7 @@ subroutine B88_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
! Local variables
integer :: iG
double precision :: alpha,beta,gamma
double precision :: alpha,beta
double precision :: r,g,x
! Output variables

5
src/eDFT/G96_gga_exchange_energy.f90
View File

@ -1,4 +1,4 @@
subroutine G96_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
subroutine G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Compute Gill's 96 GGA exchange energy
@ -8,9 +8,6 @@ subroutine G96_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
! Input variables
character(len=12),intent(in) :: DFA
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)

5
src/eDFT/G96_gga_exchange_potential.f90
View File

@ -1,4 +1,4 @@
subroutine G96_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
subroutine G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Compute Gill's GGA exchange poential
@ -7,9 +7,6 @@ subroutine G96_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho
! Input variables
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas

18
src/eDFT/LIM_RKS.f90
View File

@ -73,18 +73,16 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
!------------------------------------------------------------------------
write(*,'(A40)') '*************************************************'
write(*,'(A40)') ' EQUI-WEIGHT CALCULATION '
write(*,'(A40)') ' NON-ZERO-WEIGHT CALCULATION '
write(*,'(A40)') '*************************************************'
wLIM(1:nEns) = 1d0/dble(nEns)
do iEns=1,nEns
write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wLIM(iEns)
write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wEns(iEns)
end do
write(*,'(A40)') '*************************************************'
write(*,*)
call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,nEns,wLIM,nGrid,weight,maxSCF,thresh, &
call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thresh, &
max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,EwGICEW,F)
!------------------------------------------------------------------------
@ -92,10 +90,14 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres
!------------------------------------------------------------------------
Om(1) = 0d0
Om(2) = 2d0*(EwEW - EwZW)
OmGIC(1) = 0d0
OmGIC(2) = 2d0*(EwGICEW - EwGICZW)
Om(2) = 0d0
OmGIC(2) = 0d0
if(wEns(2) > 10d-3) then
Om(2) = (EwEW - EwZW)/wEns(2)
OmGIC(2) = (EwGICEW - EwGICZW)/wEns(2)
end if
write(*,'(A60)') '-------------------------------------------------'
write(*,'(A60)') ' LINEAR INTERPOLATION METHOD EXCITATION ENERGIES '

53
src/eDFT/RB88_gga_exchange_energy.f90 Normal file
View File

@ -0,0 +1,53 @@
subroutine RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Compute restricted version of Becke's 88 GGA exchange energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Local variables
integer :: iG
double precision :: alpha
double precision :: beta
double precision :: r,g,x
! Output variables
double precision :: Ex
! Coefficients for B88 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
beta = 0.0042d0
! Compute GGA exchange energy
Ex = 0d0
do iG=1,nGrid
r = max(0d0,0.5d0*rho(iG))
if(r > threshold) then
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
x = sqrt(g)/r**(4d0/3d0)
Ex = Ex + weight(iG)*alpha*r**(4d0/3d0) &
- weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x))
end if
end do
Ex = 2d0*Ex
end subroutine RB88_gga_exchange_energy

66
src/eDFT/RB88_gga_exchange_potential.f90 Normal file
View File

@ -0,0 +1,66 @@
subroutine RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Compute restricted Becke's GGA exchange potential
implicit none
include 'parameters.h'
! Input variables
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) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Local variables
integer :: mu,nu,iG
double precision :: alpha
double precision :: beta
double precision :: r,g,vAO,gAO
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Coefficients for B88 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
beta = 0.0042d0
! Compute GGA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,0.5d0*rho(iG))
if(r > threshold) then
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) &
+ vAO*(4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) &
+ 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0))
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gAO = weight(iG)*gAO
Fx(mu,nu) = Fx(mu,nu) - 2d0*gAO*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0)
end if
end do
end do
end do
end subroutine RB88_gga_exchange_potential

10
src/eDFT/gga_exchange_energy.f90
View File

@ -26,13 +26,19 @@ subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
case ('G96')
call G96_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
call G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Becke's 88 exchange functional
case ('RB88')
call RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Becke's 88 exchange functional
case ('B88')
call B88_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
call B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
case default

10
src/eDFT/gga_exchange_potential.f90
View File

@ -30,13 +30,19 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh
case ('G96')
call G96_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
call G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Becke's 88 exchange functional
case ('RB88')
call RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Becke's 88 exchange functional
case ('B88')
call B88_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
call B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
case default