4
1
mirror of https://github.com/pfloos/quack synced 2024-06-20 12:12:15 +02:00

removed restricted routines

This commit is contained in:
Pierre-Francois Loos 2021-11-30 09:30:09 +01:00
parent 9078568db5
commit e3b7f70025
29 changed files with 0 additions and 1642 deletions

View File

@ -1,56 +0,0 @@
subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux,occnum)
! Compute the auxiliary KS energies
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nEns
integer,intent(in) :: nO
double precision,intent(in) :: eps(nBas)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
! Local variables
integer :: iEns
! Output variables
double precision,intent(out) :: Eaux(nEns)
! Ground state density matrix
iEns = 1
Eaux(iEns) = 2d0*sum(eps(1:nO))
! Singly-excited state density matrix
iEns = 2
if(nO > 1) then
Eaux(iEns) = 2d0*sum(eps(1:nO-1))
else
Eaux(iEns) = 0d0
end if
Eaux(iEns) = Eaux(iEns) + sum(occnum(:,1,iEns))*eps(nO) + sum(occnum(:,2,iEns))*eps(nO+1)
! Doubly-excited state density matrix
iEns = 3
if(nO > 1) then
Eaux(iEns) = 2d0*sum(eps(1:nO-1))
else
Eaux(iEns) = 0d0
end if
Eaux(iEns) = Eaux(iEns) + sum(occnum(:,1,iEns))*eps(nO) + sum(occnum(:,2,iEns))*eps(nO+1)
end subroutine restricted_auxiliary_energy

View File

@ -1,65 +0,0 @@
subroutine restricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec)
! Compute the correlation part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
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) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
! Local variables
double precision :: aC
! Output variables
double precision,intent(out) :: Ec(nEns)
select case (rung)
! Hartree calculation
case(0)
Ec(:) = 0d0
! LDA functionals
case(1)
call restricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),Ec(:))
! GGA functionals
case(2)
call print_warning('!!! derivative discontinuity NYI for GGAs !!!')
stop
! Hybrid functionals
case(4)
call print_warning('!!! derivative discontinuity NYI for hybrids !!!')
stop
aC = 0.81d0
! Hartree-Fock calculation
case(666)
Ec(:) = 0d0
end select
end subroutine restricted_correlation_derivative_discontinuity

View File

@ -1,65 +0,0 @@
subroutine restricted_correlation_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,drho,Ec)
! Compute the correlation energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
logical,intent(in) :: LDA_centered
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)
double precision,intent(in) :: drho(ncart,nGrid)
! Local variables
double precision :: EcLDA
double precision :: EcGGA
double precision :: aC
! Output variables
double precision,intent(out) :: Ec
select case (rung)
! Hartree calculation
case(0)
Ec = 0d0
! LDA functionals
case(1)
call restricted_lda_correlation_energy(DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),rho(:),Ec)
! GGA functionals
case(2)
call print_warning('!!! restricted correlation energies NYI for GGAs !!!')
stop
case(4)
call print_warning('!!! restricted correlation energies NYI for Hybrids !!!')
stop
! Hartree-Fock calculation
case(666)
Ec = 0d0
end select
end subroutine restricted_correlation_energy

View File

@ -1,70 +0,0 @@
subroutine restricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec)
! Compute the correlation energy of individual states
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
logical,intent(in) :: LDA_centered
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) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Local variables
double precision :: EcLDA
double precision :: EcGGA
double precision :: aC
! Output variables
double precision,intent(out) :: Ec
select case (rung)
! Hartree calculation
case(0)
Ec = 0d0
! LDA functionals
case(1)
call restricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),rhow(:),rho(:),Ec)
! GGA functionals
case(2)
call print_warning('!!! Individual energies NYI for GGAs !!!')
stop
! Hybrid functionals
case(4)
call print_warning('!!! Individual energies NYI for hybrids !!!')
stop
aC = 0.81d0
! Hartree-Fock calculation
case(666)
Ec = 0d0
end select
end subroutine restricted_correlation_individual_energy

View File

