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 Distributed block-sparse finite localizers, retaining the covariant AO metric.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE spectral_localizer_dbcsr
13 : USE cp_dbcsr_api, ONLY: &
14 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_get, &
15 : dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
16 : dbcsr_finalize, dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
17 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
18 : dbcsr_put_block, dbcsr_release, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
19 : USE ieee_arithmetic, ONLY: ieee_is_finite
20 : USE kinds, ONLY: dp
21 : USE message_passing, ONLY: mp_para_env_type
22 : #include "./base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 : PRIVATE
26 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'spectral_localizer_dbcsr'
27 : PUBLIC :: localizer_dbcsr_build, localizer_dbcsr_dense, localizer_dbcsr_release, localizer_matrix_type
28 : PUBLIC :: localizer_dbcsr_check
29 : PUBLIC :: localizer_ao_operators
30 : PUBLIC :: localizer_dbcsr_realify, localizer_dbcsr_from_real
31 : PUBLIC :: localizer_dbcsr_gap_pencil
32 : PUBLIC :: localizer_dbcsr_3d_build
33 :
34 : TYPE localizer_matrix_type
35 : TYPE(dbcsr_type) :: real_part, imaginary_part
36 : END TYPE localizer_matrix_type
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief Build the 3D AII chiral block in DBCSR; not a Hermitian 2D localizer.
42 : !> \param h Scalar AO Hamiltonian
43 : !> \param s Scalar AO overlap
44 : !> \param positions Three Cartesian AO first moments
45 : !> \param point Cartesian position query
46 : !> \param energy Absolute energy query
47 : !> \param kappa Positive energy/length scale
48 : !> \param block Chiral block, with auxiliary index outside physical up/down spin
49 : !> \param metric Block metric I_aux tensor I_spin tensor S
50 : !> \param soc Three real antisymmetric SOC components
51 : !> \param shifted_hamiltonian Optional already energy-shifted physical-spin Hamiltonian
52 : !> \param periodic_cosines Three cos(G.r) matrices. Then positions are sin(G.r), point is
53 : !> in radians, and kappa is the periodic energy scale eta.
54 : ! **************************************************************************************************
55 6 : SUBROUTINE localizer_dbcsr_3d_build(h, s, positions, point, energy, kappa, block, metric, soc, &
56 6 : shifted_hamiltonian, periodic_cosines)
57 : TYPE(dbcsr_type), INTENT(IN) :: h, s
58 : TYPE(dbcsr_p_type), INTENT(IN) :: positions(:)
59 : REAL(dp), INTENT(IN) :: point(3), energy, kappa
60 : TYPE(localizer_matrix_type), INTENT(OUT) :: block, metric
61 : TYPE(dbcsr_p_type), INTENT(IN) :: soc(:, :)
62 : TYPE(localizer_matrix_type), INTENT(IN), OPTIONAL :: shifted_hamiltonian
63 : TYPE(dbcsr_p_type), INTENT(IN), OPTIONAL :: periodic_cosines(:)
64 :
65 : COMPLEX(dp), PARAMETER :: ii = CMPLX(0.0_dp, 1.0_dp, dp)
66 :
67 : COMPLEX(dp) :: factors(2, 2), identity(4, 4), &
68 : pauli(4, 4, 3)
69 : INTEGER :: i
70 12 : TYPE(localizer_matrix_type) :: physical(1), physical_metric
71 :
72 0 : CPASSERT(SIZE(positions) == 3)
73 6 : identity = 0.0_dp
74 6 : pauli = 0.0_dp
75 30 : DO i = 1, 4
76 30 : identity(i, i) = 1.0_dp
77 : END DO
78 18 : DO i = 1, 2
79 12 : pauli(i, i + 2, 1) = 1.0_dp
80 12 : pauli(i + 2, i, 1) = 1.0_dp
81 12 : pauli(i, i + 2, 2) = -ii
82 12 : pauli(i + 2, i, 2) = ii
83 12 : pauli(i, i, 3) = 1.0_dp
84 18 : pauli(i + 2, i + 2, 3) = -1.0_dp
85 : END DO
86 6 : CALL create_operator(s, 4, block)
87 6 : CALL create_operator(s, 4, metric)
88 6 : CALL add_real_blocks(s, identity, metric)
89 6 : factors = 0.0_dp
90 6 : factors(1, 1) = 1.0_dp
91 6 : factors(2, 2) = 1.0_dp
92 6 : IF (PRESENT(shifted_hamiltonian)) THEN
93 2 : CPASSERT(energy == 0.0_dp)
94 2 : CALL add_real_blocks(shifted_hamiltonian%real_part, factors, block)
95 14 : CALL add_real_blocks(shifted_hamiltonian%imaginary_part, ii*factors, block)
96 : ELSE
97 4 : CALL localizer_ao_operators(h, s, positions(:0), physical, physical_metric, soc)
98 4 : CALL add_real_blocks(physical(1)%real_part, factors, block)
99 28 : CALL add_real_blocks(physical(1)%imaginary_part, ii*factors, block)
100 4 : CALL localizer_dbcsr_release(physical(1))
101 4 : CALL localizer_dbcsr_release(physical_metric)
102 84 : CALL add_real_blocks(s, -energy*identity, block)
103 : END IF
104 6 : IF (PRESENT(periodic_cosines)) THEN
105 4 : CPASSERT(SIZE(periodic_cosines) == 3)
106 84 : CALL add_real_blocks(s, -3.0_dp*kappa*identity, block)
107 16 : DO i = 1, 3
108 : CALL add_real_blocks(periodic_cosines(i)%matrix, &
109 252 : kappa*(COS(point(i))*identity + ii*SIN(point(i))*pauli(:, :, i)), block)
110 : CALL add_real_blocks(positions(i)%matrix, &
111 256 : kappa*(SIN(point(i))*identity - ii*COS(point(i))*pauli(:, :, i)), block)
112 : END DO
113 : ELSE
114 8 : DO i = 1, 3
115 126 : CALL add_real_blocks(positions(i)%matrix, -ii*kappa*pauli(:, :, i), block)
116 128 : CALL add_real_blocks(s, ii*kappa*point(i)*pauli(:, :, i), block)
117 : END DO
118 : END IF
119 6 : END SUBROUTINE localizer_dbcsr_3d_build
120 :
121 : ! **************************************************************************************************
122 : !> \brief Distributed real matrix R(A)=[Re(A),-Im(A);Im(A),Re(A)].
123 : !> \param matrix Complex operator with full-storage components
124 : !> \param realified Newly initialized real DBCSR operator
125 : ! **************************************************************************************************
126 134 : SUBROUTINE localizer_dbcsr_realify(matrix, realified)
127 : TYPE(localizer_matrix_type), INTENT(IN) :: matrix
128 : TYPE(dbcsr_type), INTENT(OUT) :: realified
129 :
130 : COMPLEX(dp) :: factors(2, 2)
131 : TYPE(localizer_matrix_type) :: expanded
132 :
133 134 : CALL create_operator(matrix%real_part, 2, expanded)
134 134 : factors = RESHAPE([1.0_dp, 0.0_dp, 0.0_dp, 1.0_dp], [2, 2])
135 134 : CALL add_real_blocks(matrix%real_part, factors, expanded)
136 134 : factors = RESHAPE([0.0_dp, 1.0_dp, -1.0_dp, 0.0_dp], [2, 2])
137 134 : CALL add_real_blocks(matrix%imaginary_part, factors, expanded)
138 134 : CALL dbcsr_create(realified, template=expanded%real_part)
139 134 : CALL dbcsr_copy(realified, expanded%real_part)
140 134 : CALL localizer_dbcsr_release(expanded)
141 134 : END SUBROUTINE localizer_dbcsr_realify
142 :
143 : ! **************************************************************************************************
144 : !> \brief Extract the two left quadrants of R(A); does not average away structure errors.
145 : !> \param realified Fully stored real operator with two repeated component distributions
146 : !> \param template Original complex operator defining the physical AO block distribution
147 : !> \param matrix Newly initialized complex operator; caller must check realification residual
148 : ! **************************************************************************************************
149 126 : SUBROUTINE localizer_dbcsr_from_real(realified, template, matrix)
150 : TYPE(dbcsr_type), INTENT(INOUT) :: realified
151 : TYPE(localizer_matrix_type), INTENT(IN) :: template
152 : TYPE(localizer_matrix_type), INTENT(OUT) :: matrix
153 :
154 : INTEGER :: blocks, col, nblocks, row
155 : LOGICAL :: transposed
156 42 : REAL(dp), POINTER :: block(:, :)
157 : TYPE(dbcsr_iterator_type) :: iterator
158 :
159 42 : CALL dbcsr_get_info(template%real_part, nblkrows_total=blocks)
160 42 : CALL dbcsr_get_info(realified, nblkrows_total=nblocks)
161 42 : CPASSERT(nblocks == 2*blocks)
162 42 : CALL create_operator(template%real_part, 1, matrix)
163 42 : CALL dbcsr_iterator_start(iterator, realified)
164 2794 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
165 2752 : CALL dbcsr_iterator_next_block(iterator, row, col, block, transposed=transposed)
166 2752 : CPASSERT(.NOT. transposed)
167 2752 : IF (col > blocks) CYCLE
168 1418 : IF (row <= blocks) THEN
169 692 : CALL dbcsr_put_block(matrix%real_part, row, col, block)
170 : ELSE
171 684 : CALL dbcsr_put_block(matrix%imaginary_part, row - blocks, col, block)
172 : END IF
173 : END DO
174 42 : CALL dbcsr_iterator_stop(iterator)
175 42 : CALL dbcsr_finalize(matrix%real_part)
176 42 : CALL dbcsr_finalize(matrix%imaginary_part)
177 42 : END SUBROUTINE localizer_dbcsr_from_real
178 :
179 : ! **************************************************************************************************
180 : !> \brief Create a complex block-sparse operator with component-major scalar AO ordering.
181 : !> \param source Scalar AO template
182 : !> \param components Number of spin/auxiliary components
183 : !> \param matrix Newly initialized operator
184 : ! **************************************************************************************************
185 398 : SUBROUTINE create_operator(source, components, matrix)
186 : TYPE(dbcsr_type), INTENT(IN) :: source
187 : INTEGER, INTENT(IN) :: components
188 : TYPE(localizer_matrix_type), INTENT(OUT) :: matrix
189 :
190 : INTEGER :: blocks
191 398 : INTEGER, POINTER :: col_dist(:), expanded_cols(:), &
192 398 : expanded_rows(:), expanded_sizes(:), &
193 398 : row_dist(:), sizes(:)
194 : TYPE(dbcsr_distribution_type) :: distribution, expanded_distribution
195 :
196 398 : CALL dbcsr_get_info(source, distribution=distribution, row_blk_size=sizes)
197 398 : CALL dbcsr_distribution_get(distribution, row_dist=row_dist, col_dist=col_dist)
198 398 : blocks = SIZE(sizes)
199 1990 : ALLOCATE (expanded_sizes(blocks*components), expanded_rows(blocks*components), expanded_cols(blocks*components))
200 796 : expanded_sizes = RESHAPE(SPREAD(sizes, 2, components), [blocks*components])
201 796 : expanded_rows = RESHAPE(SPREAD(row_dist, 2, components), [blocks*components])
202 796 : expanded_cols = RESHAPE(SPREAD(col_dist, 2, components), [blocks*components])
203 : CALL dbcsr_distribution_new(expanded_distribution, template=distribution, &
204 398 : row_dist=expanded_rows, col_dist=expanded_cols)
205 : CALL dbcsr_create(matrix%real_part, 'FINITE AO RE', expanded_distribution, &
206 398 : dbcsr_type_no_symmetry, expanded_sizes, expanded_sizes)
207 398 : CALL dbcsr_create(matrix%imaginary_part, template=matrix%real_part, name='FINITE AO IM')
208 398 : CALL dbcsr_distribution_release(expanded_distribution)
209 398 : DEALLOCATE (expanded_sizes, expanded_rows, expanded_cols)
210 796 : END SUBROUTINE create_operator
211 :
212 : ! **************************************************************************************************
213 : !> \brief Shared AO Hamiltonian, Hermitian scalar observables and optional GTH physical-spin SOC.
214 : !> \param h Scalar Hamiltonian
215 : !> \param s AO metric
216 : !> \param positions Covariant scalar observables, such as first moments or trigonometric coordinates
217 : !> \param operators Newly initialized Hamiltonian followed by the supplied observables
218 : !> \param metric Newly initialized physical-spin metric
219 : !> \param soc Optional antisymmetric real GTH SOC components
220 : ! **************************************************************************************************
221 168 : SUBROUTINE localizer_ao_operators(h, s, positions, operators, metric, soc)
222 : TYPE(dbcsr_type), INTENT(IN) :: h, s
223 : TYPE(dbcsr_p_type), INTENT(IN) :: positions(:)
224 : TYPE(localizer_matrix_type), INTENT(OUT) :: operators(:), metric
225 : TYPE(dbcsr_p_type), INTENT(IN), OPTIONAL :: soc(:, :)
226 :
227 34 : COMPLEX(dp), ALLOCATABLE :: factors(:, :), identity(:, :)
228 : INTEGER :: i, nspin
229 :
230 34 : CPASSERT(SIZE(operators) == SIZE(positions) + 1)
231 34 : nspin = 1
232 34 : IF (PRESENT(soc)) nspin = 2
233 204 : ALLOCATE (identity(nspin, nspin), factors(nspin, nspin))
234 34 : identity = CMPLX(0.0_dp, 0.0_dp, dp)
235 82 : DO i = 1, nspin
236 82 : identity(i, i) = 1.0_dp
237 : END DO
238 34 : CALL create_operator(s, nspin, metric)
239 34 : CALL add_real_blocks(s, identity, metric)
240 134 : DO i = 1, SIZE(operators)
241 134 : CALL create_operator(s, nspin, operators(i))
242 : END DO
243 34 : CALL add_real_blocks(h, identity, operators(1))
244 100 : DO i = 1, SIZE(positions)
245 100 : CALL add_real_blocks(positions(i)%matrix, identity, operators(i + 1))
246 : END DO
247 34 : IF (PRESENT(soc)) THEN
248 14 : factors = CMPLX(0.0_dp, 0.0_dp, dp)
249 14 : factors(1, 2) = CMPLX(0.0_dp, 1.0_dp, dp)
250 14 : factors(2, 1) = CMPLX(0.0_dp, 1.0_dp, dp)
251 14 : CALL add_real_blocks(soc(1, 1)%matrix, factors, operators(1))
252 14 : factors(1, 2) = -1.0_dp
253 14 : factors(2, 1) = 1.0_dp
254 14 : CALL add_real_blocks(soc(2, 1)%matrix, factors, operators(1))
255 14 : factors = CMPLX(0.0_dp, 0.0_dp, dp)
256 14 : factors(1, 1) = CMPLX(0.0_dp, 1.0_dp, dp)
257 14 : factors(2, 2) = CMPLX(0.0_dp, -1.0_dp, dp)
258 14 : CALL add_real_blocks(soc(3, 1)%matrix, factors, operators(1))
259 : END IF
260 68 : END SUBROUTINE localizer_ao_operators
261 :
262 : ! **************************************************************************************************
263 : !> \brief Pair H-E*S with its negative to reuse the metric-aware localizer gap solvers.
264 : !> \param h Complex physical Hamiltonian, including spin if present
265 : !> \param s Complex physical metric with the same AO block layout
266 : !> \param energy Reference energy
267 : !> \param pencil Newly allocated diag(H-E*S,-H+E*S)
268 : !> \param metric Newly allocated diag(S,S)
269 : ! **************************************************************************************************
270 38 : SUBROUTINE localizer_dbcsr_gap_pencil(h, s, energy, pencil, metric)
271 : TYPE(localizer_matrix_type), INTENT(IN) :: h, s
272 : REAL(dp), INTENT(IN) :: energy
273 : TYPE(localizer_matrix_type), INTENT(OUT) :: pencil, metric
274 :
275 : COMPLEX(dp), PARAMETER :: imaginary = CMPLX(0.0_dp, 1.0_dp, dp)
276 :
277 : COMPLEX(dp) :: factors(2, 2)
278 :
279 38 : CALL create_operator(s%real_part, 2, pencil)
280 38 : CALL create_operator(s%real_part, 2, metric)
281 38 : factors = RESHAPE([1.0_dp, 0.0_dp, 0.0_dp, -1.0_dp], [2, 2])
282 38 : CALL add_real_blocks(h%real_part, factors, pencil)
283 266 : CALL add_real_blocks(h%imaginary_part, imaginary*factors, pencil)
284 266 : CALL add_real_blocks(s%real_part, -energy*factors, pencil)
285 266 : CALL add_real_blocks(s%imaginary_part, -energy*imaginary*factors, pencil)
286 38 : factors(2, 2) = 1.0_dp
287 38 : CALL add_real_blocks(s%real_part, factors, metric)
288 266 : CALL add_real_blocks(s%imaginary_part, imaginary*factors, metric)
289 38 : END SUBROUTINE localizer_dbcsr_gap_pencil
290 :
291 : ! **************************************************************************************************
292 : !> \brief Assemble L and B=diag(S,S), including optional physical spin and SOC, without dense AOs.
293 : !> \param h Real scalar Hamiltonian
294 : !> \param s Real AO overlap
295 : !> \param x First oriented real AO position
296 : !> \param y Second oriented real AO position
297 : !> \param point Two Cartesian coordinates
298 : !> \param energy Reference energy
299 : !> \param kappa Energy/length scale
300 : !> \param localizer Newly created complex DBCSR localizer
301 : !> \param metric Newly created complex DBCSR metric
302 : !> \param soc Real antisymmetric SOC components from V_SOC_xyz_from_pseudopotential
303 : !> \param periodic_cosines Optional two cos(G.r) matrices. Then x,y are sin(G.r), point is
304 : !> G.r0, and kappa is the energy scale eta (not the finite energy/length scale).
305 : !> \param shifted_hamiltonian Optional replacement for the physical H-E*S including SOC.
306 : !> Requires energy=0; the scalar h and SOC contributions are then not added again.
307 : !> \note Component-major ordering: auxiliary index, physical spin, then atomic AO blocks.
308 : ! **************************************************************************************************
309 78 : SUBROUTINE localizer_dbcsr_build(h, s, x, y, point, energy, kappa, localizer, metric, soc, periodic_cosines, &
310 : shifted_hamiltonian)
311 : TYPE(dbcsr_type), INTENT(IN) :: h, s, x, y
312 : REAL(KIND=dp), INTENT(IN) :: point(2), energy, kappa
313 : TYPE(localizer_matrix_type), INTENT(OUT) :: localizer, metric
314 : TYPE(dbcsr_p_type), INTENT(IN), OPTIONAL :: soc(:, :), periodic_cosines(:)
315 : TYPE(localizer_matrix_type), INTENT(IN), OPTIONAL :: shifted_hamiltonian
316 :
317 : COMPLEX(KIND=dp), PARAMETER :: imaginary = CMPLX(0.0_dp, 1.0_dp, dp)
318 :
319 : COMPLEX(dp) :: signs(2, 2)
320 78 : COMPLEX(KIND=dp), ALLOCATABLE :: factors(:, :), identity(:, :), sx(:, :), &
321 78 : sy(:, :), sz(:, :)
322 : INTEGER :: blocks, component, components, nspin, &
323 : physical_blocks
324 78 : INTEGER, POINTER :: col_dist(:), expanded_cols(:), &
325 78 : expanded_rows(:), expanded_sizes(:), &
326 78 : row_dist(:), sizes(:)
327 : TYPE(dbcsr_distribution_type) :: distribution, expanded_distribution
328 :
329 78 : nspin = 1
330 12 : IF (PRESENT(soc)) nspin = 2
331 78 : components = 2*nspin
332 78 : CALL dbcsr_get_info(s, distribution=distribution, row_blk_size=sizes)
333 78 : CALL dbcsr_distribution_get(distribution, row_dist=row_dist, col_dist=col_dist)
334 78 : blocks = SIZE(sizes)
335 390 : ALLOCATE (expanded_sizes(blocks*components), expanded_rows(blocks*components), expanded_cols(blocks*components))
336 156 : expanded_sizes = RESHAPE(SPREAD(sizes, 2, components), [blocks*components])
337 156 : expanded_rows = RESHAPE(SPREAD(row_dist, 2, components), [blocks*components])
338 156 : expanded_cols = RESHAPE(SPREAD(col_dist, 2, components), [blocks*components])
339 : CALL dbcsr_distribution_new(expanded_distribution, template=distribution, &
340 78 : row_dist=expanded_rows, col_dist=expanded_cols)
341 : CALL dbcsr_create(localizer%real_part, 'SPECTRAL LOCALIZER RE', expanded_distribution, &
342 78 : dbcsr_type_no_symmetry, expanded_sizes, expanded_sizes)
343 78 : CALL dbcsr_create(localizer%imaginary_part, template=localizer%real_part, name='SPECTRAL LOCALIZER IM')
344 78 : CALL dbcsr_create(metric%real_part, template=localizer%real_part, name='LOCALIZER METRIC RE')
345 78 : CALL dbcsr_create(metric%imaginary_part, template=localizer%real_part, name='LOCALIZER METRIC IM')
346 78 : CALL dbcsr_distribution_release(expanded_distribution)
347 78 : DEALLOCATE (expanded_sizes, expanded_rows, expanded_cols)
348 : ALLOCATE (factors(components, components), identity(components, components), &
349 936 : sx(components, components), sy(components, components), sz(components, components))
350 78 : identity = 0.0_dp
351 78 : sx = 0.0_dp
352 78 : sy = 0.0_dp
353 78 : sz = 0.0_dp
354 168 : DO component = 1, nspin
355 90 : identity(component, component) = 1.0_dp
356 90 : identity(nspin + component, nspin + component) = 1.0_dp
357 90 : sz(component, component) = 1.0_dp
358 90 : sz(nspin + component, nspin + component) = -1.0_dp
359 90 : sx(component, nspin + component) = 1.0_dp
360 90 : sx(nspin + component, component) = 1.0_dp
361 90 : sy(component, nspin + component) = -imaginary
362 168 : sy(nspin + component, component) = imaginary
363 : END DO
364 78 : CALL add_real_blocks(s, identity, metric)
365 78 : IF (PRESENT(shifted_hamiltonian)) THEN
366 18 : CPASSERT(energy == 0.0_dp)
367 18 : CALL dbcsr_get_info(shifted_hamiltonian%real_part, nblkrows_total=physical_blocks)
368 18 : CPASSERT(physical_blocks == nspin*blocks)
369 18 : signs = RESHAPE([1.0_dp, 0.0_dp, 0.0_dp, -1.0_dp], [2, 2])
370 18 : CALL add_real_blocks(shifted_hamiltonian%real_part, signs, localizer)
371 144 : CALL add_real_blocks(shifted_hamiltonian%imaginary_part, imaginary*signs, localizer)
372 : ELSE
373 60 : CALL add_real_blocks(h, sz, localizer)
374 : END IF
375 78 : IF (PRESENT(periodic_cosines)) THEN
376 10 : CPASSERT(SIZE(periodic_cosines) == 2)
377 210 : CALL add_real_blocks(s, -(energy + 2.0_dp*kappa)*sz, localizer)
378 210 : CALL add_real_blocks(periodic_cosines(1)%matrix, kappa*(COS(point(1))*sz - SIN(point(1))*sx), localizer)
379 210 : CALL add_real_blocks(periodic_cosines(2)%matrix, kappa*(COS(point(2))*sz - SIN(point(2))*sy), localizer)
380 210 : CALL add_real_blocks(x, kappa*(SIN(point(1))*sz + COS(point(1))*sx), localizer)
381 210 : CALL add_real_blocks(y, kappa*(SIN(point(2))*sz + COS(point(2))*sy), localizer)
382 : ELSE
383 504 : CALL add_real_blocks(s, -energy*sz - kappa*(point(1)*sx + point(2)*sy), localizer)
384 504 : CALL add_real_blocks(x, kappa*sx, localizer)
385 504 : CALL add_real_blocks(y, kappa*sy, localizer)
386 : END IF
387 78 : IF (PRESENT(soc) .AND. .NOT. PRESENT(shifted_hamiltonian)) THEN
388 : ! The real spherical-harmonic SOC matrices carry the missing factor i explicitly.
389 10 : factors = 0.0_dp
390 10 : factors(1, 2) = imaginary
391 10 : factors(2, 1) = imaginary
392 10 : factors(3, 4) = -imaginary
393 10 : factors(4, 3) = -imaginary
394 10 : CALL add_real_blocks(soc(1, 1)%matrix, factors, localizer)
395 10 : factors = 0.0_dp
396 10 : factors(1, 2) = -1.0_dp
397 10 : factors(2, 1) = 1.0_dp
398 10 : factors(3, 4) = 1.0_dp
399 10 : factors(4, 3) = -1.0_dp
400 10 : CALL add_real_blocks(soc(2, 1)%matrix, factors, localizer)
401 10 : factors = 0.0_dp
402 10 : factors(1, 1) = imaginary
403 10 : factors(2, 2) = -imaginary
404 10 : factors(3, 3) = -imaginary
405 10 : factors(4, 4) = imaginary
406 10 : CALL add_real_blocks(soc(3, 1)%matrix, factors, localizer)
407 : END IF
408 246 : END SUBROUTINE localizer_dbcsr_build
409 :
410 : ! **************************************************************************************************
411 : !> \brief Add a small Kronecker factor times distributed real AO blocks, preserving their union.
412 : !> \param source Real symmetric or antisymmetric AO matrix
413 : !> \param factors Small component matrix
414 : !> \param destination Complex block-sparse destination
415 : ! **************************************************************************************************
416 1192 : SUBROUTINE add_real_blocks(source, factors, destination)
417 : TYPE(dbcsr_type), INTENT(IN) :: source
418 : COMPLEX(KIND=dp), INTENT(IN) :: factors(:, :)
419 : TYPE(localizer_matrix_type), INTENT(INOUT) :: destination
420 :
421 1192 : COMPLEX(KIND=dp), ALLOCATABLE :: values(:, :)
422 : INTEGER :: blocks, col, i, j, row
423 : LOGICAL :: transposed
424 1192 : REAL(KIND=dp), POINTER :: block(:, :)
425 : TYPE(dbcsr_iterator_type) :: iterator
426 : TYPE(dbcsr_type) :: full
427 :
428 1192 : CALL dbcsr_get_info(source, nblkrows_total=blocks)
429 1192 : CALL dbcsr_create(full, template=source, matrix_type=dbcsr_type_no_symmetry)
430 1192 : CALL dbcsr_desymmetrize(source, full)
431 1192 : CALL dbcsr_iterator_start(iterator, full)
432 15783 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
433 14591 : CALL dbcsr_iterator_next_block(iterator, row, col, block, transposed=transposed)
434 14591 : CPASSERT(.NOT. transposed)
435 58364 : ALLOCATE (values(SIZE(block, 1), SIZE(block, 2)))
436 49272 : DO j = 1, SIZE(factors, 2)
437 141569 : DO i = 1, SIZE(factors, 1)
438 92297 : IF (ABS(factors(i, j)) == 0.0_dp) CYCLE
439 3836851 : values(:, :) = factors(i, j)*block
440 37433 : IF (REAL(factors(i, j), dp) /= 0.0_dp) THEN
441 : CALL dbcsr_put_block(destination%real_part, row + (i - 1)*blocks, col + (j - 1)*blocks, &
442 2761147 : REAL(values, dp), summation=.TRUE.)
443 : END IF
444 72114 : IF (AIMAG(factors(i, j)) /= 0.0_dp) THEN
445 : CALL dbcsr_put_block(destination%imaginary_part, row + (i - 1)*blocks, col + (j - 1)*blocks, &
446 1274208 : AIMAG(values), summation=.TRUE.)
447 : END IF
448 : END DO
449 : END DO
450 30374 : DEALLOCATE (values)
451 : END DO
452 1192 : CALL dbcsr_iterator_stop(iterator)
453 1192 : CALL dbcsr_release(full)
454 1192 : CALL dbcsr_finalize(destination%real_part)
455 1192 : CALL dbcsr_finalize(destination%imaginary_part)
456 2384 : END SUBROUTINE add_real_blocks
457 :
458 : ! **************************************************************************************************
459 : !> \brief Collect the size-limited dense oracle; sparse solvers must not call this routine.
460 : !> \param matrix Fully stored complex block-sparse matrix
461 : !> \param para_env Matrix communicator
462 : !> \param dense Replicated matrix
463 : ! **************************************************************************************************
464 434 : SUBROUTINE localizer_dbcsr_dense(matrix, para_env, dense)
465 : TYPE(localizer_matrix_type), INTENT(INOUT) :: matrix
466 : TYPE(mp_para_env_type), POINTER :: para_env
467 : COMPLEX(KIND=dp), ALLOCATABLE, INTENT(OUT) :: dense(:, :)
468 :
469 : INTEGER :: col, n, row
470 : LOGICAL :: transposed
471 434 : REAL(KIND=dp), POINTER :: block(:, :)
472 : TYPE(dbcsr_iterator_type) :: iterator
473 :
474 434 : CALL dbcsr_get_info(matrix%real_part, nfullrows_total=n)
475 1736 : ALLOCATE (dense(n, n))
476 434 : dense = 0.0_dp
477 434 : CALL dbcsr_iterator_start(iterator, matrix%real_part)
478 10851 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
479 : CALL dbcsr_iterator_next_block(iterator, block=block, row_offset=row, col_offset=col, &
480 10417 : transposed=transposed)
481 10417 : CPASSERT(.NOT. transposed)
482 1064685 : dense(row:row + SIZE(block, 1) - 1, col:col + SIZE(block, 2) - 1) = block
483 : END DO
484 434 : CALL dbcsr_iterator_stop(iterator)
485 434 : CALL dbcsr_iterator_start(iterator, matrix%imaginary_part)
486 8582 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
487 : CALL dbcsr_iterator_next_block(iterator, block=block, row_offset=row, col_offset=col, &
488 8148 : transposed=transposed)
489 8148 : CPASSERT(.NOT. transposed)
490 : dense(row:row + SIZE(block, 1) - 1, col:col + SIZE(block, 2) - 1) = &
491 779774 : dense(row:row + SIZE(block, 1) - 1, col:col + SIZE(block, 2) - 1) + CMPLX(0.0_dp, 1.0_dp, dp)*block
492 : END DO
493 434 : CALL dbcsr_iterator_stop(iterator)
494 434 : CALL para_env%sum(dense)
495 1302 : END SUBROUTINE localizer_dbcsr_dense
496 :
497 : ! **************************************************************************************************
498 : !> \brief Check both components before a sparse solver discards the upper triangle.
499 : !> \param matrix Fully stored complex matrix
500 : !> \param para_env Matrix communicator
501 : !> \param tolerance Relative Hermiticity tolerance
502 : !> \param valid Whether all entries are finite and the matrix is Hermitian
503 : !> \param relative_error Largest component transpose residual divided by the largest component entry;
504 : !> HUGE for nonfinite entries or residuals
505 : ! **************************************************************************************************
506 414 : SUBROUTINE localizer_dbcsr_check(matrix, para_env, tolerance, valid, relative_error)
507 : TYPE(localizer_matrix_type), INTENT(INOUT) :: matrix
508 : TYPE(mp_para_env_type), POINTER :: para_env
509 : REAL(KIND=dp), INTENT(IN) :: tolerance
510 : LOGICAL, INTENT(OUT) :: valid
511 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: relative_error
512 :
513 : LOGICAL :: valid_real, valid_imaginary
514 : REAL(KIND=dp) :: real_norm, imaginary_norm, real_error, imaginary_error
515 :
516 414 : CALL check_component(matrix%real_part, -1.0_dp, real_norm, real_error, valid_real)
517 414 : CALL check_component(matrix%imaginary_part, 1.0_dp, imaginary_norm, imaginary_error, valid_imaginary)
518 414 : valid = valid_real .AND. valid_imaginary
519 : valid = valid .AND. MAX(real_error, imaginary_error) <= &
520 414 : tolerance*MAX(TINY(1.0_dp), real_norm, imaginary_norm)
521 414 : IF (PRESENT(relative_error)) THEN
522 124 : relative_error = HUGE(1.0_dp)
523 124 : IF (valid_real .AND. valid_imaginary) THEN
524 124 : relative_error = MAX(real_error, imaginary_error)/MAX(TINY(1.0_dp), real_norm, imaginary_norm)
525 : END IF
526 : END IF
527 :
528 : CONTAINS
529 :
530 : ! **************************************************************************************************
531 : !> \brief Check real symmetry or imaginary antisymmetry using a distributed transpose.
532 : !> \param component Fully stored real matrix
533 : !> \param sign Sign multiplying the transpose
534 : !> \param norm Largest absolute entry
535 : !> \param error Largest transpose residual
536 : !> \param ok Whether all entries and residuals are finite
537 : ! **************************************************************************************************
538 2484 : SUBROUTINE check_component(component, sign, norm, error, ok)
539 : TYPE(dbcsr_type), INTENT(INOUT) :: component
540 : REAL(KIND=dp), INTENT(IN) :: sign
541 : REAL(KIND=dp), INTENT(OUT) :: norm, error
542 : LOGICAL, INTENT(OUT) :: ok
543 :
544 : INTEGER :: invalid
545 828 : REAL(KIND=dp), POINTER :: block(:, :)
546 : TYPE(dbcsr_distribution_type) :: distribution
547 : TYPE(dbcsr_iterator_type) :: iterator
548 : TYPE(dbcsr_type) :: difference
549 :
550 828 : invalid = 0
551 828 : norm = 0.0_dp
552 828 : error = 0.0_dp
553 828 : CALL dbcsr_iterator_start(iterator, component)
554 14918 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
555 14090 : CALL dbcsr_iterator_next_block(iterator, block=block)
556 527194 : IF (.NOT. ALL(ieee_is_finite(block))) THEN
557 0 : invalid = 1
558 : ELSE
559 526366 : norm = MAX(norm, MAXVAL(ABS(block)))
560 : END IF
561 : END DO
562 828 : CALL dbcsr_iterator_stop(iterator)
563 828 : CALL para_env%max(invalid)
564 828 : ok = .FALSE.
565 828 : IF (invalid /= 0) RETURN
566 828 : CALL para_env%max(norm)
567 828 : CALL dbcsr_get_info(component, distribution=distribution)
568 828 : CALL dbcsr_transposed(difference, component, shallow_data_copy=.FALSE., use_distribution=distribution)
569 828 : CALL dbcsr_add(difference, component, sign, 1.0_dp)
570 828 : error = 0.0_dp
571 828 : CALL dbcsr_iterator_start(iterator, difference)
572 14918 : DO WHILE (dbcsr_iterator_blocks_left(iterator))
573 14090 : CALL dbcsr_iterator_next_block(iterator, block=block)
574 527194 : IF (.NOT. ALL(ieee_is_finite(block))) THEN
575 0 : invalid = 1
576 : ELSE
577 526366 : error = MAX(error, MAXVAL(ABS(block)))
578 : END IF
579 : END DO
580 828 : CALL dbcsr_iterator_stop(iterator)
581 828 : CALL dbcsr_release(difference)
582 828 : CALL para_env%max(invalid)
583 828 : CALL para_env%max(error)
584 828 : ok = invalid == 0
585 828 : END SUBROUTINE check_component
586 :
587 : END SUBROUTINE localizer_dbcsr_check
588 :
589 : ! **************************************************************************************************
590 : !> \brief Release both real components of a complex sparse operator.
591 : !> \param matrix Sparse operator
592 : ! **************************************************************************************************
593 804 : SUBROUTINE localizer_dbcsr_release(matrix)
594 : TYPE(localizer_matrix_type), INTENT(INOUT) :: matrix
595 :
596 804 : CALL dbcsr_release(matrix%real_part)
597 804 : CALL dbcsr_release(matrix%imaginary_part)
598 804 : END SUBROUTINE localizer_dbcsr_release
599 :
600 0 : END MODULE spectral_localizer_dbcsr
|
