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 Post-SCF finite-system spectral localizers and analytic AO position export.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE qs_spectral_localizer
13 : USE bibliography, ONLY: Dixon2023,&
14 : Doll2025,&
15 : Loring2015,&
16 : Wong2026
17 : USE cell_types, ONLY: cell_type,&
18 : use_perd_none,&
19 : use_perd_xy,&
20 : use_perd_xyz,&
21 : use_perd_xz,&
22 : use_perd_yz
23 : USE cp_control_types, ONLY: dft_control_type
24 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
25 : dbcsr_get_info,&
26 : dbcsr_p_type,&
27 : dbcsr_set
28 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
29 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
30 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
31 : USE cp_log_handling, ONLY: cp_get_default_logger,&
32 : cp_logger_get_default_io_unit,&
33 : cp_logger_type
34 : USE cp_output_handling, ONLY: cp_p_file,&
35 : cp_print_key_finished_output,&
36 : cp_print_key_should_output,&
37 : cp_print_key_unit_nr
38 : USE ieee_arithmetic, ONLY: ieee_is_finite
39 : USE input_constants, ONLY: do_method_gapw,&
40 : do_method_gpw
41 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
42 : section_vals_type,&
43 : section_vals_val_get
44 : USE kinds, ONLY: dp
45 : USE mathconstants, ONLY: twopi
46 : USE message_passing, ONLY: mp_para_env_type
47 : USE qs_environment_types, ONLY: get_qs_env,&
48 : qs_environment_type
49 : USE qs_finite_ao, ONLY: finite_ao_positions,&
50 : finite_ao_soc
51 : USE qs_localizer_torus, ONLY: localizer_torus_operators
52 : USE qs_moments, ONLY: build_berry_moment_matrix
53 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
54 : USE reference_manager, ONLY: cite_reference
55 : USE spectral_localizer, ONLY: localizer_3d_matrix_evaluate,&
56 : localizer_matrix_evaluate
57 : USE spectral_localizer_bloch, ONLY: localizer_flatten_bloch
58 : USE spectral_localizer_dbcsr, ONLY: localizer_ao_operators,&
59 : localizer_dbcsr_3d_build,&
60 : localizer_dbcsr_build,&
61 : localizer_dbcsr_dense,&
62 : localizer_dbcsr_release,&
63 : localizer_matrix_type
64 : USE spectral_localizer_flatten, ONLY: localizer_flatten_checked
65 : #include "./base/base_uses.f90"
66 :
67 : IMPLICIT NONE
68 : PRIVATE
69 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_spectral_localizer'
70 : PUBLIC :: qs_scf_post_localizer
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief Export finite AO positions and optionally evaluate localizers after SCF.
76 : !> \param qs_env Converged Quickstep environment
77 : ! **************************************************************************************************
78 14147 : SUBROUTINE qs_scf_post_localizer(qs_env)
79 : TYPE(qs_environment_type), POINTER :: qs_env
80 :
81 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_localizer'
82 : CHARACTER(len=1), PARAMETER :: axes(3) = ['X', 'Y', 'Z']
83 : INTEGER, PARAMETER :: plane_periodicity(3) = [use_perd_xy, use_perd_xz, use_perd_yz]
84 :
85 14147 : COMPLEX(KIND=dp), ALLOCATABLE :: local_dense(:, :), metric_dense(:, :)
86 : INTEGER :: after, axis, directions(3), flatten_handle, formulation, handle, ie, ik, &
87 : invariant, ip, iw, max_ao, n, nao, ndim, nposition, output_unit, plane, poisson_periodic, &
88 : spin_channel, status, topology
89 14147 : INTEGER, POINTER :: grid(:)
90 : LOGICAL :: active, aii, do_kpoints, do_soc, explicit_scale, flatten, omit_headers, &
91 : periodic_localizer, print_position, print_soc, time_reversal, torus
92 : REAL(dp) :: build_energy, electronic_gap(2), eps_electronic_gap, eps_flattening, &
93 : flatten_residuals(4), flatten_scale, local_point(3), physical_ratio
94 : REAL(KIND=dp) :: eps_gap, eps_metric, gap, &
95 : matrix_tolerance, metric_ratio
96 28294 : REAL(KIND=dp), POINTER :: energies(:), kappas(:), point(:)
97 : TYPE(cell_type), POINTER :: cell
98 : TYPE(cp_logger_type), POINTER :: logger
99 14147 : TYPE(dbcsr_p_type), POINTER :: cosines(:), image_matrices(:, :), &
100 14147 : matrix_ks(:), matrix_s(:), &
101 14147 : positions(:), soc_matrices(:, :)
102 : TYPE(dft_control_type), POINTER :: dft_control
103 28294 : TYPE(localizer_matrix_type) :: localizer, metric, physical_h(1), &
104 : physical_metric
105 : TYPE(localizer_matrix_type), POINTER :: flat_hamiltonian
106 : TYPE(mp_para_env_type), POINTER :: para_env
107 14147 : TYPE(neighbor_list_set_p_type), POINTER :: sab_orb(:)
108 : TYPE(section_vals_type), POINTER :: input, section
109 :
110 14147 : CALL timeset(routineN, handle)
111 14147 : NULLIFY (input, section, cell, dft_control, para_env, sab_orb, matrix_ks, matrix_s, &
112 14147 : positions, soc_matrices, energies, kappas, point, cosines, flat_hamiltonian)
113 14147 : CALL get_qs_env(qs_env, input=input)
114 14147 : section => section_vals_get_subs_vals(input, 'PROPERTIES%SPECTRAL_LOCALIZER')
115 14147 : CALL section_vals_val_get(section, '_SECTION_PARAMETERS_', l_val=active)
116 14147 : logger => cp_get_default_logger()
117 14147 : output_unit = cp_logger_get_default_io_unit(logger)
118 : print_position = BTEST(cp_print_key_should_output(logger%iter_info, input, &
119 14147 : 'DFT%PRINT%AO_MATRICES/POSITION'), cp_p_file)
120 : print_soc = BTEST(cp_print_key_should_output(logger%iter_info, input, &
121 14147 : 'DFT%PRINT%AO_MATRICES/SOC'), cp_p_file)
122 14147 : IF (.NOT. active .AND. .NOT. print_position .AND. .NOT. print_soc) THEN
123 14119 : CALL timestop(handle)
124 14119 : RETURN
125 : END IF
126 : CALL get_qs_env(qs_env, cell=cell, do_kpoints=do_kpoints, dft_control=dft_control, &
127 28 : matrix_s_kp=image_matrices, para_env=para_env, sab_orb=sab_orb)
128 28 : matrix_s => image_matrices(:, 1)
129 28 : CALL section_vals_val_get(section, 'FORMULATION', i_val=formulation)
130 28 : periodic_localizer = active .AND. formulation == 1
131 72 : IF ((print_position .OR. print_soc) .AND. (ANY(cell%perd /= 0) .OR. do_kpoints)) THEN
132 0 : CPABORT('AO_MATRICES POSITION/SOC requires finite systems without KPOINTS.')
133 : END IF
134 28 : IF (do_kpoints .AND. .NOT. periodic_localizer) THEN
135 0 : CPABORT('Explicit KPOINTS require FORMULATION PERIODIC for localizers.')
136 : END IF
137 72 : IF (.NOT. periodic_localizer .AND. ANY(cell%perd /= 0)) THEN
138 0 : CPABORT('AO position/localizer requires PERIODIC NONE without KPOINTS.')
139 : END IF
140 28 : CALL section_vals_val_get(section, 'DIMENSION', i_val=ndim)
141 28 : IF (ndim /= 2 .AND. ndim /= 3) CPABORT('SPECTRAL_LOCALIZER DIMENSION must be 2 or 3.')
142 28 : CALL section_vals_val_get(section, 'PLANE', i_val=plane, explicit=explicit_scale)
143 28 : IF (active .AND. ndim == 3 .AND. explicit_scale) THEN
144 0 : CPABORT('PLANE is only valid for DIMENSION 2.')
145 : END IF
146 28 : CALL section_vals_val_get(section, 'MP_GRID', i_vals=grid, explicit=explicit_scale)
147 112 : IF (ANY(grid < 1)) CPABORT('SPECTRAL_LOCALIZER MP_GRID must be positive.')
148 28 : IF (.NOT. periodic_localizer .AND. explicit_scale) THEN
149 0 : CPABORT('SPECTRAL_LOCALIZER MP_GRID requires FORMULATION PERIODIC.')
150 : END IF
151 112 : IF (ANY(grid /= 1 .AND. cell%perd == 0)) THEN
152 0 : CPABORT('Localizer MP_GRID must be one in nonperiodic directions.')
153 : END IF
154 76 : torus = periodic_localizer .AND. (do_kpoints .OR. ANY(grid /= 1))
155 28 : directions = [1, 2, 3]
156 28 : IF (ndim == 2) THEN
157 : SELECT CASE (plane)
158 : CASE (1)
159 42 : directions(:2) = [1, 2]
160 : CASE (2)
161 18 : directions(:2) = [1, 3]
162 : CASE (3)
163 26 : directions(:2) = [2, 3]
164 : END SELECT
165 : END IF
166 28 : CALL section_vals_val_get(input, 'DFT%POISSON%PERIODIC', i_val=poisson_periodic)
167 28 : IF (periodic_localizer) THEN
168 14 : IF (ndim == 3) THEN
169 16 : IF (ANY(cell%perd == 0) .OR. poisson_periodic /= use_perd_xyz) THEN
170 0 : CPABORT('PERIODIC DIMENSION 3 requires CELL and POISSON PERIODIC XYZ.')
171 : END IF
172 : ELSE
173 60 : IF (COUNT(cell%perd /= 0) /= 2 .OR. ANY(cell%perd(directions(:2)) == 0)) THEN
174 0 : CPABORT('PERIODIC localizer requires two periodic CELL axes matching PLANE.')
175 : END IF
176 10 : IF (poisson_periodic /= plane_periodicity(plane)) THEN
177 0 : CPABORT('PERIODIC localizer requires matching CELL and POISSON periodicity.')
178 : END IF
179 : END IF
180 14 : CALL section_vals_val_get(section, 'KAPPA', explicit=explicit_scale)
181 14 : IF (explicit_scale) CPABORT('Use ETA, not KAPPA, for FORMULATION PERIODIC.')
182 14 : ELSE IF (active .AND. poisson_periodic /= use_perd_none) THEN
183 0 : CPABORT('SPECTRAL_LOCALIZER requires POISSON PERIODIC NONE.')
184 : END IF
185 28 : IF (active .AND. .NOT. periodic_localizer) THEN
186 14 : CALL section_vals_val_get(section, 'ETA', explicit=explicit_scale)
187 14 : IF (explicit_scale) CPABORT('Use KAPPA, not ETA, for FORMULATION FINITE.')
188 : END IF
189 28 : SELECT CASE (dft_control%qs_control%method_id)
190 : CASE (do_method_gpw, do_method_gapw)
191 : CASE DEFAULT
192 28 : CPABORT('AO position/localizer currently supports GPW and GAPW only.')
193 : END SELECT
194 28 : IF (active) THEN
195 28 : IF (.NOT. qs_env%scf_convergence_available .OR. .NOT. qs_env%scf_converged) THEN
196 0 : CPABORT('SPECTRAL_LOCALIZER requires a converged SCF calculation.')
197 : END IF
198 28 : CALL cite_reference(Dixon2023)
199 28 : IF (periodic_localizer) CALL cite_reference(Doll2025)
200 28 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
201 28 : CALL section_vals_val_get(section, 'MAX_AO', i_val=max_ao)
202 28 : IF (max_ao < 1 .OR. nao > max_ao) THEN
203 0 : CPABORT('SPECTRAL_LOCALIZER exceeds the MAX_AO limit.')
204 : END IF
205 112 : DO axis = 1, 3
206 84 : IF (grid(axis) > max_ao/nao) THEN
207 0 : CPABORT('SPECTRAL_LOCALIZER torus exceeds the MAX_AO limit.')
208 : END IF
209 112 : nao = nao*grid(axis)
210 : END DO
211 : END IF
212 :
213 28 : IF (periodic_localizer .AND. .NOT. torus) THEN
214 2 : CALL periodic_ao_moments(qs_env, directions(:ndim), cosines, positions)
215 26 : ELSE IF (.NOT. periodic_localizer .AND. (active .OR. print_position)) THEN
216 14 : CALL finite_ao_positions(qs_env, positions)
217 : END IF
218 28 : IF (print_position) THEN
219 10 : CALL section_vals_val_get(input, 'DFT%PRINT%AO_MATRICES%NDIGITS', i_val=after)
220 10 : CALL section_vals_val_get(input, 'DFT%PRINT%AO_MATRICES%OMIT_HEADERS', l_val=omit_headers)
221 10 : after = MIN(MAX(after, 1), 16)
222 40 : DO axis = 1, 3
223 : iw = cp_print_key_unit_nr(logger, input, 'DFT%PRINT%AO_MATRICES/POSITION', &
224 30 : extension='.Log', middle_name=axes(axis))
225 30 : IF (iw > 0 .AND. .NOT. omit_headers) THEN
226 15 : WRITE (iw, '(A)') 'AO POSITION '//axes(axis)//' [bohr]; coordinate origin; covariant AO matrix'
227 : END IF
228 : CALL cp_dbcsr_write_sparse_matrix(positions(axis)%matrix, 4, after, qs_env, para_env, &
229 30 : output_unit=iw, omit_headers=omit_headers)
230 40 : CALL cp_print_key_finished_output(iw, logger, input, 'DFT%PRINT%AO_MATRICES/POSITION')
231 : END DO
232 : END IF
233 28 : IF (print_soc) THEN
234 0 : CALL finite_ao_soc(qs_env, soc_matrices)
235 0 : CALL section_vals_val_get(input, 'DFT%PRINT%AO_MATRICES%NDIGITS', i_val=after)
236 0 : CALL section_vals_val_get(input, 'DFT%PRINT%AO_MATRICES%OMIT_HEADERS', l_val=omit_headers)
237 0 : after = MIN(MAX(after, 1), 16)
238 0 : DO axis = 1, 3
239 : iw = cp_print_key_unit_nr(logger, input, 'DFT%PRINT%AO_MATRICES/SOC', &
240 0 : extension='.Log', middle_name='SOC_'//axes(axis))
241 0 : IF (iw > 0 .AND. .NOT. omit_headers) THEN
242 0 : WRITE (iw, '(A)') 'AO SOC '//axes(axis)//' [hartree]; real antisymmetric component; multiply by i'
243 : END IF
244 : CALL cp_dbcsr_write_sparse_matrix(soc_matrices(axis, 1)%matrix, 4, after, qs_env, para_env, &
245 0 : output_unit=iw, omit_headers=omit_headers)
246 0 : CALL cp_print_key_finished_output(iw, logger, input, 'DFT%PRINT%AO_MATRICES/SOC')
247 : END DO
248 : END IF
249 28 : IF (.NOT. active) THEN
250 0 : IF (ASSOCIATED(positions)) CALL dbcsr_deallocate_matrix_set(positions)
251 0 : IF (ASSOCIATED(soc_matrices)) CALL dbcsr_deallocate_matrix_set(soc_matrices)
252 0 : CALL timestop(handle)
253 0 : RETURN
254 : END IF
255 :
256 28 : CALL section_vals_val_get(section, 'INVARIANT', i_val=topology)
257 28 : aii = topology == 1
258 28 : IF (ndim == 3 .AND. .NOT. aii) CPABORT('DIMENSION 3 requires INVARIANT Z2.')
259 28 : IF (aii .AND. ndim == 2) CALL cite_reference(Wong2026)
260 28 : IF (ndim == 3) CALL cite_reference(Loring2015)
261 28 : CALL section_vals_val_get(section, 'SOC', l_val=do_soc)
262 28 : CALL section_vals_val_get(section, 'TIME_REVERSAL', l_val=time_reversal)
263 28 : CALL section_vals_val_get(section, 'SPIN_CHANNEL', i_val=spin_channel)
264 28 : IF (aii .AND. (.NOT. do_soc .OR. .NOT. time_reversal)) THEN
265 0 : CPABORT('SPECTRAL_LOCALIZER INVARIANT Z2 requires SOC T and TIME_REVERSAL T.')
266 : END IF
267 28 : CALL get_qs_env(qs_env, matrix_ks_kp=image_matrices)
268 28 : matrix_ks => image_matrices(:, 1)
269 28 : IF (spin_channel < 1 .OR. spin_channel > SIZE(matrix_ks)) THEN
270 0 : CPABORT('SPECTRAL_LOCALIZER SPIN_CHANNEL is unavailable.')
271 : END IF
272 28 : IF (do_soc .AND. (SIZE(matrix_ks) /= 1 .OR. spin_channel /= 1)) THEN
273 0 : CPABORT('SPECTRAL_LOCALIZER SOC requires a restricted SCF and SPIN_CHANNEL 1.')
274 : END IF
275 28 : CALL section_vals_val_get(section, 'ENERGY', r_vals=energies)
276 28 : IF (periodic_localizer) THEN
277 14 : CALL section_vals_val_get(section, 'ETA', r_vals=kappas)
278 : ELSE
279 14 : CALL section_vals_val_get(section, 'KAPPA', r_vals=kappas)
280 : END IF
281 100 : IF (.NOT. ALL(ieee_is_finite(energies)) .OR. .NOT. ALL(ieee_is_finite(kappas))) THEN
282 0 : CPABORT('SPECTRAL_LOCALIZER requires finite ENERGY and KAPPA/ETA.')
283 : END IF
284 64 : IF (ANY(kappas <= 0.0_dp)) THEN
285 0 : CPABORT('SPECTRAL_LOCALIZER requires positive KAPPA/ETA.')
286 : END IF
287 28 : CALL section_vals_val_get(section, 'EPS_GAP', r_val=eps_gap)
288 28 : CALL section_vals_val_get(section, 'EPS_METRIC', r_val=eps_metric)
289 28 : CALL section_vals_val_get(section, 'MATRIX_TOLERANCE', r_val=matrix_tolerance)
290 112 : IF (.NOT. ALL(ieee_is_finite([eps_gap, eps_metric, matrix_tolerance]))) THEN
291 0 : CPABORT('SPECTRAL_LOCALIZER requires finite numerical tolerances.')
292 : END IF
293 : IF (eps_gap < 0.0_dp .OR. eps_metric <= 0.0_dp .OR. eps_metric >= 1.0_dp .OR. &
294 28 : matrix_tolerance <= 0.0_dp .OR. matrix_tolerance >= 1.0_dp) THEN
295 0 : CPABORT('SPECTRAL_LOCALIZER: invalid EPS_GAP, EPS_METRIC or MATRIX_TOLERANCE.')
296 : END IF
297 28 : CALL section_vals_val_get(section, 'SPECTRAL_FLATTENING', l_val=flatten)
298 28 : CALL section_vals_val_get(section, 'FLATTENING_SCALE', r_val=flatten_scale)
299 28 : CALL section_vals_val_get(section, 'EPS_FLATTENING', r_val=eps_flattening)
300 28 : CALL section_vals_val_get(section, 'EPS_ELECTRONIC_GAP', r_val=eps_electronic_gap)
301 28 : IF (flatten) THEN
302 24 : IF (.NOT. ALL(ieee_is_finite([flatten_scale, eps_flattening, eps_electronic_gap]))) THEN
303 0 : CPABORT('SPECTRAL_FLATTENING requires finite scales and tolerances.')
304 : END IF
305 6 : IF (flatten_scale <= 0.0_dp .OR. eps_flattening <= 0.0_dp .OR. eps_flattening >= 1.0_dp .OR. &
306 : eps_electronic_gap < 0.0_dp) THEN
307 0 : CPABORT('SPECTRAL_FLATTENING: invalid scale or tolerance.')
308 : END IF
309 : END IF
310 28 : CALL section_vals_val_get(section, 'POSITION', n_rep_val=nposition)
311 28 : n = nao
312 28 : IF (do_soc) n = 2*nao
313 28 : IF (torus) THEN
314 12 : NULLIFY (matrix_ks, matrix_s)
315 : CALL localizer_torus_operators(qs_env, grid, directions(:ndim), do_soc, matrix_ks, matrix_s, &
316 12 : cosines, positions, soc_matrices)
317 16 : ELSE IF (do_soc) THEN
318 6 : IF (.NOT. ASSOCIATED(soc_matrices)) CALL finite_ao_soc(qs_env, soc_matrices)
319 : END IF
320 28 : IF (output_unit > 0) THEN
321 14 : IF (periodic_localizer) THEN
322 7 : WRITE (output_unit, '(/,T2,A)') 'SPECTRAL_LOCALIZER| Periodic Born-von-Karman DBCSR localizer'
323 28 : WRITE (output_unit, '(T2,A,3I6)') 'SPECTRAL_LOCALIZER| Full analysis mesh: ', grid
324 7 : WRITE (output_unit, '(T2,A)') 'SPECTRAL_LOCALIZER| Analytic cos(G.r)/sin(G.r); no finite-coordinate substitution.'
325 7 : WRITE (output_unit, '(T2,A)') 'SPECTRAL_LOCALIZER| Finite AO torus diagnostic; converge basis, supercell and ETA.'
326 : ELSE
327 7 : WRITE (output_unit, '(/,T2,A)') 'SPECTRAL_LOCALIZER| Finite-system DBCSR localizer'
328 : END IF
329 14 : IF (ndim == 3) THEN
330 3 : WRITE (output_unit, '(T2,A)') 'SPECTRAL_LOCALIZER| LAPACK real LU/QR sign; independent generalized chiral gap.'
331 : ELSE
332 11 : WRITE (output_unit, '(T2,A)') 'SPECTRAL_LOCALIZER| LAPACK LDL inertia with eigenvalue validation.'
333 : END IF
334 25 : WRITE (output_unit, '(T2,A,I0)') 'SPECTRAL_LOCALIZER| Localizer dimension: ', MERGE(4*n, 2*n, ndim == 3)
335 14 : IF (ndim == 3) THEN
336 3 : WRITE (output_unit, '(T2,A)') 'SPECTRAL_LOCALIZER| Three-dimensional AII; chiral determinant sign.'
337 : ELSE
338 33 : WRITE (output_unit, '(T2,A,2A)') 'SPECTRAL_LOCALIZER| Oriented plane: ', axes(directions(:2))
339 : END IF
340 14 : IF (.NOT. periodic_localizer) THEN
341 7 : WRITE (output_unit, '(T2,A)') 'SPECTRAL_LOCALIZER| Check basis, volume and KAPPA convergence.'
342 : END IF
343 14 : IF (.NOT. do_soc) WRITE (output_unit, '(T2,A)') &
344 5 : 'SPECTRAL_LOCALIZER| Real scalar Hamiltonian: no nonzero class A Chern index is expected.'
345 : END IF
346 28 : IF (flatten) THEN
347 6 : ALLOCATE (flat_hamiltonian)
348 6 : IF (do_soc) THEN
349 : CALL localizer_ao_operators(matrix_ks(spin_channel)%matrix, matrix_s(1)%matrix, positions(:0), &
350 4 : physical_h, physical_metric, soc_matrices)
351 : ELSE
352 : CALL localizer_ao_operators(matrix_ks(spin_channel)%matrix, matrix_s(1)%matrix, positions(:0), &
353 2 : physical_h, physical_metric)
354 : END IF
355 6 : IF (output_unit > 0) THEN
356 3 : WRITE (output_unit, '(T2,A,ES23.14)') 'SPECTRAL_LOCALIZER| Flattening scale [hartree]: ', flatten_scale
357 3 : WRITE (output_unit, '(T2,A)') 'SPECTRAL_LOCALIZER| Gaps and indices refer to the flattened Hamiltonian.'
358 3 : IF (periodic_localizer) THEN
359 2 : WRITE (output_unit, '(T2,A)') 'SPECTRAL_LOCALIZER| Flattening in complete primitive Bloch blocks.'
360 : WRITE (output_unit, '(T2,A)') &
361 2 : 'SPECTRAL_LOCALIZER| Residuals: translation, metric orthonormality, eigenproblem, Hermiticity.'
362 : END IF
363 : END IF
364 : END IF
365 64 : DO ie = 1, SIZE(energies)
366 36 : build_energy = energies(ie)
367 36 : IF (flatten) THEN
368 8 : CALL timeset('localizer_spectral_flattening', flatten_handle)
369 8 : IF (periodic_localizer) THEN
370 : CALL localizer_flatten_bloch(physical_h(1), physical_metric, grid, MERGE(2, 1, do_soc), &
371 : energies(ie), flatten_scale, eps_flattening, eps_metric, eps_electronic_gap, &
372 : para_env, flat_hamiltonian, electronic_gap(1), physical_ratio, &
373 4 : flatten_residuals, status)
374 4 : electronic_gap(2) = electronic_gap(1)
375 : ELSE
376 : CALL localizer_flatten_checked(physical_h(1), physical_metric, energies(ie), flatten_scale, eps_flattening, &
377 : eps_metric, eps_electronic_gap, para_env, flat_hamiltonian, &
378 4 : electronic_gap, physical_ratio, flatten_residuals, status)
379 : END IF
380 8 : CALL timestop(flatten_handle)
381 8 : IF (output_unit > 0) THEN
382 4 : WRITE (output_unit, '(T2,A,ES23.14)') 'SPECTRAL_LOCALIZER| Flattening query [hartree]: ', energies(ie)
383 : WRITE (output_unit, '(T2,A,2ES23.14)') &
384 4 : 'SPECTRAL_LOCALIZER| Original electronic gap bracket [hartree]: ', electronic_gap
385 4 : WRITE (output_unit, '(T2,A,ES23.14)') 'SPECTRAL_LOCALIZER| Original metric ratio bound: ', physical_ratio
386 4 : WRITE (output_unit, '(T2,A,4ES14.5)') 'SPECTRAL_LOCALIZER| Flattening residuals: ', flatten_residuals
387 : END IF
388 0 : SELECT CASE (status)
389 : CASE (0)
390 : CASE (1)
391 0 : CPABORT('SPECTRAL_FLATTENING: electronic gap unresolved; no index.')
392 : CASE (-2)
393 0 : CPABORT('SPECTRAL_FLATTENING: invalid original AO metric.')
394 : CASE (-6)
395 0 : CPABORT('SPECTRAL_FLATTENING: sign iteration or structure failed.')
396 : CASE (-7)
397 0 : CPABORT('SPECTRAL_FLATTENING: broken torus translations.')
398 : CASE DEFAULT
399 8 : CPABORT('SPECTRAL_FLATTENING: input or gap factorization failed.')
400 : END SELECT
401 8 : build_energy = 0.0_dp
402 : END IF
403 88 : DO ip = 1, nposition
404 52 : CALL section_vals_val_get(section, 'POSITION', i_rep_val=ip, r_vals=point)
405 208 : IF (.NOT. ALL(ieee_is_finite(point))) THEN
406 0 : CPABORT('SPECTRAL_LOCALIZER requires finite POSITION.')
407 : END IF
408 162 : local_point(:ndim) = point(directions(:ndim))
409 52 : IF (periodic_localizer) THEN
410 46 : DO axis = 1, ndim
411 : local_point(axis) = twopi*DOT_PRODUCT(cell%h_inv(directions(axis), :), point)/ &
412 142 : REAL(grid(directions(axis)), dp)
413 : END DO
414 : END IF
415 172 : DO ik = 1, SIZE(kappas)
416 84 : IF (ndim == 3) THEN
417 : CALL localizer_dbcsr_3d_build(matrix_ks(spin_channel)%matrix, matrix_s(1)%matrix, positions, &
418 : local_point, build_energy, kappas(ik), localizer, metric, soc_matrices, &
419 6 : shifted_hamiltonian=flat_hamiltonian, periodic_cosines=cosines)
420 78 : ELSE IF (do_soc) THEN
421 : CALL localizer_dbcsr_build(matrix_ks(spin_channel)%matrix, matrix_s(1)%matrix, &
422 : positions(directions(1))%matrix, positions(directions(2))%matrix, &
423 : local_point(:2), build_energy, kappas(ik), localizer, metric, soc_matrices, cosines, &
424 12 : shifted_hamiltonian=flat_hamiltonian)
425 : ELSE
426 : CALL localizer_dbcsr_build(matrix_ks(spin_channel)%matrix, matrix_s(1)%matrix, &
427 : positions(directions(1))%matrix, positions(directions(2))%matrix, &
428 : local_point(:2), build_energy, kappas(ik), localizer, metric, &
429 66 : periodic_cosines=cosines, shifted_hamiltonian=flat_hamiltonian)
430 : END IF
431 84 : CALL localizer_dbcsr_dense(localizer, para_env, local_dense)
432 84 : CALL localizer_dbcsr_dense(metric, para_env, metric_dense)
433 84 : IF (para_env%is_source()) THEN
434 42 : IF (ndim == 3) THEN
435 : CALL localizer_3d_matrix_evaluate(local_dense, metric_dense(:n, :n), &
436 : matrix_tolerance, eps_metric, eps_gap, invariant, gap, &
437 3 : metric_ratio, status)
438 : ELSE
439 : CALL localizer_matrix_evaluate(local_dense, metric_dense(:n, :n), aii, &
440 39 : matrix_tolerance, eps_metric, eps_gap, invariant, gap, metric_ratio, status)
441 : END IF
442 42 : IF (status < 0 .AND. output_unit > 0 .AND. ndim == 2) THEN
443 : WRITE (output_unit, '(T2,A,I0,2ES14.5)') &
444 0 : 'SPECTRAL_LOCALIZER| Failure status and L/B Hermiticity residuals: ', status, &
445 0 : MAXVAL(ABS(local_dense - CONJG(TRANSPOSE(local_dense)))), &
446 0 : MAXVAL(ABS(metric_dense - CONJG(TRANSPOSE(metric_dense))))
447 : END IF
448 : END IF
449 84 : DEALLOCATE (local_dense, metric_dense)
450 84 : CALL localizer_dbcsr_release(localizer)
451 84 : CALL localizer_dbcsr_release(metric)
452 84 : CALL para_env%bcast(status, para_env%source)
453 84 : IF (status < 0 .OR. status > 1) THEN
454 0 : SELECT CASE (status)
455 : CASE (-2)
456 0 : CPABORT('SPECTRAL_LOCALIZER: invalid AO metric; check basis/EPS_METRIC.')
457 : CASE (-3)
458 0 : CPABORT('SPECTRAL_LOCALIZER: broken physical TR/AII structure.')
459 : CASE DEFAULT
460 0 : CPABORT('SPECTRAL_LOCALIZER: invalid data or factorization failure.')
461 : END SELECT
462 : END IF
463 84 : CALL para_env%bcast(invariant, para_env%source)
464 84 : CALL para_env%bcast(gap, para_env%source)
465 84 : CALL para_env%bcast(metric_ratio, para_env%source)
466 136 : IF (output_unit > 0) THEN
467 168 : WRITE (output_unit, '(T2,A,3ES23.14)') 'SPECTRAL_LOCALIZER| Position [bohr]: ', point
468 42 : WRITE (output_unit, '(T2,A,ES23.14)') 'SPECTRAL_LOCALIZER| Energy [hartree]: ', energies(ie)
469 42 : IF (periodic_localizer) THEN
470 7 : WRITE (output_unit, '(T2,A,ES23.14)') 'SPECTRAL_LOCALIZER| Eta [hartree]: ', kappas(ik)
471 : ELSE
472 35 : WRITE (output_unit, '(T2,A,ES23.14)') 'SPECTRAL_LOCALIZER| Kappa [hartree/bohr]: ', kappas(ik)
473 : END IF
474 42 : WRITE (output_unit, '(T2,A,ES23.14)') 'SPECTRAL_LOCALIZER| Metric ratio: ', metric_ratio
475 42 : WRITE (output_unit, '(T2,A,ES23.14)') 'SPECTRAL_LOCALIZER| Gap [hartree]: ', gap
476 42 : IF (status == 1) THEN
477 8 : WRITE (output_unit, '(T2,A)') 'SPECTRAL_LOCALIZER| UNRESOLVED: gap or factorization below accuracy; no index.'
478 34 : ELSE IF (aii) THEN
479 9 : WRITE (output_unit, '(T2,A,I0)') 'SPECTRAL_LOCALIZER| Z2 index: ', invariant
480 : ELSE
481 25 : WRITE (output_unit, '(T2,A,I0)') 'SPECTRAL_LOCALIZER| Chern index: ', invariant
482 : END IF
483 : END IF
484 : END DO
485 : END DO
486 64 : IF (flatten) CALL localizer_dbcsr_release(flat_hamiltonian)
487 : END DO
488 28 : IF (flatten) THEN
489 6 : CALL localizer_dbcsr_release(physical_h(1))
490 6 : CALL localizer_dbcsr_release(physical_metric)
491 6 : DEALLOCATE (flat_hamiltonian)
492 : END IF
493 28 : CALL dbcsr_deallocate_matrix_set(positions)
494 28 : IF (ASSOCIATED(cosines)) CALL dbcsr_deallocate_matrix_set(cosines)
495 28 : IF (ASSOCIATED(soc_matrices)) CALL dbcsr_deallocate_matrix_set(soc_matrices)
496 28 : IF (torus) THEN
497 12 : CALL dbcsr_deallocate_matrix_set(matrix_ks)
498 12 : CALL dbcsr_deallocate_matrix_set(matrix_s)
499 : END IF
500 28 : CALL timestop(handle)
501 42749 : END SUBROUTINE qs_scf_post_localizer
502 :
503 : ! **************************************************************************************************
504 : !> \brief Analytic torus coordinates in the reciprocal frame of a Gamma supercell.
505 : !> \param qs_env Quickstep environment
506 : !> \param directions Two or three periodic lattice directions
507 : !> \param cosines Newly allocated cosine matrices, indexed by oriented direction
508 : !> \param sines Newly allocated sine matrices, indexed by lattice direction
509 : ! **************************************************************************************************
510 2 : SUBROUTINE periodic_ao_moments(qs_env, directions, cosines, sines)
511 : TYPE(qs_environment_type), POINTER :: qs_env
512 : INTEGER, INTENT(IN) :: directions(:)
513 : TYPE(dbcsr_p_type), POINTER :: cosines(:), sines(:)
514 :
515 : INTEGER :: j
516 : TYPE(cell_type), POINTER :: cell
517 2 : TYPE(dbcsr_p_type), POINTER :: matrix_s(:)
518 2 : TYPE(neighbor_list_set_p_type), POINTER :: sab_orb(:)
519 :
520 2 : CALL get_qs_env(qs_env, cell=cell, matrix_s=matrix_s, sab_orb=sab_orb)
521 16 : ALLOCATE (cosines(SIZE(directions)), sines(3))
522 8 : DO j = 1, 3
523 6 : ALLOCATE (sines(j)%matrix)
524 6 : CALL dbcsr_create(sines(j)%matrix, template=matrix_s(1)%matrix, name='TORUS SIN')
525 6 : CALL cp_dbcsr_alloc_block_from_nbl(sines(j)%matrix, sab_orb)
526 8 : CALL dbcsr_set(sines(j)%matrix, 0.0_dp)
527 : END DO
528 6 : DO j = 1, SIZE(directions)
529 4 : ALLOCATE (cosines(j)%matrix)
530 4 : CALL dbcsr_create(cosines(j)%matrix, template=matrix_s(1)%matrix, name='TORUS COS')
531 4 : CALL cp_dbcsr_alloc_block_from_nbl(cosines(j)%matrix, sab_orb)
532 4 : CALL dbcsr_set(cosines(j)%matrix, 0.0_dp)
533 : CALL build_berry_moment_matrix(qs_env, cosines(j)%matrix, sines(directions(j))%matrix, &
534 18 : twopi*cell%h_inv(directions(j), :))
535 : END DO
536 2 : END SUBROUTINE periodic_ao_moments
537 :
538 : END MODULE qs_spectral_localizer
|
