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            2 : PROGRAM low_rank_preconditioner_unittest
       9            2 :    USE kinds,                           ONLY: dp
      10              :    USE low_rank_preconditioner_model,   ONLY: apply_low_rank_dense,&
      11              :                                               low_rank_inverse_weight,&
      12              :                                               low_rank_select_rank
      13              : 
      14              :    IMPLICIT NONE
      15              : 
      16              :    INTEGER, PARAMETER                                  :: nao = 5, nocc = 3, rank = 2
      17              :    REAL(KIND=dp), PARAMETER                            :: tolerance = 2.0E-14_dp
      18              :    INTEGER                                             :: i, selected_rank
      19              :    REAL(KIND=dp)                                       :: angle, energy_gap, mu, window
      20              :    REAL(KIND=dp), DIMENSION(7)                         :: full_eigenvalues
      21              :    REAL(KIND=dp), DIMENSION(rank)                      :: eigenvalues
      22              :    REAL(KIND=dp), DIMENSION(nao, nao)                  :: identity, overlap_inverse
      23              :    REAL(KIND=dp), DIMENSION(nao, rank)                 :: vectors
      24              :    REAL(KIND=dp), DIMENSION(nao, nocc)                 :: gradient, gradient_rotated, RESULT, &
      25              :                                                           result_rotated
      26              :    REAL(KIND=dp), DIMENSION(nocc, nocc)                :: rotation
      27              : 
      28            2 :    identity = 0.0_dp
      29           12 :    DO i = 1, nao
      30           12 :       identity(i, i) = 1.0_dp
      31              :    END DO
      32              : 
      33            2 :    overlap_inverse = identity
      34            2 :    overlap_inverse(1, 1) = 1.4_dp
      35            2 :    overlap_inverse(2, 2) = 0.8_dp
      36            2 :    overlap_inverse(1, 2) = 0.1_dp
      37            2 :    overlap_inverse(2, 1) = 0.1_dp
      38              : 
      39            2 :    vectors = 0.0_dp
      40            2 :    vectors(3, 1) = 1.0_dp
      41            2 :    vectors(4, 2) = 0.8_dp
      42            2 :    vectors(5, 2) = 0.6_dp
      43            2 :    eigenvalues = [0.35_dp, 0.8_dp]
      44            2 :    mu = -0.2_dp
      45            2 :    energy_gap = 0.08_dp
      46            2 :    window = 1.0_dp
      47              : 
      48              :    gradient = RESHAPE([0.2_dp, -0.1_dp, 0.4_dp, 0.3_dp, -0.2_dp, &
      49              :                        0.6_dp, 0.5_dp, -0.3_dp, 0.1_dp, 0.7_dp, &
      50            2 :                        -0.4_dp, 0.8_dp, 0.2_dp, -0.5_dp, 0.9_dp], [nao, nocc])
      51            2 :    angle = 0.37_dp
      52            2 :    rotation = 0.0_dp
      53            2 :    rotation(1, 1) = COS(angle)
      54            2 :    rotation(1, 2) = -SIN(angle)
      55            2 :    rotation(2, 1) = SIN(angle)
      56            2 :    rotation(2, 2) = COS(angle)
      57            2 :    rotation(3, 3) = 1.0_dp
      58          116 :    gradient_rotated = MATMUL(gradient, rotation)
      59              : 
      60              :    CALL apply_low_rank_dense(overlap_inverse, vectors, eigenvalues, mu, energy_gap, window, &
      61            2 :                              gradient, RESULT)
      62              :    CALL apply_low_rank_dense(overlap_inverse, vectors, eigenvalues, mu, energy_gap, window, &
      63            2 :                              gradient_rotated, result_rotated)
      64          152 :    IF (MAXVAL(ABS(result_rotated - MATMUL(RESULT, rotation))) > tolerance) THEN
      65            0 :       ERROR STOP "Low-rank preconditioner is not right-covariant under occupied rotations"
      66              :    END IF
      67              : 
      68            2 :    full_eigenvalues = [-0.8_dp, -0.2_dp, 0.2_dp, 0.4_dp, 0.7_dp, 0.7_dp, 1.4_dp]
      69            2 :    selected_rank = low_rank_select_rank(full_eigenvalues, 2, 3, 1.0E-8_dp)
      70            2 :    IF (selected_rank /= 2) ERROR STOP "Rank cap split a degenerate complementary manifold"
      71            2 :    selected_rank = low_rank_select_rank(full_eigenvalues, 2, 5, 1.0E-8_dp)
      72            2 :    IF (selected_rank /= 5) ERROR STOP "Full complementary space did not retain its rank"
      73            2 :    IF (ABS(low_rank_inverse_weight(1.5_dp, 0.0_dp, 0.08_dp) - 2.0_dp/3.0_dp) > tolerance) THEN
      74            0 :       ERROR STOP "Common-reference inverse weight differs from the spectral model"
      75              :    END IF
      76              : 
      77            2 : END PROGRAM low_rank_preconditioner_unittest