(git:08825b9)
Loading...
Searching...
No Matches
spectral_localizer_dbcsr.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
38CONTAINS
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 SUBROUTINE localizer_dbcsr_3d_build(h, s, positions, point, energy, kappa, block, metric, soc, &
56 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 TYPE(localizer_matrix_type) :: physical(1), physical_metric
71
72 cpassert(SIZE(positions) == 3)
73 identity = 0.0_dp
74 pauli = 0.0_dp
75 DO i = 1, 4
76 identity(i, i) = 1.0_dp
77 END DO
78 DO i = 1, 2
79 pauli(i, i + 2, 1) = 1.0_dp
80 pauli(i + 2, i, 1) = 1.0_dp
81 pauli(i, i + 2, 2) = -ii
82 pauli(i + 2, i, 2) = ii
83 pauli(i, i, 3) = 1.0_dp
84 pauli(i + 2, i + 2, 3) = -1.0_dp
85 END DO
86 CALL create_operator(s, 4, block)
87 CALL create_operator(s, 4, metric)
88 CALL add_real_blocks(s, identity, metric)
89 factors = 0.0_dp
90 factors(1, 1) = 1.0_dp
91 factors(2, 2) = 1.0_dp
92 IF (PRESENT(shifted_hamiltonian)) THEN
93 cpassert(energy == 0.0_dp)
94 CALL add_real_blocks(shifted_hamiltonian%real_part, factors, block)
95 CALL add_real_blocks(shifted_hamiltonian%imaginary_part, ii*factors, block)
96 ELSE
97 CALL localizer_ao_operators(h, s, positions(:0), physical, physical_metric, soc)
98 CALL add_real_blocks(physical(1)%real_part, factors, block)
99 CALL add_real_blocks(physical(1)%imaginary_part, ii*factors, block)
100 CALL localizer_dbcsr_release(physical(1))
101 CALL localizer_dbcsr_release(physical_metric)
102 CALL add_real_blocks(s, -energy*identity, block)
103 END IF
104 IF (PRESENT(periodic_cosines)) THEN
105 cpassert(SIZE(periodic_cosines) == 3)
106 CALL add_real_blocks(s, -3.0_dp*kappa*identity, block)
107 DO i = 1, 3
108 CALL add_real_blocks(periodic_cosines(i)%matrix, &
109 kappa*(cos(point(i))*identity + ii*sin(point(i))*pauli(:, :, i)), block)
110 CALL add_real_blocks(positions(i)%matrix, &
111 kappa*(sin(point(i))*identity - ii*cos(point(i))*pauli(:, :, i)), block)
112 END DO
113 ELSE
114 DO i = 1, 3
115 CALL add_real_blocks(positions(i)%matrix, -ii*kappa*pauli(:, :, i), block)
116 CALL add_real_blocks(s, ii*kappa*point(i)*pauli(:, :, i), block)
117 END DO
118 END IF
119 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 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 CALL create_operator(matrix%real_part, 2, expanded)
134 factors = reshape([1.0_dp, 0.0_dp, 0.0_dp, 1.0_dp], [2, 2])
135 CALL add_real_blocks(matrix%real_part, factors, expanded)
136 factors = reshape([0.0_dp, 1.0_dp, -1.0_dp, 0.0_dp], [2, 2])
137 CALL add_real_blocks(matrix%imaginary_part, factors, expanded)
138 CALL dbcsr_create(realified, template=expanded%real_part)
139 CALL dbcsr_copy(realified, expanded%real_part)
140 CALL localizer_dbcsr_release(expanded)
141 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 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 REAL(dp), POINTER :: block(:, :)
157 TYPE(dbcsr_iterator_type) :: iterator
158
159 CALL dbcsr_get_info(template%real_part, nblkrows_total=blocks)
160 CALL dbcsr_get_info(realified, nblkrows_total=nblocks)
161 cpassert(nblocks == 2*blocks)
162 CALL create_operator(template%real_part, 1, matrix)
163 CALL dbcsr_iterator_start(iterator, realified)
164 DO WHILE (dbcsr_iterator_blocks_left(iterator))
165 CALL dbcsr_iterator_next_block(iterator, row, col, block, transposed=transposed)
166 cpassert(.NOT. transposed)
167 IF (col > blocks) cycle
168 IF (row <= blocks) THEN
169 CALL dbcsr_put_block(matrix%real_part, row, col, block)
170 ELSE
171 CALL dbcsr_put_block(matrix%imaginary_part, row - blocks, col, block)
172 END IF
173 END DO
174 CALL dbcsr_iterator_stop(iterator)
175 CALL dbcsr_finalize(matrix%real_part)
176 CALL dbcsr_finalize(matrix%imaginary_part)
177 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 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 INTEGER, POINTER :: col_dist(:), expanded_cols(:), &
192 expanded_rows(:), expanded_sizes(:), &
193 row_dist(:), sizes(:)
194 TYPE(dbcsr_distribution_type) :: distribution, expanded_distribution
195
196 CALL dbcsr_get_info(source, distribution=distribution, row_blk_size=sizes)
197 CALL dbcsr_distribution_get(distribution, row_dist=row_dist, col_dist=col_dist)
198 blocks = SIZE(sizes)
199 ALLOCATE (expanded_sizes(blocks*components), expanded_rows(blocks*components), expanded_cols(blocks*components))
200 expanded_sizes = reshape(spread(sizes, 2, components), [blocks*components])
201 expanded_rows = reshape(spread(row_dist, 2, components), [blocks*components])
202 expanded_cols = reshape(spread(col_dist, 2, components), [blocks*components])
203 CALL dbcsr_distribution_new(expanded_distribution, template=distribution, &
204 row_dist=expanded_rows, col_dist=expanded_cols)
205 CALL dbcsr_create(matrix%real_part, 'FINITE AO RE', expanded_distribution, &
206 dbcsr_type_no_symmetry, expanded_sizes, expanded_sizes)
207 CALL dbcsr_create(matrix%imaginary_part, template=matrix%real_part, name='FINITE AO IM')
208 CALL dbcsr_distribution_release(expanded_distribution)
209 DEALLOCATE (expanded_sizes, expanded_rows, expanded_cols)
210 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 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 COMPLEX(dp), ALLOCATABLE :: factors(:, :), identity(:, :)
228 INTEGER :: i, nspin
229
230 cpassert(SIZE(operators) == SIZE(positions) + 1)
231 nspin = 1
232 IF (PRESENT(soc)) nspin = 2
233 ALLOCATE (identity(nspin, nspin), factors(nspin, nspin))
234 identity = cmplx(0.0_dp, 0.0_dp, dp)
235 DO i = 1, nspin
236 identity(i, i) = 1.0_dp
237 END DO
238 CALL create_operator(s, nspin, metric)
239 CALL add_real_blocks(s, identity, metric)
240 DO i = 1, SIZE(operators)
241 CALL create_operator(s, nspin, operators(i))
242 END DO
243 CALL add_real_blocks(h, identity, operators(1))
244 DO i = 1, SIZE(positions)
245 CALL add_real_blocks(positions(i)%matrix, identity, operators(i + 1))
246 END DO
247 IF (PRESENT(soc)) THEN
248 factors = cmplx(0.0_dp, 0.0_dp, dp)
249 factors(1, 2) = cmplx(0.0_dp, 1.0_dp, dp)
250 factors(2, 1) = cmplx(0.0_dp, 1.0_dp, dp)
251 CALL add_real_blocks(soc(1, 1)%matrix, factors, operators(1))
252 factors(1, 2) = -1.0_dp
253 factors(2, 1) = 1.0_dp
254 CALL add_real_blocks(soc(2, 1)%matrix, factors, operators(1))
255 factors = cmplx(0.0_dp, 0.0_dp, dp)
256 factors(1, 1) = cmplx(0.0_dp, 1.0_dp, dp)
257 factors(2, 2) = cmplx(0.0_dp, -1.0_dp, dp)
258 CALL add_real_blocks(soc(3, 1)%matrix, factors, operators(1))
259 END IF
260 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 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 CALL create_operator(s%real_part, 2, pencil)
280 CALL create_operator(s%real_part, 2, metric)
281 factors = reshape([1.0_dp, 0.0_dp, 0.0_dp, -1.0_dp], [2, 2])
282 CALL add_real_blocks(h%real_part, factors, pencil)
283 CALL add_real_blocks(h%imaginary_part, imaginary*factors, pencil)
284 CALL add_real_blocks(s%real_part, -energy*factors, pencil)
285 CALL add_real_blocks(s%imaginary_part, -energy*imaginary*factors, pencil)
286 factors(2, 2) = 1.0_dp
287 CALL add_real_blocks(s%real_part, factors, metric)
288 CALL add_real_blocks(s%imaginary_part, imaginary*factors, metric)
289 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 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 COMPLEX(KIND=dp), ALLOCATABLE :: factors(:, :), identity(:, :), sx(:, :), &
321 sy(:, :), sz(:, :)
322 INTEGER :: blocks, component, components, nspin, &
323 physical_blocks
324 INTEGER, POINTER :: col_dist(:), expanded_cols(:), &
325 expanded_rows(:), expanded_sizes(:), &
326 row_dist(:), sizes(:)
327 TYPE(dbcsr_distribution_type) :: distribution, expanded_distribution
328
329 nspin = 1
330 IF (PRESENT(soc)) nspin = 2
331 components = 2*nspin
332 CALL dbcsr_get_info(s, distribution=distribution, row_blk_size=sizes)
333 CALL dbcsr_distribution_get(distribution, row_dist=row_dist, col_dist=col_dist)
334 blocks = SIZE(sizes)
335 ALLOCATE (expanded_sizes(blocks*components), expanded_rows(blocks*components), expanded_cols(blocks*components))
336 expanded_sizes = reshape(spread(sizes, 2, components), [blocks*components])
337 expanded_rows = reshape(spread(row_dist, 2, components), [blocks*components])
338 expanded_cols = reshape(spread(col_dist, 2, components), [blocks*components])
339 CALL dbcsr_distribution_new(expanded_distribution, template=distribution, &
340 row_dist=expanded_rows, col_dist=expanded_cols)
341 CALL dbcsr_create(localizer%real_part, 'SPECTRAL LOCALIZER RE', expanded_distribution, &
342 dbcsr_type_no_symmetry, expanded_sizes, expanded_sizes)
343 CALL dbcsr_create(localizer%imaginary_part, template=localizer%real_part, name='SPECTRAL LOCALIZER IM')
344 CALL dbcsr_create(metric%real_part, template=localizer%real_part, name='LOCALIZER METRIC RE')
345 CALL dbcsr_create(metric%imaginary_part, template=localizer%real_part, name='LOCALIZER METRIC IM')
346 CALL dbcsr_distribution_release(expanded_distribution)
347 DEALLOCATE (expanded_sizes, expanded_rows, expanded_cols)
348 ALLOCATE (factors(components, components), identity(components, components), &
349 sx(components, components), sy(components, components), sz(components, components))
350 identity = 0.0_dp
351 sx = 0.0_dp
352 sy = 0.0_dp
353 sz = 0.0_dp
354 DO component = 1, nspin
355 identity(component, component) = 1.0_dp
356 identity(nspin + component, nspin + component) = 1.0_dp
357 sz(component, component) = 1.0_dp
358 sz(nspin + component, nspin + component) = -1.0_dp
359 sx(component, nspin + component) = 1.0_dp
360 sx(nspin + component, component) = 1.0_dp
361 sy(component, nspin + component) = -imaginary
362 sy(nspin + component, component) = imaginary
363 END DO
364 CALL add_real_blocks(s, identity, metric)
365 IF (PRESENT(shifted_hamiltonian)) THEN
366 cpassert(energy == 0.0_dp)
367 CALL dbcsr_get_info(shifted_hamiltonian%real_part, nblkrows_total=physical_blocks)
368 cpassert(physical_blocks == nspin*blocks)
369 signs = reshape([1.0_dp, 0.0_dp, 0.0_dp, -1.0_dp], [2, 2])
370 CALL add_real_blocks(shifted_hamiltonian%real_part, signs, localizer)
371 CALL add_real_blocks(shifted_hamiltonian%imaginary_part, imaginary*signs, localizer)
372 ELSE
373 CALL add_real_blocks(h, sz, localizer)
374 END IF
375 IF (PRESENT(periodic_cosines)) THEN
376 cpassert(SIZE(periodic_cosines) == 2)
377 CALL add_real_blocks(s, -(energy + 2.0_dp*kappa)*sz, localizer)
378 CALL add_real_blocks(periodic_cosines(1)%matrix, kappa*(cos(point(1))*sz - sin(point(1))*sx), localizer)
379 CALL add_real_blocks(periodic_cosines(2)%matrix, kappa*(cos(point(2))*sz - sin(point(2))*sy), localizer)
380 CALL add_real_blocks(x, kappa*(sin(point(1))*sz + cos(point(1))*sx), localizer)
381 CALL add_real_blocks(y, kappa*(sin(point(2))*sz + cos(point(2))*sy), localizer)
382 ELSE
383 CALL add_real_blocks(s, -energy*sz - kappa*(point(1)*sx + point(2)*sy), localizer)
384 CALL add_real_blocks(x, kappa*sx, localizer)
385 CALL add_real_blocks(y, kappa*sy, localizer)
386 END IF
387 IF (PRESENT(soc) .AND. .NOT. PRESENT(shifted_hamiltonian)) THEN
388 ! The real spherical-harmonic SOC matrices carry the missing factor i explicitly.
389 factors = 0.0_dp
390 factors(1, 2) = imaginary
391 factors(2, 1) = imaginary
392 factors(3, 4) = -imaginary
393 factors(4, 3) = -imaginary
394 CALL add_real_blocks(soc(1, 1)%matrix, factors, localizer)
395 factors = 0.0_dp
396 factors(1, 2) = -1.0_dp
397 factors(2, 1) = 1.0_dp
398 factors(3, 4) = 1.0_dp
399 factors(4, 3) = -1.0_dp
400 CALL add_real_blocks(soc(2, 1)%matrix, factors, localizer)
401 factors = 0.0_dp
402 factors(1, 1) = imaginary
403 factors(2, 2) = -imaginary
404 factors(3, 3) = -imaginary
405 factors(4, 4) = imaginary
406 CALL add_real_blocks(soc(3, 1)%matrix, factors, localizer)
407 END IF
408 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 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 COMPLEX(KIND=dp), ALLOCATABLE :: values(:, :)
422 INTEGER :: blocks, col, i, j, row
423 LOGICAL :: transposed
424 REAL(kind=dp), POINTER :: block(:, :)
425 TYPE(dbcsr_iterator_type) :: iterator
426 TYPE(dbcsr_type) :: full
427
428 CALL dbcsr_get_info(source, nblkrows_total=blocks)
429 CALL dbcsr_create(full, template=source, matrix_type=dbcsr_type_no_symmetry)
430 CALL dbcsr_desymmetrize(source, full)
431 CALL dbcsr_iterator_start(iterator, full)
432 DO WHILE (dbcsr_iterator_blocks_left(iterator))
433 CALL dbcsr_iterator_next_block(iterator, row, col, block, transposed=transposed)
434 cpassert(.NOT. transposed)
435 ALLOCATE (values(SIZE(block, 1), SIZE(block, 2)))
436 DO j = 1, SIZE(factors, 2)
437 DO i = 1, SIZE(factors, 1)
438 IF (abs(factors(i, j)) == 0.0_dp) cycle
439 values(:, :) = factors(i, j)*block
440 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 REAL(values, dp), summation=.true.)
443 END IF
444 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 aimag(values), summation=.true.)
447 END IF
448 END DO
449 END DO
450 DEALLOCATE (values)
451 END DO
452 CALL dbcsr_iterator_stop(iterator)
453 CALL dbcsr_release(full)
454 CALL dbcsr_finalize(destination%real_part)
455 CALL dbcsr_finalize(destination%imaginary_part)
456 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 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 REAL(kind=dp), POINTER :: block(:, :)
472 TYPE(dbcsr_iterator_type) :: iterator
473
474 CALL dbcsr_get_info(matrix%real_part, nfullrows_total=n)
475 ALLOCATE (dense(n, n))
476 dense = 0.0_dp
477 CALL dbcsr_iterator_start(iterator, matrix%real_part)
478 DO WHILE (dbcsr_iterator_blocks_left(iterator))
479 CALL dbcsr_iterator_next_block(iterator, block=block, row_offset=row, col_offset=col, &
480 transposed=transposed)
481 cpassert(.NOT. transposed)
482 dense(row:row + SIZE(block, 1) - 1, col:col + SIZE(block, 2) - 1) = block
483 END DO
484 CALL dbcsr_iterator_stop(iterator)
485 CALL dbcsr_iterator_start(iterator, matrix%imaginary_part)
486 DO WHILE (dbcsr_iterator_blocks_left(iterator))
487 CALL dbcsr_iterator_next_block(iterator, block=block, row_offset=row, col_offset=col, &
488 transposed=transposed)
489 cpassert(.NOT. transposed)
490 dense(row:row + SIZE(block, 1) - 1, col:col + SIZE(block, 2) - 1) = &
491 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 CALL dbcsr_iterator_stop(iterator)
494 CALL para_env%sum(dense)
495 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 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 CALL check_component(matrix%real_part, -1.0_dp, real_norm, real_error, valid_real)
517 CALL check_component(matrix%imaginary_part, 1.0_dp, imaginary_norm, imaginary_error, valid_imaginary)
518 valid = valid_real .AND. valid_imaginary
519 valid = valid .AND. max(real_error, imaginary_error) <= &
520 tolerance*max(tiny(1.0_dp), real_norm, imaginary_norm)
521 IF (PRESENT(relative_error)) THEN
522 relative_error = huge(1.0_dp)
523 IF (valid_real .AND. valid_imaginary) THEN
524 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 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 REAL(kind=dp), POINTER :: block(:, :)
546 TYPE(dbcsr_distribution_type) :: distribution
547 TYPE(dbcsr_iterator_type) :: iterator
548 TYPE(dbcsr_type) :: difference
549
550 invalid = 0
551 norm = 0.0_dp
552 error = 0.0_dp
553 CALL dbcsr_iterator_start(iterator, component)
554 DO WHILE (dbcsr_iterator_blocks_left(iterator))
555 CALL dbcsr_iterator_next_block(iterator, block=block)
556 IF (.NOT. all(ieee_is_finite(block))) THEN
557 invalid = 1
558 ELSE
559 norm = max(norm, maxval(abs(block)))
560 END IF
561 END DO
562 CALL dbcsr_iterator_stop(iterator)
563 CALL para_env%max(invalid)
564 ok = .false.
565 IF (invalid /= 0) RETURN
566 CALL para_env%max(norm)
567 CALL dbcsr_get_info(component, distribution=distribution)
568 CALL dbcsr_transposed(difference, component, shallow_data_copy=.false., use_distribution=distribution)
569 CALL dbcsr_add(difference, component, sign, 1.0_dp)
570 error = 0.0_dp
571 CALL dbcsr_iterator_start(iterator, difference)
572 DO WHILE (dbcsr_iterator_blocks_left(iterator))
573 CALL dbcsr_iterator_next_block(iterator, block=block)
574 IF (.NOT. all(ieee_is_finite(block))) THEN
575 invalid = 1
576 ELSE
577 error = max(error, maxval(abs(block)))
578 END IF
579 END DO
580 CALL dbcsr_iterator_stop(iterator)
581 CALL dbcsr_release(difference)
582 CALL para_env%max(invalid)
583 CALL para_env%max(error)
584 ok = invalid == 0
585 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 SUBROUTINE localizer_dbcsr_release(matrix)
594 TYPE(localizer_matrix_type), INTENT(INOUT) :: matrix
595
596 CALL dbcsr_release(matrix%real_part)
597 CALL dbcsr_release(matrix%imaginary_part)
598 END SUBROUTINE localizer_dbcsr_release
599
600END MODULE spectral_localizer_dbcsr
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_transposed
subroutine, public dbcsr_transposed(transposed, normal, shallow_data_copy, transpose_distribution, use_distribution)
...
Definition cp_dbcsr_api.F:1242
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:643
cp_dbcsr_api::dbcsr_distribution_new
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
Definition cp_dbcsr_api.F:616
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:572
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
spectral_localizer_dbcsr
Distributed block-sparse finite localizers, retaining the covariant AO metric.
Definition spectral_localizer_dbcsr.F:12
spectral_localizer_dbcsr::localizer_dbcsr_from_real
subroutine, public localizer_dbcsr_from_real(realified, template, matrix)
Extract the two left quadrants of R(A); does not average away structure errors.
Definition spectral_localizer_dbcsr.F:150
spectral_localizer_dbcsr::localizer_ao_operators
subroutine, public localizer_ao_operators(h, s, positions, operators, metric, soc)
Shared AO Hamiltonian, Hermitian scalar observables and optional GTH physical-spin SOC.
Definition spectral_localizer_dbcsr.F:222
spectral_localizer_dbcsr::localizer_dbcsr_dense
subroutine, public localizer_dbcsr_dense(matrix, para_env, dense)
Collect the size-limited dense oracle; sparse solvers must not call this routine.
Definition spectral_localizer_dbcsr.F:465
spectral_localizer_dbcsr::localizer_dbcsr_gap_pencil
subroutine, public localizer_dbcsr_gap_pencil(h, s, energy, pencil, metric)
Pair H-E*S with its negative to reuse the metric-aware localizer gap solvers.
Definition spectral_localizer_dbcsr.F:271
spectral_localizer_dbcsr::localizer_dbcsr_release
subroutine, public localizer_dbcsr_release(matrix)
Release both real components of a complex sparse operator.
Definition spectral_localizer_dbcsr.F:594
spectral_localizer_dbcsr::localizer_dbcsr_check
subroutine, public localizer_dbcsr_check(matrix, para_env, tolerance, valid, relative_error)
Check both components before a sparse solver discards the upper triangle.
Definition spectral_localizer_dbcsr.F:507
spectral_localizer_dbcsr::localizer_dbcsr_realify
subroutine, public localizer_dbcsr_realify(matrix, realified)
Distributed real matrix R(A)=[Re(A),-Im(A);Im(A),Re(A)].
Definition spectral_localizer_dbcsr.F:127
spectral_localizer_dbcsr::localizer_dbcsr_build
subroutine, public localizer_dbcsr_build(h, s, x, y, point, energy, kappa, localizer, metric, soc, periodic_cosines, shifted_hamiltonian)
Assemble L and B=diag(S,S), including optional physical spin and SOC, without dense AOs.
Definition spectral_localizer_dbcsr.F:311
spectral_localizer_dbcsr::localizer_dbcsr_3d_build
subroutine, public localizer_dbcsr_3d_build(h, s, positions, point, energy, kappa, block, metric, soc, shifted_hamiltonian, periodic_cosines)
Build the 3D AII chiral block in DBCSR; not a Hermitian 2D localizer.
Definition spectral_localizer_dbcsr.F:57
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
spectral_localizer_dbcsr::localizer_matrix_type
Definition spectral_localizer_dbcsr.F:34