@ -1,69 +0,0 @@
subroutine restricted_correlation_potential(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc)
! Compute the correlation potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
logical,intent(in) :: LDA_centered
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) :: dAO(ncart,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Local variables
double precision,allocatable :: FcLDA(:,:)
! Output variables
double precision,intent(out) :: Fc(nBas,nBas)
! Memory allocation
select case (rung)
! Hartree calculation
case(0)
Fc(:,:) = 0d0
! LDA functionals
case(1)
call restricted_lda_correlation_potential(DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:))
! GGA functionals
case(2)
call print_warning('!!! restricted correlation potentials NYI for GGAs !!!')
stop
! Hybrid functionals
case(4)
call print_warning('!!! restricted correlation potentials NYI for Hybrids !!!')
stop
! Hartree-Fock calculation
case(666)
Fc(:,:) = 0d0
end select
end subroutine restricted_correlation_potential

View File

@ -1,56 +0,0 @@
subroutine restricted_density_matrix(nBas,nEns,nO,c,P,occnum)
! Calculate density matrices
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nEns
integer,intent(in) :: nO
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: occnum(nBas,nspin,nEns)
! Local variables
integer :: iEns
! Output variables
double precision,intent(out) :: P(nBas,nBas,nEns)
! Ground state density matrix
iEns = 1
P(:,:,iEns) = 2.d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Doubly-excited state density matrix
iEns = 2
if(nO > 1) then
P(:,:,iEns) = 2d0*matmul(c(:,1:nO-1),transpose(c(:,1:nO-1)))
else
P(:,:,iEns) = 0d0
end if
P(:,:,iEns) = P(:,:,iEns) + sum(occnum(:,1,iEns))*matmul(c(:,nO:nO),transpose(c(:,nO:nO))) &
+ sum(occnum(:,2,iEns))*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
! Doubly-excited state density matrix
iEns = 3
if(nO > 1) then
P(:,:,iEns) = 2d0*matmul(c(:,1:nO-1),transpose(c(:,1:nO-1)))
else
P(:,:,iEns) = 0d0
end if
P(:,:,iEns) = P(:,:,iEns) + sum(occnum(:,1,iEns))*matmul(c(:,nO:nO),transpose(c(:,nO:nO))) &
+ sum(occnum(:,2,iEns))*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
end subroutine restricted_density_matrix

View File

@ -1,39 +0,0 @@
subroutine restricted_elda_correlation_energy(aMFL,nGrid,weight,rho,Ec)
! Compute the restricted LDA correlation energy of 2-glomium for various states
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: aMFL(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: iG
double precision :: r,e
! Output variables
double precision,intent(out) :: Ec
! Compute eLDA correlation energy
Ec = 0d0
do iG=1,nGrid
r = max(0d0,rho(iG))
if(r > threshold) then
e = aMFL(1)/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0))
Ec = Ec + weight(iG)*e*r
end if
end do
end subroutine restricted_elda_correlation_energy

View File

@ -1,50 +0,0 @@
subroutine restricted_elda_correlation_individual_energy(aMFL,nGrid,weight,rhow,rho,Ec)
! Compute LDA correlation individual energy of 2-glomium for various states
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: aMFL(3)
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 :: r
double precision :: rI
double precision :: ec_p,dFcdr
! Output variables
double precision,intent(out) :: Ec
! Compute eLDA correlation potential
Ec = 0d0
do iG=1,nGrid
r = max(0d0,rhow(iG))
rI = max(0d0,rho(iG))
if(r > threshold .or. rI > 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)
Ec = Ec + weight(iG)*(ec_p*rI + dFcdr*r*rI - dFcdr*r*r)
end if
end do
end subroutine restricted_elda_correlation_individual_energy

View File

@ -1,54 +0,0 @@
subroutine restricted_elda_correlation_potential(aMFL,nGrid,weight,nBas,AO,rho,Fc)
! Compute LDA correlation energy of 2-glomium for various states
implicit none
include 'parameters.h'
! Input variables
double precision,intent(in) :: aMFL(3)
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,ec_p
double precision :: dFcdr
! Output variables
double precision,intent(out) :: Fc(nBas,nBas)
! Compute eLDA correlation potential
Fc(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = max(0d0,rho(iG))
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)
dFcdr = ec_p + dFcdr*r
Fc(mu,nu) = Fc(mu,nu) + weight(iG)*AO(mu,iG)*AO(nu,iG)*dFcdr
endif
end do
end do
end do
end subroutine restricted_elda_correlation_potential

