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