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 Test mesh-wide rank validation and unitary covariance of AO trial selection.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 6 : PROGRAM wannier90_projections_unittest
13 : #if defined (__HAS_IEEE_EXCEPTIONS)
14 : USE ieee_exceptions, ONLY: IEEE_ALL, IEEE_DIVIDE_BY_ZERO, IEEE_INVALID, IEEE_OVERFLOW, &
15 : ieee_get_halting_mode, ieee_set_halting_mode, ieee_support_halting
16 : #endif
17 2 : USE kinds, ONLY: dp
18 : USE message_passing, ONLY: mp_comm_type, &
19 : mp_para_env_release, &
20 : mp_para_env_type, &
21 : mp_world_finalize, &
22 : mp_world_init
23 : USE wannier90_projections, ONLY: wannier90_hybrid_candidates, &
24 : wannier90_window_mask, &
25 : wannier90_projection_rank, &
26 : wannier90_select_projections
27 : #include "./base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 :
31 2 : COMPLEX(KIND=dp), ALLOCATABLE :: candidates(:, :, :)
32 2 : COMPLEX(KIND=dp), ALLOCATABLE :: hybrids(:, :, :), rotated_hybrids(:, :, :)
33 : COMPLEX(KIND=dp) :: quartet(2, 4)
34 : COMPLEX(KIND=dp) :: amn(2, 2, 2), reference(2, 2, 2), rotation(2, 2, 2)
35 : COMPLEX(KIND=dp) :: repaired(2, 2, 3)
36 : COMPLEX(KIND=dp) :: mixed_reference(4, 3, 2)
37 : LOGICAL :: allowed(4)
38 : INTEGER :: anchor, ik, nfail, pivots(2), reference_pivots(2)
39 : REAL(KIND=dp) :: angle, min_svalue
40 : TYPE(mp_comm_type) :: comm
41 : TYPE(mp_para_env_type), POINTER :: para_env
42 : #if defined (__HAS_IEEE_EXCEPTIONS)
43 : LOGICAL :: saved_halt(SIZE(IEEE_ALL)), expected_halt(SIZE(IEEE_ALL)), actual_halt(SIZE(IEEE_ALL))
44 : #endif
45 :
46 2 : CALL mp_world_init(comm)
47 2 : ALLOCATE (para_env)
48 2 : CALL para_env%from_dup(comm)
49 2 : nfail = 0
50 : #if defined (__HAS_IEEE_EXCEPTIONS)
51 : CALL ieee_get_halting_mode(IEEE_ALL, saved_halt)
52 : IF (ieee_support_halting(IEEE_INVALID)) CALL ieee_set_halting_mode(IEEE_INVALID, .TRUE.)
53 : IF (ieee_support_halting(IEEE_DIVIDE_BY_ZERO)) CALL ieee_set_halting_mode(IEEE_DIVIDE_BY_ZERO, .TRUE.)
54 : IF (ieee_support_halting(IEEE_OVERFLOW)) CALL ieee_set_halting_mode(IEEE_OVERFLOW, .TRUE.)
55 : CALL ieee_get_halting_mode(IEEE_ALL, expected_halt)
56 : #endif
57 2 : IF (para_env%is_source()) THEN
58 1 : ALLOCATE (candidates(2, 3, 2))
59 1 : candidates = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
60 : ! The strongest stacked column vanishes at the second k-point.
61 1 : candidates(1, 1, 1) = 5.0_dp
62 3 : candidates(2, 2, :) = 2.0_dp
63 3 : candidates(1, 3, :) = 0.4_dp
64 : ELSE
65 1 : ALLOCATE (candidates(0, 0, 0))
66 : END IF
67 2 : CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue, anchor_kpoint=0)
68 2 : IF (min_svalue > 1.0e-14_dp) nfail = nfail + 1
69 2 : CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
70 2 : IF (ABS(min_svalue - 0.4_dp) > 1.0e-14_dp) nfail = nfail + 1
71 6 : IF (.NOT. ALL(pivots == [2, 3])) nfail = nfail + 1
72 2 : reference = amn
73 2 : reference_pivots = pivots
74 :
75 6 : DO ik = 1, 2
76 4 : angle = 0.31_dp*REAL(ik, KIND=dp)
77 4 : rotation(1, 1, ik) = COS(angle)
78 4 : rotation(2, 2, ik) = COS(angle)
79 4 : rotation(1, 2, ik) = CMPLX(0.0_dp, SIN(angle), KIND=dp)
80 4 : rotation(2, 1, ik) = CMPLX(0.0_dp, SIN(angle), KIND=dp)
81 108 : IF (para_env%is_source()) candidates(:, :, ik) = MATMUL(rotation(:, :, ik), candidates(:, :, ik))
82 : END DO
83 2 : CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
84 2 : IF (ABS(min_svalue - 0.4_dp) > 1.0e-14_dp) nfail = nfail + 1
85 6 : IF (.NOT. ALL(pivots == reference_pivots)) nfail = nfail + 1
86 6 : DO ik = 1, 2
87 86 : IF (MAXVAL(ABS(amn(:, :, ik) - MATMUL(rotation(:, :, ik), reference(:, :, ik)))) > 1.0e-14_dp) THEN
88 0 : nfail = nfail + 1
89 : END IF
90 : END DO
91 2 : IF (para_env%is_source()) candidates(:, :, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
92 2 : CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
93 2 : IF (min_svalue /= 0.0_dp) nfail = nfail + 1
94 2 : DEALLOCATE (candidates)
95 2 : IF (para_env%is_source()) THEN
96 1 : ALLOCATE (candidates(2, 5, 2))
97 3 : DO ik = 1, 2
98 13 : DO anchor = 1, 5
99 10 : candidates(1, anchor, ik) = CMPLX(0.1_dp*anchor, 0.2_dp*ik, KIND=dp)
100 12 : candidates(2, anchor, ik) = CMPLX(-0.3_dp*anchor, 0.05_dp*ik*anchor, KIND=dp)
101 : END DO
102 : END DO
103 : ELSE
104 1 : ALLOCATE (candidates(0, 0, 0))
105 : END IF
106 2 : CALL wannier90_hybrid_candidates(candidates, RESHAPE([1, 2, 3, 4], [4, 1]), hybrids)
107 2 : IF (para_env%is_source()) THEN
108 4 : IF (.NOT. ALL(SHAPE(hybrids) == [2, 5, 2])) nfail = nfail + 1
109 7 : IF (MAXVAL(ABS(hybrids(:, 1, :) - candidates(:, 5, :))) > 1.0e-14_dp) nfail = nfail + 1
110 3 : DO ik = 1, 2
111 : ! A full tetrahedral quartet is a unitary change of the four AO trials.
112 26 : quartet = hybrids(:, 2:5, ik)
113 14 : IF (MAXVAL(ABS(MATMUL(quartet, CONJG(TRANSPOSE(quartet))) - &
114 8 : MATMUL(candidates(:, 1:4, ik), CONJG(TRANSPOSE(candidates(:, 1:4, ik)))))) &
115 108 : > 1.0e-14_dp) nfail = nfail + 1
116 173 : candidates(:, :, ik) = MATMUL(rotation(:, :, ik), candidates(:, :, ik))
117 : END DO
118 : END IF
119 2 : CALL wannier90_hybrid_candidates(candidates, RESHAPE([1, 2, 3, 4], [4, 1]), rotated_hybrids)
120 2 : IF (para_env%is_source()) THEN
121 3 : DO ik = 1, 2
122 38 : IF (MAXVAL(ABS(rotated_hybrids(:, :, ik) - MATMUL(rotation(:, :, ik), hybrids(:, :, ik)))) &
123 115 : > 1.0e-14_dp) nfail = nfail + 1
124 : END DO
125 : ELSE
126 7 : IF (SIZE(hybrids) /= 0 .OR. SIZE(rotated_hybrids) /= 0) nfail = nfail + 1
127 : END IF
128 2 : DEALLOCATE (hybrids, rotated_hybrids)
129 2 : CALL wannier90_hybrid_candidates(candidates, RESHAPE([1, 2, 3, 4, 5, 2, 3, 4], [4, 2]), hybrids)
130 2 : IF (para_env%is_source()) THEN
131 4 : IF (.NOT. ALL(SHAPE(hybrids) == [2, 8, 2])) nfail = nfail + 1
132 3 : DO ik = 1, 2
133 22 : IF (MAXVAL(ABS(SUM(hybrids(:, 1:4, ik), DIM=2)*0.5_dp - candidates(:, 1, ik))) &
134 0 : > 1.0e-14_dp) nfail = nfail + 1
135 22 : IF (MAXVAL(ABS(SUM(hybrids(:, 5:8, ik), DIM=2)*0.5_dp - candidates(:, 5, ik))) &
136 1 : > 1.0e-14_dp) nfail = nfail + 1
137 : END DO
138 : END IF
139 2 : DEALLOCATE (hybrids)
140 2 : CALL wannier90_hybrid_candidates(candidates, RESHAPE([1, 2, 3, 4], [4, 0]), hybrids)
141 2 : IF (para_env%is_source()) THEN
142 33 : IF (MAXVAL(ABS(hybrids - candidates)) > 1.0e-14_dp) nfail = nfail + 1
143 : END IF
144 2 : DEALLOCATE (hybrids, candidates)
145 2 : IF (para_env%is_source()) THEN
146 1 : ALLOCATE (candidates(2, 4, 3))
147 1 : candidates = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
148 : ! Every QR anchor chooses at least one column that vanishes elsewhere.
149 : ! Exchanging both of those columns recovers the fixed full-rank pair.
150 1 : candidates(1, 1, 1) = 10.0_dp
151 1 : candidates(1, 1, 3) = 10.0_dp
152 1 : candidates(2, 2, 2) = 10.0_dp
153 1 : candidates(2, 2, 3) = 10.0_dp
154 4 : candidates(1, 3, :) = 1.0_dp
155 4 : candidates(2, 4, :) = 1.0_dp
156 : ELSE
157 1 : ALLOCATE (candidates(0, 0, 0))
158 : END IF
159 10 : DO anchor = 0, 3
160 8 : CALL wannier90_select_projections(candidates, para_env, repaired, pivots, min_svalue, anchor)
161 10 : IF (min_svalue > 1.0e-14_dp) nfail = nfail + 1
162 : END DO
163 2 : CALL wannier90_select_projections(candidates, para_env, repaired, pivots, min_svalue)
164 2 : IF (ABS(min_svalue - 1.0_dp) > 1.0e-14_dp) nfail = nfail + 1
165 6 : IF (.NOT. ALL(pivots == [3, 4])) nfail = nfail + 1
166 2 : DEALLOCATE (candidates)
167 2 : IF (para_env%is_source()) THEN
168 1 : ALLOCATE (candidates(2, 3, 2))
169 1 : candidates = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
170 3 : candidates(1, 1, :) = 1.0_dp
171 3 : candidates(1, 2, :) = 1.0_dp + 1.0e-14_dp
172 3 : candidates(2, 3, :) = 1.0_dp
173 : ELSE
174 1 : ALLOCATE (candidates(0, 0, 0))
175 : END IF
176 2 : CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
177 6 : IF (.NOT. ALL(pivots == [1, 3])) nfail = nfail + 1
178 2 : IF (para_env%is_source()) THEN
179 3 : DO ik = 1, 2
180 105 : candidates(:, :, ik) = MATMUL(rotation(:, :, ik), candidates(:, :, ik))
181 : END DO
182 : END IF
183 2 : CALL wannier90_select_projections(candidates, para_env, amn, pivots, min_svalue)
184 6 : IF (.NOT. ALL(pivots == [1, 3])) nfail = nfail + 1
185 2 : DEALLOCATE (candidates)
186 6 : DO ik = 1, 2
187 18 : DO anchor = 1, 3
188 12 : mixed_reference(1, anchor, ik) = CMPLX(0.1_dp*anchor, 0.2_dp*ik, KIND=dp)
189 12 : mixed_reference(2, anchor, ik) = CMPLX(-0.3_dp*anchor, 0.05_dp*ik*anchor, KIND=dp)
190 12 : mixed_reference(3, anchor, ik) = CMPLX(0.2_dp*anchor**2, 0.3_dp*ik, KIND=dp)
191 16 : mixed_reference(4, anchor, ik) = CMPLX(0.7_dp, -0.1_dp*ik*anchor, KIND=dp)
192 : END DO
193 : END DO
194 2 : min_svalue = wannier90_projection_rank(mixed_reference)
195 : #if defined (__HAS_IEEE_EXCEPTIONS)
196 : CALL ieee_get_halting_mode(IEEE_ALL, actual_halt)
197 : IF (ANY(actual_halt .NEQV. expected_halt)) nfail = nfail + 1
198 : ! Also check that a caller with traps disabled is left unchanged.
199 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
200 : IF (ABS(wannier90_projection_rank(mixed_reference) - min_svalue) > 1.0e-14_dp) nfail = nfail + 1
201 : CALL ieee_get_halting_mode(IEEE_ALL, actual_halt)
202 : IF (ANY(actual_halt)) nfail = nfail + 1
203 : CALL ieee_set_halting_mode(IEEE_ALL, expected_halt)
204 : #endif
205 : CALL wannier90_window_mask([-2.0_dp, -1.0_dp, 1.0_dp, 2.0_dp], 2, &
206 : [-1.5_dp, 1.5_dp, 0.0_dp, 0.0_dp], &
207 2 : [.TRUE., .TRUE., .FALSE., .FALSE.], allowed)
208 10 : IF (ANY(allowed .NEQV. [.FALSE., .TRUE., .TRUE., .FALSE.])) nfail = nfail + 1
209 : CALL wannier90_window_mask([-2.0_dp, -1.0_dp, 1.0_dp, 2.0_dp], 2, &
210 : [-3.0_dp, 3.0_dp, 0.0_dp, 0.0_dp], &
211 2 : [.TRUE., .TRUE., .FALSE., .TRUE.], allowed)
212 10 : IF (ANY(allowed .NEQV. [.TRUE., .TRUE., .FALSE., .FALSE.])) nfail = nfail + 1
213 : #if defined (__HAS_IEEE_EXCEPTIONS)
214 : CALL ieee_get_halting_mode(IEEE_ALL, actual_halt)
215 : IF (ANY(actual_halt .NEQV. expected_halt)) nfail = nfail + 1
216 : CALL ieee_set_halting_mode(IEEE_ALL, saved_halt)
217 : #endif
218 2 : CALL comm%sum(nfail)
219 2 : CPASSERT(nfail == 0)
220 2 : CALL mp_para_env_release(para_env)
221 2 : CALL mp_world_finalize()
222 2 : END PROGRAM wannier90_projections_unittest
|
