diff --git a/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 index d8eab4c..455f29f 100644 --- a/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 @@ -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 do iEns=1,nEns - call restricted_elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rhow(:),dEcdw(iEns)) - end do EcDD(:) = 0d0 diff --git a/src/eDFT/RMFL20_lda_correlation_energy.f90 b/src/eDFT/RMFL20_lda_correlation_energy.f90 index 534ca8c..53791c4 100644 --- a/src/eDFT/RMFL20_lda_correlation_energy.f90 +++ b/src/eDFT/RMFL20_lda_correlation_energy.f90 @@ -50,17 +50,10 @@ subroutine RMFL20_lda_correlation_energy(nEns,wEns,nGrid,weight,rho,Ec) call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),EcLDA) - if(LDA_centered) then - do iEns=1,nEns - EceLDA(iEns) = EceLDA(iEns) + EcLDA - EceLDA(1) - end do - end if + if(LDA_centered) EceLDA(:) = EceLDA(:) + EcLDA - EceLDA(1) ! Weight-denpendent functional for ensembles - Ec = 0d0 - do iEns=1,nEns - Ec = Ec + wEns(iEns)*EceLDA(iEns) - end do + Ec = dot_product(wEns(:),EceLDA(:)) end subroutine RMFL20_lda_correlation_energy diff --git a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 index 9b0d37d..ea93e68 100644 --- a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 @@ -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) - if(LDA_centered) then - do iEns=1,nEns - EceLDA(iEns) = EceLDA(iEns) + EcLDA - EceLDA(1) - end do - end if + if(LDA_centered) EceLDA(:) = EceLDA(:) + EcLDA - EceLDA(1) ! Weight-denpendent functional for ensembles - Ec = 0d0 - do iEns=1,nEns - Ec = Ec + wEns(iEns)*EceLDA(iEns) - enddo + Ec = dot_product(wEns(:),EceLDA(:)) end subroutine RMFL20_lda_correlation_individual_energy diff --git a/src/eDFT/RMFL20_lda_correlation_potential.f90 b/src/eDFT/RMFL20_lda_correlation_potential.f90 index efc53b3..24f7a30 100644 --- a/src/eDFT/RMFL20_lda_correlation_potential.f90 +++ b/src/eDFT/RMFL20_lda_correlation_potential.f90 @@ -44,7 +44,7 @@ include 'parameters.h' ! Compute correlation energy for ground, singly-excited and doubly-excited states 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 ! LDA-centered functional @@ -55,7 +55,7 @@ include 'parameters.h' do iEns=1,nEns FceLDA(:,:,iEns) = FceLDA(:,:,iEns) + FcLDA(:,:) - FceLDA(:,:,1) end do - end if + end if ! Weight-denpendent functional for ensembles diff --git a/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 index 6c57c4a..3a8414c 100644 --- a/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 @@ -16,7 +16,8 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r ! Local variables integer :: iEns,jEns - double precision :: Cx(nEns) + integer :: iG + double precision :: r double precision :: dExdw(nEns) double precision,external :: Kronecker_delta @@ -24,16 +25,20 @@ subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,r 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 - do iEns=1,nEns - call restricted_elda_exchange_energy(nEns,Cx(iEns),nGrid,weight(:),rhow(:),dExdw(iEns)) - end do + dExdw(:) = 0d0 + + 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 ExDD(:) = 0d0 diff --git a/src/eDFT/RMFL20_lda_exchange_energy.f90 b/src/eDFT/RMFL20_lda_exchange_energy.f90 index 092f8a9..555c257 100644 --- a/src/eDFT/RMFL20_lda_exchange_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_energy.f90 @@ -33,16 +33,18 @@ subroutine RMFL20_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 end if -! Cxw = CxLDA + (Cx1 - Cx0)*wEns(2)*(cos(2d0*pi*wEns(2)) + 1d0) - ! Compute LDA exchange energy Ex = 0d0 + do iG=1,nGrid + r = max(0d0,rho(iG)) + if(r > threshold) then Ex = Ex + weight(iG)*Cxw*r**(4d0/3d0) endif + enddo end subroutine RMFL20_lda_exchange_energy diff --git a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 index f35a235..4feab64 100644 --- a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 @@ -34,8 +34,6 @@ subroutine RMFL20_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 end if -! Cxw = CxLDA + (Cx1 - Cx0)*wEns(2)*(cos(2d0*pi*wEns(2)) + 1d0) - ! Compute LDA exchange matrix in the AO basis Ex = 0d0 diff --git a/src/eDFT/RMFL20_lda_exchange_potential.f90 b/src/eDFT/RMFL20_lda_exchange_potential.f90 index 040f217..989573b 100644 --- a/src/eDFT/RMFL20_lda_exchange_potential.f90 +++ b/src/eDFT/RMFL20_lda_exchange_potential.f90 @@ -34,11 +34,10 @@ subroutine RMFL20_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 end if -! Cxw = CxLDA + (Cx1 - Cx0)*wEns(2)*(cos(2d0*pi*wEns(2)) + 1d0) - ! Compute LDA exchange matrix in the AO basis Fx(:,:) = 0d0 + do mu=1,nBas do nu=1,nBas do iG=1,nGrid diff --git a/src/eDFT/allocate_grid.f90 b/src/eDFT/allocate_grid.f90 index 8dd7330..22ab6cf 100644 --- a/src/eDFT/allocate_grid.f90 +++ b/src/eDFT/allocate_grid.f90 @@ -1,5 +1,5 @@ -subroutine allocate_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent, & - radial_precision,nRad,nAng,nGrid) +subroutine allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent, & + radial_precision,nAng,nGrid) ! 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 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) double precision,intent(in) :: min_exponent(nNuc,maxL+1) double precision,intent(in) :: max_exponent(nNuc) double precision :: radial_precision - integer,intent(in) :: nRad integer,intent(in) :: nAng ! Local variables @@ -33,9 +27,7 @@ subroutine allocate_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_m integer :: min_num_angular_points integer :: max_num_angular_points - integer :: num_points - integer :: center_index type(c_ptr) :: context ! Output variables diff --git a/src/eDFT/build_grid.f90 b/src/eDFT/build_grid.f90 index 44622c6..c255231 100644 --- a/src/eDFT/build_grid.f90 +++ b/src/eDFT/build_grid.f90 @@ -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) ! 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) :: 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) double precision,intent(in) :: min_exponent(nNuc,maxL+1) double precision,intent(in) :: max_exponent(nNuc) diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 842b0cb..464bc01 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -136,14 +136,13 @@ program eDFT call read_grid(SGn,radial_precision,nRad,nAng,nGrid) ! nGrid = nRad*nAng - call allocate_grid(nNuc,ZNuc,rNuc,nShell,TotAngMomShell,ExpShell,max_ang_mom,min_exponent,max_exponent, & - radial_precision,nRad,nAng,nGrid) + call allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent,radial_precision,nAng,nGrid) allocate(root(ncart,nGrid),weight(nGrid)) ! 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) !------------------------------------------------------------------------ diff --git a/src/eDFT/elda_correlation_individual_energy.f90 b/src/eDFT/elda_correlation_individual_energy.f90 index d0b86b6..58ca58f 100644 --- a/src/eDFT/elda_correlation_individual_energy.f90 +++ b/src/eDFT/elda_correlation_individual_energy.f90 @@ -31,8 +31,8 @@ subroutine elda_correlation_individual_energy(nEns,aLF,nGrid,weight,rhow,rho,Ec) do iG=1,nGrid - ra = max(0d0,rho(iG,1)) - rb = max(0d0,rho(iG,2)) + ra = max(0d0,rhow(iG,1)) + rb = max(0d0,rhow(iG,2)) raI = max(0d0,rho(iG,1)) rbI = max(0d0,rho(iG,2)) diff --git a/src/eDFT/restricted_elda_correlation_energy.f90 b/src/eDFT/restricted_elda_correlation_energy.f90 index 06bda0a..9643911 100644 --- a/src/eDFT/restricted_elda_correlation_energy.f90 +++ b/src/eDFT/restricted_elda_correlation_energy.f90 @@ -31,11 +31,8 @@ subroutine restricted_elda_correlation_energy(aMFL,nGrid,weight,rho,Ec) 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 diff --git a/src/eDFT/restricted_elda_correlation_individual_energy.f90 b/src/eDFT/restricted_elda_correlation_individual_energy.f90 index 3cf1b78..cf052ca 100644 --- a/src/eDFT/restricted_elda_correlation_individual_energy.f90 +++ b/src/eDFT/restricted_elda_correlation_individual_energy.f90 @@ -31,7 +31,7 @@ subroutine restricted_elda_correlation_individual_energy(nEns,aLF,nGrid,weight,r do iG=1,nGrid - r = max(0d0,rho(iG)) + r = max(0d0,rhow(iG)) rI = max(0d0,rho(iG)) if(r > threshold .and. rI > threshold) then diff --git a/src/eDFT/restricted_elda_correlation_potential.f90 b/src/eDFT/restricted_elda_correlation_potential.f90 index e914a65..999ec70 100644 --- a/src/eDFT/restricted_elda_correlation_potential.f90 +++ b/src/eDFT/restricted_elda_correlation_potential.f90 @@ -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 @@ -7,7 +7,6 @@ subroutine restricted_elda_correlation_potential(nEns,aMFL,nGrid,weight,nBas,AO, ! Input variables - integer,intent(in) :: nEns double precision,intent(in) :: aMFL(3) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) diff --git a/src/eDFT/restricted_elda_exchange_energy.f90 b/src/eDFT/restricted_elda_exchange_energy.f90 deleted file mode 100644 index f56ce0e..0000000 --- a/src/eDFT/restricted_elda_exchange_energy.f90 +++ /dev/null @@ -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