(git:98357aa)
Loading...
Searching...
No Matches
greenx_interface.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 Interface to the Greenx library
10!> \par History
11!> 07.2025 Refactored from RPA and BSE modules [Frederick Stein]
12! **************************************************************************************************
13MODULE 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
44CONTAINS
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 SUBROUTINE greenx_refine_pade(e_min, e_max, x_eval, number_of_simulation_steps, number_of_pade_points, &
60 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 COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: omega_complex, &
72 moments_ft_complex
73 COMPLEX(kind=dp), DIMENSION(:, :), ALLOCATABLE :: moments_eval_complex
74
75 ! Report Padé refinement
76 IF (bse_unit > 0) WRITE (bse_unit, '(A10,A27,E23.8E3,E20.8E3)') &
77 " PADE_FT| ", "Evaluation grid bounds [eV]", e_min, e_max
78 ALLOCATE (omega_complex(number_of_simulation_steps + 2))
79 ALLOCATE (moments_ft_complex(number_of_simulation_steps + 2))
80 ALLOCATE (moments_eval_complex(3, number_of_pade_points))
81 omega_complex(:) = cmplx(omega_series(:), 0.0, kind=dp)
82 DO i = 1, 3
83 moments_ft_complex(:) = cmplx(ft_full_series(2*i - 1, :), &
84 ft_full_series(2*i, :), &
85 kind=dp)
86 ! Copy the fitting parameters
87 ! TODO : Optional direct setting of parameters?
88 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 file_form="FORMATTED", file_position="REWIND")
93 IF (ft_unit > 0) THEN
94 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 REAL(x_eval(i)), REAL(moments_eval_complex(1, i)), aimag(moments_eval_complex(1, i)), &
97 REAL(moments_eval_complex(2, i)), aimag(moments_eval_complex(2, i)), &
98 REAL(moments_eval_complex(3, i)), aimag(moments_eval_complex(3, i))
99 END DO
100 END IF
101 CALL cp_print_key_finished_output(ft_unit, logger, ft_section)
102 DEALLOCATE (omega_complex)
103 DEALLOCATE (moments_ft_complex)
104 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 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 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 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 file_form="FORMATTED", file_position="REWIND")
145 ! Printing for both the stdout and separate file
146 IF (pol_unit > 0) THEN
147 IF (pol_unit == bse_unit) THEN
148 ! Print the stdout preline
149 WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|"
150 ELSE
151 ! Print also the energy in atomic units
152 WRITE (pol_unit, '(A1,A19)', advance="no") "#", "omega [a.u.]"
153 END IF
154 ! Common - print the energy in eV
155 WRITE (pol_unit, '(A20)', advance="no") "Energy [eV]"
156 ! Print a header for each polarizability element
157 DO k = 1, n_elems - 1
158 WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)', advance="no") &
159 "Real pol.", pol_elements(k, 1), pol_elements(k, 2), &
160 "Imag pol.", pol_elements(k, 1), pol_elements(k, 2)
161 END DO
162 WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)') &
163 "Real pol.", pol_elements(n_elems, 1), pol_elements(n_elems, 2), &
164 "Imag pol.", pol_elements(n_elems, 1), pol_elements(n_elems, 2)
165 DO i = 1, SIZE(x_eval)
166 IF (pol_unit == bse_unit) THEN
167 ! Print the stdout preline
168 WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|"
169 ELSE
170 ! omega in a.u.
171 WRITE (pol_unit, '(E20.8E3)', advance="no") real(x_eval(i), kind=dp)
172 END IF
173 ! Common values
174 WRITE (pol_unit, '(E20.8E3)', advance="no") real(x_eval(i), kind=dp)*evolt
175 DO k = 1, n_elems - 1
176 WRITE (pol_unit, '(E20.8E3,E20.8E3)', advance="no") &
177 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 REAL(polarizability_refined(n_elems, i)), aimag(polarizability_refined(n_elems, i))
182 END DO
183 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 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 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 CALL timeset(routinen, handle)
228
229 ! Get the sizes from arrays
230 max_fit = SIZE(x_fit)
231 n_eval = SIZE(x_eval)
232
233 ! Search for the fit start and end indices
234 fit_start = -1
235 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 IF (fit_e_max < 0) fit_end = max_fit
240 DO i = 1, max_fit
241 IF (fit_start == -1 .AND. real(x_fit(i)) >= fit_e_min) fit_start = i
242 IF (fit_end == -1 .AND. real(x_fit(i)) > fit_e_max) fit_end = i - 1
243 IF (fit_start > 0 .AND. fit_end > 0) EXIT
244 END DO
245 IF (fit_start == -1) fit_start = 1
246 IF (fit_end == -1) fit_end = max_fit
247 n_fit = fit_end - fit_start + 1
248
249 n_pade = n_fit/2
250 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 IF (n_pade < 1) THEN
255 cpwarn(é"FT deck too short for Pad; raise STEPS or disable &FT.")
256 y_eval(1:n_eval) = cmplx(0.0, 0.0, kind=dp)
257 CALL timestop(handle)
258 RETURN
259 END IF
260
261 ! Warn about a large number of Padé parameters
262 IF (n_pade > 1000) THEN
263 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 logger => cp_get_default_logger()
269 unit_nr = cp_logger_get_default_io_unit(logger)
270 IF (unit_nr > 0) THEN
271 WRITE (unit=unit_nr, fmt="(T3,A,T45,I6,I8,2F11.4)") &
272 "GREENX FT_PADE| n_pade, n_fit, window [eV]", n_pade, n_fit, &
273 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 enforce_symmetry="conjugate")
280
281 ! Check whetner the splice is needed or not
282 y_eval(1:n_eval) = evaluate_thiele_pade_at(pade_params, x_eval)
283
284 CALL free_params(pade_params)
285 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 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 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 regularization=regularization_minimax)
342 ! Factor 4 is hard-coded in the RPA weights in the internal CP2K minimax routines
343 wj(:) = wj(:)*4.0_dp
344 IF (ierr == 0) THEN
345 IF (unit_nr > 0) THEN
346 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
347 "GREENX MINIMAX_INFO| Number of integration points:", num_integ_points
348 WRITE (unit=unit_nr, fmt="(T3,A,T61,3F20.7)") &
349 "GREENX MINIMAX_INFO| Gap (Emin):", emin
350 WRITE (unit=unit_nr, fmt="(T3,A,T61,3F20.7)") &
351 "GREENX MINIMAX_INFO| Maximum eigenvalue difference (Emax):", emax
352 WRITE (unit=unit_nr, fmt="(T3,A,T61,3F20.4)") &
353 "GREENX MINIMAX_INFO| Energy range (Emax/Emin):", emax/emin
354 WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
355 "GREENX MINIMAX_INFO| Frequency grid (scaled):", "Weights", "Abscissas"
356 DO gi = 1, num_integ_points
357 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 "GREENX MINIMAX_INFO| Time grid (scaled):", "Weights", "Abscissas"
361 DO gi = 1, num_integ_points
362 WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") tau_wj(gi), tau_tj(gi)
363 END DO
364 CALL m_flush(unit_nr)
365 END IF
366 ELSE
367 IF (unit_nr > 0) THEN
368 WRITE (unit=unit_nr, fmt="(T3,A,T75)") &
369 "GREENX MINIMAX_INFO| Grid not available, use internal CP2K grids"
370 CALL m_flush(unit_nr)
371 END IF
372 IF (ALLOCATED(tau_tj)) THEN
373 DEALLOCATE (tau_tj)
374 END IF
375 IF (ALLOCATED(tau_wj)) THEN
376 DEALLOCATE (tau_wj)
377 END IF
378 IF (ALLOCATED(tj)) THEN
379 DEALLOCATE (tj)
380 END IF
381 IF (ALLOCATED(wj)) THEN
382 DEALLOCATE (wj)
383 END IF
384 IF (ALLOCATED(weights_cos_tf_t_to_w)) THEN
385 DEALLOCATE (weights_cos_tf_t_to_w)
386 END IF
387 IF (ALLOCATED(weights_cos_tf_w_to_t)) THEN
388 DEALLOCATE (weights_cos_tf_w_to_t)
389 END IF
390 IF (ALLOCATED(weights_sin_tf_t_to_w)) THEN
391 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 END SUBROUTINE greenx_get_minimax_grid
411
412END MODULE greenx_interface
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::low_print_level
integer, parameter, public low_print_level
Definition cp_output_handling.F:73
cp_output_handling::medium_print_level
integer, parameter, public medium_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition cp_output_handling.F:754
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
greenx_interface
Interface to the Greenx library.
Definition greenx_interface.F:13
greenx_interface::greenx_refine_ft
subroutine, public greenx_refine_ft(fit_e_min, fit_e_max, x_fit, y_fit, x_eval, y_eval, n_pade_opt)
Refines the FT grid using Padé approximants.
Definition greenx_interface.F:205
greenx_interface::greenx_refine_pade
subroutine, public greenx_refine_pade(e_min, e_max, x_eval, number_of_simulation_steps, number_of_pade_points, logger, ft_section, bse_unit, omega_series, ft_full_series)
Refines Pade approximants using GreenX, skips this step if GreenX is not available.
Definition greenx_interface.F:61
greenx_interface::greenx_output_polarizability
subroutine, public greenx_output_polarizability(logger, pol_section, bse_unit, pol_elements, x_eval, polarizability_refined)
Outputs the isotropic polarizability tensor element alpha _ ij = mu_i(omega)/E_j(omega),...
Definition greenx_interface.F:131
greenx_interface::greenx_get_minimax_grid
subroutine, public greenx_get_minimax_grid(unit_nr, num_integ_points, emin, emax, tau_tj, tau_wj, regularization_minimax, tj, wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, ierr)
...
Definition greenx_interface.F:319
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127