View File

@ -1,65 +0,0 @@
subroutine restricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD,Cx_choice)
! Compute the exchange part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
! Output variables
double precision,intent(out) :: ExDD(nEns)
select case (rung)
! Hartree calculation
case(0)
ExDD(:) = 0d0
! LDA functionals
case(1)
call restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),&
rhow(:),ExDD(:),Cx_choice)
! GGA functionals
case(2)
call restricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:))
! Hybrid functionals
case(4)
call print_warning('!!! exchange part of derivative discontinuity NYI for hybrids !!!')
stop
! Hartree-Fock calculation
case(666)
ExDD(:) = 0d0
end select
end subroutine restricted_exchange_derivative_discontinuity

View File

@ -1,85 +0,0 @@
subroutine restricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF,rho,drho,Ex,Cx_choice)
! Compute the exchange energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: FxHF(nBas,nBas)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
double precision :: ExLDA,ExGGA,ExHF
double precision :: cX,aX,aC
! Output variables
double precision,intent(out) :: Ex
select case (rung)
! Hartree calculation
case(0)
Ex = 0d0
! LDA functionals
case(1)
call restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,ExLDA,Cx_choice)
Ex = ExLDA
! GGA functionals
case(2)
call restricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
Ex = ExGGA
! Hybrid functionals
case(4)
cX = 0.20d0
aX = 0.72d0
aC = 0.81d0
call restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,ExLDA,Cx_choice)
call restricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA)
call restricted_fock_exchange_energy(nBas,P,FxHF,ExHF)
Ex = ExLDA &
+ cX*(ExHF - ExLDA) &
+ aX*(ExGGA - ExLDA)
! Hartree-Fock calculation
case(666)
call restricted_fock_exchange_energy(nBas,P,FxHF,ExHF)
Ex = ExHF
end select
end subroutine restricted_exchange_energy

View File

@ -1,81 +0,0 @@
subroutine restricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
ERI,Pw,P,rhow,drhow,rho,drho,Ex,Cx_choice)
! Compute the exchange individual energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Pw(nBas,nBas)
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
integer,intent(in) :: Cx_choice
! Local variables
double precision :: ExLDA
double precision :: ExGGA
double precision :: ExHF
! Output variables
double precision,intent(out) :: Ex
select case (rung)
! Hartree calculation
case(0)
Ex = 0d0
! LDA functionals
case(1)
call restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,ExLDA,Cx_choice)
Ex = ExLDA
! GGA functionals
case(2)
call restricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,ExGGA)
Ex = ExGGA
! Hybrid functionals
case(4)
call print_warning('!!! Individual energies NYI for Hybrids !!!')
stop
! Hartree-Fock calculation
case(666)
call restricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,ExHF)
Ex = ExHF
end select
end subroutine restricted_exchange_individual_energy

View File

@ -1,87 +0,0 @@
subroutine restricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P, &
ERI,AO,dAO,rho,drho,Fx,FxHF,Cx_choice)
! Compute the exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
character(len=12),intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,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)
integer,intent(in) :: Cx_choice
! Local variables
double precision,allocatable :: FxLDA(:,:),FxGGA(:,:)
double precision :: cX,aX
! Output variables
double precision,intent(out) :: Fx(nBas,nBas),FxHF(nBas,nBas)
! Memory allocation
select case (rung)
! Hartree calculation
case(0)
Fx(:,:) = 0d0
! LDA functionals
case(1)
call restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
! GGA functionals
case(2)
call restricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Hybrid functionals
case(4)
allocate(FxLDA(nBas,nBas),FxGGA(nBas,nBas))
cX = 0.20d0
aX = 0.72d0
call restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,FxLDA,Cx_choice)
call restricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
call restricted_fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = FxLDA(:,:) &
+ cX*(FxHF(:,:) - FxLDA(:,:)) &
+ aX*(FxGGA(:,:) - FxLDA(:,:))
! Hartree-Fock calculation
case(666)
call restricted_fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = FxHF(:,:)
end select
end subroutine restricted_exchange_potential

View File

