diff --git a/src/eDFT/GOK_RKS.f90 b/src/eDFT/GOK_RKS.f90 deleted file mode 100644 index 7aa201d..0000000 --- a/src/eDFT/GOK_RKS.f90 +++ /dev/null @@ -1,351 +0,0 @@ -subroutine GOK_RKS(restart,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,c,occnum,Cx_choice) - -! Perform restricted Kohn-Sham calculation for ensembles - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: restart - 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) - 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) :: maxSCF - integer,intent(in) :: max_diis - integer,intent(in) :: guess_type - double precision,intent(in) :: thresh - integer,intent(in) :: nBas - double precision,intent(in) :: AO(nBas,nGrid) - double precision,intent(in) :: dAO(ncart,nBas,nGrid) - - integer,intent(in) :: nO - integer,intent(in) :: nV - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ENuc - - double precision,intent(inout):: c(nBas,nBas) - double precision,intent(in) :: occnum(nBas,nspin,nEns) - integer,intent(in) :: Cx_choice - -! Local variables - - integer :: xc_rung - integer :: nSCF,nBasSq - integer :: n_diis - double precision :: conv - double precision :: rcond - double precision :: ET - double precision :: EV - double precision :: EJ - double precision :: Ex - double precision :: Ec - - double precision,allocatable :: eps(:) - double precision,allocatable :: cp(:,:) - double precision,allocatable :: J(:,:) - double precision,allocatable :: F(:,:) - double precision,allocatable :: Fp(:,:) - double precision,allocatable :: Fx(:,:) - double precision,allocatable :: FxHF(:,:) - double precision,allocatable :: Fc(:,:) - double precision,allocatable :: err(:,:) - double precision,allocatable :: err_diis(:,:) - double precision,allocatable :: F_diis(:,:) - double precision,external :: trace_matrix - double precision,external :: electron_number - - double precision,allocatable :: Pw(:,:) - double precision,allocatable :: rhow(:) - double precision,allocatable :: drhow(:,:) - double precision :: nEl - - double precision,allocatable :: P(:,:,:) - double precision,allocatable :: rho(:,:) - double precision,allocatable :: drho(:,:,:) - - double precision :: E(nEns) - double precision :: Om(nEns) - - integer :: iEns - -! Output variables - - double precision,intent(out) :: Ew - -! Hello world - - write(*,*) - write(*,*)'************************************************' - write(*,*)'* Restricted Kohn-Sham calculation *' - write(*,*)'* *** for ensembles *** *' - write(*,*)'************************************************' - write(*,*) - -! Useful stuff - - nBasSq = nBas*nBas - -!------------------------------------------------------------------------ -! Rung of Jacob's ladder -!------------------------------------------------------------------------ - -! Select rung for exchange - - write(*,*) - write(*,*) '*******************************************************************' - write(*,*) '* EXCHANGE RUNG *' - write(*,*) '*******************************************************************' - - call select_rung(x_rung,x_DFA) - -! Select rung for correlation - - write(*,*) - write(*,*) '*******************************************************************' - write(*,*) '* CORRELATION RUNG *' - write(*,*) '*******************************************************************' - - call select_rung(c_rung,c_DFA) - -! Overall rung - - xc_rung = max(x_rung,c_rung) - -! Memory allocation - - allocate(eps(nBas),cp(nBas,nBas),J(nBas,nBas), & - F(nBas,nBas),Fp(nBas,nBas),Fx(nBas,nBas), & - FxHF(nBas,nBas),Fc(nBas,nBas),err(nBas,nBas), & - Pw(nBas,nBas),rhow(nGrid),drhow(ncart,nGrid), & - err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis), & - P(nBas,nBas,nEns),rho(nGrid,nEns),drho(ncart,nGrid,nEns)) - -! Guess coefficients and eigenvalues - - if(.not. restart) then - -! call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,Fx,X,cp,F,Fp,eps,c,P) - - if(guess_type == 1) then - - cp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) - call diagonalize_matrix(nBas,cp(:,:),eps(:)) - c(:,:) = matmul(X(:,:),cp(:,:)) - - else - - print*,'Wrong guess option' - stop - - end if - - end if - -! Initialization - - nSCF = 0 - conv = 1d0 - - nEl = 0d0 - Ex = 0d0 - Ec = 0d0 - - Fx(:,:) = 0d0 - FxHF(:,:) = 0d0 - Fc(:,:) = 0d0 - - n_diis = 0 - F_diis(:,:) = 0d0 - err_diis(:,:) = 0d0 - -!------------------------------------------------------------------------ -! Main SCF loop -!------------------------------------------------------------------------ - - write(*,*) - write(*,*)'------------------------------------------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','E(KS)','|','Ex(KS)','|','Ec(KS)','|','Conv','|','nEl','|' - write(*,*)'------------------------------------------------------------------------------------------' - - do while(conv > thresh .and. nSCF < maxSCF) - -! Increment - - nSCF = nSCF + 1 - -!------------------------------------------------------------------------ -! Compute density matrix -!------------------------------------------------------------------------ - - call restricted_density_matrix(nBas,nEns,nO,c(:,:),P(:,:,:),occnum) - -! Weight-dependent density matrix - - Pw(:,:) = 0d0 - do iEns=1,nEns - Pw(:,:) = Pw(:,:) + wEns(iEns)*P(:,:,iEns) - end do - -!------------------------------------------------------------------------ -! Compute one-electron density and its gradient if necessary -!------------------------------------------------------------------------ - - do iEns=1,nEns - call density(nGrid,nBas,P(:,:,iEns),AO(:,:),rho(:,iEns)) - end do - -! Weight-dependent one-electron density - - rhow(:) = 0d0 - do iEns=1,nEns - rhow(:) = rhow(:) + wEns(iEns)*rho(:,iEns) - end do - - if(xc_rung > 1) then - -! Compute gradient of the one-electron density - - do iEns=1,nEns - call gradient_density(nGrid,nBas,P(:,:,iEns),AO(:,:),dAO(:,:,:),drho(:,:,iEns)) - end do - -! Weight-dependent one-electron density gradient - - drhow(:,:) = 0d0 - do iEns=1,nEns - drhow(:,:) = drhow(:,:) + wEns(iEns)*drho(:,:,iEns) - end do - - end if - -! Build Coulomb repulsion - - call hartree_coulomb(nBas,Pw(:,:),ERI(:,:,:,:),J(:,:)) - -! Compute exchange potential - - 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(:,:),Cx_choice) - -! Compute correlation potential - - call restricted_correlation_potential(c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & - nBas,AO(:,:),dAO(:,:,:),rhow(:),drhow(:,:),Fc(:,:)) - -! Build Fock operator - - F(:,:) = Hc(:,:) + J(:,:) + Fx(:,:) + Fc(:,:) - -! Check convergence - - err(:,:) = matmul(F(:,:),matmul(Pw(:,:),S(:,:))) - matmul(matmul(S(:,:),Pw(:,:)),F(:,:)) - - conv = maxval(abs(err(:,:))) - -! DIIS extrapolation - - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis(:,:),F_diis(:,:),err(:,:),F(:,:)) - -! Reset DIIS if required - - if(abs(rcond) < 1d-15) n_diis = 0 - -! Transform Fock matrix in orthogonal basis - - Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:))) - -! Diagonalize Fock matrix to get eigenvectors and eigenvalues - - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp(:,:),eps(:)) - -! Back-transform eigenvectors in non-orthogonal basis - - c(:,:) = matmul(X(:,:),cp(:,:)) - -!------------------------------------------------------------------------ -! Compute KS energy -!------------------------------------------------------------------------ - -! Kinetic energy - - ET = trace_matrix(nBas,matmul(Pw(:,:),T(:,:))) - -! Potential energy - - EV = trace_matrix(nBas,matmul(Pw(:,:),V(:,:))) - -! Coulomb energy - - EJ = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:),J(:,:))) - -! Exchange energy - - call restricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:),nBas, & - Pw(:,:),FxHF(:,:),rhow(:),drhow(:,:),Ex,Cx_choice) - -! Correlation energy - - call restricted_correlation_energy(c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),Ec) - -! Total energy - - Ew = ET + EV + EJ + Ex + Ec - -! Check the grid accuracy by computing the number of electrons - - nEl = electron_number(nGrid,weight(:),rhow(:)) - -! Dump results - - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',Ew + ENuc,'|',Ex,'|',Ec,'|',conv,'|',nEl,'|' - - end do - write(*,*)'------------------------------------------------------------------------------------------' -!------------------------------------------------------------------------ -! End of SCF loop -!------------------------------------------------------------------------ - -! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) - -! stop - - end if - -! Compute final KS energy - - call print_RKS(nBas,nO,eps(:),c(:,:),ENuc,ET,EV,EJ,Ex,Ec,Ew) - -!------------------------------------------------------------------------ -! Compute individual energies from ensemble energy -!------------------------------------------------------------------------ - - call restricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),aCC_w1,aCC_w2,nGrid,weight(:), & - nBas,nO,nV,T(:,:),V(:,:),ERI(:,:,:,:),ENuc,eps(:),Pw(:,:),rhow(:),drhow(:,:), & - J(:,:),P(:,:,:),rho(:,:),drho(:,:,:),Ew,E(:),Om(:),occnum,Cx_choice) - -end subroutine GOK_RKS diff --git a/src/eDFT/GOK_UKS.f90 b/src/eDFT/GOK_UKS.f90 deleted file mode 100644 index 9e2f82b..0000000 --- a/src/eDFT/GOK_UKS.f90 +++ /dev/null @@ -1,375 +0,0 @@ -subroutine GOK_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,aCC_w1,aCC_w2,maxSCF,thresh,max_diis,guess_type, & - nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew,occnum,Cx_choice) - -! Perform unrestricted Kohn-Sham calculation for ensembles - - 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) - double precision,intent(in) :: aCC_w1(3) - double precision,intent(in) :: aCC_w2(3) - integer,intent(in) :: maxSCF,max_diis,guess_type - double precision,intent(in) :: thresh - 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) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ENuc - double precision,intent(in) :: occnum(nspin,2,nEns) - integer,intent(in) :: Cx_choice - -! Local variables - - integer :: xc_rung - logical :: LDA_centered = .false. - integer :: nSCF,nBasSq - integer :: n_diis - double precision :: conv - double precision :: rcond(nspin) - double precision :: ET(nspin) - double precision :: EV(nspin) - double precision :: EJ(nsp) - double precision :: Ex(nspin) - double precision :: Ec(nsp) - double precision :: Ew - - double precision,allocatable :: eps(:,:) - double precision,allocatable :: c(:,:,:) - double precision,allocatable :: cp(:,:,:) - double precision,allocatable :: J(:,:,:) - double precision,allocatable :: F(:,:,:) - double precision,allocatable :: Fp(:,:,:) - double precision,allocatable :: Fx(:,:,:) - double precision,allocatable :: FxHF(:,:,:) - double precision,allocatable :: Fc(:,:,:) - double precision,allocatable :: err(:,:,:) - double precision,allocatable :: err_diis(:,:,:) - double precision,allocatable :: F_diis(:,:,:) - double precision,external :: trace_matrix - double precision,external :: electron_number - - double precision,allocatable :: Pw(:,:,:) - double precision,allocatable :: rhow(:,:) - double precision,allocatable :: drhow(:,:,:) - double precision :: nEl(nspin) - - double precision,allocatable :: P(:,:,:,:) - double precision,allocatable :: rho(:,:,:) - double precision,allocatable :: drho(:,:,:,:) - - double precision :: E(nEns) - double precision :: Om(nEns) - - integer :: ispin,iEns - -! Hello world - - write(*,*) - write(*,*)'************************************************' - write(*,*)'* Unrestricted Kohn-Sham calculation *' - write(*,*)'* *** for ensembles *** *' - write(*,*)'************************************************' - write(*,*) - -! Useful stuff - - nBasSq = nBas*nBas - -!------------------------------------------------------------------------ -! Rung of Jacob's ladder -!------------------------------------------------------------------------ - -! Select rung for exchange - - write(*,*) - write(*,*) '*******************************************************************' - write(*,*) '* Exchange rung *' - write(*,*) '*******************************************************************' - - call select_rung(x_rung,x_DFA) - -! Select rung for correlation - - write(*,*) - write(*,*) '*******************************************************************' - write(*,*) '* Correlation rung *' - write(*,*) '*******************************************************************' - - call select_rung(c_rung,c_DFA) - -! Overall rung - - xc_rung = max(x_rung,c_rung) - -! Memory allocation - - allocate(eps(nBas,nspin),c(nBas,nBas,nspin),cp(nBas,nBas,nspin), & - J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas,nspin), & - Fx(nBas,nBas,nspin),FxHF(nBas,nBas,nspin),Fc(nBas,nBas,nspin),err(nBas,nBas,nspin), & - Pw(nBas,nBas,nspin),rhow(nGrid,nspin),drhow(ncart,nGrid,nspin), & - err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin), & - P(nBas,nBas,nspin,nEns),rho(nGrid,nspin,nEns),drho(ncart,nGrid,nspin,nEns)) - -! Guess coefficients and eigenvalues - - if(guess_type == 1) then - - do ispin=1,nspin - cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) - call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin)) - c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin)) - end do - - else if(guess_type == 2) then - - do ispin=1,nspin - call random_number(F(:,:,ispin)) - end do - else - - print*,'Wrong guess option' - stop - - end if - -! Initialization - - nSCF = 0 - conv = 1d0 - - nEl(:) = 0d0 - Ex(:) = 0d0 - Ec(:) = 0d0 - - Fx(:,:,:) = 0d0 - FxHF(:,:,:) = 0d0 - Fc(:,:,:) = 0d0 - - n_diis = 0 - F_diis(:,:,:) = 0d0 - err_diis(:,:,:) = 0d0 - -!------------------------------------------------------------------------ -! Main SCF loop -!------------------------------------------------------------------------ - - write(*,*) - write(*,*)'------------------------------------------------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','E(KS)','|','Ex(KS)','|','Ec(KS)','|','Conv','|','nEl','|' - write(*,*)'------------------------------------------------------------------------------------------' - - do while(conv > thresh .and. nSCF < maxSCF) - -! Increment - - nSCF = nSCF + 1 - -!------------------------------------------------------------------------ -! Compute density matrix -!------------------------------------------------------------------------ - - call unrestricted_density_matrix(nBas,nEns,nO(:),c(:,:,:),P(:,:,:,:),occnum) - -! Weight-dependent density matrix - - Pw(:,:,:) = 0d0 - do iEns=1,nEns - Pw(:,:,:) = Pw(:,:,:) + wEns(iEns)*P(:,:,:,iEns) - end do - -!------------------------------------------------------------------------ -! Compute one-electron density and its gradient if necessary -!------------------------------------------------------------------------ - - do ispin=1,nspin - do iEns=1,nEns - call density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),rho(:,ispin,iEns)) - end do - end do - -! Weight-dependent one-electron density - - rhow(:,:) = 0d0 - do iEns=1,nEns - rhow(:,:) = rhow(:,:) + wEns(iEns)*rho(:,:,iEns) - end do - - if(xc_rung > 1) then - -! Ground state density - - do ispin=1,nspin - do iEns=1,nEns - call gradient_density(nGrid,nBas,P(:,:,ispin,iEns),AO(:,:),dAO(:,:,:),drho(:,:,ispin,iEns)) - end do - end do - -! Weight-dependent one-electron density - - drhow(:,:,:) = 0d0 - do iEns=1,nEns - drhow(:,:,:) = drhow(:,:,:) + wEns(iEns)*drho(:,:,:,iEns) - end do - - end if - -! Build Coulomb repulsion - - do ispin=1,nspin - call hartree_coulomb(nBas,Pw(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) - end do - -! Compute exchange potential - - do ispin=1,nspin - 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),Cx_choice) - end do - -! Compute correlation potential - - 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 - F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + Fx(:,:,ispin) + Fc(:,:,ispin) - end do - -! Check convergence - - do ispin=1,nspin - err(:,:,ispin) = matmul(F(:,:,ispin),matmul(Pw(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),Pw(:,:,ispin)),F(:,:,ispin)) - end do - - conv = maxval(abs(err(:,:,:))) - -! DIIS extrapolation - - n_diis = min(n_diis+1,max_diis) - do ispin=1,nspin - call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis, & - err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin)) - end do - -! Reset DIIS if required - - if(minval(rcond(:)) < 1d-15) n_diis = 0 - -! Transform Fock matrix in orthogonal basis - - do ispin=1,nspin - Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:))) - end do - -! Diagonalize Fock matrix to get eigenvectors and eigenvalues - - cp(:,:,:) = Fp(:,:,:) - do ispin=1,nspin - call diagonalize_matrix(nBas,cp(:,:,ispin),eps(:,ispin)) - end do - -! Back-transform eigenvectors in non-orthogonal basis - - do ispin=1,nspin - c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin)) - end do - -!------------------------------------------------------------------------ -! Compute KS energy -!------------------------------------------------------------------------ - -! Kinetic energy - - do ispin=1,nspin - ET(ispin) = trace_matrix(nBas,matmul(Pw(:,:,ispin),T(:,:))) - end do - -! Potential energy - - do ispin=1,nspin - EV(ispin) = trace_matrix(nBas,matmul(Pw(:,:,ispin),V(:,:))) - end do - -! Coulomb energy - - EJ(1) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) - EJ(2) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & - + 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) - EJ(3) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) - -! Exchange energy - - do ispin=1,nspin - 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),Cx_choice) - end do - -! Correlation energy - - call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) - -! Total energy - - Ew = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:)) + sum(Ec(:)) - -! Check the grid accuracy by computing the number of electrons - - do ispin=1,nspin - nEl(ispin) = electron_number(nGrid,weight,rhow(:,ispin)) - end do - -! Dump results - - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',Ew + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',conv,'|',sum(nEl(:)),'|' - - end do - write(*,*)'------------------------------------------------------------------------------------------' -!------------------------------------------------------------------------ -! End of SCF loop -!------------------------------------------------------------------------ - -! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*) - - stop - - end if - -! Compute final KS energy - - call print_UKS(nBas,nO,eps,c,ENuc,ET,EV,EJ,Ex,Ec,Ew) - -!------------------------------------------------------------------------ -! Compute individual energies from ensemble energy -!------------------------------------------------------------------------ - - - call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,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,occnum,Cx_choice) - -end subroutine GOK_UKS diff --git a/src/eDFT/RB88_gga_exchange_energy.f90 b/src/eDFT/RB88_gga_exchange_energy.f90 deleted file mode 100644 index 596710d..0000000 --- a/src/eDFT/RB88_gga_exchange_energy.f90 +++ /dev/null @@ -1,51 +0,0 @@ -subroutine RB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) - -! Compute restricted version of Becke's 88 GGA exchange energy - - implicit none - - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nGrid - double precision,intent(in) :: weight(nGrid) - double precision,intent(in) :: rho(nGrid) - double precision,intent(in) :: drho(ncart,nGrid) - -! Local variables - - integer :: iG - double precision :: beta - double precision :: r,g,x - -! Output variables - - double precision :: Ex - -! Coefficients for B88 GGA exchange functional - - beta = 0.0042d0 - -! Compute GGA exchange energy - - Ex = 0d0 - - do iG=1,nGrid - - r = max(0d0,0.5d0*rho(iG)) - - if(r > threshold) then - g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2) - x = sqrt(g)/r**(4d0/3d0) - - Ex = Ex + weight(iG)*CxLDA*r**(4d0/3d0) & - - weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x)) - - end if - - end do - - Ex = 2d0*Ex - -end subroutine RB88_gga_exchange_energy diff --git a/src/eDFT/RB88_gga_exchange_individual_energy.f90 b/src/eDFT/RB88_gga_exchange_individual_energy.f90 deleted file mode 100644 index 56c64bf..0000000 --- a/src/eDFT/RB88_gga_exchange_individual_energy.f90 +++ /dev/null @@ -1,59 +0,0 @@ -subroutine RB88_gga_exchange_individual_energy(nGrid,weight,rhow,drhow,rho,drho,Ex) - -! Compute restricted Becke's GGA indivudal 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) :: drhow(ncart,nGrid) - double precision,intent(in) :: rho(nGrid) - double precision,intent(in) :: drho(ncart,nGrid) - -! Local variables - - integer :: iG - double precision :: beta - double precision :: r,rI,g,x - double precision :: ex_p,dexdr_p - -! Output variables - - double precision,intent(out) :: Ex - -! Coefficients for B88 GGA exchange functional - - beta = 0.0042d0 - -! Compute GGA exchange matrix in the AO basis - - Ex = 0d0 - - do iG=1,nGrid - - r = max(0d0,0.5d0*rhow(iG)) - rI = max(0d0,0.5d0*rho(iG)) - - if(r > threshold .or. rI > threshold) then - - g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2) - x = sqrt(g)/r**(4d0/3d0) - - dexdr_p = 4d0/3d0*r**(1d0/3d0)*(CxLDA - beta*g**(3d0/4d0)/r**2) & - + 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0) & - - 2d0*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0) - - ex_p = CxLDA*r**(4d0/3d0) & - - weight(iG)*beta*x**2*r**(4d0/3d0)/(1d0 + 6d0*beta*x*asinh(x)) - - Ex = Ex + weight(iG)*(ex_p*rI + dexdr_p*r*rI - dexdr_p*r*r) - - end if - - end do - -end subroutine RB88_gga_exchange_individual_energy diff --git a/src/eDFT/RB88_gga_exchange_potential.f90 b/src/eDFT/RB88_gga_exchange_potential.f90 deleted file mode 100644 index 6829644..0000000 --- a/src/eDFT/RB88_gga_exchange_potential.f90 +++ /dev/null @@ -1,64 +0,0 @@ -subroutine RB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) - -! Compute restricted Becke's GGA exchange potential - - implicit none - include 'parameters.h' - -! Input variables - - 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) - double precision,intent(in) :: rho(nGrid) - double precision,intent(in) :: drho(ncart,nGrid) - -! Local variables - - integer :: mu,nu,iG - double precision :: beta - double precision :: r,g,vAO,gAO - -! Output variables - - double precision,intent(out) :: Fx(nBas,nBas) - -! Coefficients for B88 GGA exchange functional - - beta = 0.0042d0 - -! Compute GGA exchange matrix in the AO basis - - Fx(:,:) = 0d0 - - do mu=1,nBas - do nu=1,nBas - do iG=1,nGrid - - r = max(0d0,0.5d0*rho(iG)) - - if(r > threshold) then - - g = 0.25d0*(drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2) - vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) - Fx(mu,nu) = Fx(mu,nu) & - + vAO*(4d0/3d0*r**(1d0/3d0)*(CxLDA - beta*g**(3d0/4d0)/r**2) & - + 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0)) - - gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) & - + drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) & - + drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG)) - - gAO = weight(iG)*gAO - - Fx(mu,nu) = Fx(mu,nu) - 2d0*gAO*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0) - - end if - - end do - end do - end do - -end subroutine RB88_gga_exchange_potential diff --git a/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 deleted file mode 100644 index 14d1083..0000000 --- a/src/eDFT/RCC_lda_exchange_derivative_discontinuity.f90 +++ /dev/null @@ -1,123 +0,0 @@ -subroutine RCC_lda_exchange_derivative_discontinuity(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,ExDD,Cx_choice) - -! Compute the restricted version of the curvature-corrected exchange ensemble derivative - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nEns - double precision,intent(in) :: wEns(nEns) - 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) - integer,intent(in) :: Cx_choice - -! Local variables - - integer :: iEns,jEns - integer :: iG - double precision :: r - double precision,allocatable :: dExdw(:) - double precision,external :: Kronecker_delta - - double precision :: a1,b1,c1,w1 - double precision :: a2,b2,c2,w2 - double precision :: dCxdw1,dCxdw2 - -! Output variables - - double precision,intent(out) :: ExDD(nEns) - -! Memory allocation - - allocate(dExdw(nEns)) - -! Single excitation parameters - -! a1 = 0.0d0 -! b1 = 0.0d0 -! c1 = 0.0d0 - -! Parameters for H2 at equilibrium - -! a2 = +0.5751782560799208d0 -! b2 = -0.021108186591137282d0 -! c2 = -0.36718902716347124d0 - -! Parameters for stretch H2 - -! a2 = + 0.01922622507087411d0 -! b2 = - 0.01799647558018601d0 -! c2 = - 0.022945430666782573d0 - -! Parameters for He - -! a2 = 1.9125735895875828d0 -! b2 = 2.715266992840757d0 -! c2 = 2.1634223380633086d0 - - a1 = aCC_w1(1) - b1 = aCC_w1(2) - c1 = aCC_w1(3) - - a2 = aCC_w2(1) - b2 = aCC_w2(2) - c2 = aCC_w2(3) - - w1 = wEns(2) - w2 = wEns(3) - - select case (Cx_choice) - case(1) - dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) - dCxdw2 = 0.d0 - case(2) - dCxdw1 = 0.d0 - dCxdw2 =(0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) - case(3) - dCxdw1 = (0.5d0*b1 + (2d0*a1 + 0.5d0*c1)*(w1 - 0.5d0) - (1d0 - w1)*w1*(3d0*b1 + 4d0*c1*(w1 - 0.5d0))) & - * (1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2)) - - dCxdw2 = (1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2)) & - * (0.5d0*b2 + (2d0*a2 + 0.5d0*c2)*(w2 - 0.5d0) - (1d0 - w2)*w2*(3d0*b2 + 4d0*c2*(w2 - 0.5d0))) - - case default - dCxdw1 = 0.d0 - dCxdw2 = 0.d0 - - end select - - dCxdw1 = CxLDA*dCxdw1 - dCxdw2 = CxLDA*dCxdw2 - - dExdw(:) = 0d0 - - do iG=1,nGrid - - r = max(0d0,rhow(iG)) - - if(r > threshold) then - - dExdw(1) = 0d0 - dExdw(2) = dExdw(2) + weight(iG)*dCxdw1*r**(4d0/3d0) - dExdw(3) = dExdw(3) + weight(iG)*dCxdw2*r**(4d0/3d0) - - end if - - end do - - ExDD(:) = 0d0 - - do iEns=1,nEns - do jEns=2,nEns - - ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) - - end do - end do - -end subroutine RCC_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/RCC_lda_exchange_energy.f90 b/src/eDFT/RCC_lda_exchange_energy.f90 deleted file mode 100644 index 5f55d95..0000000 --- a/src/eDFT/RCC_lda_exchange_energy.f90 +++ /dev/null @@ -1,101 +0,0 @@ -subroutine RCC_lda_exchange_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rho,Ex,Cx_choice) - -! Compute the restricted version of the curvature-corrected exchange functional - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nEns - double precision,intent(in) :: wEns(nEns) - 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) - integer,intent(in) :: Cx_choice - -! Local variables - - integer :: iG - double precision :: r - - double precision :: a1,b1,c1,w1 - double precision :: a2,b2,c2,w2 - double precision :: Fx1,Fx2,Cx - -! Output variables - - double precision :: Ex - -! Single excitation parameter - -! a1 = 0.0d0 -! b1 = 0.0d0 -! c1 = 0.0d0 - -! Parameters for H2 at equilibrium - -! a2 = +0.5751782560799208d0 -! b2 = -0.021108186591137282d0 -! c2 = -0.36718902716347124d0 - -! Parameters for stretch H2 - -! a2 = + 0.01922622507087411d0 -! b2 = - 0.01799647558018601d0 -! c2 = - 0.022945430666782573d0 - -! Parameters for He - -! a2 = 1.9125735895875828d0 -! b2 = 2.715266992840757d0 -! c2 = 2.1634223380633086d0 - - a1 = aCC_w1(1) - b1 = aCC_w1(2) - c1 = aCC_w1(3) - - a2 = aCC_w2(1) - b2 = aCC_w2(2) - c2 = aCC_w2(3) - - - 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) - - select case (Cx_choice) - - case(1) - Cx = CxLDA*Fx1 - - case(2) - Cx = CxLDA*Fx2 - - case(3) - Cx = CxLDA*Fx2*Fx1 - - case default - Cx = CxLDA - - end select - -! Compute GIC-LDA exchange energy - - Ex = 0d0 - - do iG=1,nGrid - - r = max(0d0,rho(iG)) - - if(r > threshold) then - Ex = Ex + weight(iG)*Cx*r**(4d0/3d0) - endif - - enddo - -end subroutine RCC_lda_exchange_energy diff --git a/src/eDFT/RCC_lda_exchange_individual_energy.f90 b/src/eDFT/RCC_lda_exchange_individual_energy.f90 deleted file mode 100644 index ade4f56..0000000 --- a/src/eDFT/RCC_lda_exchange_individual_energy.f90 +++ /dev/null @@ -1,111 +0,0 @@ -subroutine RCC_lda_exchange_individual_energy(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,rhow,rho,Ex,Cx_choice) - -! Compute the restricted version of the curvature-corrected exchange functional - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nEns - double precision,intent(in) :: wEns(nEns) - 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) - integer,intent(in) :: Cx_choice - -! Local variables - - integer :: iG - double precision :: r,rI - double precision :: e_p,dedr - - double precision :: a1,b1,c1,w1 - double precision :: a2,b2,c2,w2 - double precision :: Fx1,Fx2,Cx - -! Output variables - - double precision,intent(out) :: Ex - -! Single excitation parameter - -! a1 = 0.0d0 -! b1 = 0.0d0 -! c1 = 0.0d0 - -! Parameters for H2 at equilibrium - -! a2 = +0.5751782560799208d0 -! b2 = -0.021108186591137282d0 -! c2 = -0.36718902716347124d0 - -! Parameters for stretch H2 - -! a2 = + 0.01922622507087411d0 -! b2 = - 0.01799647558018601d0 -! c2 = - 0.022945430666782573d0 - -! Parameters for He - -! a2 = 1.9125735895875828d0 -! b2 = 2.715266992840757d0 -! c2 = 2.1634223380633086d0 - - - - a1 = aCC_w1(1) - b1 = aCC_w1(2) - c1 = aCC_w1(3) - - - - a2 = aCC_w2(1) - b2 = aCC_w2(2) - c2 = aCC_w2(3) - - - 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) - - select case (Cx_choice) - - case(1) - Cx = CxLDA*Fx1 - - case(2) - Cx = CxLDA*Fx2 - - case(3) - Cx = CxLDA*Fx2*Fx1 - - case default - Cx = CxLDA - - end select - -! 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 .or. rI > threshold) then - - e_p = Cx*r**(1d0/3d0) - dedr = 1d0/3d0*Cx*r**(-2d0/3d0) - Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r) - - endif - - enddo - -end subroutine RCC_lda_exchange_individual_energy diff --git a/src/eDFT/RCC_lda_exchange_potential.f90 b/src/eDFT/RCC_lda_exchange_potential.f90 deleted file mode 100644 index 519dfa9..0000000 --- a/src/eDFT/RCC_lda_exchange_potential.f90 +++ /dev/null @@ -1,114 +0,0 @@ -subroutine RCC_lda_exchange_potential(nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,AO,rho,Fx,Cx_choice) - -! Compute the restricted version of the curvature-corrected exchange potential - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nEns - double precision,intent(in) :: wEns(nEns) - 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) - integer,intent(in) :: Cx_choice - -! Local variables - - integer :: mu,nu,iG - double precision :: r,vAO - - double precision :: a1,b1,c1,w1 - double precision :: a2,b2,c2,w2 - double precision :: Fx1,Fx2,Cx - -! Output variables - - double precision,intent(out) :: Fx(nBas,nBas) - -! Single excitation parameter - -! a1 = 0.0d0 -! b1 = 0.0d0 -! c1 = 0.0d0 - -! Parameters for H2 at equilibrium - -! a2 = +0.5751782560799208d0 -! b2 = -0.021108186591137282d0 -! c2 = -0.36718902716347124d0 - -! Parameters for stretch H2 - -! a2 = + 0.01922622507087411d0 -! b2 = - 0.01799647558018601d0 -! c2 = - 0.022945430666782573d0 - -! Parameters for He - -! a2 = 1.9125735895875828d0 -! b2 = 2.715266992840757d0 -! c2 = 2.1634223380633086d0 - -! Parameters for He N -> N-1 - - a1 = aCC_w1(1) - b1 = aCC_w1(2) - c1 = aCC_w1(3) - -! Parameters for He N -> N+1 - - a2 = aCC_w2(1) - b2 = aCC_w2(2) - c2 = aCC_w2(3) - - 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) - - select case (Cx_choice) - - case(1) - Cx = CxLDA*Fx1 - - case(2) - Cx = CxLDA*Fx2 - - case(3) - Cx = CxLDA*Fx2*Fx1 - - case default - Cx = CxLDA - - end select - - -! Compute LDA exchange matrix in the AO basis - - Fx(:,:) = 0d0 - - do mu=1,nBas - do nu=1,nBas - do iG=1,nGrid - - r = max(0d0,rho(iG)) - - if(r > threshold) then - - vAO = weight(iG)*AO(mu,iG)*AO(nu,iG) - Fx(mu,nu) = Fx(mu,nu) + vAO*4d0/3d0*Cx*r**(1d0/3d0) - - endif - - enddo - enddo - enddo - -end subroutine RCC_lda_exchange_potential diff --git a/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 deleted file mode 100644 index 5f3e4a3..0000000 --- a/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 +++ /dev/null @@ -1,61 +0,0 @@ -subroutine RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,EcDD) - -! Compute the restricted version of the 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) - -! Local variables - - integer :: iEns,jEns - double precision,allocatable :: aMFL(:,:) - double precision :: dEcdw(nEns) - double precision,external :: Kronecker_delta - -! Output variables - - double precision,intent(out) :: EcDD(nEns) - -! Allocation - - allocate(aMFL(3,nEns)) - -! Parameters for weight-dependent LDA correlation functional - - aMFL(1,1) = -0.0238184d0 - aMFL(2,1) = +0.00540994d0 - aMFL(3,1) = +0.0830766d0 - - aMFL(1,2) = -0.0282814d0 - aMFL(2,2) = +0.00273925d0 - aMFL(3,2) = +0.0664914d0 - - aMFL(1,3) = -0.0144633d0 - aMFL(2,3) = -0.0506019d0 - aMFL(3,3) = +0.0331417d0 - -! Compute correlation energy for ground, singly-excited and doubly-excited states - - do iEns=1,nEns - call restricted_elda_correlation_energy(aMFL(1:3,iEns),nGrid,weight(:),rhow(:),dEcdw(iEns)) - end do - - EcDD(:) = 0d0 - - do iEns=1,nEns - do jEns=1,nEns - - EcDD(iEns) = EcDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*(dEcdw(jEns) - dEcdw(1)) - - end do - end do - -end subroutine RMFL20_lda_correlation_derivative_discontinuity diff --git a/src/eDFT/RMFL20_lda_correlation_energy.f90 b/src/eDFT/RMFL20_lda_correlation_energy.f90 deleted file mode 100644 index 731378f..0000000 --- a/src/eDFT/RMFL20_lda_correlation_energy.f90 +++ /dev/null @@ -1,69 +0,0 @@ -subroutine RMFL20_lda_correlation_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ec) - -! Compute the restricted version of the Marut-Fromager-Loos weight-dependent correlation functional -! The RMFL20 is a two-state, single-weight correlation functional for spin-unpolarized systems - - implicit none - include 'parameters.h' - -! Input variables - - 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) - double precision,intent(in) :: rho(nGrid) - -! Local variables - - integer :: iEns - double precision :: EcLDA - double precision,allocatable :: aMFL(:,:) - double precision,allocatable :: EceLDA(:) - -! Output variables - - double precision :: Ec - -! Allocation - - allocate(aMFL(3,nEns),EceLDA(nEns)) - -! Parameters for weight-dependent LDA correlation functional - - aMFL(1,1) = -0.0238184d0 - aMFL(2,1) = +0.00540994d0 - aMFL(3,1) = +0.0830766d0 - - aMFL(1,2) = -0.0282814d0 - aMFL(2,2) = +0.00273925d0 - aMFL(3,2) = +0.0664914d0 - - aMFL(1,3) = -0.0144633d0 - aMFL(2,3) = -0.0506019d0 - aMFL(3,3) = +0.0331417d0 - -! Compute correlation energy for ground and doubly-excited states - - do iEns=1,nEns - call restricted_elda_correlation_energy(aMFL(1:3,iEns),nGrid,weight(:),rho(:),EceLDA(iEns)) - end do - -! LDA-centered functional - - if(LDA_centered) then - - call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),EcLDA) - - EceLDA(3) = EcLDA + wEns(3)*(EceLDA(3) - EceLDA(1)) - EceLDA(2) = EcLDA + wEns(2)*(EceLDA(2) - EceLDA(1)) - EceLDA(1) = EcLDA - - end if - -! Weight-denpendent functional for ensembles - - Ec = dot_product(wEns(:),EceLDA(:)) - -end subroutine RMFL20_lda_correlation_energy diff --git a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 deleted file mode 100644 index 93a6170..0000000 --- a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 +++ /dev/null @@ -1,69 +0,0 @@ -subroutine RMFL20_lda_correlation_individual_energy(LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ec) - -! Compute eLDA correlation energy - - implicit none - include 'parameters.h' - -! Input variables - - 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) - double precision,intent(in) :: rhow(nGrid) - double precision,intent(in) :: rho(nGrid) - -! Local variables - - integer :: iEns - double precision :: EcLDA - double precision,allocatable :: aMFL(:,:) - double precision,allocatable :: EceLDA(:) - -! Output variables - - double precision :: Ec - -! Allocation - - allocate(aMFL(3,nEns),EceLDA(nEns)) - -! Parameters for weight-dependent LDA correlation functional - - aMFL(1,1) = -0.0238184d0 - aMFL(2,1) = +0.00540994d0 - aMFL(3,1) = +0.0830766d0 - - aMFL(1,2) = -0.0282814d0 - aMFL(2,2) = +0.00273925d0 - aMFL(3,2) = +0.0664914d0 - - aMFL(1,3) = -0.0144633d0 - aMFL(2,3) = -0.0506019d0 - aMFL(3,3) = +0.0331417d0 - -! Compute correlation energy for ground- and doubly-excited states - - do iEns=1,nEns - call restricted_elda_correlation_individual_energy(aMFL(1:3,iEns),nGrid,weight(:),rhow(:),rho(:),EceLDA(iEns)) - end do - -! LDA-centered functional - - if(LDA_centered) then - - call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA) - - EceLDA(3) = EcLDA + wEns(3)*(EceLDA(3) - EceLDA(1)) - EceLDA(2) = EcLDA + wEns(2)*(EceLDA(2) - EceLDA(1)) - EceLDA(1) = EcLDA - - end if - -! Weight-denpendent functional for ensembles - - Ec = dot_product(wEns(:),EceLDA(:)) - -end subroutine RMFL20_lda_correlation_individual_energy diff --git a/src/eDFT/RMFL20_lda_correlation_potential.f90 b/src/eDFT/RMFL20_lda_correlation_potential.f90 deleted file mode 100644 index e1ffffc..0000000 --- a/src/eDFT/RMFL20_lda_correlation_potential.f90 +++ /dev/null @@ -1,73 +0,0 @@ -subroutine RMFL20_lda_correlation_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) - -! Compute Marut-Fromager-Loos weight-dependent LDA correlation potential - - implicit none - include 'parameters.h' - -! Input variables - - 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) :: rho(nGrid) - -! Local variables - - integer :: iEns - double precision,allocatable :: aMFL(:,:) - double precision,allocatable :: FcLDA(:,:) - double precision,allocatable :: FceLDA(:,:,:) - -! Output variables - - double precision,intent(out) :: Fc(nBas,nBas) - -! Allocation - - allocate(aMFL(3,nEns),FcLDA(nBas,nBas),FceLDA(nBas,nBas,nEns)) - -! Parameters for weight-dependent LDA correlation functional - - aMFL(1,1) = -0.0238184d0 - aMFL(2,1) = +0.00540994d0 - aMFL(3,1) = +0.0830766d0 - - aMFL(1,2) = -0.0282814d0 - aMFL(2,2) = +0.00273925d0 - aMFL(3,2) = +0.0664914d0 - - aMFL(1,3) = -0.0144633d0 - aMFL(2,3) = -0.0506019d0 - aMFL(3,3) = +0.0331417d0 - -! Compute correlation energy for ground- and doubly-excited states - - do iEns=1,nEns - call restricted_elda_correlation_potential(aMFL(1:3,iEns),nGrid,weight(:),nBas,AO(:,:),rho(:),FceLDA(:,:,iEns)) - end do - -! LDA-centered functional - - if(LDA_centered) then - - call RVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:),FcLDA(:,:)) - - FceLDA(:,:,3) = FcLDA(:,:) + wEns(3)*(FceLDA(:,:,3) - FceLDA(:,:,1)) - FceLDA(:,:,2) = FcLDA(:,:) + wEns(2)*(FceLDA(:,:,2) - FceLDA(:,:,1)) - FceLDA(:,:,1) = FcLDA(:,:) - - end if - -! Weight-denpendent functional for ensembles - - Fc(:,:) = 0d0 - do iEns=1,nEns - Fc(:,:) = Fc(:,:) + wEns(iEns)*FceLDA(:,:,iEns) - enddo - -end subroutine RMFL20_lda_correlation_potential diff --git a/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 deleted file mode 100644 index d3620c8..0000000 --- a/src/eDFT/RMFL20_lda_exchange_derivative_discontinuity.f90 +++ /dev/null @@ -1,58 +0,0 @@ -subroutine RMFL20_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rhow,ExDD) - -! Compute the restricted version of the eLDA exchange 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) - -! Local variables - - integer :: iEns,jEns - integer :: iG - double precision :: r - double precision :: dExdw(nEns) - double precision,external :: Kronecker_delta - - double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) - double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) - -! Output variables - - double precision,intent(out) :: ExDD(nEns) - -! Compute correlation energy for ground- and doubly-excited states - - dExdw(:) = 0d0 - - do iG=1,nGrid - - r = max(0d0,rhow(iG)) - - if(r > threshold) then - - dExdw(1) = 0d0 - dExdw(2) = dExdw(2) + weight(iG)*(Cx1 - Cx0)*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 RMFL20_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/RMFL20_lda_exchange_energy.f90 b/src/eDFT/RMFL20_lda_exchange_energy.f90 deleted file mode 100644 index 3d50004..0000000 --- a/src/eDFT/RMFL20_lda_exchange_energy.f90 +++ /dev/null @@ -1,53 +0,0 @@ -subroutine RMFL20_lda_exchange_energy(LDA_centered,nEns,wEns,nGrid,weight,rho,Ex) - -! Compute the restricted version of the Marut-Fromager-Loos weight-dependent exchange functional -! The RMFL20 is a two-state, single-weight exchange functional - - implicit none - include 'parameters.h' - -! Input variables - - 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) - double precision,intent(in) :: rho(nGrid) - - double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) - double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) - -! Local variables - - integer :: iG - double precision :: Cxw - double precision :: r - -! Output variables - - double precision :: Ex - -! Weight-denepdent Cx coefficient - - if(LDA_centered) then - Cxw = CxLDA + (Cx1 - Cx0)*wEns(2) - else - Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 - end if - -! Compute LDA exchange energy - - Ex = 0d0 - - do iG=1,nGrid - - r = max(0d0,rho(iG)) - - if(r > threshold) then - Ex = Ex + weight(iG)*Cxw*r**(4d0/3d0) - endif - - enddo - -end subroutine RMFL20_lda_exchange_energy diff --git a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 b/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 deleted file mode 100644 index 49a1299..0000000 --- a/src/eDFT/RMFL20_lda_exchange_individual_energy.f90 +++ /dev/null @@ -1,58 +0,0 @@ -subroutine RMFL20_lda_exchange_individual_energy(LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,Ex) - -! Compute the restricted version of the Marut-Fromager-Loos 2020 weight-dependent exchange functional - - implicit none - include 'parameters.h' - -! Input variables - - 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) - double precision,intent(in) :: rhow(nGrid) - double precision,intent(in) :: rho(nGrid) - -! Local variables - - integer :: iG - double precision :: Cxw - double precision :: r,rI - double precision :: e_p,dedr - - double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) - double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) - -! Output variables - - double precision,intent(out) :: Ex - -! Weight-dependent Cx coefficient for RMFL20 exchange functional - - if(LDA_centered) then - Cxw = CxLDA + (Cx1 - Cx0)*wEns(2) - else - Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 - end if - -! 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 .or. rI > threshold) then - - e_p = Cxw*r**(1d0/3d0) - dedr = 1d0/3d0*Cxw*r**(-2d0/3d0) - Ex = Ex + weight(iG)*(e_p*rI + dedr*r*rI - dedr*r*r) - - endif - - enddo - -end subroutine RMFL20_lda_exchange_individual_energy diff --git a/src/eDFT/RMFL20_lda_exchange_potential.f90 b/src/eDFT/RMFL20_lda_exchange_potential.f90 deleted file mode 100644 index 4d0ca82..0000000 --- a/src/eDFT/RMFL20_lda_exchange_potential.f90 +++ /dev/null @@ -1,61 +0,0 @@ -subroutine RMFL20_lda_exchange_potential(LDA_centered,nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) - -! Compute the restricted version of the weight-dependent MFL20 exchange potential - - implicit none - include 'parameters.h' - -! Input variables - - 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) :: rho(nGrid) - -! Local variables - - integer :: mu,nu,iG - double precision :: Cxw - double precision :: r,vAO - - double precision,parameter :: Cx0 = - (4d0/3d0)*(1d0/pi)**(1d0/3d0) - double precision,parameter :: Cx1 = - (176d0/105d0)*(1d0/pi)**(1d0/3d0) - -! Output variables - - double precision,intent(out) :: Fx(nBas,nBas) - -! Weight-dependent Cx coefficient for RMFL20 exchange functional - - if(LDA_centered) then - Cxw = CxLDA + (Cx1 - Cx0)*wEns(2) - else - Cxw = wEns(1)*Cx0 + wEns(2)*Cx1 - end if - -! 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*Cxw*r**(1d0/3d0) - - endif - - enddo - enddo - enddo - -end subroutine RMFL20_lda_exchange_potential diff --git a/src/eDFT/RS51_lda_exchange_energy.f90 b/src/eDFT/RS51_lda_exchange_energy.f90 deleted file mode 100644 index b12cf3f..0000000 --- a/src/eDFT/RS51_lda_exchange_energy.f90 +++ /dev/null @@ -1,38 +0,0 @@ -subroutine RS51_lda_exchange_energy(nGrid,weight,rho,Ex) - -! Compute restriced version of Slater's LDA exchange energy - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nGrid - double precision,intent(in) :: weight(nGrid) - double precision,intent(in) :: rho(nGrid) - -! Local variables - - integer :: iG - double precision :: r - -! Output variables - - double precision :: Ex - -! Compute LDA exchange energy - - Ex = 0d0 - do iG=1,nGrid - - r = max(0d0,rho(iG)) - - if(r > threshold) then - - Ex = Ex + weight(iG)*CxLDA*r**(4d0/3d0) - - endif - - enddo - -end subroutine RS51_lda_exchange_energy diff --git a/src/eDFT/RS51_lda_exchange_individual_energy.f90 b/src/eDFT/RS51_lda_exchange_individual_energy.f90 deleted file mode 100644 index f2cdfad..0000000 --- a/src/eDFT/RS51_lda_exchange_individual_energy.f90 +++ /dev/null @@ -1,44 +0,0 @@ -subroutine RS51_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 - double precision :: e,dedr - -! Output variables - - double precision,intent(out) :: Ex - -! 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 .or. rI > threshold) then - - e = CxLDA*r**(1d0/3d0) - dedr = 1d0/3d0*CxLDA*r**(-2d0/3d0) - - Ex = Ex + weight(iG)*(e*rI + dedr*r*rI - dedr*r*r) - - endif - - enddo - -end subroutine RS51_lda_exchange_individual_energy diff --git a/src/eDFT/RS51_lda_exchange_potential.f90 b/src/eDFT/RS51_lda_exchange_potential.f90 deleted file mode 100644 index 2270a47..0000000 --- a/src/eDFT/RS51_lda_exchange_potential.f90 +++ /dev/null @@ -1,45 +0,0 @@ -subroutine RS51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) - -! Compute the restricted version of Slater's LDA exchange potential - - implicit none - include 'parameters.h' - -! Input variables - - 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 - -! Output variables - - double precision,intent(out) :: Fx(nBas,nBas) - -! 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*CxLDA*r**(1d0/3d0) - - endif - - enddo - enddo - enddo - -end subroutine RS51_lda_exchange_potential diff --git a/src/eDFT/RVWN5_lda_correlation_energy.f90 b/src/eDFT/RVWN5_lda_correlation_energy.f90 deleted file mode 100644 index 4e3d3db..0000000 --- a/src/eDFT/RVWN5_lda_correlation_energy.f90 +++ /dev/null @@ -1,62 +0,0 @@ -subroutine RVWN5_lda_correlation_energy(nGrid,weight,rho,Ec) - -! Compute the restricted VWN5 LDA correlation energy - - implicit none - - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nGrid - double precision,intent(in) :: weight(nGrid) - double precision,intent(in) :: rho(nGrid) - -! Local variables - - integer :: iG - double precision :: r,rs,x - double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p - - double precision :: ec_p - -! Output variables - - double precision :: Ec - -! Parameters of the functional - - a_p = +0.0621814D0/2D0 - x0_p = -0.10498d0 - b_p = +3.72744d0 - c_p = +12.9352d0 - -! Initialization - - Ec = 0d0 - - do iG=1,nGrid - - r = max(0d0,rho(iG)) - - if(r > threshold) then - - rs = (4d0*pi*r/3d0)**(-1d0/3d0) - x = sqrt(rs) - - x_p = x*x + b_p*x + c_p - - xx0_p = x0_p*x0_p + b_p*x0_p + c_p - - q_p = sqrt(4d0*c_p - b_p*b_p) - - ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & - - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) - - Ec = Ec + weight(iG)*ec_p*r - - end if - - end do - -end subroutine RVWN5_lda_correlation_energy diff --git a/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 deleted file mode 100644 index d6fb0a1..0000000 --- a/src/eDFT/RVWN5_lda_correlation_individual_energy.f90 +++ /dev/null @@ -1,73 +0,0 @@ -subroutine RVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) - -! Compute the restricted VWN5 LDA correlation individual energies - - 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,rs,x - double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p - double precision :: dxdrs,dxdx_p,decdx_p - double precision :: drsdr,decdr_p - double precision :: ec_p - -! Output variables - - double precision :: Ec - -! Parameters of the functional - - a_p = +0.0621814D0/2D0 - x0_p = -0.10498d0 - b_p = +3.72744d0 - c_p = +12.9352d0 - -! Initialization - - Ec = 0d0 - - do iG=1,nGrid - - r = max(0d0,rhow(iG)) - rI = max(0d0,rho(iG)) - - if(r > threshold .or. rI > threshold) then - - rs = (4d0*pi*r/3d0)**(-1d0/3d0) - x = sqrt(rs) - - x_p = x*x + b_p*x + c_p - xx0_p = x0_p*x0_p + b_p*x0_p + c_p - q_p = sqrt(4d0*c_p - b_p*b_p) - - ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & - - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) - - drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) - dxdrs = 0.5d0/sqrt(rs) - - dxdx_p = 2d0*x + b_p - - decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & - - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) ) - - decdr_p = drsdr*dxdrs*decdx_p - - Ec = Ec + weight(iG)*(ec_p*rI + decdr_p*r*rI - decdr_p*r*r) - - end if - - end do - -end subroutine RVWN5_lda_correlation_individual_energy diff --git a/src/eDFT/RVWN5_lda_correlation_potential.f90 b/src/eDFT/RVWN5_lda_correlation_potential.f90 deleted file mode 100644 index 918f18a..0000000 --- a/src/eDFT/RVWN5_lda_correlation_potential.f90 +++ /dev/null @@ -1,78 +0,0 @@ -subroutine RVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) - -! Compute the restricted VWN5 LDA correlation potential - - implicit none - - include 'parameters.h' - -! Input variables - - 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,rs,x - double precision :: a_p,x0_p,xx0_p,b_p,c_p,x_p,q_p - double precision :: dxdrs,dxdx_p,decdx_p - double precision :: drsdr,decdr_p - - double precision :: ec_p - -! Output variables - - double precision :: Fc(nBas,nBas) - -! Parameters of the functional - - a_p = +0.0621814d0/2d0 - x0_p = -0.10498d0 - b_p = +3.72744d0 - c_p = +12.9352d0 - -! Initialization - - Fc(:,:) = 0d0 - - do mu=1,nBas - do nu=1,nBas - do iG=1,nGrid - - r = max(0d0,rho(iG)) - - if(r > threshold) then - - rs = (4d0*pi*r/3d0)**(-1d0/3d0) - x = sqrt(rs) - - x_p = x*x + b_p*x + c_p - xx0_p = x0_p*x0_p + b_p*x0_p + c_p - q_p = sqrt(4d0*c_p - b_p*b_p) - - ec_p = a_p*( log(x**2/x_p) + 2d0*b_p/q_p*atan(q_p/(2d0*x + b_p)) & - - b_p*x0_p/xx0_p*( log((x - x0_p)**2/x_p) + 2d0*(b_p + 2d0*x0_p)/q_p*atan(q_p/(2d0*x + b_p)) ) ) - - drsdr = - (36d0*pi)**(-1d0/3d0)*r**(-4d0/3d0) - dxdrs = 0.5d0/sqrt(rs) - - dxdx_p = 2d0*x + b_p - - decdx_p = a_p*( 2d0/x - 4d0*b_p/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p & - - b_p*x0_p/xx0_p*( 2/(x-x0_p) - 4d0*(b_p+2d0*x0_p)/( (b_p+2d0*x)**2 + q_p**2) - dxdx_p/x_p ) ) - - decdr_p = drsdr*dxdrs*decdx_p - - Fc(mu,nu) = Fc(mu,nu) + weight(iG)*AO(mu,iG)*AO(nu,iG)*(ec_p + decdr_p*r) - - end if - - end do - end do - end do - -end subroutine RVWN5_lda_correlation_potential