PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-07-23 11:17:38 +02:00
Code Releases Activity

fix bug in LDA shift

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-31 22:15:45 +02:00
parent 95e10e2ee8
commit 196ac67c26
21 changed files with 102 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/dft
View File

@ -6,7 +6,7 @@
# GGA = 2: RB88
# Hybrid = 4
# Hartree-Fock = 666
1 RMFL20
1 RS51
# correlation rung:
# Hartree = 0
# LDA = 1: RVWN5,RMFL20
@ -15,7 +15,7 @@
# Hartree-Fock = 666
1 RMFL20
# quadrature grid SG-n
3
1
# Number of states in ensemble (nEns)
2
# Ensemble weights: wEns(1),...,wEns(nEns-1)

5
src/eDFT/B88_gga_exchange_potential.f90
View File

@ -1,4 +1,4 @@
subroutine B88_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
subroutine B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Compute Becke's GGA exchange potential
@ -7,9 +7,6 @@ subroutine B88_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

14
src/eDFT/GOK_RKS.f90
View File

@ -236,7 +236,7 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS
call restricted_correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:), &
nBas,AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fc(:,:))
! print*,'Done with restricted_correlation_potential'
! Build Fock operator
F(:,:) = Hc(:,:) + J(:,:) + Fx(:,:) + Fc(:,:)
@ -256,11 +256,11 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS
if(abs(rcond) < 1d-15) n_diis = 0
! Transform Fock matrix in orthogonal basis
! Transform Fock matrix in orthogonal basis
Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:)))
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp(:,:),eps(:))
@ -335,9 +335,9 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxS
! Compute individual energies from ensemble energy
!------------------------------------------------------------------------
call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas, &
AO(:,:),dAO(:,:,:),nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc, &
eps(:),Pw(:,:),rhow(:),drhow(:,:),J(:,:),Fx(:,:),FxHF(:,:), &
Fc(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,EwGIC,E(:),Om(:))
call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:), &
nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc, &
eps(:),Pw(:,:),rhow(:),drhow(:,:),J(:,:),P(:,:,:), &
rho(:,:),drho(:,:,:),Ew,EwGIC,E(:),Om(:))
end subroutine GOK_RKS

2
src/eDFT/MFL20_lda_correlation_derivative_discontinuity.f90
View File

@ -46,7 +46,7 @@ subroutine MFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight
do iEns=1,nEns
call elda_correlation_energy(nEns,aMFL(:,iEns),nGrid,weight(:),rhow(:,:),dEc(:,iEns))
call elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rhow(:,:),dEc(:,iEns))
end do

2
src/eDFT/MFL20_lda_correlation_individual_energy.f90
View File

@ -48,7 +48,7 @@ subroutine MFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,r
do iEns=1,nEns
call elda_correlation_individual_energy(nEns,aMFL(:,iEns),nGrid,weight(:),rhow(:,:),rho(:,:),EceLDA(:,iEns))
call elda_correlation_individual_energy(aMFL(:,iEns),nGrid,weight(:),rhow(:,:),rho(:,:),EceLDA(:,iEns))
end do

2
src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90
View File

@ -41,7 +41,7 @@ subroutine RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weigh
! Compute correlation energy for ground, singly-excited and doubly-excited states
do iEns=1,nEns
call restricted_elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rhow(:),dEcdw(iEns))
call restricted_elda_correlation_energy(aMFL(1:3,iEns),nGrid,weight(:),rhow(:),dEcdw(iEns))
end do
EcDD(:) = 0d0

11
src/eDFT/RMFL20_lda_correlation_energy.f90
View File

@ -43,14 +43,19 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec)
! Compute correlation energy for ground and doubly-excited states
do iEns=1,nEns
call restricted_elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rho(:),EceLDA(iEns))
call restricted_elda_correlation_energy(aMFL(1:3,iEns),nGrid,weight(:),rho(:),EceLDA(iEns))
end do
! LDA-centered functional
call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),EcLDA)
if(LDA_centered) then
if(LDA_centered) EceLDA(:) = EceLDA(:) + EcLDA - EceLDA(1)
call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),EcLDA)
EceLDA(2) = EcLDA + wEns(2)*(EceLDA(2) - EceLDA(1))
EceLDA(1) = EcLDA
end if
! Weight-denpendent functional for ensembles

12
src/eDFT/RMFL20_lda_correlation_individual_energy.f90
View File

@ -43,14 +43,20 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,
! Compute correlation energy for ground- and doubly-excited states
do iEns=1,nEns
call restricted_elda_correlation_individual_energy(nEns,aMFL(:,iEns),nGrid,weight(:),rhow(:),rho(:),EceLDA(iEns))
call restricted_elda_correlation_individual_energy(aMFL(1:3,iEns),nGrid,weight(:),rhow(:),rho(:),EceLDA(iEns))
end do
! LDA-centered functional
call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA)
if(LDA_centered) EceLDA(:) = EceLDA(:) + EcLDA - EceLDA(1)
if(LDA_centered) then
call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA)
EceLDA(2) = EcLDA + wEns(2)*(EceLDA(2) - EceLDA(1))
EceLDA(1) = EcLDA
end if
! Weight-denpendent functional for ensembles

