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