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 : USE mathconstants, ONLY: sqrthalf
14 :
15 : IMPLICIT NONE
16 :
17 : INTEGER, PARAMETER :: nao = 5, nocc = 3, rank = 2
18 : REAL(KIND=dp), PARAMETER :: tolerance = 2.0E-14_dp
19 : INTEGER :: i, selected_rank
20 : REAL(KIND=dp) :: angle, common_shift, energy_gap, mu, window
21 : REAL(KIND=dp), DIMENSION(7) :: full_eigenvalues
22 : REAL(KIND=dp), DIMENSION(rank) :: eigenvalues
23 : REAL(KIND=dp), DIMENSION(nao, nao) :: complement_operator, &
24 : complement_operator_rotated, hamiltonian, &
25 : identity, overlap_inverse, projector, &
26 : projector_rotated
27 : REAL(KIND=dp), DIMENSION(nao, 2) :: occupied, occupied_rotated, vectors
28 : REAL(KIND=dp), DIMENSION(nao, nocc) :: gradient, gradient_rotated, RESULT, &
29 : result_rotated
30 : REAL(KIND=dp), DIMENSION(nocc, nocc) :: rotation
31 : REAL(KIND=dp), DIMENSION(2, 2) :: occupied_rotation
32 :
33 2 : identity = 0.0_dp
34 12 : DO i = 1, nao
35 12 : identity(i, i) = 1.0_dp
36 : END DO
37 :
38 2 : overlap_inverse = identity
39 2 : overlap_inverse(1, 1) = 1.4_dp
40 2 : overlap_inverse(2, 2) = 0.8_dp
41 2 : overlap_inverse(1, 2) = 0.1_dp
42 2 : overlap_inverse(2, 1) = 0.1_dp
43 :
44 2 : vectors = 0.0_dp
45 2 : vectors(3, 1) = 1.0_dp
46 2 : vectors(4, 2) = 0.8_dp
47 2 : vectors(5, 2) = 0.6_dp
48 2 : eigenvalues = [0.35_dp, 0.8_dp]
49 2 : mu = -0.2_dp
50 2 : energy_gap = 0.08_dp
51 2 : window = 1.0_dp
52 :
53 : gradient = RESHAPE([0.2_dp, -0.1_dp, 0.4_dp, 0.3_dp, -0.2_dp, &
54 : 0.6_dp, 0.5_dp, -0.3_dp, 0.1_dp, 0.7_dp, &
55 2 : -0.4_dp, 0.8_dp, 0.2_dp, -0.5_dp, 0.9_dp], [nao, nocc])
56 2 : angle = 0.37_dp
57 2 : rotation = 0.0_dp
58 2 : rotation(1, 1) = COS(angle)
59 2 : rotation(1, 2) = -SIN(angle)
60 2 : rotation(2, 1) = SIN(angle)
61 2 : rotation(2, 2) = COS(angle)
62 2 : rotation(3, 3) = 1.0_dp
63 116 : gradient_rotated = MATMUL(gradient, rotation)
64 :
65 : CALL apply_low_rank_dense(overlap_inverse, vectors, eigenvalues, mu, energy_gap, window, &
66 2 : gradient, RESULT)
67 : CALL apply_low_rank_dense(overlap_inverse, vectors, eigenvalues, mu, energy_gap, window, &
68 2 : gradient_rotated, result_rotated)
69 152 : IF (MAXVAL(ABS(result_rotated - MATMUL(RESULT, rotation))) > tolerance) THEN
70 0 : ERROR STOP "Low-rank preconditioner is not right-covariant under occupied rotations"
71 : END IF
72 :
73 : ! Verify that the projected complementary operator depends only on the occupied subspace.
74 2 : occupied = 0.0_dp
75 2 : occupied(1, 1) = sqrthalf
76 2 : occupied(2, 1) = sqrthalf
77 2 : occupied(3, 2) = sqrthalf
78 2 : occupied(4, 2) = sqrthalf
79 10 : occupied_rotation = RESHAPE([COS(angle), SIN(angle), -SIN(angle), COS(angle)], [2, 2])
80 54 : occupied_rotated = MATMUL(occupied, occupied_rotation)
81 :
82 : hamiltonian = RESHAPE([1.0_dp, 0.2_dp, 0.1_dp, 0.0_dp, 0.3_dp, &
83 : 0.2_dp, 1.4_dp, 0.0_dp, 0.2_dp, 0.1_dp, &
84 : 0.1_dp, 0.0_dp, 2.0_dp, 0.4_dp, 0.2_dp, &
85 : 0.0_dp, 0.2_dp, 0.4_dp, 2.3_dp, 0.1_dp, &
86 2 : 0.3_dp, 0.1_dp, 0.2_dp, 0.1_dp, 3.0_dp], [nao, nao])
87 192 : projector = identity - MATMUL(occupied, TRANSPOSE(occupied))
88 192 : projector_rotated = identity - MATMUL(occupied_rotated, TRANSPOSE(occupied_rotated))
89 2 : common_shift = -10.2_dp
90 : complement_operator = MATMUL(TRANSPOSE(projector), MATMUL(hamiltonian, projector)) + &
91 812 : common_shift*MATMUL(occupied, TRANSPOSE(occupied))
92 : complement_operator_rotated = &
93 : MATMUL(TRANSPOSE(projector_rotated), MATMUL(hamiltonian, projector_rotated)) + &
94 812 : common_shift*MATMUL(occupied_rotated, TRANSPOSE(occupied_rotated))
95 62 : IF (MAXVAL(ABS(complement_operator_rotated - complement_operator)) > tolerance) THEN
96 0 : ERROR STOP "Complementary-state construction depends on the occupied orbital gauge"
97 : END IF
98 :
99 2 : full_eigenvalues = [-0.8_dp, -0.2_dp, 0.2_dp, 0.4_dp, 0.7_dp, 0.7_dp, 1.4_dp]
100 2 : selected_rank = low_rank_select_rank(full_eigenvalues, 2, 3, 1.0E-8_dp)
101 2 : IF (selected_rank /= 2) ERROR STOP "Rank cap split a degenerate complementary manifold"
102 2 : selected_rank = low_rank_select_rank(full_eigenvalues, 2, 5, 1.0E-8_dp)
103 2 : IF (selected_rank /= 5) ERROR STOP "Full complementary space did not retain its rank"
104 2 : IF (ABS(low_rank_inverse_weight(1.5_dp, 0.0_dp, 0.08_dp) - 2.0_dp/3.0_dp) > tolerance) THEN
105 0 : ERROR STOP "Common-reference inverse weight differs from the spectral model"
106 : END IF
107 :
108 2 : END PROGRAM low_rank_preconditioner_unittest
|