17
src/eDFT/RMFL20_lda_correlation_potential.f90
View File

@ -3,7 +3,7 @@ subroutine RMFL20_lda_correlation_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,F
! Compute Marut-Fromager-Loos weight-dependent LDA correlation potential
implicit none
include 'parameters.h'
include 'parameters.h'
! Input variables
@ -41,20 +41,21 @@ include 'parameters.h'
aMFL(2,2) = -0.0506019d0
aMFL(3,2) = +0.0331417d0
! Compute correlation energy for ground, singly-excited and doubly-excited states
! Compute correlation energy for ground- and doubly-excited states
do iEns=1,nEns
call restricted_elda_correlation_potential(aMFL(:,iEns),nGrid,weight,nBas,AO,rho,FceLDA(:,:,iEns))
call restricted_elda_correlation_potential(aMFL(1:3,iEns),nGrid,weight(:),nBas,AO(:,:),rho(:),FceLDA(:,:,iEns))
end do
! LDA-centered functional
call RVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,FcLDA)
if(LDA_centered) then
do iEns=1,nEns
FceLDA(:,:,iEns) = FceLDA(:,:,iEns) + FcLDA(:,:) - FceLDA(:,:,1)
end do
call RVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:),FcLDA(:,:))
FceLDA(:,:,2) = FcLDA(:,:) + wEns(2)*(FceLDA(:,:,2) - FceLDA(:,:,1))
FceLDA(:,:,1) = FcLDA(:,:)
end if
! Weight-denpendent functional for ensembles

2
src/eDFT/RVWN5_lda_correlation_potential.f90
View File

@ -30,7 +30,7 @@ subroutine RVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc)
! Parameters of the functional
a_p = +0.0621814D0/2D0
a_p = +0.0621814d0/2d0
x0_p = -0.10498d0
b_p = +3.72744d0
c_p = +12.9352d0

5
src/eDFT/elda_correlation_energy.f90
View File

@ -1,4 +1,4 @@
subroutine elda_correlation_energy(nEns,aLF,nGrid,weight,rho,Ec)
subroutine elda_correlation_energy(aLF,nGrid,weight,rho,Ec)
! Compute LDA correlation energy of 2-glomium for various states
@ -7,8 +7,7 @@ subroutine elda_correlation_energy(nEns,aLF,nGrid,weight,rho,Ec)
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: aLF(nEns)
double precision,intent(in) :: aLF(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid,nspin)

5
src/eDFT/elda_correlation_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine elda_correlation_individual_energy(nEns,aLF,nGrid,weight,rhow,rho,Ec)
subroutine elda_correlation_individual_energy(aLF,nGrid,weight,rhow,rho,Ec)
! Compute LDA correlation individual energy of 2-glomium for various states
@ -7,8 +7,7 @@ subroutine elda_correlation_individual_energy(nEns,aLF,nGrid,weight,rhow,rho,Ec)
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: aLF(nEns)
double precision,intent(in) :: aLF(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid,nspin)

5
src/eDFT/elda_correlation_potential.f90
View File

@ -1,4 +1,4 @@
subroutine elda_correlation_potential(nEns,aLF,nGrid,weight,nBas,AO,rho,Fc)
subroutine elda_correlation_potential(aLF,nGrid,weight,nBas,AO,rho,Fc)
! Compute LDA correlation energy of 2-glomium for various states
@ -7,8 +7,7 @@ subroutine elda_correlation_potential(nEns,aLF,nGrid,weight,nBas,AO,rho,Fc)
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: aLF(nEns)
double precision,intent(in) :: aLF(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas

6
src/eDFT/exchange_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas, &
ERI,P,FxHF,rhow,drhow,rho,drho,Ex)
ERI,P,rhow,drhow,rho,drho,Ex)
! Compute the exchange individual energy
@ -17,7 +17,6 @@ subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas, &
integer,intent(in) :: nBas
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: FxHF(nBas,nBas)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
@ -66,8 +65,7 @@ subroutine exchange_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,nBas, &
case(666)
call fock_exchange_potential(nBas,P,ERI,FxHF)
call fock_exchange_energy(nBas,P,FxHF,ExHF)
call fock_exchange_individual_energy(nBas,P,ERI,ExHF)
Ex = ExHF

41
src/eDFT/fock_exchange_individual_energy.f90 Normal file
View File

@ -0,0 +1,41 @@
subroutine fock_exchange_individual_energy(nBas,P,ERI,Ex)
! Compute the Fock exchange potential
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu,la,si
double precision,allocatable :: Fx(:,:)
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: Ex
! Compute HF exchange matrix
allocate(Fx(nBas,nBas))
Fx(:,:) = 0d0
do si=1,nBas
do la=1,nBas
do nu=1,nBas
do mu=1,nBas
Fx(mu,nu) = Fx(mu,nu) + P(la,si)*ERI(mu,la,si,nu)
enddo
enddo
enddo
enddo
Ex = -0.25d0*trace_matrix(nBas,matmul(P,Fx))
end subroutine fock_exchange_individual_energy

