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 Finite or periodic quadratic pseudospectrum analysis after SCF.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE qs_quadratic_pseudospectrum
13 : USE bibliography, ONLY: Bairnsfather2026,&
14 : Cerjan2023
15 : USE cell_types, ONLY: cell_transform_input_cartesian,&
16 : cell_transform_input_reciprocal,&
17 : cell_type,&
18 : use_perd_none
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
21 : dbcsr_p_type
22 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_get_default_io_unit,&
25 : cp_logger_type
26 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
27 : cp_print_key_unit_nr
28 : USE ieee_arithmetic, ONLY: ieee_is_finite
29 : USE input_constants, ONLY: do_method_gapw,&
30 : do_method_gpw
31 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
32 : section_vals_type,&
33 : section_vals_val_get
34 : USE kinds, ONLY: dp
35 : USE mathconstants, ONLY: twopi
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE qs_ao_translations, ONLY: ao_translation_matrices
38 : USE qs_environment_types, ONLY: get_qs_env,&
39 : qs_environment_type
40 : USE qs_finite_ao, ONLY: finite_ao_positions,&
41 : finite_ao_soc
42 : USE qs_localizer_torus, ONLY: localizer_torus_operators,&
43 : torus_translation_operators
44 : USE quadratic_pseudospectrum, ONLY: quadratic_dense,&
45 : quadratic_result_type
46 : USE reference_manager, ONLY: cite_reference
47 : USE spectral_localizer_dbcsr, ONLY: localizer_ao_operators,&
48 : localizer_dbcsr_check,&
49 : localizer_dbcsr_dense,&
50 : localizer_dbcsr_release,&
51 : localizer_matrix_type
52 : #include "./base/base_uses.f90"
53 :
54 : IMPLICIT NONE
55 : PRIVATE
56 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_quadratic_pseudospectrum'
57 : PUBLIC :: qs_scf_post_quadratic
58 : CONTAINS
59 : ! **************************************************************************************************
60 : !> \brief Evaluate scans of quadratic gaps and localized joint approximate eigenstates.
61 : !> \param qs_env Converged Quickstep environment
62 : ! **************************************************************************************************
63 14147 : SUBROUTINE qs_scf_post_quadratic(qs_env)
64 : TYPE(qs_environment_type), POINTER :: qs_env
65 :
66 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_quadratic'
67 :
68 : CHARACTER(len=256) :: matrix_error
69 14147 : COMPLEX(dp), ALLOCATABLE :: complex_shift(:), dense_metric(:, :), &
70 14147 : dense_operators(:, :, :), &
71 14147 : temporary(:, :)
72 : INTEGER :: ao, axis, cell_periodic, formulation, handle, i, ie, ik, ip, iw, j, max_ao, &
73 : max_iter, max_space, n, nao, ndim, nobs, npositions, nstates, nt, output_unit, periodic, &
74 : spin, status
75 14147 : INTEGER, ALLOCATABLE :: directions(:)
76 14147 : INTEGER, POINTER :: grid(:)
77 : LOGICAL :: active, do_kpoints, do_soc, explicit_grid, explicit_value, periodic_analysis, &
78 : torus_analysis, translation_analysis, valid
79 14147 : LOGICAL, ALLOCATABLE :: unitary(:)
80 : REAL(dp) :: eps_eigen, eps_metric, &
81 : hermiticity_error, &
82 : max_hermiticity_error, query(3)
83 14147 : REAL(dp), ALLOCATABLE :: angles(:), radii(:), shift(:), &
84 14147 : translations(:, :), weights(:)
85 14147 : REAL(dp), POINTER :: energies(:), kappas(:), point(:)
86 : TYPE(cell_type), POINTER :: cell
87 : TYPE(cp_logger_type), POINTER :: logger
88 14147 : TYPE(dbcsr_p_type), POINTER :: cosines(:), images_h(:, :), &
89 14147 : images_s(:, :), matrix_ks(:), &
90 14147 : matrix_s(:), positions(:), sines(:), &
91 14147 : soc(:, :)
92 : TYPE(dft_control_type), POINTER :: dft_control
93 : TYPE(localizer_matrix_type) :: metric
94 14147 : TYPE(localizer_matrix_type), ALLOCATABLE :: operators(:)
95 : TYPE(mp_para_env_type), POINTER :: para_env
96 14147 : TYPE(quadratic_result_type) :: result
97 : TYPE(section_vals_type), POINTER :: input, section
98 :
99 14147 : CALL timeset(routineN, handle)
100 14147 : CALL get_qs_env(qs_env, input=input)
101 14147 : section => section_vals_get_subs_vals(input, 'PROPERTIES%QUADRATIC_PSEUDOSPECTRUM')
102 14147 : CALL section_vals_val_get(section, '_SECTION_PARAMETERS_', l_val=active)
103 14147 : IF (.NOT. active) THEN
104 14131 : CALL timestop(handle)
105 14131 : RETURN
106 : END IF
107 16 : IF (.NOT. qs_env%scf_convergence_available .OR. .NOT. qs_env%scf_converged) THEN
108 0 : CPABORT('QUADRATIC_PSEUDOSPECTRUM requires converged SCF.')
109 : END IF
110 : CALL get_qs_env(qs_env, matrix_s_kp=images_s, matrix_ks_kp=images_h, para_env=para_env, &
111 16 : cell=cell, dft_control=dft_control, do_kpoints=do_kpoints)
112 16 : matrix_s => images_s(:, 1)
113 16 : matrix_ks => images_h(:, 1)
114 16 : SELECT CASE (dft_control%qs_control%method_id)
115 : CASE (do_method_gpw, do_method_gapw)
116 : CASE DEFAULT
117 16 : CPABORT('QUADRATIC_PSEUDOSPECTRUM supports GPW/GAPW only.')
118 : END SELECT
119 16 : CALL section_vals_val_get(section, 'FORMULATION', i_val=formulation)
120 16 : periodic_analysis = formulation == 1
121 16 : translation_analysis = formulation == 2
122 34 : torus_analysis = periodic_analysis .OR. (translation_analysis .AND. ANY(cell%perd /= 0))
123 16 : CALL section_vals_val_get(section, 'MP_GRID', i_vals=grid, explicit=explicit_grid)
124 64 : IF (ANY(grid < 1)) CPABORT('Quadratic MP_GRID must be positive.')
125 16 : IF (.NOT. torus_analysis .AND. explicit_grid) THEN
126 0 : CPABORT('Quadratic MP_GRID requires a periodic property analysis.')
127 : END IF
128 64 : IF (ANY(grid /= 1 .AND. cell%perd == 0)) THEN
129 0 : CPABORT('Quadratic MP_GRID must be one in open directions.')
130 : END IF
131 16 : CALL section_vals_val_get(input, 'DFT%POISSON%PERIODIC', i_val=periodic)
132 16 : IF (torus_analysis) THEN
133 10 : CALL section_vals_val_get(input, 'SUBSYS%CELL%PERIODIC', i_val=cell_periodic)
134 40 : IF (cell_periodic /= periodic .OR. COUNT(cell%perd /= 0) == 0) THEN
135 0 : CPABORT('Quadratic torus needs matching, nonempty CELL/POISSON PERIODIC.')
136 : END IF
137 40 : ndim = COUNT(cell%perd /= 0)
138 30 : ALLOCATE (directions(ndim))
139 40 : directions(:) = PACK([1, 2, 3], cell%perd /= 0)
140 20 : IF (periodic_analysis) THEN
141 4 : nobs = 1 + 2*ndim
142 16 : ALLOCATE (radii(ndim), angles(ndim))
143 14 : DO i = 1, ndim
144 10 : axis = directions(i)
145 44 : radii(i) = REAL(grid(axis), dp)/(twopi*SQRT(SUM(cell%h_inv(axis, :)**2)))
146 : END DO
147 : END IF
148 : ELSE
149 24 : IF (periodic /= use_perd_none .OR. ANY(cell%perd /= 0) .OR. do_kpoints) THEN
150 0 : CPABORT('Finite quadratic analysis requires PERIODIC NONE without KPOINTS.')
151 : END IF
152 : nobs = 4
153 : END IF
154 16 : IF (translation_analysis) THEN
155 8 : CALL section_vals_val_get(section, 'TRANSLATION_VECTOR', n_rep_val=nt)
156 8 : IF (nt < 1) THEN
157 0 : CPABORT('Quadratic TRANSLATION requires at least one TRANSLATION_VECTOR.')
158 : END IF
159 8 : nobs = 1 + nt
160 56 : ALLOCATE (translations(3, nt), complex_shift(nobs), unitary(nobs))
161 26 : unitary = .TRUE.
162 8 : unitary(1) = .FALSE.
163 18 : DO i = 1, nt
164 10 : CALL section_vals_val_get(section, 'TRANSLATION_VECTOR', i_rep_val=i, r_vals=point)
165 40 : IF (.NOT. ALL(ieee_is_finite(point))) THEN
166 0 : CPABORT('Quadratic TRANSLATION_VECTOR must be finite.')
167 : END IF
168 40 : translations(:, i) = point
169 18 : CALL cell_transform_input_cartesian(cell, translations(:, i))
170 : END DO
171 8 : CALL section_vals_val_get(section, 'POSITION', explicit=explicit_value)
172 8 : IF (explicit_value) THEN
173 0 : CPABORT('Quadratic TRANSLATION uses WAVE_VECTOR, not POSITION.')
174 : END IF
175 8 : CALL section_vals_val_get(section, 'KAPPA', explicit=explicit_value)
176 24 : IF (explicit_value) THEN
177 0 : CPABORT('Quadratic TRANSLATION uses TRANSLATION_SCALE, not KAPPA.')
178 : END IF
179 : ELSE
180 8 : CALL section_vals_val_get(section, 'TRANSLATION_VECTOR', explicit=explicit_value)
181 8 : IF (explicit_value) THEN
182 0 : CPABORT('TRANSLATION_VECTOR requires quadratic FORMULATION TRANSLATION.')
183 : END IF
184 8 : CALL section_vals_val_get(section, 'WAVE_VECTOR', explicit=explicit_value)
185 8 : IF (explicit_value) THEN
186 0 : CPABORT('WAVE_VECTOR requires quadratic FORMULATION TRANSLATION.')
187 : END IF
188 8 : CALL section_vals_val_get(section, 'TRANSLATION_SCALE', explicit=explicit_value)
189 8 : IF (explicit_value) THEN
190 0 : CPABORT('TRANSLATION_SCALE requires quadratic FORMULATION TRANSLATION.')
191 : END IF
192 : END IF
193 154 : ALLOCATE (operators(nobs), shift(nobs), weights(nobs))
194 16 : CALL section_vals_val_get(section, 'SPIN_CHANNEL', i_val=spin)
195 16 : IF (spin < 1 .OR. spin > SIZE(matrix_ks)) THEN
196 0 : CPABORT('Quadratic SPIN_CHANNEL is unavailable.')
197 : END IF
198 16 : CALL section_vals_val_get(section, 'SOC', l_val=do_soc)
199 16 : IF (do_soc .AND. spin /= 1) CPABORT('Quadratic SOC requires SPIN_CHANNEL 1.')
200 16 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
201 16 : CALL section_vals_val_get(section, 'MAX_AO', i_val=max_ao)
202 16 : IF (nao > max_ao .OR. max_ao < 1) CPABORT('QUADRATIC_PSEUDOSPECTRUM exceeds MAX_AO.')
203 64 : DO axis = 1, 3
204 48 : IF (grid(axis) > max_ao/nao) CPABORT('Quadratic full AO torus exceeds MAX_AO.')
205 64 : nao = nao*grid(axis)
206 : END DO
207 16 : n = nao
208 16 : IF (do_soc) n = 2*nao
209 16 : CALL section_vals_val_get(section, 'NSTATES', i_val=nstates)
210 16 : max_space = n
211 16 : max_iter = 1
212 16 : CALL section_vals_val_get(section, 'EPS_EIGEN', r_val=eps_eigen)
213 16 : CALL section_vals_val_get(section, 'EPS_METRIC', r_val=eps_metric)
214 16 : CALL section_vals_val_get(section, 'ENERGY', r_vals=energies)
215 16 : CALL section_vals_val_get(section, 'KAPPA', r_vals=kappas)
216 16 : CALL section_vals_val_get(section, 'POSITION', n_rep_val=npositions)
217 16 : IF (translation_analysis) THEN
218 8 : CALL section_vals_val_get(section, 'TRANSLATION_SCALE', r_vals=kappas)
219 8 : CALL section_vals_val_get(section, 'WAVE_VECTOR', n_rep_val=npositions)
220 : END IF
221 16 : IF (nstates < 1 .OR. nstates > n .OR. max_iter < 1) THEN
222 0 : CPABORT('Invalid quadratic state/iteration limits.')
223 : END IF
224 116 : IF (.NOT. ALL(ieee_is_finite([energies, kappas, eps_eigen, eps_metric]))) THEN
225 0 : CPABORT('Quadratic analysis requires finite parameters.')
226 : END IF
227 32 : IF (ANY(kappas <= 0.0_dp) .OR. eps_eigen <= 0.0_dp .OR. eps_metric <= 0.0_dp .OR. eps_metric >= 1.0_dp) THEN
228 0 : CPABORT('Invalid quadratic scales or tolerances.')
229 : END IF
230 16 : CALL cite_reference(Cerjan2023)
231 16 : IF (translation_analysis) CALL cite_reference(Bairnsfather2026)
232 16 : NULLIFY (positions, soc)
233 16 : IF (periodic_analysis) THEN
234 4 : CALL localizer_torus_operators(qs_env, grid, directions, do_soc, matrix_ks, matrix_s, cosines, sines, soc)
235 32 : ALLOCATE (positions(2*ndim))
236 14 : DO i = 1, ndim
237 10 : positions(2*i - 1)%matrix => cosines(i)%matrix
238 14 : positions(2*i)%matrix => sines(directions(i))%matrix
239 : END DO
240 12 : ELSE IF (translation_analysis) THEN
241 8 : IF (torus_analysis) THEN
242 6 : CALL torus_translation_operators(qs_env, grid, translations, do_soc, matrix_ks, matrix_s, positions, soc)
243 : ELSE
244 2 : CALL ao_translation_matrices(qs_env, translations, positions)
245 : END IF
246 : ELSE
247 4 : CALL finite_ao_positions(qs_env, positions)
248 : END IF
249 16 : IF (do_soc .AND. .NOT. torus_analysis) THEN
250 2 : CALL finite_ao_soc(qs_env, soc)
251 : END IF
252 16 : IF (do_soc) THEN
253 4 : CALL localizer_ao_operators(matrix_ks(spin)%matrix, matrix_s(1)%matrix, positions, operators, metric, soc)
254 4 : CALL dbcsr_deallocate_matrix_set(soc)
255 : ELSE
256 12 : CALL localizer_ao_operators(matrix_ks(spin)%matrix, matrix_s(1)%matrix, positions, operators, metric)
257 : END IF
258 16 : IF (periodic_analysis) THEN
259 4 : DEALLOCATE (positions)
260 4 : CALL dbcsr_deallocate_matrix_set(cosines)
261 4 : CALL dbcsr_deallocate_matrix_set(sines)
262 : ELSE
263 12 : CALL dbcsr_deallocate_matrix_set(positions)
264 : END IF
265 16 : IF (torus_analysis) THEN
266 10 : CALL dbcsr_deallocate_matrix_set(matrix_ks)
267 10 : CALL dbcsr_deallocate_matrix_set(matrix_s)
268 : END IF
269 16 : CALL localizer_dbcsr_check(metric, para_env, 1.e-10_dp, valid, max_hermiticity_error)
270 16 : IF (.NOT. valid) CPABORT('Quadratic AO metric is nonfinite or not Hermitian.')
271 74 : DO i = 1, nobs
272 58 : IF (translation_analysis .AND. i > 1) CYCLE
273 48 : CALL localizer_dbcsr_check(operators(i), para_env, 1.e-10_dp, valid, hermiticity_error)
274 48 : max_hermiticity_error = MAX(max_hermiticity_error, hermiticity_error)
275 64 : IF (.NOT. valid) THEN
276 0 : WRITE (matrix_error, '(A,I0,A,ES12.4,A)') 'Quadratic AO operator ', i, &
277 0 : ' is nonfinite or not Hermitian (relative error ', hermiticity_error, &
278 0 : '). Check QS%EPS_PPNL/EPS_PGF_ORB and integration thresholds; no symmetrization is applied.'
279 0 : CPABORT(TRIM(matrix_error))
280 : END IF
281 : END DO
282 80 : ALLOCATE (dense_operators(n, n, nobs))
283 16 : CALL localizer_dbcsr_dense(metric, para_env, dense_metric)
284 74 : DO i = 1, nobs
285 58 : CALL localizer_dbcsr_dense(operators(i), para_env, temporary)
286 45034 : dense_operators(:, :, i) = temporary
287 74 : DEALLOCATE (temporary)
288 : END DO
289 16 : logger => cp_get_default_logger()
290 16 : output_unit = cp_logger_get_default_io_unit(logger)
291 16 : iw = cp_print_key_unit_nr(logger, section, 'STATES', extension='.states')
292 16 : IF (iw > 0) WRITE (iw, '(A)') '# Quadratic eigenstates: S-normalized AO coefficients, spin-major for SOC'
293 31 : IF (iw > 0 .AND. torus_analysis) WRITE (iw, '(A,3I8)') '# Cell-major AO torus within each spin block; grid ', grid
294 16 : IF (output_unit > 0) THEN
295 8 : WRITE (output_unit, '(/,T2,A,I0)') 'QUADRATIC_PSEUDOSPECTRUM| AO/spinor dimension: ', n
296 8 : WRITE (output_unit, '(T2,A,ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Relative Hermiticity error: ', max_hermiticity_error
297 8 : WRITE (output_unit, '(T2,A)') 'QUADRATIC_PSEUDOSPECTRUM| Dense finite-basis reference'
298 8 : WRITE (output_unit, '(T2,A)') 'QUADRATIC_PSEUDOSPECTRUM| Joint approximate states, not a topology index.'
299 8 : IF (translation_analysis) THEN
300 4 : WRITE (output_unit, '(T2,A)') 'QUADRATIC_PSEUDOSPECTRUM| Continuum translations U(a) psi(r) = psi(r+a).'
301 4 : WRITE (output_unit, '(T2,A)') 'QUADRATIC_PSEUDOSPECTRUM| Full translation norms include loss outside the AO span.'
302 4 : WRITE (output_unit, '(T2,A)') 'QUADRATIC_PSEUDOSPECTRUM| Energy residual still uses the projected AO Hamiltonian.'
303 9 : DO i = 1, nt
304 5 : WRITE (output_unit, '(T2,A,I0,A,3ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Translation ', i, &
305 14 : ' internal [bohr]: ', translations(:, i)
306 : END DO
307 : END IF
308 8 : IF (torus_analysis) THEN
309 20 : WRITE (output_unit, '(T2,A,3I6)') 'QUADRATIC_PSEUDOSPECTRUM| Full analysis mesh: ', grid
310 5 : WRITE (output_unit, '(T2,A,3I4)') 'QUADRATIC_PSEUDOSPECTRUM| Periodic lattice directions: ', directions
311 5 : IF (translation_analysis) THEN
312 3 : WRITE (output_unit, '(T2,A)') 'QUADRATIC_PSEUDOSPECTRUM| Periodized translations on the full AO torus.'
313 : ELSE
314 2 : WRITE (output_unit, '(T2,A)') 'QUADRATIC_PSEUDOSPECTRUM| Torus chord metric; open directions are not localized.'
315 : END IF
316 : END IF
317 : END IF
318 36 : DO ip = 1, npositions
319 20 : IF (translation_analysis) THEN
320 10 : CALL section_vals_val_get(section, 'WAVE_VECTOR', i_rep_val=ip, r_vals=point)
321 : ELSE
322 10 : CALL section_vals_val_get(section, 'POSITION', i_rep_val=ip, r_vals=point)
323 : END IF
324 80 : IF (.NOT. ALL(ieee_is_finite(point))) THEN
325 0 : CPABORT('Quadratic POSITION or WAVE_VECTOR must be finite.')
326 : END IF
327 20 : IF (translation_analysis) THEN
328 40 : query = point
329 10 : CALL cell_transform_input_reciprocal(cell, query)
330 24 : DO i = 1, nt
331 66 : complex_shift(i + 1) = EXP(CMPLX(0.0_dp, DOT_PRODUCT(query, translations(:, i)), dp))
332 : END DO
333 10 : ELSE IF (periodic_analysis) THEN
334 16 : query = point
335 4 : CALL cell_transform_input_cartesian(cell, query)
336 14 : DO i = 1, ndim
337 10 : axis = directions(i)
338 40 : angles(i) = twopi*MODULO(DOT_PRODUCT(cell%h_inv(axis, :), query)/REAL(grid(axis), dp), 1.0_dp)
339 10 : shift(2*i) = COS(angles(i))
340 14 : shift(2*i + 1) = SIN(angles(i))
341 : END DO
342 : ELSE
343 24 : shift(2:) = point
344 : END IF
345 56 : DO ik = 1, SIZE(kappas)
346 64 : DO ie = 1, SIZE(energies)
347 24 : shift(1) = energies(ie)
348 112 : weights = kappas(ik)
349 24 : weights(1) = 1.0_dp
350 24 : IF (periodic_analysis) THEN
351 14 : DO i = 1, ndim
352 34 : weights(2*i:2*i + 1) = kappas(ik)*radii(i)
353 : END DO
354 : END IF
355 24 : IF (para_env%is_source()) THEN
356 12 : IF (translation_analysis) THEN
357 5 : complex_shift(1) = energies(ie)
358 : CALL quadratic_dense(dense_operators, dense_metric, complex_shift, weights, nstates, .FALSE., max_space, &
359 5 : max_iter, eps_eigen, eps_metric, RESULT, status, unitary)
360 : ELSE
361 : CALL quadratic_dense(dense_operators, dense_metric, shift, weights, nstates, .FALSE., max_space, &
362 7 : max_iter, eps_eigen, eps_metric, RESULT, status)
363 : END IF
364 : END IF
365 24 : CALL para_env%bcast(status, para_env%source)
366 24 : IF (status /= 0) THEN
367 0 : IF (status == 1) THEN
368 0 : IF (output_unit > 0) THEN
369 : WRITE (output_unit, '(T2,A,I0,A,ES23.14)') &
370 0 : 'QUADRATIC_PSEUDOSPECTRUM| Unconverged after ', result%iterations, &
371 0 : ' iterations; maximum residual: ', MAXVAL(result%residuals)
372 : END IF
373 0 : CPABORT('Quadratic eigensolver unconverged.')
374 : END IF
375 0 : CPABORT('Quadratic analysis failed numerical checks.')
376 : END IF
377 24 : IF (output_unit > 0) THEN
378 12 : IF (translation_analysis) THEN
379 20 : WRITE (output_unit, '(T2,A,3ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Wave vector [bohr^-1]: ', point
380 5 : WRITE (output_unit, '(T2,A,ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Translation scale [hartree]: ', kappas(ik)
381 : ELSE
382 28 : WRITE (output_unit, '(T2,A,3ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Position [bohr]: ', point
383 7 : WRITE (output_unit, '(T2,A,ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Kappa [hartree/bohr]: ', kappas(ik)
384 : END IF
385 12 : WRITE (output_unit, '(T2,A,ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Energy [hartree]: ', energies(ie)
386 12 : WRITE (output_unit, '(T2,A,I0)') 'QUADRATIC_PSEUDOSPECTRUM| Ritz iterations: ', result%iterations
387 41 : DO j = 1, nstates
388 29 : WRITE (output_unit, '(T2,A,I0)') 'QUADRATIC_PSEUDOSPECTRUM| State: ', j
389 29 : WRITE (output_unit, '(T2,A,ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Gap [hartree]: ', &
390 58 : SQRT(result%eigenvalues(j))
391 29 : WRITE (output_unit, '(T2,A,ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Energy residual [hartree]: ', &
392 58 : result%errors(1, j)
393 29 : IF (translation_analysis) THEN
394 30 : DO i = 2, nobs
395 17 : WRITE (output_unit, '(T2,A,I0,A,3ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Translation residual ', i - 1, &
396 34 : ' full/projected/leakage: ', result%errors(i, j), result%projected_errors(i, j), result%leakage(i, j)
397 17 : WRITE (output_unit, '(T2,A,I0,A,2ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Translation expectation ', i - 1, &
398 47 : ' real/imag: ', result%complex_means(i, j)
399 : END DO
400 16 : ELSE IF (periodic_analysis) THEN
401 14 : DO i = 1, ndim
402 10 : WRITE (output_unit, '(T2,A,I1,A,ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Chord residual axis ', &
403 40 : directions(i), ' [bohr]: ', radii(i)*SQRT(SUM(result%errors(2*i:2*i + 1, j)**2))
404 10 : WRITE (output_unit, '(T2,A,I1,A,2ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Mean cos/sin axis ', &
405 24 : directions(i), ': ', result%means(2*i:2*i + 1, j)
406 : END DO
407 : ELSE
408 12 : WRITE (output_unit, '(T2,A,3ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Position residuals [bohr]: ', &
409 24 : result%errors(2:4, j)
410 : END IF
411 29 : WRITE (output_unit, '(T2,A,ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Eigenpair residual [hartree^2]: ', &
412 58 : result%residuals(j)
413 29 : WRITE (output_unit, '(T2,A,ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Mean energy [hartree]: ', result%means(1, j)
414 41 : IF (.NOT. periodic_analysis .AND. .NOT. translation_analysis) THEN
415 12 : WRITE (output_unit, '(T2,A,3ES23.14)') 'QUADRATIC_PSEUDOSPECTRUM| Mean position [bohr]: ', &
416 24 : result%means(2:4, j)
417 : END IF
418 : END DO
419 : END IF
420 44 : IF (iw > 0) THEN
421 12 : IF (translation_analysis) THEN
422 20 : WRITE (iw, '(A,3I8,A,5ES24.15)') '# Query ', ip, ie, ik, ' E/kx/ky/kz/scale ', energies(ie), point, kappas(ik)
423 : ELSE
424 28 : WRITE (iw, '(A,3I8,A,5ES24.15)') '# Query ', ip, ie, ik, ' E/x/y/z/kappa ', energies(ie), point, kappas(ik)
425 : END IF
426 41 : DO j = 1, nstates
427 29 : WRITE (iw, '(A,I8,A,ES24.15)') '# State ', j, ' lambda [hartree^2] ', result%eigenvalues(j)
428 579 : DO ao = 1, n
429 567 : WRITE (iw, '(I10,2ES25.16)') ao, REAL(result%states(ao, j), dp), AIMAG(result%states(ao, j))
430 : END DO
431 : END DO
432 : END IF
433 : END DO
434 : END DO
435 : END DO
436 16 : CALL cp_print_key_finished_output(iw, logger, section, 'STATES')
437 74 : DO i = 1, nobs
438 74 : CALL localizer_dbcsr_release(operators(i))
439 : END DO
440 16 : CALL localizer_dbcsr_release(metric)
441 16 : CALL timestop(handle)
442 42569 : END SUBROUTINE qs_scf_post_quadratic
443 : END MODULE qs_quadratic_pseudospectrum
|
