From 1753f4190dbd1e007df223ab537b9d5f7ed4058c Mon Sep 17 00:00:00 2001 From: Clotilde Marut Date: Thu, 2 Jul 2020 14:27:38 +0200 Subject: [PATCH] clean up unrestricted --- input/dft | 2 +- src/eDFT/GOK_UKS.f90 | 12 +- ...a_correlation_derivative_discontinuity.f90 | 63 ----- ...FL20_lda_correlation_individual_energy.f90 | 81 ------ ...gy.f90 => UC16_lda_correlation_energy.f90} | 6 +- ...f90 => UC16_lda_correlation_potential.f90} | 6 +- ..._lda_exchange_derivative_discontinuity.f90 | 121 +++++++++ src/eDFT/UCC_lda_exchange_energy.f90 | 94 +++++++ .../UCC_lda_exchange_individual_energy.f90 | 97 +++++++ src/eDFT/UCC_lda_exchange_potential.f90 | 103 ++++++++ ...nergy.f90 => US51_lda_exchange_energy.f90} | 20 +- .../US51_lda_exchange_individual_energy.f90 | 47 ++++ ...al.f90 => US51_lda_exchange_potential.f90} | 4 +- ...y.f90 => UVWN5_lda_correlation_energy.f90} | 6 +- ...90 => UVWN5_lda_correlation_potential.f90} | 6 +- ...gy.f90 => UW38_lda_correlation_energy.f90} | 6 +- ...W38_lda_correlation_individual_energy.f90} | 6 +- ...f90 => UW38_lda_correlation_potential.f90} | 6 +- ...WN5_lda_correlation_individual_energy.f90} | 4 +- src/eDFT/eDFT_UKS.f90 | 10 +- src/eDFT/individual_energy.f90 | 149 ----------- .../lda_exchange_derivative_discontinuity.f90 | 6 +- src/eDFT/lda_exchange_energy.f90 | 8 +- src/eDFT/lda_exchange_individual_energy.f90 | 8 + src/eDFT/lda_exchange_potential.f90 | 8 +- .../print_unrestricted_individual_energy.f90 | 204 +++++++++++++++ src/eDFT/restricted_individual_energy.f90 | 6 +- src/eDFT/unrestricted_auxiliary_energy.f90 | 52 ++++ ..._correlation_derivative_discontinuity.f90} | 6 +- ...90 => unrestricted_correlation_energy.f90} | 16 +- ...tricted_correlation_individual_energy.f90} | 7 +- ...=> unrestricted_correlation_potential.f90} | 24 +- ...> unrestricted_gga_correlation_energy.f90} | 6 +- ...nrestricted_gga_correlation_potential.f90} | 6 +- src/eDFT/unrestricted_individual_energy.f90 | 240 ++++++++++++++++++ ..._correlation_derivative_discontinuity.f90} | 14 +- ...> unrestricted_lda_correlation_energy.f90} | 18 +- ...ted_lda_correlation_individual_energy.f90} | 21 +- ...nrestricted_lda_correlation_potential.f90} | 16 +- 39 files changed, 1098 insertions(+), 417 deletions(-) delete mode 100644 src/eDFT/MFL20_lda_correlation_derivative_discontinuity.f90 delete mode 100644 src/eDFT/MFL20_lda_correlation_individual_energy.f90 rename src/eDFT/{C16_lda_correlation_energy.f90 => UC16_lda_correlation_energy.f90} (91%) rename src/eDFT/{C16_lda_correlation_potential.f90 => UC16_lda_correlation_potential.f90} (95%) create mode 100644 src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 create mode 100644 src/eDFT/UCC_lda_exchange_energy.f90 create mode 100644 src/eDFT/UCC_lda_exchange_individual_energy.f90 create mode 100644 src/eDFT/UCC_lda_exchange_potential.f90 rename src/eDFT/{S51_lda_exchange_energy.f90 => US51_lda_exchange_energy.f90} (53%) create mode 100644 src/eDFT/US51_lda_exchange_individual_energy.f90 rename src/eDFT/{S51_lda_exchange_potential.f90 => US51_lda_exchange_potential.f90} (89%) rename src/eDFT/{VWN5_lda_correlation_energy.f90 => UVWN5_lda_correlation_energy.f90} (95%) rename src/eDFT/{VWN5_lda_correlation_potential.f90 => UVWN5_lda_correlation_potential.f90} (97%) rename src/eDFT/{W38_lda_correlation_energy.f90 => UW38_lda_correlation_energy.f90} (82%) rename src/eDFT/{W38_lda_correlation_individual_energy.f90 => UW38_lda_correlation_individual_energy.f90} (85%) rename src/eDFT/{W38_lda_correlation_potential.f90 => UW38_lda_correlation_potential.f90} (88%) rename src/eDFT/{VWN5_lda_correlation_individual_energy.f90 => UWN5_lda_correlation_individual_energy.f90} (98%) delete mode 100644 src/eDFT/individual_energy.f90 create mode 100644 src/eDFT/print_unrestricted_individual_energy.f90 create mode 100644 src/eDFT/unrestricted_auxiliary_energy.f90 rename src/eDFT/{correlation_derivative_discontinuity.f90 => unrestricted_correlation_derivative_discontinuity.f90} (79%) rename src/eDFT/{correlation_energy.f90 => unrestricted_correlation_energy.f90} (58%) rename src/eDFT/{correlation_individual_energy.f90 => unrestricted_correlation_individual_energy.f90} (78%) rename src/eDFT/{correlation_potential.f90 => unrestricted_correlation_potential.f90} (64%) rename src/eDFT/{gga_correlation_energy.f90 => unrestricted_gga_correlation_energy.f90} (82%) rename src/eDFT/{gga_correlation_potential.f90 => unrestricted_gga_correlation_potential.f90} (77%) create mode 100644 src/eDFT/unrestricted_individual_energy.f90 rename src/eDFT/{lda_correlation_derivative_discontinuity.f90 => unrestricted_lda_correlation_derivative_discontinuity.f90} (69%) rename src/eDFT/{lda_correlation_energy.f90 => unrestricted_lda_correlation_energy.f90} (64%) rename src/eDFT/{lda_correlation_individual_energy.f90 => unrestricted_lda_correlation_individual_energy.f90} (53%) rename src/eDFT/{lda_correlation_potential.f90 => unrestricted_lda_correlation_potential.f90} (69%) diff --git a/input/dft b/input/dft index f3deb40..10aa014 100644 --- a/input/dft +++ b/input/dft @@ -19,6 +19,6 @@ # Number of states in ensemble (nEns) 3 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 1.0 0.0 + 0.333 0.333 # GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type 32 0.00001 T 5 1 1 diff --git a/src/eDFT/GOK_UKS.f90 b/src/eDFT/GOK_UKS.f90 index 41556ad..0ab3741 100644 --- a/src/eDFT/GOK_UKS.f90 +++ b/src/eDFT/GOK_UKS.f90 @@ -252,7 +252,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Compute correlation potential - call correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rhow(:,:),drhow(:,:,:),Fc(:,:,:)) + call unrestricted_correlation_potential(c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rhow,drhow,Fc) ! Build Fock operator do ispin=1,nspin @@ -310,7 +310,7 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Correlation energy - call correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),Ec) + call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),Ec) ! Total energy @@ -355,9 +355,9 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres ! Compute individual energies from ensemble energy !------------------------------------------------------------------------ - call individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas, & - AO(:,:),dAO(:,:,:),nO(:),nV(:),T(:,:),V(:,:),ERI(:,:,:,:),ENuc, & - Pw(:,:,:),rhow(:,:),drhow(:,:,:),J(:,:,:),Fx(:,:,:),FxHF(:,:,:), & - Fc(:,:,:),P(:,:,:,:),rho(:,:,:),drho(:,:,:,:),E(:),Om(:)) + call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas, & + AO(:,:),dAO(:,:,:),nO(:),nV(:),T(:,:),V(:,:),ERI(:,:,:,:),ENuc, & + Pw(:,:,:),rhow(:,:),drhow(:,:,:),J(:,:,:),Fx(:,:,:),FxHF(:,:,:), & + Fc(:,:,:),P(:,:,:,:),rho(:,:,:),drho(:,:,:,:),E(:),Om(:)) end subroutine GOK_UKS diff --git a/src/eDFT/MFL20_lda_correlation_derivative_discontinuity.f90 b/src/eDFT/MFL20_lda_correlation_derivative_discontinuity.f90 deleted file mode 100644 index b445413..0000000 --- a/src/eDFT/MFL20_lda_correlation_derivative_discontinuity.f90 +++ /dev/null @@ -1,63 +0,0 @@ -subroutine MFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,Ec) - -! Compute eLDA correlation part of the derivative discontinuity - - 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,nspin) - -! Local variables - - integer :: iEns,jEns - double precision,allocatable :: aMFL(:,:) - double precision :: dEc(nsp,nEns) - double precision,external :: Kronecker_delta - -! Output variables - - double precision,intent(out) :: Ec(nsp,nEns) - -! Allocation - - allocate(aMFL(3,nEns)) - -! Parameters for weight-dependent LDA correlation functional - - aMFL(1,1) = -0.0238184d0 - aMFL(2,1) = +0.00575719d0 - aMFL(3,1) = +0.0830576d0 - - aMFL(1,2) = -0.0282814d0 - aMFL(2,2) = +0.00340758d0 - aMFL(3,2) = +0.0663967d0 - - aMFL(1,3) = -0.0144633d0 - aMFL(2,3) = -0.0504501d0 - aMFL(3,3) = +0.0331287d0 - -! Compute correlation energy for ground, singly-excited and doubly-excited states - - do iEns=1,nEns - - call elda_correlation_energy(aMFL(:,iEns),nGrid,weight(:),rhow(:,:),dEc(:,iEns)) - - end do - - Ec(:,:) = 0d0 - - do iEns=1,nEns - do jEns=1,nEns - - Ec(:,iEns) = Ec(:,iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*(dEc(:,jEns) - dEc(:,1)) - - end do - end do - -end subroutine MFL20_lda_correlation_derivative_discontinuity diff --git a/src/eDFT/MFL20_lda_correlation_individual_energy.f90 b/src/eDFT/MFL20_lda_correlation_individual_energy.f90 deleted file mode 100644 index 00bb982..0000000 --- a/src/eDFT/MFL20_lda_correlation_individual_energy.f90 +++ /dev/null @@ -1,81 +0,0 @@ -subroutine MFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ec) - -! Compute eLDA correlation energy - - 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,nspin) - double precision,intent(in) :: rho(nGrid,nspin) - -! Local variables - - logical :: LDA_centered = .true. - integer :: iEns,isp - double precision :: EcLDA(nsp) - double precision,allocatable :: aMFL(:,:) - double precision,allocatable :: EceLDA(:,:) - -! Output variables - - double precision :: Ec(nsp) - -! Allocation - - allocate(aMFL(3,nEns),EceLDA(nsp,nEns)) - -! Parameters for weight-dependent LDA correlation functional - - aMFL(1,1) = -0.0238184d0 - aMFL(2,1) = +0.00575719d0 - aMFL(3,1) = +0.0830576d0 - - aMFL(1,2) = -0.0282814d0 - aMFL(2,2) = +0.00340758d0 - aMFL(3,2) = +0.0663967d0 - - aMFL(1,3) = -0.0144633d0 - aMFL(2,3) = -0.0504501d0 - aMFL(3,3) = +0.0331287d0 - -! Compute correlation energy for ground, singly-excited and doubly-excited states - - do iEns=1,nEns - - call elda_correlation_individual_energy(aMFL(:,iEns),nGrid,weight(:),rhow(:,:),rho(:,:),EceLDA(:,iEns)) - - end do - -! LDA-centered functional - - EcLDA(:) = 0d0 - - if(LDA_centered) then - - call W38_lda_correlation_individual_energy(nGrid,weight(:),rhow(:,:),rho(:,:),EcLDA(:)) - - do iEns=1,nEns - do isp=1,nsp - EceLDA(isp,iEns) = EceLDA(isp,iEns) + EcLDA(isp) - EceLDA(isp,1) - end do - end do - - end if - -! Weight-denpendent functional for ensembles - - Ec(:) = 0d0 - - do iEns=1,nEns - do isp=1,nsp - Ec(isp) = Ec(isp) + wEns(iEns)*EceLDA(isp,iEns) - enddo - enddo - -end subroutine MFL20_lda_correlation_individual_energy diff --git a/src/eDFT/C16_lda_correlation_energy.f90 b/src/eDFT/UC16_lda_correlation_energy.f90 similarity index 91% rename from src/eDFT/C16_lda_correlation_energy.f90 rename to src/eDFT/UC16_lda_correlation_energy.f90 index 1ea0d89..d069e4f 100644 --- a/src/eDFT/C16_lda_correlation_energy.f90 +++ b/src/eDFT/UC16_lda_correlation_energy.f90 @@ -1,6 +1,6 @@ -subroutine C16_lda_correlation_energy(nGrid,weight,rho,Ec) +subroutine UC16_lda_correlation_energy(nGrid,weight,rho,Ec) -! Compute Chachiyo's LDA correlation energy +! Compute unrestricted Chachiyo's LDA correlation energy implicit none include 'parameters.h' @@ -90,4 +90,4 @@ subroutine C16_lda_correlation_energy(nGrid,weight,rho,Ec) Ec(2) = Ec(2) - Ec(1) - Ec(3) -end subroutine C16_lda_correlation_energy +end subroutine UC16_lda_correlation_energy diff --git a/src/eDFT/C16_lda_correlation_potential.f90 b/src/eDFT/UC16_lda_correlation_potential.f90 similarity index 95% rename from src/eDFT/C16_lda_correlation_potential.f90 rename to src/eDFT/UC16_lda_correlation_potential.f90 index 60d2c8d..43fa599 100644 --- a/src/eDFT/C16_lda_correlation_potential.f90 +++ b/src/eDFT/UC16_lda_correlation_potential.f90 @@ -1,6 +1,6 @@ -subroutine C16_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) +subroutine UC16_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) -! Compute LDA correlation potential +! Compute unrestricted LDA correlation potential implicit none include 'parameters.h' @@ -128,4 +128,4 @@ include 'parameters.h' enddo enddo -end subroutine C16_lda_correlation_potential +end subroutine UC16_lda_correlation_potential diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 new file mode 100644 index 0000000..9941626 --- /dev/null +++ b/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 @@ -0,0 +1,121 @@ +subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,ExDD) + +! Compute the unrestricted 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,alpha + 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 + +! Parameters for He N -> N-1 + + a1 = 0.420243d0 + b1 = 0.0700561d0 + c1 = -0.288301d0 + +! Parameters for He N -> N+1 + + a2 = 0.135068d0 + b2 = -0.00774769d0 + c2 = -0.0278205d0 + +! Cx coefficient for unrestricted Slater LDA exchange + + alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + + w1 = wEns(2) + w2 = wEns(3) + +! Double weight-dependency + + 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))) + +! Single weight-dependency + +! dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) +! dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) + + dCxdw1 = alpha*dCxdw1 + dCxdw2 = alpha*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 UCC_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/UCC_lda_exchange_energy.f90 b/src/eDFT/UCC_lda_exchange_energy.f90 new file mode 100644 index 0000000..73507e6 --- /dev/null +++ b/src/eDFT/UCC_lda_exchange_energy.f90 @@ -0,0 +1,94 @@ +subroutine UCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) + +! Compute the unrestricted 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,alpha + + 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 + +! Parameters for He N -> N-1 + + a1 = 0.420243d0 + b1 = 0.0700561d0 + c1 = -0.288301d0 + +! Parameters for He N -> N+1 + + a2 = 0.135068d0 + b2 = -0.00774769d0 + c2 = -0.0278205d0 + +! Cx coefficient for unrestricted Slater LDA exchange + + alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + +! Fx1 for states N and N-1 +! Fx2 for states N and N+1 + + 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 = alpha*Fx2*Fx1 + +! 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 UCC_lda_exchange_energy diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/UCC_lda_exchange_individual_energy.f90 new file mode 100644 index 0000000..6c52e6b --- /dev/null +++ b/src/eDFT/UCC_lda_exchange_individual_energy.f90 @@ -0,0 +1,97 @@ +subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,Ex) + +! Compute the unrestricted 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,alpha + 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 + +! Parameters for He N -> N-1 + + a1 = 0.420243d0 + b1 = 0.0700561d0 + c1 = -0.288301d0 + +! Parameters for He N -> N+1 + + a2 = 0.135068d0 + b2 = -0.00774769d0 + c2 = -0.0278205d0 + +! Cx coefficient for unrestricted Slater LDA exchange + + alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + + 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 = alpha*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 UCC_lda_exchange_individual_energy diff --git a/src/eDFT/UCC_lda_exchange_potential.f90 b/src/eDFT/UCC_lda_exchange_potential.f90 new file mode 100644 index 0000000..4e73d74 --- /dev/null +++ b/src/eDFT/UCC_lda_exchange_potential.f90 @@ -0,0 +1,103 @@ +subroutine UCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + +! Compute the unrestricted 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,alpha + + 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 + +! Parameters for He N -> N-1 + + a1 = 0.420243d0 + b1 = 0.0700561d0 + c1 = -0.288301d0 + +! Parameters for He N -> N+1 + + a2 = 0.135068d0 + b2 = -0.00774769d0 + c2 = -0.0278205d0 + +! Cx coefficient for unrestricted Slater LDA exchange + + alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + +! Fx1 for states N and N-1 +! Fx2 for states N and N+1 + + 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 = alpha*Fx2*Fx1 + +! 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 UCC_lda_exchange_potential diff --git a/src/eDFT/S51_lda_exchange_energy.f90 b/src/eDFT/US51_lda_exchange_energy.f90 similarity index 53% rename from src/eDFT/S51_lda_exchange_energy.f90 rename to src/eDFT/US51_lda_exchange_energy.f90 index a273bf5..920e059 100644 --- a/src/eDFT/S51_lda_exchange_energy.f90 +++ b/src/eDFT/US51_lda_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine S51_lda_exchange_energy(nGrid,weight,rho,Ex) +subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex,wEns,nEns) ! Compute Slater's LDA exchange energy @@ -11,19 +11,33 @@ subroutine S51_lda_exchange_energy(nGrid,weight,rho,Ex) double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rho(nGrid) + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + ! Local variables integer :: iG - double precision :: alpha,r + double precision :: alpha,r,alphaw,a2,b2,c2,a1,b1,c1 ! Output variables double precision :: Ex +! Cxw2 parameters for He N->N+1 +! a2 = 0.135068d0 +! b2 = -0.00774769d0 +! c2 = -0.0278205d0 + +! Cxw1 parameters for He N->N-1 +! a1 = 0.420243d0 +! b1 = 0.0700561d0 +! c1 = -0.288301d0 + ! Cx coefficient for Slater LDA exchange alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) +! alphaw = alpha*(1d0 - wEns(2)*(1d0 - wEns(2))*(a1 + b1*(wEns(2) - 0.5d0) + c1*(wEns(2) - 0.5d0)**2)) ! Compute LDA exchange energy Ex = 0d0 @@ -39,4 +53,4 @@ subroutine S51_lda_exchange_energy(nGrid,weight,rho,Ex) enddo -end subroutine S51_lda_exchange_energy +end subroutine US51_lda_exchange_energy diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/US51_lda_exchange_individual_energy.f90 new file mode 100644 index 0000000..74b1749 --- /dev/null +++ b/src/eDFT/US51_lda_exchange_individual_energy.f90 @@ -0,0 +1,47 @@ +subroutine US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex) + +! Compute the restricted version of Slater's LDA exchange individual energy + + implicit none + include 'parameters.h' + +! Input variables + + 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,alpha + double precision :: e,dedr + +! Output variables + + double precision,intent(out) :: Ex + +! Compute LDA exchange matrix in the AO basis + + alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0) + + + Ex = 0d0 + do iG=1,nGrid + + r = max(0d0,rhow(iG)) + rI = max(0d0,rho(iG)) + + if(r > threshold .or. rI > threshold) then + + e = alpha*r**(1d0/3d0) + dedr = 1d0/3d0*alpha*r**(-2d0/3d0) + + Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) + + endif + + enddo + +end subroutine US51_lda_exchange_individual_energy diff --git a/src/eDFT/S51_lda_exchange_potential.f90 b/src/eDFT/US51_lda_exchange_potential.f90 similarity index 89% rename from src/eDFT/S51_lda_exchange_potential.f90 rename to src/eDFT/US51_lda_exchange_potential.f90 index 802332c..4066379 100644 --- a/src/eDFT/S51_lda_exchange_potential.f90 +++ b/src/eDFT/US51_lda_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) +subroutine US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) ! Compute Slater's LDA exchange potential @@ -47,4 +47,4 @@ subroutine S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) enddo enddo -end subroutine S51_lda_exchange_potential +end subroutine US51_lda_exchange_potential diff --git a/src/eDFT/VWN5_lda_correlation_energy.f90 b/src/eDFT/UVWN5_lda_correlation_energy.f90 similarity index 95% rename from src/eDFT/VWN5_lda_correlation_energy.f90 rename to src/eDFT/UVWN5_lda_correlation_energy.f90 index 0cff57e..07e8d24 100644 --- a/src/eDFT/VWN5_lda_correlation_energy.f90 +++ b/src/eDFT/UVWN5_lda_correlation_energy.f90 @@ -1,6 +1,6 @@ -subroutine VWN5_lda_correlation_energy(nGrid,weight,rho,Ec) +subroutine UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec) -! Compute VWN5 LDA correlation energy +! Compute unrestricted VWN5 LDA correlation energy implicit none @@ -134,4 +134,4 @@ subroutine VWN5_lda_correlation_energy(nGrid,weight,rho,Ec) Ec(2) = Ec(2) - Ec(1) - Ec(3) -end subroutine VWN5_lda_correlation_energy +end subroutine UVWN5_lda_correlation_energy diff --git a/src/eDFT/VWN5_lda_correlation_potential.f90 b/src/eDFT/UVWN5_lda_correlation_potential.f90 similarity index 97% rename from src/eDFT/VWN5_lda_correlation_potential.f90 rename to src/eDFT/UVWN5_lda_correlation_potential.f90 index 2a11aaf..d36e8a1 100644 --- a/src/eDFT/VWN5_lda_correlation_potential.f90 +++ b/src/eDFT/UVWN5_lda_correlation_potential.f90 @@ -1,6 +1,6 @@ -subroutine VWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) +subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) -! Compute VWN5 LDA correlation potential +! Compute unrestricted VWN5 LDA correlation potential implicit none @@ -199,4 +199,4 @@ subroutine VWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) end do end do -end subroutine VWN5_lda_correlation_potential +end subroutine UVWN5_lda_correlation_potential diff --git a/src/eDFT/W38_lda_correlation_energy.f90 b/src/eDFT/UW38_lda_correlation_energy.f90 similarity index 82% rename from src/eDFT/W38_lda_correlation_energy.f90 rename to src/eDFT/UW38_lda_correlation_energy.f90 index a13e806..c6f2905 100644 --- a/src/eDFT/W38_lda_correlation_energy.f90 +++ b/src/eDFT/UW38_lda_correlation_energy.f90 @@ -1,6 +1,6 @@ -subroutine W38_lda_correlation_energy(nGrid,weight,rho,Ec) +subroutine UW38_lda_correlation_energy(nGrid,weight,rho,Ec) -! Compute Wigner's LDA correlation energy +! Compute the unrestricted version of the Wigner's LDA correlation energy implicit none include 'parameters.h' @@ -49,4 +49,4 @@ subroutine W38_lda_correlation_energy(nGrid,weight,rho,Ec) Ec(2) = -4d0*a*Ec(2) -end subroutine W38_lda_correlation_energy +end subroutine UW38_lda_correlation_energy diff --git a/src/eDFT/W38_lda_correlation_individual_energy.f90 b/src/eDFT/UW38_lda_correlation_individual_energy.f90 similarity index 85% rename from src/eDFT/W38_lda_correlation_individual_energy.f90 rename to src/eDFT/UW38_lda_correlation_individual_energy.f90 index 31f7bd8..dc00815 100644 --- a/src/eDFT/W38_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UW38_lda_correlation_individual_energy.f90 @@ -1,6 +1,6 @@ -subroutine W38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) +subroutine UW38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) -! Compute Wigner's LDA individual energy +! Compute the unrestricted version of the Wigner's LDA individual energy implicit none @@ -59,4 +59,4 @@ subroutine W38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) Ec(2) = -4d0*a*Ec(2) -end subroutine W38_lda_correlation_individual_energy +end subroutine UW38_lda_correlation_individual_energy diff --git a/src/eDFT/W38_lda_correlation_potential.f90 b/src/eDFT/UW38_lda_correlation_potential.f90 similarity index 88% rename from src/eDFT/W38_lda_correlation_potential.f90 rename to src/eDFT/UW38_lda_correlation_potential.f90 index 755974b..8c1c200 100644 --- a/src/eDFT/W38_lda_correlation_potential.f90 +++ b/src/eDFT/UW38_lda_correlation_potential.f90 @@ -1,6 +1,6 @@ -subroutine W38_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) +subroutine UW38_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) -! Compute Wigner's LDA correlation potential +! Compute the unrestricted version of the Wigner's LDA correlation potential implicit none include 'parameters.h' @@ -73,4 +73,4 @@ include 'parameters.h' Fc(:,:,:) = -4d0*a*Fc(:,:,:) -end subroutine W38_lda_correlation_potential +end subroutine UW38_lda_correlation_potential diff --git a/src/eDFT/VWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UWN5_lda_correlation_individual_energy.f90 similarity index 98% rename from src/eDFT/VWN5_lda_correlation_individual_energy.f90 rename to src/eDFT/UWN5_lda_correlation_individual_energy.f90 index 513354c..03b185f 100644 --- a/src/eDFT/VWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UWN5_lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine VWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) +subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) ! Compute VWN5 LDA correlation potential @@ -201,4 +201,4 @@ subroutine VWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) end do -end subroutine VWN5_lda_correlation_individual_energy +end subroutine UVWN5_lda_correlation_individual_energy diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index 8e89843..791e68f 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -242,7 +242,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thre ! Compute correlation potential - call correlation_potential(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rhow(:,:),drhow(:,:,:),Fc(:,:,:)) + call unrestricted_correlation_potential(c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rhow,drhow,Fc) ! Build Fock operator do ispin=1,nspin @@ -319,7 +319,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thre ! Correlation energy - call correlation_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),Ec) + call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) ! Total energy @@ -364,9 +364,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thre ! Compute individual energies from ensemble energy !------------------------------------------------------------------------ -! call individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),nBas, & -! AO(:,:),dAO(:,:,:),nO(:),nV(:),T(:,:),V(:,:),ERI(:,:,:,:),ENuc, & -! Pw(:,:,:),rhow(:,:),drhow(:,:,:),J(:,:,:),Fx(:,:,:),FxHF(:,:,:), & -! Fc(:,:,:),P(:,:,:,:),rho(:,:,:),drho(:,:,:,:),E(:),Om(:)) + call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,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,E,Om) end subroutine eDFT_UKS diff --git a/src/eDFT/individual_energy.f90 b/src/eDFT/individual_energy.f90 deleted file mode 100644 index cefa131..0000000 --- a/src/eDFT/individual_energy.f90 +++ /dev/null @@ -1,149 +0,0 @@ -subroutine individual_energy(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO, & - nO,nV,T,V,ERI,ENuc,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,E,Om) - -! 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 - 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) - - integer,intent(in) :: nO(nspin),nV(nspin) - 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) :: Pw(nBas,nBas) - double precision,intent(in) :: rhow(nGrid,nspin) - double precision,intent(in) :: drhow(ncart,nGrid,nspin) - - double precision,intent(in) :: P(nBas,nBas,nspin,nEns) - double precision,intent(in) :: rho(nGrid,nspin,nEns) - double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns) - - double precision,intent(in) :: J(nBas,nBas,nspin) - double precision,intent(in) :: Fx(nBas,nBas,nspin) - double precision,intent(in) :: FxHF(nBas,nBas,nspin) - double precision,intent(in) :: Fc(nBas,nBas,nspin) - -! Local variables - - double precision :: ET(nspin,nEns) - double precision :: EV(nspin,nEns) - double precision :: EJ(nsp,nEns) - double precision :: Ex(nspin,nEns) - double precision :: Ec(nsp,nEns) - double precision :: EcDD(nsp,nEns) - - double precision,external :: trace_matrix - - integer :: ispin,iEns - -! Output variables - - double precision,intent(out) :: E(nEns) - double precision,intent(out) :: Om(nEns) - -!------------------------------------------------------------------------ -! Kinetic energy -!------------------------------------------------------------------------ - - do ispin=1,nspin - do iEns=1,nEns - ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) - end do - end do - -!------------------------------------------------------------------------ -! Potential energy -!------------------------------------------------------------------------ - - do iEns=1,nEns - do ispin=1,nspin - EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) - end do - end do - -!------------------------------------------------------------------------ -! Hartree energy -!------------------------------------------------------------------------ - - do iEns=1,nEns - - do ispin=1,nspin - call hartree_coulomb(nBas,P(:,:,ispin,iEns),ERI,J(:,:,ispin)) - end do - - EJ(1,iEns) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) - EJ(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) - EJ(3,iEns) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) - - end do - -!------------------------------------------------------------------------ -! Exchange energy -!------------------------------------------------------------------------ - - do iEns=1,nEns - do ispin=1,nspin - - call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,ispin,iEns),ERI(:,:,:,:), & - AO(:,:),dAO(:,:,:),rho(:,ispin,iEns),drho(:,:,ispin,iEns),Fx(:,:,ispin),FxHF(:,:,ispin)) - call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,ispin,iEns),FxHF(:,:,ispin), & - rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns)) - - end do - end do - -!------------------------------------------------------------------------ -! Correlation energy -!------------------------------------------------------------------------ - - do iEns=1,nEns - - call correlation_individual_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:), & - rho(:,:,iEns),drho(:,:,:,iEns),Ec(:,iEns)) - - end do - -!------------------------------------------------------------------------ -! Compute derivative discontinuities -!------------------------------------------------------------------------ - - call correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),EcDD(:,:)) - -!------------------------------------------------------------------------ -! Total energy -!------------------------------------------------------------------------ - - do iEns=1,nEns - E(iEns) = ENuc + sum(ET(:,iEns)) + sum(EV(:,iEns)) + sum(EJ(:,iEns)) & - + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + sum(EcDD(:,iEns)) - end do - -!------------------------------------------------------------------------ -! Excitation energies -!------------------------------------------------------------------------ - - do iEns=1,nEns - Om(iEns) = E(iEns) - E(1) - end do - -!------------------------------------------------------------------------ -! Dump results -!------------------------------------------------------------------------ - - call print_individual_energy(nEns,EJ,Ex,Ec,EcDD,E,Om) - -end subroutine individual_energy diff --git a/src/eDFT/lda_exchange_derivative_discontinuity.f90 b/src/eDFT/lda_exchange_derivative_discontinuity.f90 index 75296f7..ed0ce0c 100644 --- a/src/eDFT/lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/lda_exchange_derivative_discontinuity.f90 @@ -25,7 +25,7 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow select case (DFA) - case ('S51') + case ('US51') ExDD(:) = 0d0 @@ -41,6 +41,10 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow call RCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) + case ('UCC') + + call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) + case default call print_warning('!!! LDA exchange derivative discontinuity not available !!!') diff --git a/src/eDFT/lda_exchange_energy.f90 b/src/eDFT/lda_exchange_energy.f90 index 5b652b4..e513b53 100644 --- a/src/eDFT/lda_exchange_energy.f90 +++ b/src/eDFT/lda_exchange_energy.f90 @@ -23,9 +23,9 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) select case (DFA) - case ('S51') + case ('US51') - call S51_lda_exchange_energy(nGrid,weight,rho,Ex) + call US51_lda_exchange_energy(nGrid,weight,rho,Ex,wEns,nEns) case ('RS51') @@ -39,6 +39,10 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) call RCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) + case ('UCC') + + call UCC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) + case default call print_warning('!!! LDA exchange functional not available !!!') diff --git a/src/eDFT/lda_exchange_individual_energy.f90 b/src/eDFT/lda_exchange_individual_energy.f90 index bf8caa6..e9ea3cc 100644 --- a/src/eDFT/lda_exchange_individual_energy.f90 +++ b/src/eDFT/lda_exchange_individual_energy.f90 @@ -24,6 +24,10 @@ subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weigh select case (DFA) + case ('US51') + + call US51_lda_exchange_individual_energy(nGrid,weight(:),rhow(:),rho(:),Ex) + case ('RS51') call RS51_lda_exchange_individual_energy(nGrid,weight(:),rhow(:),rho(:),Ex) @@ -36,6 +40,10 @@ subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weigh call RCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) + case ('UCC') + + call UCC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:),rho(:),Ex) + case default call print_warning('!!! LDA exchange individual energy not available !!!') diff --git a/src/eDFT/lda_exchange_potential.f90 b/src/eDFT/lda_exchange_potential.f90 index 38d77ab..14d46b8 100644 --- a/src/eDFT/lda_exchange_potential.f90 +++ b/src/eDFT/lda_exchange_potential.f90 @@ -30,9 +30,9 @@ subroutine lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,A call RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) - case ('S51') + case ('US51') - call S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) + call US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) case ('RMFL20') @@ -42,6 +42,10 @@ subroutine lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,A call RCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + case ('UCC') + + call UCC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + case default call print_warning('!!! LDA exchange functional not available !!!') diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90 new file mode 100644 index 0000000..371d7f8 --- /dev/null +++ b/src/eDFT/print_unrestricted_individual_energy.f90 @@ -0,0 +1,204 @@ +subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Eaux,ExDD,EcDD,ExcDD,E, & + Om,Omx,Omc,Omxc,Omaux,OmxDD,OmcDD,OmxcDD) + +! Print individual energies for eDFT calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nEns + double precision,intent(in) :: ENuc + double precision,intent(in) :: Ew + double precision,intent(in) :: ET(nspin,nEns) + double precision,intent(in) :: EV(nspin,nEns) + double precision,intent(in) :: EJ(nsp,nEns) + double precision,intent(in) :: Ex(nspin,nEns), Ec(nsp,nEns), Exc(nEns) + double precision,intent(in) :: Eaux(nspin,nEns) + double precision,intent(in) :: ExDD(nspin,nEns), EcDD(nsp,nEns), ExcDD(nspin,nEns) + double precision,intent(in) :: Omx(nEns), Omc(nEns), Omxc(nEns) + double precision,intent(in) :: Omaux(nEns) + double precision,intent(in) :: OmxDD(nEns),OmcDD(nEns),OmxcDD(nEns) + double precision,intent(in) :: E(nEns) + double precision,intent(in) :: Om(nEns) + +! Local variables + + integer :: iEns + +!------------------------------------------------------------------------ +! Ensemble energies +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' ENSEMBLE ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A44,F16.10,A3)') ' Ensemble energy: ',Ew + ENuc,' au' + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Kinetic energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL KINETIC ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Kinetic energy state ',iEns,': ',sum(ET(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Potential energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL POTENTIAL ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Potential energy state ',iEns,': ',sum(EV(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Hartree energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL HARTREE ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Hartree energy state ',iEns,': ',sum(EJ(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Exchange energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL EXCHANGE ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Exchange energy state ',iEns,': ',sum(Ex(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Correlation energy +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' INDIVIDUAL CORRELATION ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Correlation energy state ',iEns,': ',sum(Ec(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Auxiliary energies +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' AUXILIARY KS ENERGIES' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') 'Auxiliary KS energy state ',iEns,': ',sum(Eaux(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Compute derivative discontinuities +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' ENSEMBLE DERIVATIVE CONTRIBUTIONS' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,*) + write(*,'(A40,I2,A2,F16.10,A3)') ' x ensemble derivative state ',iEns,': ',sum(ExDD(:,iEns)), ' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' c ensemble derivative state ',iEns,': ',sum(EcDD(:,iEns)), ' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' xc ensemble derivative state ',iEns,': ',sum(ExcDD(:,iEns)),' au' + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +!------------------------------------------------------------------------ +! Total and Excitation energies +!------------------------------------------------------------------------ + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' EXCITATION ENERGIES FROM AUXILIARY ENERGIES ' + write(*,'(A60)') '-------------------------------------------------' + + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Omaux(iEns)+OmxcDD(iEns),' au' + write(*,*) + write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns),' au' + write(*,*) + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',(Omaux(iEns)+OmxcDD(iEns))*HaToeV,' eV' + write(*,*) + write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV' + write(*,*) + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' EXCITATION ENERGIES FROM INDIVIDUAL ENERGIES ' + write(*,'(A60)') '-------------------------------------------------' + do iEns=1,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' + end do + write(*,'(A60)') '-------------------------------------------------' + + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns), ' au' + write(*,*) + write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns), ' au' + write(*,*) + write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns),' au' + write(*,*) + end do + write(*,'(A60)') '-------------------------------------------------' + + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV, ' eV' + write(*,*) + write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' c energy contribution : ',Omc(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' xc energy contribution : ',Omxc(iEns)*HaToeV, ' eV' + write(*,*) + write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV' + write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV' + write(*,*) + end do + write(*,'(A60)') '-------------------------------------------------' + write(*,*) + +end subroutine print_unrestricted_individual_energy diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index 48c6bfd..29ce420 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -75,7 +75,7 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n end do !------------------------------------------------------------------------ -! Individua Hartree energy +! Individual Hartree energy !------------------------------------------------------------------------ do iEns=1,nEns @@ -89,8 +89,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n !------------------------------------------------------------------------ do iEns=1,nEns - call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,ERI(:,:,:,:), & - Pw(:,:),P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns)) + call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),nBas,ERI(:,:,:,:), & + Pw(:,:),P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns)) end do !------------------------------------------------------------------------ diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/unrestricted_auxiliary_energy.f90 new file mode 100644 index 0000000..b4b477d --- /dev/null +++ b/src/eDFT/unrestricted_auxiliary_energy.f90 @@ -0,0 +1,52 @@ +subroutine unrestricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux) + +! Compute the auxiliary KS energies + + implicit none + + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nEns + integer,intent(in) :: nO(nspin) + double precision,intent(in) :: eps(nBas,nspin) + +! Local variables + + integer :: iEns,ispin + +! Output variables + + double precision,intent(out) :: Eaux(nspin,nEns) + +! N-electron ground state + + iEns = 1 + do ispin=1,nspin + Eaux(ispin,iEns) = sum(eps(1:nO(ispin),ispin)) + end do + +! N-1-electron ground state + + iEns = 2 + + Eaux(1,iEns) = sum(eps(1:nO(1),1)) + + if(nO(2) > 1) then + Eaux(2,iEns) = sum(eps(1:nO(2)-1,2)) + else + Eaux(2,iEns) = 0d0 + end if + + +! N+1 ground state + + iEns = 3 + + Eaux(1,iEns) = sum(eps(1:nO(1)+1,1)) + Eaux(2,iEns) = sum(eps(1:nO(2),2)) + + +end subroutine unrestricted_auxiliary_energy diff --git a/src/eDFT/correlation_derivative_discontinuity.f90 b/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 similarity index 79% rename from src/eDFT/correlation_derivative_discontinuity.f90 rename to src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 index 4cb46aa..e04ff0d 100644 --- a/src/eDFT/correlation_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) +subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) ! Compute the correlation part of the derivative discontinuity @@ -36,7 +36,7 @@ subroutine correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight, case(1) - call lda_correlation_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),Ec(:,:)) + call unrestricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! GGA functionals @@ -62,4 +62,4 @@ subroutine correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight, end select -end subroutine correlation_derivative_discontinuity +end subroutine unrestricted_correlation_derivative_discontinuity diff --git a/src/eDFT/correlation_energy.f90 b/src/eDFT/unrestricted_correlation_energy.f90 similarity index 58% rename from src/eDFT/correlation_energy.f90 rename to src/eDFT/unrestricted_correlation_energy.f90 index 9326ba9..9a49cde 100644 --- a/src/eDFT/correlation_energy.f90 +++ b/src/eDFT/unrestricted_correlation_energy.f90 @@ -1,6 +1,6 @@ -subroutine correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) +subroutine unrestricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) -! Compute the correlation energy +! Compute the unrestricted version of the correlation energy implicit none include 'parameters.h' @@ -14,7 +14,7 @@ subroutine correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rho(nGrid,nspin) - double precision,intent(in) :: drho(3,nGrid,nspin) + double precision,intent(in) :: drho(ncart,nGrid,nspin) ! Local variables @@ -38,13 +38,13 @@ subroutine correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) case(1) - call lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),Ec(:)) + call unrestricted_lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),Ec(:)) ! GGA functionals case(2) - call gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),Ec(:)) + call unrestricted_gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),Ec(:)) ! Hybrid functionals @@ -52,8 +52,8 @@ subroutine correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) aC = 0.81d0 - call lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),EcLDA(:)) - call gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),EcGGA(:)) + call unrestricted_lda_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),EcLDA(:)) + call unrestricted_gga_correlation_energy(DFA,nEns,wEns(:),nGrid,weight(:),rho(:,:),drho(:,:,:),EcGGA(:)) Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:)) @@ -65,4 +65,4 @@ subroutine correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) end select -end subroutine correlation_energy +end subroutine unrestricted_correlation_energy diff --git a/src/eDFT/correlation_individual_energy.f90 b/src/eDFT/unrestricted_correlation_individual_energy.f90 similarity index 78% rename from src/eDFT/correlation_individual_energy.f90 rename to src/eDFT/unrestricted_correlation_individual_energy.f90 index 4fd9fbb..a0431f4 100644 --- a/src/eDFT/correlation_individual_energy.f90 +++ b/src/eDFT/unrestricted_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine correlation_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec) +subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ec) ! Compute the correlation energy of individual states @@ -9,6 +9,7 @@ subroutine correlation_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,dr 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 @@ -40,7 +41,7 @@ subroutine correlation_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,dr case(1) - call lda_correlation_individual_energy(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),rho(:,:),Ec(:)) + call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ec) ! GGA functionals @@ -66,4 +67,4 @@ subroutine correlation_individual_energy(rung,DFA,nEns,wEns,nGrid,weight,rhow,dr end select -end subroutine correlation_individual_energy +end subroutine unrestricted_correlation_individual_energy diff --git a/src/eDFT/correlation_potential.f90 b/src/eDFT/unrestricted_correlation_potential.f90 similarity index 64% rename from src/eDFT/correlation_potential.f90 rename to src/eDFT/unrestricted_correlation_potential.f90 index d79f47b..e847b28 100644 --- a/src/eDFT/correlation_potential.f90 +++ b/src/eDFT/unrestricted_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) +subroutine unrestricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Compute the correlation potential @@ -35,41 +35,41 @@ subroutine correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho ! Hartree calculation - case(0) + case(0) Fc(:,:,:) = 0d0 ! LDA functionals - case(1) + case(1) - call lda_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) + call unrestricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) ! GGA functionals - case(2) + case(2) - call gga_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rho(:,:),drho(:,:,:),Fc(:,:,:)) + call unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Hybrid functionals - case(4) + case(4) allocate(FcLDA(nBas,nBas,nspin),FcGGA(nBas,nBas,nspin)) aC = 0.81d0 - call lda_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),rho(:,:),FcLDA(:,:,:)) - call gga_correlation_potential(DFA,nEns,wEns(:),nGrid,weight(:),nBas,AO(:,:),dAO(:,:,:),rho(:,:),drho(:,:,:),FcGGA(:,:,:)) + call unrestricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,FcLDA) + call unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FcGGA) - Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:)) + Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:)) ! Hartree-Fock calculation - case(666) + case(666) Fc(:,:,:) = 0d0 end select -end subroutine correlation_potential +end subroutine unrestricted_correlation_potential diff --git a/src/eDFT/gga_correlation_energy.f90 b/src/eDFT/unrestricted_gga_correlation_energy.f90 similarity index 82% rename from src/eDFT/gga_correlation_energy.f90 rename to src/eDFT/unrestricted_gga_correlation_energy.f90 index 3a4c10d..1ba444a 100644 --- a/src/eDFT/gga_correlation_energy.f90 +++ b/src/eDFT/unrestricted_gga_correlation_energy.f90 @@ -1,6 +1,6 @@ -subroutine gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) +subroutine unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) -! Compute GGA correlation energy +! Compute unrstricted GGA correlation energy implicit none include 'parameters.h' @@ -39,4 +39,4 @@ subroutine gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) enddo -end subroutine gga_correlation_energy +end subroutine unrestricted_gga_correlation_energy diff --git a/src/eDFT/gga_correlation_potential.f90 b/src/eDFT/unrestricted_gga_correlation_potential.f90 similarity index 77% rename from src/eDFT/gga_correlation_potential.f90 rename to src/eDFT/unrestricted_gga_correlation_potential.f90 index 1839e07..0f49fb7 100644 --- a/src/eDFT/gga_correlation_potential.f90 +++ b/src/eDFT/unrestricted_gga_correlation_potential.f90 @@ -1,6 +1,6 @@ -subroutine gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) +subroutine unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) -! Compute GGA correlation potential +! Compute unrestricted GGA correlation potential implicit none include 'parameters.h' @@ -28,4 +28,4 @@ subroutine gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho, ! Compute GGA correlation matrix in the AO basis -end subroutine gga_correlation_potential +end subroutine unrestricted_gga_correlation_potential diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 new file mode 100644 index 0000000..42e4d01 --- /dev/null +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -0,0 +1,240 @@ +subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,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,E,Om) + +! Compute unrestricted 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) + 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(nspin),nV(nspin) + 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,nspin) !!!!! + double precision,intent(in) :: Pw(nBas,nBas,nspin) !!!!! + double precision,intent(in) :: rhow(nGrid,nspin) + double precision,intent(in) :: drhow(ncart,nGrid,nspin) + + double precision,intent(in) :: P(nBas,nBas,nspin,nEns) + double precision,intent(in) :: rho(nGrid,nspin,nEns) + double precision,intent(in) :: drho(ncart,nGrid,nspin,nEns) + + double precision,intent(in) :: J(nBas,nBas,nspin) + double precision,intent(in) :: Fx(nBas,nBas,nspin) + double precision,intent(in) :: FxHF(nBas,nBas,nspin) + double precision,intent(in) :: Fc(nBas,nBas,nspin) + double precision :: Ew + + +! Local variables + + double precision :: ET(nspin,nEns) + double precision :: EV(nspin,nEns) + double precision :: EJ(nsp,nEns) + double precision :: Ex(nspin,nEns) + double precision :: Ec(nsp,nEns) + double precision :: Exc(nEns) !!!!! + double precision :: Eaux(nspin,nEns) !!!!! + + double precision :: ExDD(nspin,nEns) !!!!! + double precision :: EcDD(nsp,nEns) + double precision :: ExcDD(nspin,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 :: ispin,iEns + +! Output variables + + double precision,intent(out) :: E(nEns) + double precision,intent(out) :: Om(nEns) + +!------------------------------------------------------------------------ +! Kinetic energy +!------------------------------------------------------------------------ + + do ispin=1,nspin + do iEns=1,nEns + ET(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),T(:,:))) + end do + end do + +!------------------------------------------------------------------------ +! Potential energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + do ispin=1,nspin + EV(ispin,iEns) = trace_matrix(nBas,matmul(P(:,:,ispin,iEns),V(:,:))) + end do + end do + +!------------------------------------------------------------------------ +! Individual Hartree energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + + do ispin=1,nspin + call hartree_coulomb(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) + end do + + EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) + + EJ(2,iEns) = 2.0d0*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) + + EJ(3,iEns) = trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,2))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) + + end do + +!------------------------------------------------------------------------ +! Individual exchange energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + do ispin=1,nspin + call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nGrid,weight,nBas,ERI, & + Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin), & + rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns)) + end do + end do + +!------------------------------------------------------------------------ +! Individual correlation energy +!------------------------------------------------------------------------ + + do iEns=1,nEns + call unrestricted_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 unrestricted_auxiliary_energy(nBas,nEns,nO(:),eps(:,:),Eaux(:,:)) + +!------------------------------------------------------------------------ +! Compute derivative discontinuities +!------------------------------------------------------------------------ + do iEns=1,nEns + do ispin=1,nspin + + call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:), & + rhow(:,ispin),drhow(:,:,ispin),ExDD(ispin,iEns)) + + call restricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:), & + rhow(:,ispin),drhow(:,:,ispin),EcDD(:,iEns)) +! EcDD(ispin,:) = 0.d0 !!!!!!!!! + ExcDD(ispin,iEns) = ExDD(ispin,iEns) +sum(EcDD(:,iEns)) + end do + end do + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + +!------------------------------------------------------------------------ +! Exchange energy +!------------------------------------------------------------------------ + +! do iEns=1,nEns +! do ispin=1,nspin + +! call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,ispin,iEns),ERI(:,:,:,:), & +! AO(:,:),dAO(:,:,:),rho(:,ispin,iEns),drho(:,:,ispin,iEns),Fx(:,:,ispin),FxHF(:,:,ispin)) +! call exchange_energy(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),nBas,P(:,:,ispin,iEns),FxHF(:,:,ispin), & +! rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns)) + +! end do +! end do + +!------------------------------------------------------------------------ +! Correlation energy +!------------------------------------------------------------------------ + +! do iEns=1,nEns + +! call correlation_individual_energy(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:), & +! rho(:,:,iEns),drho(:,:,:,iEns),Ec(:,iEns)) + +! end do + +!------------------------------------------------------------------------ +! Compute derivative discontinuities +!------------------------------------------------------------------------ + +! call correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:,:),drhow(:,:,:),EcDD(:,:)) + +!------------------------------------------------------------------------ +! Total energy +!------------------------------------------------------------------------ + +! do iEns=1,nEns +! E(iEns) = ENuc + sum(ET(:,iEns)) + sum(EV(:,iEns)) + sum(EJ(:,iEns)) & +! + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + sum(EcDD(:,iEns)) +! end do + + do iEns=1,nEns + Exc(iEns) = sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + E(iEns) =sum( ET(:,iEns)) + sum( EV(:,iEns)) + sum(EJ(:,iEns)) & + + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) + sum(ExcDD(:,iEns)) + end do + + + + + +!------------------------------------------------------------------------ +! Excitation energies +!------------------------------------------------------------------------ + + do iEns=1,nEns + Om(iEns) = E(iEns) - E(1) + + Omx(iEns) = sum(Ex(:,iEns)) -sum(Ex(:,1)) + Omc(iEns) = sum(Ec(:,iEns)) -sum(Ec(:,1)) + Omxc(iEns) = Exc(iEns) - Exc(1) + + Omaux(iEns) = sum(Eaux(:,iEns)) -sum(Eaux(:,1)) + + OmxDD(iEns) = sum(ExDD(:,iEns)) - sum(ExDD(:,1)) + OmcDD(iEns) = sum(EcDD(:,iEns)) - sum(EcDD(:,1)) + OmxcDD(iEns) = sum(ExcDD(:,iEns)) - sum(ExcDD(:,1)) + + + end do + +!------------------------------------------------------------------------ +! Dump results +!------------------------------------------------------------------------ + + call print_unrestricted_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 unrestricted_individual_energy diff --git a/src/eDFT/lda_correlation_derivative_discontinuity.f90 b/src/eDFT/unrestricted_lda_correlation_derivative_discontinuity.f90 similarity index 69% rename from src/eDFT/lda_correlation_derivative_discontinuity.f90 rename to src/eDFT/unrestricted_lda_correlation_derivative_discontinuity.f90 index bc5dd88..f7c2f83 100644 --- a/src/eDFT/lda_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/unrestricted_lda_correlation_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) +subroutine unrestricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! Compute the correlation LDA part of the derivative discontinuity @@ -28,22 +28,16 @@ subroutine lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,r ! Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678 - case ('W38') + case ('RW38') Ec(:,:) = 0d0 ! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200 - case ('VWN5') + case ('RVWN5') Ec(:,:) = 0d0 -! Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation) - - case ('MFL20') - - call MFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:,:),Ec(:,:)) - case default call print_warning('!!! LDA correlation functional not available !!!') @@ -51,4 +45,4 @@ subroutine lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,r end select -end subroutine lda_correlation_derivative_discontinuity +end subroutine unrestricted_lda_correlation_derivative_discontinuity diff --git a/src/eDFT/lda_correlation_energy.f90 b/src/eDFT/unrestricted_lda_correlation_energy.f90 similarity index 64% rename from src/eDFT/lda_correlation_energy.f90 rename to src/eDFT/unrestricted_lda_correlation_energy.f90 index 31008fc..abefb8d 100644 --- a/src/eDFT/lda_correlation_energy.f90 +++ b/src/eDFT/unrestricted_lda_correlation_energy.f90 @@ -1,6 +1,6 @@ -subroutine lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) +subroutine unrestricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) -! Select LDA correlation functional +! Select the unrestrited version of the LDA correlation functional implicit none include 'parameters.h' @@ -28,21 +28,21 @@ subroutine lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) Ec(:) = 0d0 - case ('W38') + case ('UW38') - call W38_lda_correlation_energy(nGrid,weight(:),rho(:,:),Ec(:)) + call UW38_lda_correlation_energy(nGrid,weight,rho,Ec) ! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200 - case ('VWN5') + case ('UVWN5') - call VWN5_lda_correlation_energy(nGrid,weight(:),rho(:,:),Ec(:)) + call UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec) ! Chachiyo's LDA correlation functional: Chachiyo, JCP 145 (2016) 021101 - case ('C16') + case ('UC16') - call C16_lda_correlation_energy(nGrid,weight(:),rho(:,:),Ec(:)) + call UC16_lda_correlation_energy(nGrid,weight,rho,Ec) case default @@ -51,4 +51,4 @@ subroutine lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) end select -end subroutine lda_correlation_energy +end subroutine unrestricted_lda_correlation_energy diff --git a/src/eDFT/lda_correlation_individual_energy.f90 b/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 similarity index 53% rename from src/eDFT/lda_correlation_individual_energy.f90 rename to src/eDFT/unrestricted_lda_correlation_individual_energy.f90 index 4703864..2183a07 100644 --- a/src/eDFT/lda_correlation_individual_energy.f90 +++ b/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho,Ec) +subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ec) ! Compute LDA correlation energy for individual states @@ -7,6 +7,7 @@ subroutine lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho ! Input variables + logical,intent(in) :: LDA_centered character(len=12),intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) @@ -23,23 +24,11 @@ subroutine lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho select case (DFA) -! Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678 - - case ('W38') - - call W38_lda_correlation_individual_energy(nGrid,weight(:),rhow(:,:),rho(:,:),Ec(:)) - ! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200 - case ('VWN5') + case ('UVWN5') - call VWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:,:),rho(:,:),Ec(:)) - -! Loos-Fromager weight-dependent correlation functional: Loos and Fromager (in preparation) - - case ('MFL20') - - call MFL20_lda_correlation_individual_energy(nEns,wEns,nGrid,weight(:),rhow(:,:),rho(:,:),Ec(:)) + call UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) case default @@ -48,4 +37,4 @@ subroutine lda_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,rho end select -end subroutine lda_correlation_individual_energy +end subroutine unrestricted_lda_correlation_individual_energy diff --git a/src/eDFT/lda_correlation_potential.f90 b/src/eDFT/unrestricted_lda_correlation_potential.f90 similarity index 69% rename from src/eDFT/lda_correlation_potential.f90 rename to src/eDFT/unrestricted_lda_correlation_potential.f90 index b2d8d2a..fc0fd62 100644 --- a/src/eDFT/lda_correlation_potential.f90 +++ b/src/eDFT/unrestricted_lda_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) +subroutine unrestricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) ! Select LDA correlation potential @@ -32,21 +32,21 @@ include 'parameters.h' ! Wigner's LDA correlation functional: Wigner, Trans. Faraday Soc. 34 (1938) 678 - case ('W38') + case ('UW38') - call W38_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) + call UW38_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) ! Vosko, Wilk and Nusair's functional V: Can. J. Phys. 58 (1980) 1200 - case ('VWN5') + case ('UVWN5') - call VWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) + call UVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) ! Chachiyo's LDA correlation functional: Chachiyo, JCP 145 (2016) 021101 - case ('C16') + case ('UC16') - call C16_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) + call UC16_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) case default @@ -55,4 +55,4 @@ include 'parameters.h' end select -end subroutine lda_correlation_potential +end subroutine unrestricted_lda_correlation_potential