low_rank_preconditioner_unittest.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
PROGRAM low_rank_preconditioner_unittest
   USE kinds,                           ONLY: dp
   USE low_rank_preconditioner_model,   ONLY: apply_low_rank_dense,&
                                              low_rank_inverse_weight,&
                                              low_rank_select_rank
   USE mathconstants,                   ONLY: sqrthalf
 
   IMPLICIT NONE
 
   INTEGER, PARAMETER                                  :: nao = 5, nocc = 3, rank = 2
   REAL(kind=dp), PARAMETER                            :: tolerance = 2.0e-14_dp
   INTEGER                                             :: i, selected_rank
   REAL(kind=dp)                                       :: angle, common_shift, energy_gap, mu, window
   REAL(kind=dp), DIMENSION(7)                         :: full_eigenvalues
   REAL(kind=dp), DIMENSION(rank)                      :: eigenvalues
   REAL(kind=dp), DIMENSION(nao, nao)                  :: complement_operator, &
                                                          complement_operator_rotated, hamiltonian, &
                                                          identity, overlap_inverse, projector, &
                                                          projector_rotated
   REAL(kind=dp), DIMENSION(nao, 2)                    :: occupied, occupied_rotated, vectors
   REAL(kind=dp), DIMENSION(nao, nocc)                 :: gradient, gradient_rotated, result, &
                                                          result_rotated
   REAL(kind=dp), DIMENSION(nocc, nocc)                :: rotation
   REAL(kind=dp), DIMENSION(2, 2)                      :: occupied_rotation
 
   identity = 0.0_dp
   DO i = 1, nao
      identity(i, i) = 1.0_dp
   END DO
 
   overlap_inverse = identity
   overlap_inverse(1, 1) = 1.4_dp
   overlap_inverse(2, 2) = 0.8_dp
   overlap_inverse(1, 2) = 0.1_dp
   overlap_inverse(2, 1) = 0.1_dp
 
   vectors = 0.0_dp
   vectors(3, 1) = 1.0_dp
   vectors(4, 2) = 0.8_dp
   vectors(5, 2) = 0.6_dp
   eigenvalues = [0.35_dp, 0.8_dp]
   mu = -0.2_dp
   energy_gap = 0.08_dp
   window = 1.0_dp
 
   gradient = reshape([0.2_dp, -0.1_dp, 0.4_dp, 0.3_dp, -0.2_dp, &
                       0.6_dp, 0.5_dp, -0.3_dp, 0.1_dp, 0.7_dp, &
                       -0.4_dp, 0.8_dp, 0.2_dp, -0.5_dp, 0.9_dp], [nao, nocc])
   angle = 0.37_dp
   rotation = 0.0_dp
   rotation(1, 1) = cos(angle)
   rotation(1, 2) = -sin(angle)
   rotation(2, 1) = sin(angle)
   rotation(2, 2) = cos(angle)
   rotation(3, 3) = 1.0_dp
   gradient_rotated = matmul(gradient, rotation)
 
   CALL apply_low_rank_dense(overlap_inverse, vectors, eigenvalues, mu, energy_gap, window, &
                             gradient, result)
   CALL apply_low_rank_dense(overlap_inverse, vectors, eigenvalues, mu, energy_gap, window, &
                             gradient_rotated, result_rotated)
   IF (maxval(abs(result_rotated - matmul(result, rotation))) > tolerance) THEN
      error stop "Low-rank preconditioner is not right-covariant under occupied rotations"
   END IF
 
   ! Verify that the projected complementary operator depends only on the occupied subspace.
   occupied = 0.0_dp
   occupied(1, 1) = sqrthalf
   occupied(2, 1) = sqrthalf
   occupied(3, 2) = sqrthalf
   occupied(4, 2) = sqrthalf
   occupied_rotation = reshape([cos(angle), sin(angle), -sin(angle), cos(angle)], [2, 2])
   occupied_rotated = matmul(occupied, occupied_rotation)
 
   hamiltonian = reshape([1.0_dp, 0.2_dp, 0.1_dp, 0.0_dp, 0.3_dp, &
                          0.2_dp, 1.4_dp, 0.0_dp, 0.2_dp, 0.1_dp, &
                          0.1_dp, 0.0_dp, 2.0_dp, 0.4_dp, 0.2_dp, &
                          0.0_dp, 0.2_dp, 0.4_dp, 2.3_dp, 0.1_dp, &
                          0.3_dp, 0.1_dp, 0.2_dp, 0.1_dp, 3.0_dp], [nao, nao])
   projector = identity - matmul(occupied, transpose(occupied))
   projector_rotated = identity - matmul(occupied_rotated, transpose(occupied_rotated))
   common_shift = -10.2_dp
   complement_operator = matmul(transpose(projector), matmul(hamiltonian, projector)) + &
                         common_shift*matmul(occupied, transpose(occupied))
   complement_operator_rotated = &
      matmul(transpose(projector_rotated), matmul(hamiltonian, projector_rotated)) + &
      common_shift*matmul(occupied_rotated, transpose(occupied_rotated))
   IF (maxval(abs(complement_operator_rotated - complement_operator)) > tolerance) THEN
      error stop "Complementary-state construction depends on the occupied orbital gauge"
   END IF
 
   full_eigenvalues = [-0.8_dp, -0.2_dp, 0.2_dp, 0.4_dp, 0.7_dp, 0.7_dp, 1.4_dp]
   selected_rank = low_rank_select_rank(full_eigenvalues, 2, 3, 1.0e-8_dp)
   IF (selected_rank /= 2) error stop "Rank cap split a degenerate complementary manifold"
   selected_rank = low_rank_select_rank(full_eigenvalues, 2, 5, 1.0e-8_dp)
   IF (selected_rank /= 5) error stop "Full complementary space did not retain its rank"
   IF (abs(low_rank_inverse_weight(1.5_dp, 0.0_dp, 0.08_dp) - 2.0_dp/3.0_dp) > tolerance) THEN
      error stop "Common-reference inverse weight differs from the spectral model"
   END IF
 
END PROGRAM low_rank_preconditioner_unittest