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 :: apply_low_rank_dense, 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        27428 :    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        27428 :       delta = eigenvalue - reference
      36        27428 :       weight = 1.0_dp/MAX(energy_gap, delta)
      37              : 
      38        27428 :    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           60 :    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           60 :       rank_used = 0
      59           60 :       ncomplement = SIZE(eigenvalues) - nocc
      60           60 :       IF (nocc < 0 .OR. ncomplement <= 0 .OR. max_rank <= 0) RETURN
      61              : 
      62           16 :       rank_used = MIN(max_rank, ncomplement)
      63              : 
      64              :       ! Exclude the complete boundary cluster if the hard rank cap would split it.
      65           18 :       DO WHILE (rank_used > 0 .AND. rank_used < ncomplement)
      66              :          scale = MAX(1.0_dp, ABS(eigenvalues(nocc + rank_used)), &
      67            8 :                      ABS(eigenvalues(nocc + rank_used + 1)))
      68            8 :          IF (ABS(eigenvalues(nocc + rank_used + 1) - eigenvalues(nocc + rank_used)) > &
      69              :              degeneracy_tolerance*scale) EXIT
      70            8 :          rank_used = rank_used - 1
      71              :       END DO
      72              : 
      73              :    END FUNCTION low_rank_select_rank
      74              : 
      75              : ! **************************************************************************************************
      76              : !> \brief Dense reference application used to verify right covariance under orbital rotations.
      77              : !> \param overlap_inverse inverse overlap matrix
      78              : !> \param vectors retained complementary-state vectors
      79              : !> \param eigenvalues retained complementary-state eigenvalues
      80              : !> \param reference common occupied reference level
      81              : !> \param energy_gap lower bound for spectral denominators
      82              : !> \param spectral_window overlap-inverse replacement window
      83              : !> \param matrix_in input orbital-gradient matrix
      84              : !> \param matrix_out preconditioned matrix
      85              : ! **************************************************************************************************
      86            4 :    PURE SUBROUTINE apply_low_rank_dense(overlap_inverse, vectors, eigenvalues, reference, &
      87            4 :                                         energy_gap, spectral_window, matrix_in, matrix_out)
      88              : 
      89              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: overlap_inverse, vectors
      90              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: eigenvalues
      91              :       REAL(KIND=dp), INTENT(IN)                          :: reference, energy_gap, spectral_window
      92              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: matrix_in
      93              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: matrix_out
      94              : 
      95              :       INTEGER                                            :: i
      96              :       REAL(KIND=dp)                                      :: correction
      97            8 :       REAL(KIND=dp), DIMENSION(SIZE(matrix_in, 2))       :: projections
      98              : 
      99          520 :       matrix_out = MATMUL(overlap_inverse, matrix_in)/spectral_window
     100           12 :       DO i = 1, SIZE(eigenvalues)
     101              :          correction = low_rank_inverse_weight(eigenvalues(i), reference, energy_gap) - &
     102            8 :                       1.0_dp/spectral_window
     103          176 :          projections(:) = MATMUL(vectors(:, i), matrix_in)
     104              :          matrix_out = matrix_out + correction*SPREAD(vectors(:, i), 2, SIZE(matrix_in, 2))* &
     105          156 :                       SPREAD(projections, 1, SIZE(matrix_in, 1))
     106              :       END DO
     107              : 
     108            4 :    END SUBROUTINE apply_low_rank_dense
     109              : 
     110              : END MODULE low_rank_preconditioner_model