@ -1,25 +0,0 @@
subroutine restricted_fock_exchange_energy(nBas,P,Fx,Ex)
! Compute the (exact) Fock exchange energy
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: Fx(nBas,nBas)
! Local variables
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: Ex
! Compute HF exchange energy
Ex = 0.5d0*trace_matrix(nBas,matmul(P,Fx))
end subroutine restricted_fock_exchange_energy

View File

@ -1,31 +0,0 @@
subroutine restricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex)
! Compute the Fock exchange potential
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: Pw(nBas,nBas)
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
double precision,allocatable :: Fx(:,:)
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: Ex
! Compute HF exchange matrix
allocate(Fx(nBas,nBas))
call restricted_fock_exchange_potential(nBas,Pw(:,:),ERI(:,:,:,:),Fx(:,:))
Ex = trace_matrix(nBas,matmul(P(:,:),Fx(:,:))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),Fx(:,:)))
end subroutine restricted_fock_exchange_individual_energy

View File

@ -1,36 +0,0 @@
subroutine restricted_fock_exchange_potential(nBas,P,ERI,Fx)
! 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
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Compute HF exchange matrix
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
Fx(:,:) = 0.5d0*Fx(:,:)
end subroutine restricted_fock_exchange_potential

View File

@ -1,40 +0,0 @@
subroutine restricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD)
! Compute the exchange GGA part of the derivative discontinuity
implicit none
include 'parameters.h'
! 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) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
! Local variables
! Output variables
double precision,intent(out) :: ExDD(nEns)
! Select correlation functional
select case (DFA)
case ('B88')
ExDD(:) = 0d0
case default
call print_warning('!!! GGA exchange derivative discontinuity not available !!!')
stop
end select
end subroutine restricted_gga_exchange_derivative_discontinuity

View File

@ -1,36 +0,0 @@
subroutine restricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex)
! Select GGA exchange functional for energy calculation
implicit none
include 'parameters.h'
! 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)
double precision,intent(in) :: drho(3,nGrid)
! Output variables
double precision :: Ex
select case (DFA)
case ('B88')
call RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex)
case default
call print_warning('!!! GGA exchange energy not available !!!')
stop
end select
end subroutine restricted_gga_exchange_energy

View File

@ -1,39 +0,0 @@
subroutine restricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ex)
! Compute GGA exchange energy for individual states
implicit none
include 'parameters.h'
! 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) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(ncart,nGrid)
! Output variables
double precision :: Ex
! Select correlation functional
select case (DFA)
case ('B88')
call RB88_gga_exchange_individual_energy(nGrid,weight(:),rhow(:),drhow(:,:),rho(:),drho(:,:),Ex)
case default
call print_warning('!!! GGA exchange individual energy not available !!!')
stop
end select
end subroutine restricted_gga_exchange_individual_energy

View File

@ -1,40 +0,0 @@
subroutine restricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Select GGA exchange functional for potential calculation
implicit none
include 'parameters.h'
! 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
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Select GGA exchange functional
select case (DFA)
case ('B88')
call RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
case default
call print_warning('!!! GGA exchange potential not available !!!')
stop
end select
end subroutine restricted_gga_exchange_potential

View File

