PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:42 +01:00
Code Releases Activity

fix bug in RMFL20

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-03-31 12:39:26 +02:00
parent 9e8eeb4522
commit af9660b004
16 changed files with 34 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90
View File

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

11
src/eDFT/RMFL20_lda_correlation_energy.f90
View File

@ -50,17 +50,10 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec)
call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),EcLDA) call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),EcLDA)
if(LDA_centered) then if(LDA_centered) EceLDA(:) = EceLDA(:) + EcLDA - EceLDA(1)
do iEns=1,nEns
EceLDA(iEns) = EceLDA(iEns) + EcLDA - EceLDA(1)
end do
end if
! Weight-denpendent functional for ensembles ! Weight-denpendent functional for ensembles
Ec = 0d0 Ec = dot_product(wEns(:),EceLDA(:))
do iEns=1,nEns
Ec = Ec + wEns(iEns)*EceLDA(iEns)
end do
end subroutine RMFL20_lda_correlation_energy end subroutine RMFL20_lda_correlation_energy

11
src/eDFT/RMFL20_lda_correlation_individual_energy.f90
View File

@ -50,17 +50,10 @@ subroutine RMFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,
call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA) call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA)
if(LDA_centered) then if(LDA_centered) EceLDA(:) = EceLDA(:) + EcLDA - EceLDA(1)
do iEns=1,nEns
EceLDA(iEns) = EceLDA(iEns) + EcLDA - EceLDA(1)
end do
end if
! Weight-denpendent functional for ensembles ! Weight-denpendent functional for ensembles
Ec = 0d0 Ec = dot_product(wEns(:),EceLDA(:))
do iEns=1,nEns
Ec = Ec + wEns(iEns)*EceLDA(iEns)
enddo
end subroutine RMFL20_lda_correlation_individual_energy end subroutine RMFL20_lda_correlation_individual_energy

2
src/eDFT/RMFL20_lda_correlation_potential.f90
View File

@ -44,7 +44,7 @@ include 'parameters.h'
! Compute correlation energy for ground, singly-excited and doubly-excited states ! Compute correlation energy for ground, singly-excited and doubly-excited states
do iEns=1,nEns do iEns=1,nEns
call restricted_elda_correlation_potential(nEns,aMFL(:,iEns),nGrid,weight,nBas,AO,rho,FceLDA(:,:,iEns)) call restricted_elda_correlation_potential(aMFL(:,iEns),nGrid,weight,nBas,AO,rho,FceLDA(:,:,iEns))
end do end do
! LDA-centered functional ! LDA-centered functional

21
src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90
View File

@ -16,7 +16,8 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r
! Local variables ! Local variables
integer :: iEns,jEns integer :: iEns,jEns
double precision :: Cx(nEns) integer :: iG
double precision :: r
double precision :: dExdw(nEns) double precision :: dExdw(nEns)
double precision,external :: Kronecker_delta double precision,external :: Kronecker_delta
@ -24,15 +25,19 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r
double precision,intent(out) :: ExDD(nEns) double precision,intent(out) :: ExDD(nEns)
! Weight-dependent Cx coefficient for RMFL20 exchange functional
Cx(1) = Cx0
Cx(2) = Cx1
! Compute correlation energy for ground- and doubly-excited states ! Compute correlation energy for ground- and doubly-excited states
do iEns=1,nEns dExdw(:) = 0d0
call restricted_elda_exchange_energy(nEns,Cx(iEns),nGrid,weight(:),rhow(:),dExdw(iEns))
do iG=1,nGrid
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)
end if
end do end do
ExDD(:) = 0d0 ExDD(:) = 0d0

6
src/eDFT/RMFL20_lda_exchange_energy.f90
View File

@ -33,16 +33,18 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex)
Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 Cxw = wEns(1)*Cx0 + wEns(2)*Cx1
end if end if
! Cxw = CxLDA + (Cx1 - Cx0)*wEns(2)*(cos(2d0*pi*wEns(2)) + 1d0)
! Compute LDA exchange energy ! Compute LDA exchange energy
Ex = 0d0 Ex = 0d0
do iG=1,nGrid do iG=1,nGrid
r = max(0d0,rho(iG)) r = max(0d0,rho(iG))
if(r > threshold) then if(r > threshold) then
Ex = Ex + weight(iG)*Cxw*r**(4d0/3d0) Ex = Ex + weight(iG)*Cxw*r**(4d0/3d0)
endif endif
enddo enddo
end subroutine RMFL20_lda_exchange_energy end subroutine RMFL20_lda_exchange_energy

2
src/eDFT/RMFL20_lda_exchange_individual_energy.f90
View File

@ -34,8 +34,6 @@ subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho
Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 Cxw = wEns(1)*Cx0 + wEns(2)*Cx1
end if end if
! Cxw = CxLDA + (Cx1 - Cx0)*wEns(2)*(cos(2d0*pi*wEns(2)) + 1d0)
! Compute LDA exchange matrix in the AO basis ! Compute LDA exchange matrix in the AO basis
Ex = 0d0 Ex = 0d0

3
src/eDFT/RMFL20_lda_exchange_potential.f90
View File

@ -34,11 +34,10 @@ subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx)
Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 Cxw = wEns(1)*Cx0 + wEns(2)*Cx1
end if end if
! Cxw = CxLDA + (Cx1 - Cx0)*wEns(2)*(cos(2d0*pi*wEns(2)) + 1d0)
! Compute LDA exchange matrix in the AO basis ! Compute LDA exchange matrix in the AO basis
Fx(:,:) = 0d0 Fx(:,:) = 0d0
do mu=1,nBas do mu=1,nBas
do nu=1,nBas do nu=1,nBas
do iG=1,nGrid do iG=1,nGrid