1
src/eDFT/restricted_elda_correlation_energy.f90
View File

@ -25,7 +25,6 @@ subroutine restricted_elda_correlation_energy(aMFL,nGrid,weight,rho,Ec)
! Compute eLDA correlation energy
Ec = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))

11
src/eDFT/restricted_elda_correlation_individual_energy.f90
View File

@ -1,4 +1,4 @@
subroutine restricted_elda_correlation_individual_energy(nEns,aLF,nGrid,weight,rhow,rho,Ec)
subroutine restricted_elda_correlation_individual_energy(aMFL,nGrid,weight,rhow,rho,Ec)
! Compute LDA correlation individual energy of 2-glomium for various states
@ -7,8 +7,7 @@ subroutine restricted_elda_correlation_individual_energy(nEns,aLF,nGrid,weight,r
! Input variables
integer,intent(in) :: nEns
double precision,intent(in) :: aLF(nEns)
double precision,intent(in) :: aMFL(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
@ -36,10 +35,10 @@ subroutine restricted_elda_correlation_individual_energy(nEns,aLF,nGrid,weight,r
if(r > threshold .and. rI > threshold) then
ec_p = aLF(1)/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
ec_p = aMFL(1)/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0))
dFcdr = aLF(2)*r**(-1d0/6d0) + 2d0*aLF(3)*r**(-1d0/3d0)
dFcdr = dFcdr/(1d0 + aLF(2)*r**(-1d0/6d0) + aLF(3)*r**(-1d0/3d0))
dFcdr = aMFL(2)*r**(-1d0/6d0) + 2d0*aMFL(3)*r**(-1d0/3d0)
dFcdr = dFcdr/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0))
dFcdr = ec_p*dFcdr/(6d0*r)
Ec = Ec + weight(iG)*(ec_p*rI + dFcdr*r*rI - dFcdr*r*r)

1
src/eDFT/restricted_elda_correlation_potential.f90
View File

@ -37,6 +37,7 @@ subroutine restricted_elda_correlation_potential(aMFL,nGrid,weight,nBas,AO,rho,F
if(r > threshold) then
ec_p = aMFL(1)/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0))
dFcdr = aMFL(2)*r**(-1d0/6d0) + 2d0*aMFL(3)*r**(-1d0/3d0)
dFcdr = dFcdr/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0))
dFcdr = ec_p*dFcdr/(6d0*r)

15
src/eDFT/restricted_individual_energy.f90
View File

@ -1,5 +1,5 @@
subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO, &
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,Ew,EwGIC,E,Om)
subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas, &
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,EwGIC,E,Om)
! Compute individual energies as well as excitation energies
@ -15,8 +15,6 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
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)
integer,intent(in) :: nO
integer,intent(in) :: nV
@ -35,9 +33,6 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
double precision,intent(in) :: drho(ncart,nGrid,nEns)
double precision,intent(in) :: J(nBas,nBas)
double precision,intent(in) :: Fx(nBas,nBas)
double precision,intent(in) :: FxHF(nBas,nBas)
double precision,intent(in) :: Fc(nBas,nBas)
double precision :: Ew
@ -95,7 +90,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
do iEns=1,nEns
call exchange_individual_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,ERI(:,:,:,:), &
P(:,:,iEns),FxHF(:,:),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns))
P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns))
end do
@ -117,7 +112,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
! Compute auxiliary energies
!------------------------------------------------------------------------
call restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux)
call restricted_auxiliary_energy(nBas,nEns,nO,eps(:),Eaux(:))
!------------------------------------------------------------------------
! Compute derivative discontinuities
@ -160,7 +155,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGri
Omc(iEns) = Ec(iEns) - Ec(1)
Omxc(iEns) = Exc(iEns) - Exc(1)
Omaux(iEns) = Eaux(iENs) - Eaux(1)
Omaux(iEns) = Eaux(iEns) - Eaux(1)
OmxDD(iEns) = ExDD(iEns) - ExDD(1)
OmcDD(iEns) = EcDD(iEns) - EcDD(1)

2
src/eDFT/restricted_lda_correlation_individual_energy.f90
View File

@ -29,7 +29,7 @@ subroutine restricted_lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weig
call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),Ec)
! marut-Fromager-Loos weight-dependent correlation functional
! Marut-Fromager-Loos weight-dependent correlation functional
case ('RMFL20')

2
src/eDFT/restricted_lda_correlation_potential.f90
View File

@ -3,7 +3,7 @@ subroutine restricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,
! Select LDA correlation potential
implicit none
include 'parameters.h'
include 'parameters.h'
! Input variables