diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 index c304dab..3f657a0 100644 --- a/src/eDFT/GOK_RKS.f90 +++ b/src/eDFT/GOK_RKS.f90 @@ -236,8 +236,8 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_ ! Compute exchange potential - call exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,Pw(:,:),ERI(:,:,:,:), & - AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fx(:,:),FxHF(:,:)) + call restricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,Pw(:,:), & + ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fx(:,:),FxHF(:,:)) ! Compute correlation potential @@ -294,8 +294,8 @@ subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_ ! Exchange energy - call exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & - Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex) + call restricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & + Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex) ! Correlation energy diff --git a/src/eDFT/GOK_UKS.f90 b/src/eDFT/GOK_UKS.f90 index c0e9145..0ee5c86 100644 --- a/src/eDFT/GOK_UKS.f90 +++ b/src/eDFT/GOK_UKS.f90 @@ -248,8 +248,9 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w ! Compute exchange potential do ispin=1,nspin - call exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,nBas,Pw(:,:,ispin),ERI(:,:,:,:), & - AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin),Fx(:,:,ispin),FxHF(:,:,ispin)) + call unrestricted_exchange_potential(x_rung,x_DFA,nEns,wEns(:),nGrid,weight(:),aCC_w1,aCC_w2,nBas,Pw(:,:,ispin), & + ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin),Fx(:,:,ispin), & + FxHF(:,:,ispin)) end do ! Compute correlation potential @@ -306,8 +307,8 @@ subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w ! Exchange energy do ispin=1,nspin - call exchange_energy(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & - Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin)) + call unrestricted_exchange_energy(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & + Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin)) end do ! Correlation energy diff --git a/src/eDFT/B88_gga_exchange_energy.f90 b/src/eDFT/UB88_gga_exchange_energy.f90 similarity index 90% rename from src/eDFT/B88_gga_exchange_energy.f90 rename to src/eDFT/UB88_gga_exchange_energy.f90 index d47fa7a..2f6b14d 100644 --- a/src/eDFT/B88_gga_exchange_energy.f90 +++ b/src/eDFT/UB88_gga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) +subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Compute Becke's 88 GGA exchange energy @@ -47,4 +47,4 @@ subroutine B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) end do -end subroutine B88_gga_exchange_energy +end subroutine UB88_gga_exchange_energy diff --git a/src/eDFT/B88_gga_exchange_potential.f90 b/src/eDFT/UB88_gga_exchange_potential.f90 similarity index 93% rename from src/eDFT/B88_gga_exchange_potential.f90 rename to src/eDFT/UB88_gga_exchange_potential.f90 index 5b8902e..14244ff 100644 --- a/src/eDFT/B88_gga_exchange_potential.f90 +++ b/src/eDFT/UB88_gga_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Compute Becke's GGA exchange potential @@ -62,4 +62,4 @@ subroutine B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) end do end do -end subroutine B88_gga_exchange_potential +end subroutine UB88_gga_exchange_potential diff --git a/src/eDFT/G96_gga_exchange_energy.f90 b/src/eDFT/UG96_gga_exchange_energy.f90 similarity index 89% rename from src/eDFT/G96_gga_exchange_energy.f90 rename to src/eDFT/UG96_gga_exchange_energy.f90 index 90aa536..48f0a8c 100644 --- a/src/eDFT/G96_gga_exchange_energy.f90 +++ b/src/eDFT/UG96_gga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) +subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Compute Gill's 96 GGA exchange energy @@ -46,4 +46,4 @@ subroutine G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) end do -end subroutine G96_gga_exchange_energy +end subroutine UG96_gga_exchange_energy diff --git a/src/eDFT/G96_gga_exchange_potential.f90 b/src/eDFT/UG96_gga_exchange_potential.f90 similarity index 93% rename from src/eDFT/G96_gga_exchange_potential.f90 rename to src/eDFT/UG96_gga_exchange_potential.f90 index 2243865..2dd097b 100644 --- a/src/eDFT/G96_gga_exchange_potential.f90 +++ b/src/eDFT/UG96_gga_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Compute Gill's GGA exchange poential @@ -63,4 +63,4 @@ subroutine G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) enddo enddo -end subroutine G96_gga_exchange_potential +end subroutine UG96_gga_exchange_potential diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index e3dd564..1ba8e6c 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -238,9 +238,9 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig ! Compute exchange potential do ispin=1,nspin - call exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & - Pw(:,:,ispin),ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin), & - Fx(:,:,ispin),FxHF(:,:,ispin)) + call unrestricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & + Pw(:,:,ispin),ERI(:,:,:,:),AO(:,:),dAO(:,:,:),rhow(:,ispin),drhow(:,:,ispin), & + Fx(:,:,ispin),FxHF(:,:,ispin)) end do ! Compute correlation potential @@ -317,8 +317,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig ! Exchange energy do ispin=1,nspin - call exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & - Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin)) + call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & + Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin),Ex(ispin)) end do ! Correlation energy diff --git a/src/eDFT/huckel_guess.f90 b/src/eDFT/huckel_guess.f90 index c15b0fa..1bdf79d 100644 --- a/src/eDFT/huckel_guess.f90 +++ b/src/eDFT/huckel_guess.f90 @@ -36,7 +36,7 @@ subroutine huckel_guess(nBas,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) c(:,:) = matmul(X(:,:),cp(:,:)) call hartree_coulomb(nBas,P,ERI,J) - call fock_exchange_potential(nBas,P,ERI,K) + call restricted_fock_exchange_potential(nBas,P,ERI,K) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) diff --git a/src/eDFT/exchange_derivative_discontinuity.f90 b/src/eDFT/restricted_exchange_derivative_discontinuity.f90 similarity index 70% rename from src/eDFT/exchange_derivative_discontinuity.f90 rename to src/eDFT/restricted_exchange_derivative_discontinuity.f90 index 9f822de..5b7441b 100644 --- a/src/eDFT/exchange_derivative_discontinuity.f90 +++ b/src/eDFT/restricted_exchange_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD) +subroutine restricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD) ! Compute the exchange part of the derivative discontinuity @@ -37,13 +37,13 @@ subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nG case(1) - call lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:)) + call restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:)) ! GGA functionals case(2) - call gga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) + call restricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) ! Hybrid functionals @@ -60,4 +60,4 @@ subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nG end select -end subroutine exchange_derivative_discontinuity +end subroutine restricted_exchange_derivative_discontinuity diff --git a/src/eDFT/exchange_energy.f90 b/src/eDFT/restricted_exchange_energy.f90 similarity index 67% rename from src/eDFT/exchange_energy.f90 rename to src/eDFT/restricted_exchange_energy.f90 index 269819f..c9e1d67 100644 --- a/src/eDFT/exchange_energy.f90 +++ b/src/eDFT/restricted_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF,rho,drho,Ex) +subroutine restricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF,rho,drho,Ex) ! Compute the exchange energy @@ -43,7 +43,7 @@ subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,w case(1) - call lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,ExLDA) + call restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,ExLDA) Ex = ExLDA @@ -51,7 +51,7 @@ subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,w case(2) - call gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA) + call restricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA) Ex = ExGGA @@ -63,9 +63,9 @@ subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,w aX = 0.72d0 aC = 0.81d0 - call lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA) - call gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA) - call fock_exchange_energy(nBas,P,FxHF,ExHF) + call restricted_lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA) + call restricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA) + call restricted_fock_exchange_energy(nBas,P,FxHF,ExHF) Ex = ExLDA & + cX*(ExHF - ExLDA) & @@ -75,10 +75,10 @@ subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,w case(666) - call fock_exchange_energy(nBas,P,FxHF,ExHF) + call restricted_fock_exchange_energy(nBas,P,FxHF,ExHF) Ex = ExHF end select -end subroutine exchange_energy +end subroutine restricted_exchange_energy diff --git a/src/eDFT/exchange_individual_energy.f90 b/src/eDFT/restricted_exchange_individual_energy.f90 similarity index 72% rename from src/eDFT/exchange_individual_energy.f90 rename to src/eDFT/restricted_exchange_individual_energy.f90 index 9ee0bbc..a4aa93f 100644 --- a/src/eDFT/exchange_individual_energy.f90 +++ b/src/eDFT/restricted_exchange_individual_energy.f90 @@ -1,5 +1,5 @@ -subroutine exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & - ERI,Pw,P,rhow,drhow,rho,drho,Ex) +subroutine restricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & + ERI,Pw,P,rhow,drhow,rho,drho,Ex) ! Compute the exchange individual energy @@ -48,7 +48,7 @@ subroutine exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC case(1) - call lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,ExLDA) + call restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,ExLDA) Ex = ExLDA @@ -56,7 +56,7 @@ subroutine exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC case(2) - call gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,ExGGA) + call restricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,ExGGA) Ex = ExGGA @@ -71,10 +71,10 @@ subroutine exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC case(666) - call fock_exchange_individual_energy(nBas,Pw,P,ERI,ExHF) + call restricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,ExHF) Ex = ExHF end select -end subroutine exchange_individual_energy +end subroutine restricted_exchange_individual_energy diff --git a/src/eDFT/exchange_potential.f90 b/src/eDFT/restricted_exchange_potential.f90 similarity index 68% rename from src/eDFT/exchange_potential.f90 rename to src/eDFT/restricted_exchange_potential.f90 index b7fd6cd..81fcea0 100644 --- a/src/eDFT/exchange_potential.f90 +++ b/src/eDFT/restricted_exchange_potential.f90 @@ -1,4 +1,5 @@ -subroutine exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,ERI,AO,dAO,rho,drho,Fx,FxHF) +subroutine restricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P, & + ERI,AO,dAO,rho,drho,Fx,FxHF) ! Compute the exchange potential @@ -47,13 +48,13 @@ subroutine exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGri case(1) - call lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx) + call restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx) ! GGA functionals case(2) - call gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call restricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Hybrid functionals @@ -64,9 +65,9 @@ subroutine exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGri cX = 0.20d0 aX = 0.72d0 - call lda_exchange_potential(DFA,nGrid,weight,nBas,AO,rho,FxLDA) - call gga_exchange_potential(DFA,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) - call fock_exchange_potential(nBas,P,ERI,FxHF) + call restricted_lda_exchange_potential(DFA,nGrid,weight,nBas,AO,rho,FxLDA) + call restricted_gga_exchange_potential(DFA,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) + call restricted_fock_exchange_potential(nBas,P,ERI,FxHF) Fx(:,:) = FxLDA(:,:) & + cX*(FxHF(:,:) - FxLDA(:,:)) & @@ -76,10 +77,10 @@ subroutine exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGri case(666) - call fock_exchange_potential(nBas,P,ERI,FxHF) + call restricted_fock_exchange_potential(nBas,P,ERI,FxHF) Fx(:,:) = FxHF(:,:) end select -end subroutine exchange_potential +end subroutine restricted_exchange_potential diff --git a/src/eDFT/restricted_fock_exchange_energy.f90 b/src/eDFT/restricted_fock_exchange_energy.f90 new file mode 100644 index 0000000..1739efc --- /dev/null +++ b/src/eDFT/restricted_fock_exchange_energy.f90 @@ -0,0 +1,25 @@ +subroutine restricted_fock_exchange_energy(nBas,P,Fx,Ex) + +! Compute the (exact) Fock exchange energy + + implicit none + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: Fx(nBas,nBas) + +! Local variables + + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: Ex + +! Compute HF exchange energy + + Ex = 0.5d0*trace_matrix(nBas,matmul(P,Fx)) + +end subroutine restricted_fock_exchange_energy diff --git a/src/eDFT/fock_exchange_individual_energy.f90 b/src/eDFT/restricted_fock_exchange_individual_energy.f90 similarity index 74% rename from src/eDFT/fock_exchange_individual_energy.f90 rename to src/eDFT/restricted_fock_exchange_individual_energy.f90 index 9b6d3ae..c06eac5 100644 --- a/src/eDFT/fock_exchange_individual_energy.f90 +++ b/src/eDFT/restricted_fock_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex) +subroutine restricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex) ! Compute the Fock exchange potential @@ -24,8 +24,8 @@ subroutine fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex) allocate(Fx(nBas,nBas)) - call fock_exchange_potential(nBas,Pw(:,:),ERI(:,:,:,:),Fx(:,:)) + call restricted_fock_exchange_potential(nBas,Pw(:,:),ERI(:,:,:,:),Fx(:,:)) Ex = trace_matrix(nBas,matmul(P(:,:),Fx(:,:))) & - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),Fx(:,:))) -end subroutine fock_exchange_individual_energy +end subroutine restricted_fock_exchange_individual_energy diff --git a/src/eDFT/fock_exchange_potential.f90 b/src/eDFT/restricted_fock_exchange_potential.f90 similarity index 84% rename from src/eDFT/fock_exchange_potential.f90 rename to src/eDFT/restricted_fock_exchange_potential.f90 index 4423160..c43115c 100644 --- a/src/eDFT/fock_exchange_potential.f90 +++ b/src/eDFT/restricted_fock_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine fock_exchange_potential(nBas,P,ERI,Fx) +subroutine restricted_fock_exchange_potential(nBas,P,ERI,Fx) ! Compute the Fock exchange potential @@ -33,4 +33,4 @@ subroutine fock_exchange_potential(nBas,P,ERI,Fx) Fx(:,:) = 0.5d0*Fx(:,:) -end subroutine fock_exchange_potential +end subroutine restricted_fock_exchange_potential diff --git a/src/eDFT/gga_exchange_derivative_discontinuity.f90 b/src/eDFT/restricted_gga_exchange_derivative_discontinuity.f90 similarity index 80% rename from src/eDFT/gga_exchange_derivative_discontinuity.f90 rename to src/eDFT/restricted_gga_exchange_derivative_discontinuity.f90 index eeedc32..755c849 100644 --- a/src/eDFT/gga_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/restricted_gga_exchange_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD) +subroutine restricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD) ! Compute the exchange GGA part of the derivative discontinuity @@ -26,7 +26,7 @@ subroutine gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow select case (DFA) - case ('RB88') + case ('B88') ExDD(:) = 0d0 @@ -37,4 +37,4 @@ subroutine gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow end select -end subroutine gga_exchange_derivative_discontinuity +end subroutine restricted_gga_exchange_derivative_discontinuity diff --git a/src/eDFT/gga_exchange_energy.f90 b/src/eDFT/restricted_gga_exchange_energy.f90 similarity index 72% rename from src/eDFT/gga_exchange_energy.f90 rename to src/eDFT/restricted_gga_exchange_energy.f90 index e81a1af..4752444 100644 --- a/src/eDFT/gga_exchange_energy.f90 +++ b/src/eDFT/restricted_gga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) +subroutine restricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) ! Select GGA exchange functional for energy calculation @@ -22,17 +22,9 @@ subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) select case (DFA) - case ('G96') - - call G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) - - case ('RB88') - - call RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) - case ('B88') - call B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + call RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) case default @@ -41,4 +33,4 @@ subroutine gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) end select -end subroutine gga_exchange_energy +end subroutine restricted_gga_exchange_energy diff --git a/src/eDFT/gga_exchange_individual_energy.f90 b/src/eDFT/restricted_gga_exchange_individual_energy.f90 similarity index 83% rename from src/eDFT/gga_exchange_individual_energy.f90 rename to src/eDFT/restricted_gga_exchange_individual_energy.f90 index d233066..5cab9c6 100644 --- a/src/eDFT/gga_exchange_individual_energy.f90 +++ b/src/eDFT/restricted_gga_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ex) +subroutine restricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ex) ! Compute GGA exchange energy for individual states @@ -25,7 +25,7 @@ subroutine gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow, select case (DFA) - case ('RB88') + case ('B88') call RB88_gga_exchange_individual_energy(nGrid,weight(:),rhow(:),drhow(:,:),rho(:),drho(:,:),Ex) @@ -36,4 +36,4 @@ subroutine gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow, end select -end subroutine gga_exchange_individual_energy +end subroutine restricted_gga_exchange_individual_energy diff --git a/src/eDFT/gga_exchange_potential.f90 b/src/eDFT/restricted_gga_exchange_potential.f90 similarity index 73% rename from src/eDFT/gga_exchange_potential.f90 rename to src/eDFT/restricted_gga_exchange_potential.f90 index 2e9867d..83b8bde 100644 --- a/src/eDFT/gga_exchange_potential.f90 +++ b/src/eDFT/restricted_gga_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine restricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Select GGA exchange functional for potential calculation @@ -26,17 +26,9 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh select case (DFA) - case ('G96') - - call G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) - - case ('RB88') - - call RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) - case ('B88') - call B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) case default @@ -45,4 +37,4 @@ subroutine gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drh end select -end subroutine gga_exchange_potential +end subroutine restricted_gga_exchange_potential diff --git a/src/eDFT/restricted_individual_energy.f90 b/src/eDFT/restricted_individual_energy.f90 index 173376e..67968f0 100644 --- a/src/eDFT/restricted_individual_energy.f90 +++ b/src/eDFT/restricted_individual_energy.f90 @@ -91,8 +91,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(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas,ERI(:,:,:,:), & - Pw(:,:),P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns)) + call restricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & + ERI(:,:,:,:),Pw(:,:),P(:,:,iEns),rhow(:),drhow(:,:),rho(:,iEns),drho(:,:,iEns),Ex(iEns)) end do !------------------------------------------------------------------------ @@ -114,7 +114,8 @@ subroutine restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,n ! Compute derivative discontinuities !------------------------------------------------------------------------ - call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) + call restricted_exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:), & + rhow(:),drhow(:,:),ExDD(:)) call restricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),EcDD(:)) diff --git a/src/eDFT/lda_exchange_derivative_discontinuity.f90 b/src/eDFT/restricted_lda_exchange_derivative_discontinuity.f90 similarity index 72% rename from src/eDFT/lda_exchange_derivative_discontinuity.f90 rename to src/eDFT/restricted_lda_exchange_derivative_discontinuity.f90 index 45fd060..10eebb9 100644 --- a/src/eDFT/lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/restricted_lda_exchange_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD) +subroutine restricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD) ! Compute the exchange LDA part of the derivative discontinuity @@ -28,26 +28,18 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGr select case (DFA) - case ('US51') + case ('S51') ExDD(:) = 0d0 - case ('RS51') - - ExDD(:) = 0d0 - - case ('RMFL20') + case ('MFL20') call RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight(:),rhow(:),ExDD(:)) - case ('RCC') + case ('CC') call RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:)) - case ('UCC') - - call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:)) - case default call print_warning('!!! LDA exchange derivative discontinuity not available !!!') @@ -55,4 +47,4 @@ subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGr end select -end subroutine lda_exchange_derivative_discontinuity +end subroutine restricted_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/lda_exchange_energy.f90 b/src/eDFT/restricted_lda_exchange_energy.f90 similarity index 66% rename from src/eDFT/lda_exchange_energy.f90 rename to src/eDFT/restricted_lda_exchange_energy.f90 index b7141c6..519c0a8 100644 --- a/src/eDFT/lda_exchange_energy.f90 +++ b/src/eDFT/restricted_lda_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex) +subroutine restricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex) ! Select LDA exchange functional @@ -25,31 +25,23 @@ subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,we select case (DFA) - case ('US51') - - call US51_lda_exchange_energy(nGrid,weight,rho,Ex) - - case ('RS51') + case ('S51') call RS51_lda_exchange_energy(nGrid,weight,rho,Ex) - case ('RMFL20') + case ('MFL20') call RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) - case ('RCC') + case ('CC') call RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex) - case ('UCC') - - call UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex) - case default - call print_warning('!!! LDA exchange functional not available !!!') + call print_warning('!!! LDA restricted exchange functional not available !!!') stop end select -end subroutine lda_exchange_energy +end subroutine restricted_lda_exchange_energy diff --git a/src/eDFT/lda_exchange_individual_energy.f90 b/src/eDFT/restricted_lda_exchange_individual_energy.f90 similarity index 71% rename from src/eDFT/lda_exchange_individual_energy.f90 rename to src/eDFT/restricted_lda_exchange_individual_energy.f90 index c1bf10b..fcdde39 100644 --- a/src/eDFT/lda_exchange_individual_energy.f90 +++ b/src/eDFT/restricted_lda_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex) +subroutine restricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex) ! Compute LDA exchange energy for individual states @@ -26,26 +26,18 @@ subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_ select case (DFA) - case ('US51') - - call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex) - - case ('RS51') + case ('S51') call RS51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex) - case ('RMFL20') + case ('MFL20') call RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ex) - case ('RCC') + case ('CC') call RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex) - case ('UCC') - - call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex) - case default call print_warning('!!! LDA exchange individual energy not available !!!') @@ -53,4 +45,4 @@ subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_ end select -end subroutine lda_exchange_individual_energy +end subroutine restricted_lda_exchange_individual_energy diff --git a/src/eDFT/restricted_lda_exchange_potential.f90 b/src/eDFT/restricted_lda_exchange_potential.f90 new file mode 100644 index 0000000..9f1b5bf --- /dev/null +++ b/src/eDFT/restricted_lda_exchange_potential.f90 @@ -0,0 +1,50 @@ +subroutine restricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx) + +! Select LDA correlation potential + + implicit none + + include 'parameters.h' + +! Input variables + + logical,intent(in) :: LDA_centered + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Select exchange functional + + select case (DFA) + + case ('S51') + + call RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) + + case ('MFL20') + + call RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) + + case ('CC') + + call RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx) + + case default + + call print_warning('!!! LDA exchange functional not available !!!') + stop + + end select + +end subroutine restricted_lda_exchange_potential diff --git a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 new file mode 100644 index 0000000..04aae3a --- /dev/null +++ b/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 @@ -0,0 +1,63 @@ +subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,drhow,ExDD) + +! Compute the exchange part of the derivative discontinuity + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: rung + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + +! Local variables + + +! Output variables + + double precision,intent(out) :: ExDD(nEns) + + select case (rung) + +! Hartree calculation + + case(0) + + ExDD(:) = 0d0 + +! LDA functionals + + case(1) + + call unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:)) + +! GGA functionals + + case(2) + + call unrestricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) + +! Hybrid functionals + + case(4) + + call print_warning('!!! exchange part of derivative discontinuity NYI for hybrids !!!') + stop + +! Hartree-Fock calculation + + case(666) + + ExDD(:) = 0d0 + + end select + +end subroutine unrestricted_exchange_derivative_discontinuity diff --git a/src/eDFT/unrestricted_exchange_energy.f90 b/src/eDFT/unrestricted_exchange_energy.f90 new file mode 100644 index 0000000..f31c2aa --- /dev/null +++ b/src/eDFT/unrestricted_exchange_energy.f90 @@ -0,0 +1,84 @@ +subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P,FxHF,rho,drho,Ex) + +! Compute the exchange energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: rung + character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: FxHF(nBas,nBas) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + +! Local variables + + double precision :: ExLDA,ExGGA,ExHF + double precision :: cX,aX,aC + +! Output variables + + double precision,intent(out) :: Ex + + select case (rung) + +! Hartree calculation + + case(0) + + Ex = 0d0 + +! LDA functionals + + case(1) + + call unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,ExLDA) + + Ex = ExLDA + +! GGA functionals + + case(2) + + call unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA) + + Ex = ExGGA + +! Hybrid functionals + + case(4) + + cX = 0.20d0 + aX = 0.72d0 + aC = 0.81d0 + + call unrestricted_lda_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,ExLDA) + call unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,ExGGA) + call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF) + + Ex = ExLDA & + + cX*(ExHF - ExLDA) & + + aX*(ExGGA - ExLDA) + +! Hartree-Fock calculation + + case(666) + + call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF) + + Ex = ExHF + + end select + +end subroutine unrestricted_exchange_energy diff --git a/src/eDFT/unrestricted_exchange_individual_energy.f90 b/src/eDFT/unrestricted_exchange_individual_energy.f90 new file mode 100644 index 0000000..3a0844e --- /dev/null +++ b/src/eDFT/unrestricted_exchange_individual_energy.f90 @@ -0,0 +1,80 @@ +subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas, & + ERI,Pw,P,rhow,drhow,rho,drho,Ex) + +! Compute the exchange individual energy + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: rung + character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: Pw(nBas,nBas) + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + +! Local variables + + double precision :: ExLDA + double precision :: ExGGA + double precision :: ExHF + +! Output variables + + double precision,intent(out) :: Ex + + select case (rung) + +! Hartree calculation + + case(0) + + Ex = 0d0 + +! LDA functionals + + case(1) + + call unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,ExLDA) + + Ex = ExLDA + +! GGA functionals + + case(2) + + call unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,ExGGA) + + Ex = ExGGA + +! Hybrid functionals + + case(4) + + call print_warning('!!! Individual energies NYI for Hybrids !!!') + stop + +! Hartree-Fock calculation + + case(666) + + call unrestricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,ExHF) + + Ex = ExHF + + end select + +end subroutine unrestricted_exchange_individual_energy diff --git a/src/eDFT/unrestricted_exchange_potential.f90 b/src/eDFT/unrestricted_exchange_potential.f90 new file mode 100644 index 0000000..573a70e --- /dev/null +++ b/src/eDFT/unrestricted_exchange_potential.f90 @@ -0,0 +1,86 @@ +subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,P, & + ERI,AO,dAO,rho,drho,Fx,FxHF) + +! Compute the exchange potential + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: rung + character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: dAO(ncart,nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + +! Local variables + + double precision,allocatable :: FxLDA(:,:),FxGGA(:,:) + double precision :: cX,aX + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas),FxHF(nBas,nBas) + +! Memory allocation + + select case (rung) + +! Hartree calculation + + case(0) + + Fx(:,:) = 0d0 + +! LDA functionals + + case(1) + + call unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx) + +! GGA functionals + + case(2) + + call unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + +! Hybrid functionals + + case(4) + + allocate(FxLDA(nBas,nBas),FxGGA(nBas,nBas)) + + cX = 0.20d0 + aX = 0.72d0 + + call unrestricted_lda_exchange_potential(DFA,nGrid,weight,nBas,AO,rho,FxLDA) + call unrestricted_gga_exchange_potential(DFA,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) + call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + + Fx(:,:) = FxLDA(:,:) & + + cX*(FxHF(:,:) - FxLDA(:,:)) & + + aX*(FxGGA(:,:) - FxLDA(:,:)) + +! Hartree-Fock calculation + + case(666) + + call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + + Fx(:,:) = FxHF(:,:) + + end select + +end subroutine unrestricted_exchange_potential diff --git a/src/eDFT/unrestricted_fock_exchange_energy.f90 b/src/eDFT/unrestricted_fock_exchange_energy.f90 new file mode 100644 index 0000000..d312fa6 --- /dev/null +++ b/src/eDFT/unrestricted_fock_exchange_energy.f90 @@ -0,0 +1,25 @@ +subroutine unrestricted_fock_exchange_energy(nBas,P,Fx,Ex) + +! Compute the (exact) Fock exchange energy + + implicit none + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: Fx(nBas,nBas) + +! Local variables + + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: Ex + +! Compute HF exchange energy + + Ex = trace_matrix(nBas,matmul(P,Fx)) + +end subroutine unrestricted_fock_exchange_energy diff --git a/src/eDFT/unrestricted_fock_exchange_individual_energy.f90 b/src/eDFT/unrestricted_fock_exchange_individual_energy.f90 new file mode 100644 index 0000000..a91a587 --- /dev/null +++ b/src/eDFT/unrestricted_fock_exchange_individual_energy.f90 @@ -0,0 +1,31 @@ +subroutine unrestricted_fock_exchange_individual_energy(nBas,Pw,P,ERI,Ex) + +! Compute the Fock exchange potential + + implicit none + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(in) :: Pw(nBas,nBas) + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision,allocatable :: Fx(:,:) + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: Ex + +! Compute HF exchange matrix + + allocate(Fx(nBas,nBas)) + + call unrestricted_fock_exchange_potential(nBas,Pw(:,:),ERI(:,:,:,:),Fx(:,:)) + Ex = trace_matrix(nBas,matmul(P(:,:),Fx(:,:))) & + - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),Fx(:,:))) + +end subroutine unrestricted_fock_exchange_individual_energy diff --git a/src/eDFT/unrestricted_fock_exchange_potential.f90 b/src/eDFT/unrestricted_fock_exchange_potential.f90 new file mode 100644 index 0000000..acb0ef7 --- /dev/null +++ b/src/eDFT/unrestricted_fock_exchange_potential.f90 @@ -0,0 +1,34 @@ +subroutine unrestricted_fock_exchange_potential(nBas,P,ERI,Fx) + +! Compute the Fock exchange potential + + implicit none + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: mu,nu,la,si + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Compute HF exchange matrix + + Fx(:,:) = 0d0 + do si=1,nBas + do la=1,nBas + do nu=1,nBas + do mu=1,nBas + Fx(mu,nu) = Fx(mu,nu) - P(la,si)*ERI(mu,la,si,nu) + enddo + enddo + enddo + enddo + +end subroutine unrestricted_fock_exchange_potential diff --git a/src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 new file mode 100644 index 0000000..bb0c4a7 --- /dev/null +++ b/src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 @@ -0,0 +1,40 @@ +subroutine unrestricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD) + +! Compute the exchange GGA part of the derivative discontinuity + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + +! Local variables + + +! Output variables + + double precision,intent(out) :: ExDD(nEns) + +! Select correlation functional + + select case (DFA) + + case ('B88') + + ExDD(:) = 0d0 + + case default + + call print_warning('!!! GGA exchange derivative discontinuity not available !!!') + stop + + end select + +end subroutine unrestricted_gga_exchange_derivative_discontinuity diff --git a/src/eDFT/unrestricted_gga_exchange_energy.f90 b/src/eDFT/unrestricted_gga_exchange_energy.f90 new file mode 100644 index 0000000..b7d2679 --- /dev/null +++ b/src/eDFT/unrestricted_gga_exchange_energy.f90 @@ -0,0 +1,40 @@ +subroutine unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) + +! Select GGA exchange functional for energy calculation + + implicit none + + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(3,nGrid) + +! Output variables + + double precision :: Ex + + select case (DFA) + + case ('G96') + + call UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + + case ('B88') + + call UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + + case default + + call print_warning('!!! GGA exchange energy not available !!!') + stop + + end select + +end subroutine unrestricted_gga_exchange_energy diff --git a/src/eDFT/unrestricted_gga_exchange_individual_energy.f90 b/src/eDFT/unrestricted_gga_exchange_individual_energy.f90 new file mode 100644 index 0000000..31977c2 --- /dev/null +++ b/src/eDFT/unrestricted_gga_exchange_individual_energy.f90 @@ -0,0 +1,35 @@ +subroutine unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,Ex) + +! Compute GGA exchange energy for individual states + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: drhow(ncart,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(ncart,nGrid) + +! Output variables + + double precision :: Ex + +! Select correlation functional + + select case (DFA) + + case default + + call print_warning('!!! GGA exchange individual energy not available !!!') + stop + + end select + +end subroutine unrestricted_gga_exchange_individual_energy diff --git a/src/eDFT/unrestricted_gga_exchange_potential.f90 b/src/eDFT/unrestricted_gga_exchange_potential.f90 new file mode 100644 index 0000000..368a0e4 --- /dev/null +++ b/src/eDFT/unrestricted_gga_exchange_potential.f90 @@ -0,0 +1,44 @@ +subroutine unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + +! Select GGA exchange functional for potential calculation + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: dAO(3,nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + double precision,intent(in) :: drho(3,nGrid) + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Select GGA exchange functional + + select case (DFA) + + case ('G96') + + call UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + + case ('B88') + + call UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + + case default + + call print_warning('!!! GGA exchange potential not available !!!') + stop + + end select + +end subroutine unrestricted_gga_exchange_potential diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index 1a5a34f..d18abbe 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -135,9 +135,9 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do iEns=1,nEns do ispin=1,nspin - call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,ERI, & - Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin), & - rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns)) + call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,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 @@ -192,8 +192,8 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered do ispin=1,nspin - call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, & - rhow(:,ispin),drhow(:,:,ispin),ExDD(ispin,:)) + call unrestricted_exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight, & + rhow(:,ispin),drhow(:,:,ispin),ExDD(ispin,:)) end do call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD) diff --git a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 new file mode 100644 index 0000000..e84be1a --- /dev/null +++ b/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 @@ -0,0 +1,46 @@ +subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD) + +! Compute the exchange LDA part of the derivative discontinuity + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + +! Local variables + + +! Output variables + + double precision,intent(out) :: ExDD(nEns) + +! Select correlation functional + + select case (DFA) + + case ('S51') + + ExDD(:) = 0d0 + + case ('CC') + + call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight(:),rhow(:),ExDD(:)) + + case default + + call print_warning('!!! LDA exchange derivative discontinuity not available !!!') + stop + + end select + +end subroutine unrestricted_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/unrestricted_lda_exchange_energy.f90 b/src/eDFT/unrestricted_lda_exchange_energy.f90 new file mode 100644 index 0000000..76829b4 --- /dev/null +++ b/src/eDFT/unrestricted_lda_exchange_energy.f90 @@ -0,0 +1,43 @@ +subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex) + +! Select LDA exchange functional + + implicit none + include 'parameters.h' + +! Input variables + + character(len=12),intent(in) :: DFA + logical,intent(in) :: LDA_centered + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Output variables + + double precision,intent(out) :: Ex + +! Select correlation functional + + select case (DFA) + + case ('S51') + + call US51_lda_exchange_energy(nGrid,weight,rho,Ex) + + case ('CC') + + call UCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex) + + case default + + call print_warning('!!! LDA exchange functional not available !!!') + stop + + end select + +end subroutine unrestricted_lda_exchange_energy diff --git a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 new file mode 100644 index 0000000..56b6a4a --- /dev/null +++ b/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 @@ -0,0 +1,44 @@ +subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex) + +! Compute LDA exchange energy for individual states + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: LDA_centered + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + double precision,intent(in) :: rhow(nGrid) + double precision,intent(in) :: rho(nGrid) + +! Output variables + + double precision :: Ex + +! Select correlation functional + + select case (DFA) + + case ('S51') + + call US51_lda_exchange_individual_energy(nGrid,weight,rhow,rho,Ex) + + case ('CC') + + call UCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex) + + case default + + call print_warning('!!! LDA exchange individual energy not available !!!') + stop + + end select + +end subroutine unrestricted_lda_exchange_individual_energy diff --git a/src/eDFT/unrestricted_lda_exchange_potential.f90 b/src/eDFT/unrestricted_lda_exchange_potential.f90 new file mode 100644 index 0000000..9f47532 --- /dev/null +++ b/src/eDFT/unrestricted_lda_exchange_potential.f90 @@ -0,0 +1,46 @@ +subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx) + +! Select LDA correlation potential + + implicit none + + include 'parameters.h' + +! Input variables + + logical,intent(in) :: LDA_centered + character(len=12),intent(in) :: DFA + integer,intent(in) :: nEns + double precision,intent(in) :: wEns(nEns) + double precision,intent(in) :: aCC_w1(3) + double precision,intent(in) :: aCC_w2(3) + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: rho(nGrid) + +! Output variables + + double precision,intent(out) :: Fx(nBas,nBas) + +! Select exchange functional + + select case (DFA) + + case ('S51') + + call US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) + + case ('CC') + + call UCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx) + + case default + + call print_warning('!!! LDA exchange functional not available !!!') + stop + + end select + +end subroutine unrestricted_lda_exchange_potential