From bf7f441d4d1348df0bbbec4f58ac5110a0e3693d Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 4 Jun 2020 22:41:38 +0200 Subject: [PATCH] RCC --- ..._lda_exchange_derivative_discontinuity.f90 | 98 +++++++++++++++++++ src/eDFT/RCC_lda_exchange_energy.f90 | 75 ++++++++++++++ .../RCC_lda_exchange_individual_energy.f90 | 81 +++++++++++++++ src/eDFT/RCC_lda_exchange_potential.f90 | 84 ++++++++++++++++ 4 files changed, 338 insertions(+) create mode 100644 src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 create mode 100644 src/eDFT/RCC_lda_exchange_energy.f90 create mode 100644 src/eDFT/RCC_lda_exchange_individual_energy.f90 create mode 100644 src/eDFT/RCC_lda_exchange_potential.f90 diff --git a/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 new file mode 100644 index 0000000..538113e --- /dev/null +++ b/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 @@ -0,0 +1,98 @@ +subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,ExDD) + +! Compute the restricted version of the curvature-corrected exchange ensemble derivative + + implicit none + include 'parameters.h' + +! Input variables + + 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 + + integer :: iEns,jEns + integer :: iG + double precision :: r + double precision,allocatable :: dExdw(:) + double precision,external :: Kronecker_delta + + double precision :: a1,b1,c1,w1 + double precision :: a2,b2,c2,w2 + double precision :: dCxdw1,dCxdw2 + +! Output variables + + double precision,intent(out) :: ExDD(nEns) + +! Memory allocation + + allocate(dExdw(nEns)) + +! Single excitation parameters + + a1 = 0.0d0 + b1 = 0.0d0 + c1 = 0.0d0 + +! Parameters for H2 at equilibrium + +! a2 = +0.5751782560799208d0 +! b2 = -0.021108186591137282d0 +! c2 = -0.36718902716347124d0 + +! Parameters for stretch H2 + + a2 = + 0.01922622507087411d0 + b2 = - 0.01799647558018601d0 + c2 = - 0.022945430666782573d0 + +! Parameters for He + +! a2 = 1.9125735895875828d0 +! b2 = 2.715266992840757d0 +! c2 = 2.1634223380633086d0 + + w1 = wEns(2) + w2 = wEns(3) + + dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) & + * (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)) + + dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) & + * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) + + dCxdw1 = CxLDA*dCxdw1 + dCxdw2 = CxLDA*dCxdw2 + + dExdw(:) = 0d0 + + do iG=1,nGrid + + r = max(0d0,rhow(iG)) + + if(r > threshold) then + + dExdw(1) = 0d0 + dExdw(2) = dExdw(2) + weight(iG)*dCxdw1*r**(4d0/3d0) + dExdw(3) = dExdw(3) + weight(iG)*dCxdw2*r**(4d0/3d0) + + end if + + end do + + ExDD(:) = 0d0 + + do iEns=1,nEns + do jEns=2,nEns + + ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) + + end do + end do + +end subroutine RCC_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/RCC_lda_exchange_energy.f90 b/src/eDFT/RCC_lda_exchange_energy.f90 new file mode 100644 index 0000000..274917e --- /dev/null +++ b/src/eDFT/RCC_lda_exchange_energy.f90 @@ -0,0 +1,75 @@ +subroutine RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) + +! Compute the restricted version of the curvature-corrected exchange functional + + implicit none + include 'parameters.h' + +! Input variables + + 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) + +! Local variables + + integer :: iG + double precision :: r + + double precision :: a1,b1,c1,w1 + double precision :: a2,b2,c2,w2 + double precision :: Fx1,Fx2,Cx + +! Output variables + + double precision :: Ex + +! Single excitation parameter + + a1 = 0.0d0 + b1 = 0.0d0 + c1 = 0.0d0 + +! Parameters for H2 at equilibrium + +! a2 = +0.5751782560799208d0 +! b2 = -0.021108186591137282d0 +! c2 = -0.36718902716347124d0 + +! Parameters for stretch H2 + + a2 = + 0.01922622507087411d0 + b2 = - 0.01799647558018601d0 + c2 = - 0.022945430666782573d0 + +! Parameters for He + +! a2 = 1.9125735895875828d0 +! b2 = 2.715266992840757d0 +! c2 = 2.1634223380633086d0 + + w1 = wEns(2) + Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2) + + w2 = wEns(3) + Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) + + Cx = CxLDA*Fx1*Fx2 + +! Compute GIC-LDA 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) + endif + + enddo + +end subroutine RCC_lda_exchange_energy diff --git a/src/eDFT/RCC_lda_exchange_individual_energy.f90 b/src/eDFT/RCC_lda_exchange_individual_energy.f90 new file mode 100644 index 0000000..0da49d8 --- /dev/null +++ b/src/eDFT/RCC_lda_exchange_individual_energy.f90 @@ -0,0 +1,81 @@ +subroutine RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex) + +! Compute the restricted version of the curvature-corrected exchange functional + + implicit none + include 'parameters.h' + +! Input variables + + 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) + +! Local variables + + integer :: iG + double precision :: r,rI + double precision :: e_p,dedr + + double precision :: a1,b1,c1,w1 + double precision :: a2,b2,c2,w2 + double precision :: Fx1,Fx2,Cx + +! Output variables + + double precision,intent(out) :: Ex + +! Single excitation parameter + + a1 = 0.0d0 + b1 = 0.0d0 + c1 = 0.0d0 + +! Parameters for H2 at equilibrium + +! a2 = +0.5751782560799208d0 +! b2 = -0.021108186591137282d0 +! c2 = -0.36718902716347124d0 + +! Parameters for stretch H2 + + a2 = + 0.01922622507087411d0 + b2 = - 0.01799647558018601d0 + c2 = - 0.022945430666782573d0 + +! Parameters for He + +! a2 = 1.9125735895875828d0 +! b2 = 2.715266992840757d0 +! c2 = 2.1634223380633086d0 + + w1 = wEns(2) + Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2) + + w2 = wEns(3) + Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) + + Cx = CxLDA*Fx1*Fx2 + +! Compute LDA exchange matrix in the AO basis + + Ex = 0d0 + do iG=1,nGrid + + r = max(0d0,rhow(iG)) + rI = max(0d0,rho(iG)) + + if(r > threshold .and. rI > threshold) then + + e_p = Cx*r**(1d0/3d0) + dedr = 1d0/3d0*Cx*r**(-2d0/3d0) + Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r) + + endif + + enddo + +end subroutine RCC_lda_exchange_individual_energy diff --git a/src/eDFT/RCC_lda_exchange_potential.f90 b/src/eDFT/RCC_lda_exchange_potential.f90 new file mode 100644 index 0000000..e4ccd00 --- /dev/null +++ b/src/eDFT/RCC_lda_exchange_potential.f90 @@ -0,0 +1,84 @@ +subroutine RCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + +! Compute the restricted version of the curvature-corrected exchange potential + + implicit none + include 'parameters.h' + +! Input variables + + 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) + +! Local variables + + integer :: mu,nu,iG + double precision :: r,vAO + + double precision :: a1,b1,c1,w1 + double precision :: a2,b2,c2,w2 + double precision :: Fx1,Fx2,Cx + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Single excitation parameter + + a1 = 0.0d0 + b1 = 0.0d0 + c1 = 0.0d0 + +! Parameters for H2 at equilibrium + +! a2 = +0.5751782560799208d0 +! b2 = -0.021108186591137282d0 +! c2 = -0.36718902716347124d0 + +! Parameters for stretch H2 + + a2 = + 0.01922622507087411d0 + b2 = - 0.01799647558018601d0 + c2 = - 0.022945430666782573d0 + +! Parameters for He + +! a2 = 1.9125735895875828d0 +! b2 = 2.715266992840757d0 +! c2 = 2.1634223380633086d0 + + w1 = wEns(2) + Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2) + + w2 = wEns(3) + Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) + + Cx = CxLDA*Fx1*Fx2 + +! Compute LDA exchange matrix in the AO basis + + Fx(:,:) = 0d0 + + do mu=1,nBas + do nu=1,nBas + do iG=1,nGrid + + r = max(0d0,rho(iG)) + + if(r > threshold) then + + vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) + Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*Cx*r**(1d0/3d0) + + endif + + enddo + enddo + enddo + +end subroutine RCC_lda_exchange_potential