@ -1,166 +0,0 @@
subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, &
nO,nV,T,V,ERI,ENuc,eps,Pw,rhow,drhow,J,P,rho,drho,Ew,E,Om,occnum,Cx_choice)
! Compute individual energies as well as excitation energies
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: x_rung,c_rung
character(len=12),intent(in) :: x_DFA,c_DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
integer,intent(in) :: nO
integer,intent(in) :: nV
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
double precision,intent(in) :: eps(nBas)
double precision,intent(in) :: Pw(nBas,nBas)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: drhow(ncart,nGrid)
double precision,intent(in) :: P(nBas,nBas,nEns)
double precision,intent(in) :: rho(nGrid,nEns)
double precision,intent(in) :: drho(ncart,nGrid,nEns)
double precision,intent(in) :: J(nBas,nBas)
double precision :: Ew
double precision,intent(in) :: occnum(nBas,nspin,nEns)
integer,intent(in) :: Cx_choice
! Local variables
double precision :: ET(nEns)
double precision :: EV(nEns)
double precision :: EJ(nEns)
double precision :: Ex(nEns), Ec(nEns), Exc(nEns)
double precision :: Eaux(nEns)
double precision :: ExDD(nEns), EcDD(nEns), ExcDD(nEns)
double precision :: Omx(nEns), Omc(nEns), Omxc(nEns)
double precision :: Omaux(nEns)
double precision :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns)
double precision,external :: trace_matrix
integer :: iEns
! Output variables
double precision,intent(out) :: E(nEns)
double precision,intent(out) :: Om(nEns)
!------------------------------------------------------------------------
! Kinetic energy
!------------------------------------------------------------------------
do iEns=1,nEns
ET(iEns) = trace_matrix(nBas,matmul(P(:,:,iEns),T(:,:)))
end do
!------------------------------------------------------------------------
! Potential energy
!------------------------------------------------------------------------
do iEns=1,nEns
EV(iEns) = trace_matrix(nBas,matmul(P(:,:,iEns),V(:,:)))
end do
!------------------------------------------------------------------------
! Individual Hartree energy
!------------------------------------------------------------------------
do iEns=1,nEns
call hartree_coulomb(nBas,Pw(:,:),ERI(:,:,:,:),J(:,:))
EJ(iEns) = trace_matrix(nBas,matmul(P(:,:,iEns),J(:,:))) &
- 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),J(:,:)))
end do
!------------------------------------------------------------------------
! Individual exchange energy
!------------------------------------------------------------------------
do iEns=1,nEns
call restricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, &
ERI(:,:,:,:),Pw(:,:),P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns)&
,Ex(iEns),Cx_choice)
end do
!------------------------------------------------------------------------
! Indivudual correlation energy
!------------------------------------------------------------------------
do iEns=1,nEns
call restricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:), &
rho(:,iEns),drho(:,:,iEns),Ec(iEns))
end do
!------------------------------------------------------------------------
! Compute auxiliary energies
!------------------------------------------------------------------------
call restricted_auxiliary_energy(nBas,nEns,nO,eps(:),Eaux(:),occnum)
!------------------------------------------------------------------------
! Compute derivative discontinuities
!------------------------------------------------------------------------
call restricted_exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:), &
rhow(:),drhow(:,:),ExDD(:),Cx_choice)
call restricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),EcDD(:))
ExcDD(:) = ExDD(:) + EcDD(:)
!------------------------------------------------------------------------
! Total energy
!------------------------------------------------------------------------
do iEns=1,nEns
Exc(iEns) = Ex(iEns) + Ec(iEns)
E(iEns) = ET(iEns) + EV(iEns) + EJ(iEns) &
+ Ex(iEns) + Ec(iEns) + ExcDD(iEns)
end do
!------------------------------------------------------------------------
! Excitation energies
!------------------------------------------------------------------------
do iEns=1,nEns
Om(iEns) = E(iEns) - E(1)
Omx(iEns) = Ex(iEns) - Ex(1)
Omc(iEns) = Ec(iEns) - Ec(1)
Omxc(iEns) = Exc(iEns) - Exc(1)
Omaux(iEns) = Eaux(iEns) - Eaux(1)
OmxDD(iEns) = ExDD(iEns) - ExDD(1)
OmcDD(iEns) = EcDD(iEns) - EcDD(1)
OmxcDD(iEns) = ExcDD(iEns) - ExcDD(1)
end do
!------------------------------------------------------------------------
! Dump results
!------------------------------------------------------------------------
call print_restricted_individual_energy(nEns,ENuc,Ew,ET(:),EV(:),EJ(:),Ex(:),Ec(:),Exc(:), &
Eaux(:),ExDD(:),EcDD(:),ExcDD(:),E(:), &
Om(:),Omx(:),Omc(:),Omxc(:),Omaux,OmxDD(:),OmcDD(:),OmxcDD(:))
end subroutine restricted_individual_energy

View File

@ -1,46 +0,0 @@
subroutine restricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec)
! Compute the correlation LDA part of the derivative discontinuity
implicit none
include 'parameters.h'
! 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) :: rhow(nGrid)
! Local variables
! Output variables
double precision,intent(out) :: Ec(nEns)
! Select correlation functional
select case (DFA)
case ('W38')
Ec(:) = 0d0
case ('VWN5')
Ec(:) = 0d0
case ('MFL20')
call RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),Ec(:))
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine restricted_lda_correlation_derivative_discontinuity