12
src/eDFT/allocate_grid.f90
View File

@ -1,5 +1,5 @@
subroutine allocate_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent, & subroutine allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent, &
radial_precision,nRad,nAng,nGrid) radial_precision,nAng,nGrid)
! Allocate quadrature grid with numgrid (Radovan Bast) ! Allocate quadrature grid with numgrid (Radovan Bast)
@ -13,18 +13,12 @@ subroutine allocate_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_m
integer,intent(in) :: nNuc integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell)
double precision,intent(in) :: ExpShell(maxShell,maxK)
integer,intent(in) :: max_ang_mom(nNuc) integer,intent(in) :: max_ang_mom(nNuc)
double precision,intent(in) :: min_exponent(nNuc,maxL+1) double precision,intent(in) :: min_exponent(nNuc,maxL+1)
double precision,intent(in) :: max_exponent(nNuc) double precision,intent(in) :: max_exponent(nNuc)
double precision :: radial_precision double precision :: radial_precision
integer,intent(in) :: nRad
integer,intent(in) :: nAng integer,intent(in) :: nAng
! Local variables ! Local variables
@ -33,9 +27,7 @@ subroutine allocate_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_m
integer :: min_num_angular_points integer :: min_num_angular_points
integer :: max_num_angular_points integer :: max_num_angular_points
integer :: num_points
integer :: center_index
type(c_ptr) :: context type(c_ptr) :: context
! Output variables ! Output variables

6
src/eDFT/build_grid.f90
View File

@ -1,4 +1,4 @@
subroutine build_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent, & subroutine build_grid(nNuc,ZNuc,rNuc,max_ang_mom,min_exponent,max_exponent, &
radial_precision,nRad,nAng,nGrid,weight,root) radial_precision,nRad,nAng,nGrid,weight,root)
! Compute quadrature grid with numgrid (Radovan Bast) ! Compute quadrature grid with numgrid (Radovan Bast)
@ -15,10 +15,6 @@ subroutine build_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,
double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: rNuc(nNuc,ncart)
integer,intent(in) :: nShell
integer,intent(in) :: TotAngMomShell(maxShell)
double precision,intent(in) :: ExpShell(maxShell,maxK)
integer,intent(in) :: max_ang_mom(nNuc) integer,intent(in) :: max_ang_mom(nNuc)
double precision,intent(in) :: min_exponent(nNuc,maxL+1) double precision,intent(in) :: min_exponent(nNuc,maxL+1)
double precision,intent(in) :: max_exponent(nNuc) double precision,intent(in) :: max_exponent(nNuc)

5
src/eDFT/eDFT.f90
View File

@ -136,14 +136,13 @@ program eDFT
call read_grid(SGn,radial_precision,nRad,nAng,nGrid) call read_grid(SGn,radial_precision,nRad,nAng,nGrid)
! nGrid = nRad*nAng ! nGrid = nRad*nAng
call allocate_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent, & call allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent,radial_precision,nAng,nGrid)
radial_precision,nRad,nAng,nGrid)
allocate(root(ncart,nGrid),weight(nGrid)) allocate(root(ncart,nGrid),weight(nGrid))
! call quadrature_grid(nRad,nAng,nGrid,root,weight) ! call quadrature_grid(nRad,nAng,nGrid,root,weight)
call build_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent, & call build_grid(nNuc,ZNuc,rNuc,max_ang_mom,min_exponent,max_exponent, &
radial_precision,nRad,nAng,nGrid,weight,root) radial_precision,nRad,nAng,nGrid,weight,root)
!------------------------------------------------------------------------ !------------------------------------------------------------------------

4
src/eDFT/elda_correlation_individual_energy.f90
View File

@ -31,8 +31,8 @@ subroutine elda_correlation_individual_energy(nEns,aLF,nGrid,weight,rhow,rho,Ec)
do iG=1,nGrid do iG=1,nGrid
ra = max(0d0,rho(iG,1)) ra = max(0d0,rhow(iG,1))
rb = max(0d0,rho(iG,2)) rb = max(0d0,rhow(iG,2))
raI = max(0d0,rho(iG,1)) raI = max(0d0,rho(iG,1))
rbI = max(0d0,rho(iG,2)) rbI = max(0d0,rho(iG,2))

3
src/eDFT/restricted_elda_correlation_energy.f90
View File

@ -31,11 +31,8 @@ subroutine restricted_elda_correlation_energy(aMFL,nGrid,weight,rho,Ec)
r = max(0d0,rho(iG)) r = max(0d0,rho(iG))
if(r > threshold) then if(r > threshold) then
e = aMFL(1)/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0)) e = aMFL(1)/(1d0 + aMFL(2)*r**(-1d0/6d0) + aMFL(3)*r**(-1d0/3d0))
Ec = Ec + weight(iG)*e*r Ec = Ec + weight(iG)*e*r
end if end if
end do end do

2
src/eDFT/restricted_elda_correlation_individual_energy.f90
View File

@ -31,7 +31,7 @@ subroutine restricted_elda_correlation_individual_energy(nEns,aLF,nGrid,weight,r
do iG=1,nGrid do iG=1,nGrid
r = max(0d0,rho(iG)) r = max(0d0,rhow(iG))
rI = max(0d0,rho(iG)) rI = max(0d0,rho(iG))
if(r > threshold .and. rI > threshold) then if(r > threshold .and. rI > threshold) then

3
src/eDFT/restricted_elda_correlation_potential.f90
View File

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

42
src/eDFT/restricted_elda_exchange_energy.f90
View File

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