pfon_occupations Subroutine

public subroutine pfon_occupations(mo_energy, nbf, nelec, occ, beta_pfon, scf_type, nsmear, is_beta, nelec_a, nelec_b)

Uses

  • proc~~pfon_occupations~~UsesGraph proc~pfon_occupations pfon_occupations module~precision precision proc~pfon_occupations->module~precision iso_fortran_env iso_fortran_env module~precision->iso_fortran_env

@brief pFON Implementation in SCF Module Author: Alireza Lashkaripour Date: January 2025 Reference paper: https://doi.org/10.1063/1.478177 This subroutine incorporates the Partial Fractional Occupation Number (pFON) method into SCF calculations, ensuring smooth occupation numbers using Fermi-Dirac distribution. It dynamically adjusts temperature and beta factors to enhance SCF convergence, particularly for near-degenerate states.

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: mo_energy(nbf)
integer, intent(in) :: nbf
integer, intent(in) :: nelec
real(kind=dp), intent(inout) :: occ(nbf)
real(kind=dp), intent(in) :: beta_pfon
integer, intent(in) :: scf_type
integer, intent(in) :: nsmear
logical, intent(in), optional :: is_beta
integer, intent(in), optional :: nelec_a
integer, intent(in), optional :: nelec_b

Called by

proc~~pfon_occupations~~CalledByGraph proc~pfon_occupations pfon_occupations proc~scf_driver scf_driver proc~scf_driver->proc~pfon_occupations proc~hf_energy hf_energy proc~hf_energy->proc~scf_driver

Source Code

    subroutine pfon_occupations(mo_energy, nbf, nelec, occ, beta_pfon, scf_type, nsmear, is_beta, nelec_a, nelec_b)
        use precision, only: dp 
        implicit none 

        integer, intent(in) :: nbf 
        integer, intent(in) :: nelec, nsmear
        real(kind=dp), intent(in) :: beta_pfon
        real(kind=dp), intent(in) :: mo_energy(nbf)
        real(kind=dp), intent(inout) :: occ(nbf)
        integer, intent(in) :: scf_type ! 1,2,3 RHF,UHF,ROHF 
        logical, intent(in), optional :: is_beta
        integer, intent(in), optional :: nelec_a, nelec_b
        real(kind=dp) :: eF, sum_occ
        integer :: i, i_homo, i_lumo, i_low, i_high
        real(kind=dp) :: tmp
        logical :: is_beta_calc
        integer :: n_electrons, n_double, n_single

        is_beta_calc = .false.
        if (present(is_beta)) is_beta_calc = is_beta

        select case (scf_type)
        case(1) ! RHF 
            i_homo = max(1, nelec/2)
            n_electrons = nelec

        case(2) ! UHF 
            if (.not. present(nelec_a) .or. .not. present(nelec_b)) then
                stop 'UHF requires nelec_a and nelec_b'
            endif
            ! UHF: completely independent alpha and beta
            if (is_beta_calc) then
                i_homo = max(1, nelec_b)
                n_electrons = nelec_b
            else
                i_homo = max(1, nelec_a)
                n_electrons = nelec_a
            endif

        case(3) ! ROHF 
            if (.not. present(nelec_a) .or. .not. present(nelec_b)) then
                stop 'ROHF requires nelec_a and nelec_b'
            endif
            ! ROHF: same spatial orbitals, different occupations
            n_double = nelec_b           
            n_single = nelec_a - nelec_b 
            if (is_beta_calc) then
                i_homo = n_double
                n_electrons = nelec_b
            else
                i_homo = n_double + n_single
                n_electrons = nelec_a
            endif
        end select 

        i_lumo = i_homo + 1 
        if (i_lumo > nbf) i_lumo = nbf 

        ! Calculate Fermi level
        eF = 0.5_dp * (mo_energy(i_homo) + mo_energy(i_lumo))

        if (nsmear <= 0) then 
            do i = 1, nbf
                tmp = beta_pfon * (mo_energy(i) - eF)
                if (scf_type == 1) then  ! RHF
                    occ(i) = 2.0_dp / (1.0_dp + exp(tmp))
                else  ! UHF or ROHF
                    occ(i) = 1.0_dp / (1.0_dp + exp(tmp))
                endif
            end do 
        else
            i_low = max(1, i_homo - nsmear)
            i_high = min(nbf, i_lumo + nsmear)

            ! Special handling for ROHF
            if (scf_type == 3) then
                if (is_beta_calc) then
                    do i = 1, n_double
                        occ(i) = 1.0_dp
                    end do
                    do i = n_double + 1, nbf
                        occ(i) = 0.0_dp
                    end do
                else
                    do i = 1, n_double
                        occ(i) = 1.0_dp
                    end do
                    do i = n_double + 1, n_double + n_single
                        occ(i) = 1.0_dp
                    end do
                    do i = n_double + n_single + 1, nbf
                        occ(i) = 0.0_dp
                    end do
                endif

                ! Apply smearing only around the Fermi level
                do i = i_low, i_high
                    tmp = beta_pfon * (mo_energy(i) - eF)
                    occ(i) = occ(i) / (1.0_dp + exp(tmp))
                end do
            else
                ! RHF/UHF handling
                do i = 1, i_low - 1 
                    if (scf_type == 1) then 
                        occ(i) = 2.0_dp 
                    else 
                        occ(i) = 1.0_dp
                    endif
                end do 

                do i = i_high + 1, nbf 
                    occ(i) = 0.0_dp 
                end do 

                do i = i_low, i_high
                    tmp = beta_pfon * (mo_energy(i) - eF)
                    if (scf_type == 1) then  
                        occ(i) = 2.0_dp / (1.0_dp + exp(tmp))
                    else 
                        occ(i) = 1.0_dp / (1.0_dp + exp(tmp))
                    endif
                end do 
            endif
        endif 

        ! Normalize occupations
        sum_occ = sum(occ(1:nbf))
        if (sum_occ < 1.0e-14_dp) then
            sum_occ = 1.0_dp 
        endif 
        occ(1:nbf) = occ(1:nbf) * (real(n_electrons,dp) / sum_occ)

    end subroutine pfon_occupations