Line data    Source code

       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Small algebraic helpers for the rotationally covariant low-rank OT preconditioner.
      10              : ! **************************************************************************************************
      11              : MODULE low_rank_preconditioner_model
      12              :    USE kinds,                           ONLY: dp
      13              : 
      14              :    IMPLICIT NONE
      15              :    PRIVATE
      16              : 
      17              :    PUBLIC :: low_rank_inverse_weight, low_rank_select_rank
      18              : 
      19              : CONTAINS
      20              : 
      21              : ! **************************************************************************************************
      22              : !> \brief Spectral inverse weight relative to a common occupied reference level.
      23              : !> \param eigenvalue complementary-state eigenvalue
      24              : !> \param reference common occupied reference level
      25              : !> \param energy_gap lower bound for the denominator
      26              : !> \return inverse spectral weight
      27              : ! **************************************************************************************************
      28          207 :    PURE ELEMENTAL REAL(KIND=dp) FUNCTION low_rank_inverse_weight( &
      29              :       eigenvalue, reference, energy_gap) RESULT(weight)
      30              : 
      31              :       REAL(KIND=dp), INTENT(IN)                          :: eigenvalue, reference, energy_gap
      32              : 
      33              :       REAL(KIND=dp)                                      :: delta
      34              : 
      35          207 :       delta = eigenvalue - reference
      36          207 :       weight = 1.0_dp/MAX(energy_gap, delta)
      37              : 
      38          207 :    END FUNCTION low_rank_inverse_weight
      39              : 
      40              : ! **************************************************************************************************
      41              : !> \brief Select a bounded spectral rank without cutting a degenerate boundary manifold.
      42              : !> \param eigenvalues ordered full-space eigenvalues
      43              : !> \param nocc number of occupied states at the start of eigenvalues
      44              : !> \param max_rank hard upper bound
      45              : !> \param degeneracy_tolerance relative degeneracy threshold
      46              : !> \return selected rank
      47              : ! **************************************************************************************************
      48           50 :    PURE INTEGER FUNCTION low_rank_select_rank(eigenvalues, nocc, max_rank, &
      49              :                                               degeneracy_tolerance) RESULT(rank_used)
      50              : 
      51              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: eigenvalues
      52              :       INTEGER, INTENT(IN)                                :: nocc, max_rank
      53              :       REAL(KIND=dp), INTENT(IN)                          :: degeneracy_tolerance
      54              : 
      55              :       INTEGER                                            :: ncomplement
      56              :       REAL(KIND=dp)                                      :: scale
      57              : 
      58           50 :       rank_used = 0
      59           50 :       ncomplement = SIZE(eigenvalues) - nocc
      60           50 :       IF (nocc < 0 .OR. ncomplement <= 0 .OR. max_rank <= 0) RETURN
      61              : 
      62            6 :       rank_used = MIN(max_rank, ncomplement)
      63              : 
      64              :       ! Exclude the complete boundary cluster if the hard rank cap would split it.
      65            8 :       DO WHILE (rank_used > 0 .AND. rank_used < ncomplement)
      66              :          scale = MAX(1.0_dp, ABS(eigenvalues(nocc + rank_used)), &
      67            4 :                      ABS(eigenvalues(nocc + rank_used + 1)))
      68            4 :          IF (ABS(eigenvalues(nocc + rank_used + 1) - eigenvalues(nocc + rank_used)) > &
      69              :              degeneracy_tolerance*scale) EXIT
      70            4 :          rank_used = rank_used - 1
      71              :       END DO
      72              : 
      73              :    END FUNCTION low_rank_select_rank
      74              : 
      75              : END MODULE low_rank_preconditioner_model