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 Interface to the Greenx library
10 : !> \par History
11 : !> 07.2025 Refactored from RPA and BSE modules [Frederick Stein]
12 : ! **************************************************************************************************
13 : MODULE greenx_interface
14 : USE kinds, ONLY: dp
15 : USE cp_log_handling, ONLY: cp_logger_type, &
16 : cp_get_default_logger, &
17 : cp_logger_get_default_io_unit
18 : USE cp_output_handling, ONLY: cp_print_key_unit_nr, &
19 : cp_print_key_finished_output, &
20 : cp_print_key_generate_filename, &
21 : low_print_level, &
22 : medium_print_level
23 : USE input_section_types, ONLY: section_vals_type
24 : USE machine, ONLY: m_flush
25 : USE physcon, ONLY: evolt
26 : #if defined (__GREENX)
27 : USE gx_ac, ONLY: create_thiele_pade, &
28 : evaluate_thiele_pade_at, &
29 : free_params, &
30 : params
31 : USE gx_minimax, ONLY: gx_minimax_grid
32 : #endif
33 :
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : PRIVATE
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'greenx_interface'
41 :
42 : PUBLIC :: greenx_refine_pade, greenx_output_polarizability, greenx_refine_ft, greenx_get_minimax_grid
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief Refines Pade approximants using GreenX, skips this step if GreenX is not available
48 : !> \param e_min ...
49 : !> \param e_max ...
50 : !> \param x_eval ...
51 : !> \param number_of_simulation_steps ...
52 : !> \param number_of_pade_points ...
53 : !> \param logger ...
54 : !> \param ft_section ...
55 : !> \param bse_unit ...
56 : !> \param omega_series ...
57 : !> \param ft_full_series ...
58 : ! **************************************************************************************************
59 0 : SUBROUTINE greenx_refine_pade(e_min, e_max, x_eval, number_of_simulation_steps, number_of_pade_points, &
60 0 : logger, ft_section, bse_unit, omega_series, ft_full_series)
61 : REAL(KIND=dp), INTENT(IN) :: e_min, e_max
62 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: x_eval
63 : INTEGER, INTENT(IN) :: number_of_simulation_steps, number_of_pade_points
64 : TYPE(cp_logger_type), POINTER :: logger
65 : TYPE(section_vals_type), POINTER :: ft_section
66 : INTEGER, INTENT(IN) :: bse_unit
67 : REAL(KIND=dp), DIMENSION(number_of_simulation_steps + 2), INTENT(INOUT) :: omega_series
68 : REAL(KIND=dp), DIMENSION(6, number_of_simulation_steps + 2), INTENT(INOUT) :: ft_full_series
69 : #if defined (__GREENX)
70 : INTEGER :: i, ft_unit
71 0 : COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: omega_complex, &
72 0 : moments_ft_complex
73 0 : COMPLEX(kind=dp), DIMENSION(:, :), ALLOCATABLE :: moments_eval_complex
74 :
75 : ! Report Padé refinement
76 0 : IF (bse_unit > 0) WRITE (bse_unit, '(A10,A27,E23.8E3,E20.8E3)') &
77 0 : " PADE_FT| ", "Evaluation grid bounds [eV]", e_min, e_max
78 0 : ALLOCATE (omega_complex(number_of_simulation_steps + 2))
79 0 : ALLOCATE (moments_ft_complex(number_of_simulation_steps + 2))
80 0 : ALLOCATE (moments_eval_complex(3, number_of_pade_points))
81 0 : omega_complex(:) = CMPLX(omega_series(:), 0.0, kind=dp)
82 0 : DO i = 1, 3
83 : moments_ft_complex(:) = CMPLX(ft_full_series(2*i - 1, :), &
84 : ft_full_series(2*i, :), &
85 0 : kind=dp)
86 : ! Copy the fitting parameters
87 : ! TODO : Optional direct setting of parameters?
88 0 : CALL greenx_refine_ft(e_min, e_max, omega_complex, moments_ft_complex, x_eval, moments_eval_complex(i, :))
89 : END DO
90 : ! Write into alternative file
91 : ft_unit = cp_print_key_unit_nr(logger, ft_section, extension="_PADE.dat", &
92 0 : file_form="FORMATTED", file_position="REWIND")
93 0 : IF (ft_unit > 0) THEN
94 0 : DO i = 1, number_of_pade_points
95 : WRITE (ft_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') &
96 0 : REAL(x_eval(i)), REAL(moments_eval_complex(1, i)), AIMAG(moments_eval_complex(1, i)), &
97 0 : REAL(moments_eval_complex(2, i)), AIMAG(moments_eval_complex(2, i)), &
98 0 : REAL(moments_eval_complex(3, i)), AIMAG(moments_eval_complex(3, i))
99 : END DO
100 : END IF
101 0 : CALL cp_print_key_finished_output(ft_unit, logger, ft_section)
102 0 : DEALLOCATE (omega_complex)
103 0 : DEALLOCATE (moments_ft_complex)
104 0 : DEALLOCATE (moments_eval_complex)
105 : #else
106 : IF (bse_unit > 0) WRITE (bse_unit, '(A10,A70)') &
107 : " PADE_FT| ", "GreenX library is not available. Refinement is skipped"
108 : MARK_USED(e_min)
109 : MARK_USED(e_max)
110 : MARK_USED(x_eval)
111 : MARK_USED(number_of_simulation_steps)
112 : MARK_USED(number_of_pade_points)
113 : MARK_USED(logger)
114 : MARK_USED(ft_section)
115 : MARK_USED(omega_series)
116 : MARK_USED(ft_full_series)
117 : #endif
118 0 : END SUBROUTINE greenx_refine_pade
119 : ! **************************************************************************************************
120 : !> \brief Outputs the isotropic polarizability tensor element alpha _ ij = mu_i(omega)/E_j(omega),
121 : !> where i and j are provided by the configuration. The tensor element is energy dependent and
122 : !> has real and imaginary parts
123 : !> \param logger ...
124 : !> \param pol_section ...
125 : !> \param bse_unit ...
126 : !> \param pol_elements ...
127 : !> \param x_eval ...
128 : !> \param polarizability_refined ...
129 : ! **************************************************************************************************
130 0 : SUBROUTINE greenx_output_polarizability(logger, pol_section, bse_unit, pol_elements, x_eval, polarizability_refined)
131 : TYPE(cp_logger_type), POINTER :: logger
132 : TYPE(section_vals_type), POINTER :: pol_section
133 : INTEGER, INTENT(IN) :: bse_unit
134 : INTEGER, DIMENSION(:, :), POINTER :: pol_elements
135 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: x_eval
136 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: polarizability_refined
137 : #if defined(__GREENX)
138 : INTEGER :: pol_unit, &
139 : i, k, n_elems
140 :
141 0 : n_elems = SIZE(pol_elements, 1)
142 : ! Print out the refined polarizability to a file
143 : pol_unit = cp_print_key_unit_nr(logger, pol_section, extension="_PADE.dat", &
144 0 : file_form="FORMATTED", file_position="REWIND")
145 : ! Printing for both the stdout and separate file
146 0 : IF (pol_unit > 0) THEN
147 0 : IF (pol_unit == bse_unit) THEN
148 : ! Print the stdout preline
149 0 : WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|"
150 : ELSE
151 : ! Print also the energy in atomic units
152 0 : WRITE (pol_unit, '(A1,A19)', advance="no") "#", "omega [a.u.]"
153 : END IF
154 : ! Common - print the energy in eV
155 0 : WRITE (pol_unit, '(A20)', advance="no") "Energy [eV]"
156 : ! Print a header for each polarizability element
157 0 : DO k = 1, n_elems - 1
158 : WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)', advance="no") &
159 0 : "Real pol.", pol_elements(k, 1), pol_elements(k, 2), &
160 0 : "Imag pol.", pol_elements(k, 1), pol_elements(k, 2)
161 : END DO
162 : WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)') &
163 0 : "Real pol.", pol_elements(n_elems, 1), pol_elements(n_elems, 2), &
164 0 : "Imag pol.", pol_elements(n_elems, 1), pol_elements(n_elems, 2)
165 0 : DO i = 1, SIZE(x_eval)
166 0 : IF (pol_unit == bse_unit) THEN
167 : ! Print the stdout preline
168 0 : WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|"
169 : ELSE
170 : ! omega in a.u.
171 0 : WRITE (pol_unit, '(E20.8E3)', advance="no") REAL(x_eval(i), kind=dp)
172 : END IF
173 : ! Common values
174 0 : WRITE (pol_unit, '(E20.8E3)', advance="no") REAL(x_eval(i), kind=dp)*evolt
175 0 : DO k = 1, n_elems - 1
176 : WRITE (pol_unit, '(E20.8E3,E20.8E3)', advance="no") &
177 0 : REAL(polarizability_refined(k, i)), AIMAG(polarizability_refined(k, i))
178 : END DO
179 : ! Print the final value and advance
180 : WRITE (pol_unit, '(E20.8E3,E20.8E3)') &
181 0 : REAL(polarizability_refined(n_elems, i)), AIMAG(polarizability_refined(n_elems, i))
182 : END DO
183 0 : CALL cp_print_key_finished_output(pol_unit, logger, pol_section)
184 : END IF
185 : #else
186 : MARK_USED(logger)
187 : MARK_USED(pol_section)
188 : MARK_USED(bse_unit)
189 : MARK_USED(pol_elements)
190 : MARK_USED(x_eval)
191 : MARK_USED(polarizability_refined)
192 : #endif
193 0 : END SUBROUTINE greenx_output_polarizability
194 : ! **************************************************************************************************
195 : !> \brief Refines the FT grid using Padé approximants
196 : !> \param fit_e_min ...
197 : !> \param fit_e_max ...
198 : !> \param x_fit Input x-variables
199 : !> \param y_fit Input y-variables
200 : !> \param x_eval Refined x-variables
201 : !> \param y_eval Refined y-variables
202 : !> \param n_pade_opt ...
203 : ! **************************************************************************************************
204 6 : SUBROUTINE greenx_refine_ft(fit_e_min, fit_e_max, x_fit, y_fit, x_eval, y_eval, n_pade_opt)
205 : REAL(kind=dp) :: fit_e_min, &
206 : fit_e_max
207 : COMPLEX(kind=dp), DIMENSION(:) :: x_fit, &
208 : y_fit, &
209 : x_eval, &
210 : y_eval
211 : INTEGER, OPTIONAL :: n_pade_opt
212 : #if defined (__GREENX)
213 : CHARACTER(len=*), PARAMETER :: routineN = 'greenx_refine_ft'
214 :
215 : INTEGER :: fit_start, &
216 : fit_end, &
217 : max_fit, &
218 : n_fit, &
219 : n_pade, &
220 : n_eval, &
221 : i, &
222 : handle, &
223 : unit_nr
224 : TYPE(cp_logger_type), POINTER :: logger
225 : TYPE(params) :: pade_params
226 :
227 6 : CALL timeset(routineN, handle)
228 :
229 : ! Get the sizes from arrays
230 6 : max_fit = SIZE(x_fit)
231 6 : n_eval = SIZE(x_eval)
232 :
233 : ! Search for the fit start and end indices
234 6 : fit_start = -1
235 6 : fit_end = -1
236 : ! Search for the subset of FT points which is within energy limits given by
237 : ! the input
238 : ! Do not search when automatic request of highest energy is made
239 6 : IF (fit_e_max < 0) fit_end = max_fit
240 144 : DO i = 1, max_fit
241 144 : IF (fit_start == -1 .AND. REAL(x_fit(i)) >= fit_e_min) fit_start = i
242 144 : IF (fit_end == -1 .AND. REAL(x_fit(i)) > fit_e_max) fit_end = i - 1
243 144 : IF (fit_start > 0 .AND. fit_end > 0) EXIT
244 : END DO
245 6 : IF (fit_start == -1) fit_start = 1
246 6 : IF (fit_end == -1) fit_end = max_fit
247 6 : n_fit = fit_end - fit_start + 1
248 :
249 6 : n_pade = n_fit/2
250 6 : IF (PRESENT(n_pade_opt)) n_pade = n_pade_opt
251 :
252 : ! Too few FT points (e.g. very short propagation with &FT on) leave n_pade < 1;
253 : ! the Thiele recurrence would then divide by zero. Skip, returning zeros.
254 6 : IF (n_pade < 1) THEN
255 0 : CPWARN("FT deck too short for Padé; raise STEPS or disable &FT.")
256 0 : y_eval(1:n_eval) = CMPLX(0.0, 0.0, kind=dp)
257 0 : CALL timestop(handle)
258 : RETURN
259 : END IF
260 :
261 : ! Warn about a large number of Padé parameters
262 6 : IF (n_pade > 1000) THEN
263 0 : CPWARN("More then 1000 Padé parameters requested - may reduce with FIT_E_MIN/FIT_E_MAX.")
264 : END IF
265 :
266 : ! The Padé order is derived from the FT bins inside [FIT_E_MIN, FIT_E_MAX]; report it so a
267 : ! spectrum's fit is reconstructable from the log. Distinct from the GW AC Padé (nparam_pade).
268 6 : logger => cp_get_default_logger()
269 6 : unit_nr = cp_logger_get_default_io_unit(logger)
270 6 : IF (unit_nr > 0) THEN
271 : WRITE (UNIT=unit_nr, FMT="(T3,A,T45,I6,I8,2F11.4)") &
272 3 : "GREENX FT_PADE| n_pade, n_fit, window [eV]", n_pade, n_fit, &
273 6 : REAL(x_fit(fit_start), kind=dp)*evolt, REAL(x_fit(fit_end), kind=dp)*evolt
274 : END IF
275 :
276 : ! TODO : Symmetry mode settable?
277 : ! Here, we assume that ft corresponds to transform of real trace
278 : pade_params = create_thiele_pade(n_pade, x_fit(fit_start:fit_end), y_fit(fit_start:fit_end), &
279 6 : enforce_symmetry="conjugate")
280 :
281 : ! Check whetner the splice is needed or not
282 6000 : y_eval(1:n_eval) = evaluate_thiele_pade_at(pade_params, x_eval)
283 :
284 6 : CALL free_params(pade_params)
285 6 : CALL timestop(handle)
286 : #else
287 : ! Mark used
288 : MARK_USED(fit_e_min)
289 : MARK_USED(fit_e_max)
290 : MARK_USED(x_fit)
291 : MARK_USED(y_fit)
292 : MARK_USED(x_eval)
293 : MARK_USED(y_eval)
294 : MARK_USED(n_pade_opt)
295 : CPABORT("Calls to GreenX require CP2K to be compiled with support for GreenX.")
296 : #endif
297 6 : END SUBROUTINE greenx_refine_ft
298 :
299 : ! **************************************************************************************************
300 : !> \brief ...
301 : !> \param unit_nr ...
302 : !> \param num_integ_points ...
303 : !> \param emin ...
304 : !> \param emax ...
305 : !> \param tau_tj ...
306 : !> \param tau_wj ...
307 : !> \param regularization_minimax ...
308 : !> \param tj ...
309 : !> \param wj ...
310 : !> \param weights_cos_tf_t_to_w ...
311 : !> \param weights_cos_tf_w_to_t ...
312 : !> \param weights_sin_tf_t_to_w ...
313 : !> \param ierr ...
314 : ! **************************************************************************************************
315 206 : SUBROUTINE greenx_get_minimax_grid(unit_nr, num_integ_points, emin, emax, &
316 : tau_tj, tau_wj, regularization_minimax, &
317 : tj, wj, weights_cos_tf_t_to_w, &
318 : weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, ierr)
319 :
320 : INTEGER, INTENT(IN) :: unit_nr, num_integ_points
321 : REAL(KIND=dp), INTENT(IN) :: emin, emax
322 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
323 : INTENT(OUT) :: tau_tj, tau_wj
324 : REAL(KIND=dp), INTENT(IN) :: regularization_minimax
325 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
326 : INTENT(INOUT) :: tj, wj
327 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
328 : INTENT(OUT) :: weights_cos_tf_t_to_w, &
329 : weights_cos_tf_w_to_t, &
330 : weights_sin_tf_t_to_w
331 : INTEGER, INTENT(OUT) :: ierr
332 : #if defined (__GREENX)
333 : INTEGER :: gi
334 : REAL(KIND=dp) :: cosft_duality_error_greenx, &
335 : max_errors_greenx(3)
336 :
337 : CALL gx_minimax_grid(num_integ_points, Emin, Emax, tau_tj, tau_wj, tj, wj, &
338 : weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, &
339 : max_errors_greenx, cosft_duality_error_greenx, ierr, &
340 : bare_cos_sin_weights=.TRUE., &
341 206 : regularization=regularization_minimax)
342 : ! Factor 4 is hard-coded in the RPA weights in the internal CP2K minimax routines
343 1498 : wj(:) = wj(:)*4.0_dp
344 206 : IF (ierr == 0) THEN
345 76 : IF (unit_nr > 0) THEN
346 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
347 38 : "GREENX MINIMAX_INFO| Number of integration points:", num_integ_points
348 : WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") &
349 38 : "GREENX MINIMAX_INFO| Gap (Emin):", Emin
350 : WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") &
351 38 : "GREENX MINIMAX_INFO| Maximum eigenvalue difference (Emax):", Emax
352 : WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.4)") &
353 38 : "GREENX MINIMAX_INFO| Energy range (Emax/Emin):", Emax/Emin
354 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") &
355 38 : "GREENX MINIMAX_INFO| Frequency grid (scaled):", "Weights", "Abscissas"
356 432 : DO gi = 1, num_integ_points
357 432 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") wj(gi), tj(gi)
358 : END DO
359 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") &
360 38 : "GREENX MINIMAX_INFO| Time grid (scaled):", "Weights", "Abscissas"
361 432 : DO gi = 1, num_integ_points
362 432 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(gi), tau_tj(gi)
363 : END DO
364 38 : CALL m_flush(unit_nr)
365 : END IF
366 : ELSE
367 130 : IF (unit_nr > 0) THEN
368 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75)") &
369 65 : "GREENX MINIMAX_INFO| Grid not available, use internal CP2K grids"
370 65 : CALL m_flush(unit_nr)
371 : END IF
372 130 : IF (ALLOCATED(tau_tj)) THEN
373 130 : DEALLOCATE (tau_tj)
374 : END IF
375 130 : IF (ALLOCATED(tau_wj)) THEN
376 130 : DEALLOCATE (tau_wj)
377 : END IF
378 130 : IF (ALLOCATED(tj)) THEN
379 130 : DEALLOCATE (tj)
380 : END IF
381 130 : IF (ALLOCATED(wj)) THEN
382 130 : DEALLOCATE (wj)
383 : END IF
384 130 : IF (ALLOCATED(weights_cos_tf_t_to_w)) THEN
385 0 : DEALLOCATE (weights_cos_tf_t_to_w)
386 : END IF
387 130 : IF (ALLOCATED(weights_cos_tf_w_to_t)) THEN
388 0 : DEALLOCATE (weights_cos_tf_w_to_t)
389 : END IF
390 130 : IF (ALLOCATED(weights_sin_tf_t_to_w)) THEN
391 0 : DEALLOCATE (weights_sin_tf_t_to_w)
392 : END IF
393 : END IF
394 : #else
395 : ierr = 1
396 : MARK_USED(unit_nr)
397 : MARK_USED(num_integ_points)
398 : MARK_USED(emin)
399 : MARK_USED(emax)
400 : MARK_USED(tau_tj)
401 : MARK_USED(tau_wj)
402 : MARK_USED(regularization_minimax)
403 : MARK_USED(tj)
404 : MARK_USED(wj)
405 : MARK_USED(weights_cos_tf_t_to_w)
406 : MARK_USED(weights_cos_tf_w_to_t)
407 : MARK_USED(weights_sin_tf_t_to_w)
408 : #endif
409 :
410 206 : END SUBROUTINE greenx_get_minimax_grid
411 :
412 : END MODULE greenx_interface
|