View File

@ -1,47 +0,0 @@
subroutine restricted_lda_correlation_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ec)
! Select LDA correlation functional
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: LDA_centered
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)
! Output variables
double precision,intent(out) :: Ec
! Select correlation functional
select case (DFA)
! Hartree-Fock
case ('HF')
Ec = 0d0
case ('VWN5')
call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),Ec)
case ('MFL20')
call RMFL20_lda_correlation_energy(LDA_centered,nEns,wEns(:),nGrid,weight(:),rho(:),Ec)
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine restricted_lda_correlation_energy

View File

@ -1,46 +0,0 @@
subroutine restricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ec)
! Compute LDA correlation energy for individual states
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: LDA_centered
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) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
! Output variables
double precision :: Ec
! Select correlation functional
select case (DFA)
! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200
case ('VWN5')
call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),Ec)
! Marut-Fromager-Loos weight-dependent correlation functional
case ('MFL20')
call RMFL20_lda_correlation_individual_energy(LDA_centered,nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ec)
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine restricted_lda_correlation_individual_energy

View File

@ -1,49 +0,0 @@
subroutine restricted_lda_correlation_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc)
! Select LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: LDA_centered
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
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
! Output variables
double precision,intent(out) :: Fc(nBas,nBas)
! Select correlation functional
select case (DFA)
! Hartree-Fock
case ('HF')
Fc(:,:) = 0d0
case ('VWN5')
call RVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:))
case ('MFL20')
call RMFL20_lda_correlation_potential(LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:),Fc(:,:))
case default
call print_warning('!!! LDA correlation functional not available !!!')
stop
end select
end subroutine restricted_lda_correlation_potential

View File

@ -1,51 +0,0 @@
subroutine restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice)
! Compute the exchange LDA part of the derivative discontinuity
implicit none
include 'parameters.h'
! Input variables
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
integer,intent(in) :: Cx_choice
! Local variables
! Output variables
double precision,intent(out) :: ExDD(nEns)
! Select correlation functional
select case (DFA)
case ('S51')
ExDD(:) = 0d0
case ('MFL20')
call RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:))
case ('CC')
call RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:),Cx_choice)
case default
call print_warning('!!! LDA exchange derivative discontinuity not available !!!')
stop
end select
end subroutine restricted_lda_exchange_derivative_discontinuity

View File

@ -1,48 +0,0 @@
subroutine restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice)
! Select LDA exchange functional
implicit none
include 'parameters.h'
! Input variables
character(len=12),intent(in) :: DFA
logical,intent(in) :: LDA_centered
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Output variables
double precision,intent(out) :: Ex
! Select correlation functional
select case (DFA)
case ('S51')
call RS51_lda_exchange_energy(nGrid,weight,rho,Ex)
case ('MFL20')
call RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex)
case ('CC')
call RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice)
case default
call print_warning('!!! LDA restricted exchange functional not available !!!')
stop
end select
end subroutine restricted_lda_exchange_energy

View File

@ -1,49 +0,0 @@
subroutine restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
! Compute LDA exchange energy for individual states
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: LDA_centered
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rhow(nGrid)
double precision,intent(in) :: rho(nGrid)
integer,intent(in) :: Cx_choice
! Output variables
double precision :: Ex
! Select correlation functional
select case (DFA)
case ('S51')
call RS51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex)
case ('MFL20')
call RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ex)
case ('CC')
call RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice)
case default
call print_warning('!!! LDA exchange individual energy not available !!!')
stop
end select
end subroutine restricted_lda_exchange_individual_energy

View File

@ -1,51 +0,0 @@
subroutine restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
! Select LDA correlation potential
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: LDA_centered
character(len=12),intent(in) :: DFA
integer,intent(in) :: nEns
double precision,intent(in) :: wEns(nEns)
double precision,intent(in) :: aCC_w1(3)
double precision,intent(in) :: aCC_w2(3)
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)
integer,intent(in) :: Cx_choice
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Select exchange functional
select case (DFA)
case ('S51')
call RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
case ('MFL20')
call RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
case ('CC')
call RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice)
case default
call print_warning('!!! LDA exchange functional not available !!!')
stop
end select
end subroutine restricted_lda_exchange_potential