Line data Source code
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 : ! **************************************************************************************************
13 : MODULE 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 :
43 : CONTAINS
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 0 : SUBROUTINE greenx_refine_pade(e_min, e_max, x_eval, number_of_simulation_steps, number_of_pade_points, &
59 0 : 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 0 : COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: omega_complex, &
71 0 : moments_ft_complex
72 0 : COMPLEX(kind=dp), DIMENSION(:, :), ALLOCATABLE :: moments_eval_complex
73 :
74 : ! Report Padé refinement
75 0 : IF (bse_unit > 0) WRITE (bse_unit, '(A10,A27,E23.8E3,E20.8E3)') &
76 0 : " PADE_FT| ", "Evaluation grid bounds [eV]", e_min, e_max
77 0 : ALLOCATE (omega_complex(number_of_simulation_steps + 2))
78 0 : ALLOCATE (moments_ft_complex(number_of_simulation_steps + 2))
79 0 : ALLOCATE (moments_eval_complex(3, number_of_pade_points))
80 0 : omega_complex(:) = CMPLX(omega_series(:), 0.0, kind=dp)
81 0 : DO i = 1, 3
82 : moments_ft_complex(:) = CMPLX(ft_full_series(2*i - 1, :), &
83 : ft_full_series(2*i, :), &
84 0 : kind=dp)
85 : ! Copy the fitting parameters
86 : ! TODO : Optional direct setting of parameters?
87 0 : 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 0 : file_form="FORMATTED", file_position="REWIND")
92 0 : IF (ft_unit > 0) THEN
93 0 : 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 0 : REAL(x_eval(i)), REAL(moments_eval_complex(1, i)), AIMAG(moments_eval_complex(1, i)), &
96 0 : REAL(moments_eval_complex(2, i)), AIMAG(moments_eval_complex(2, i)), &
97 0 : REAL(moments_eval_complex(3, i)), AIMAG(moments_eval_complex(3, i))
98 : END DO
99 : END IF
100 0 : CALL cp_print_key_finished_output(ft_unit, logger, ft_section)
101 0 : DEALLOCATE (omega_complex)
102 0 : DEALLOCATE (moments_ft_complex)
103 0 : 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 0 : 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 0 : 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 0 : 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 0 : file_form="FORMATTED", file_position="REWIND")
144 : ! Printing for both the stdout and separate file
145 0 : IF (pol_unit > 0) THEN
146 0 : IF (pol_unit == bse_unit) THEN
147 : ! Print the stdout preline
148 0 : WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|"
149 : ELSE
150 : ! Print also the energy in atomic units
151 0 : WRITE (pol_unit, '(A1,A19)', advance="no") "#", "omega [a.u.]"
152 : END IF
153 : ! Common - print the energy in eV
154 0 : WRITE (pol_unit, '(A20)', advance="no") "Energy [eV]"
155 : ! Print a header for each polarizability element
156 0 : DO k = 1, n_elems - 1
157 : WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)', advance="no") &
158 0 : "Real pol.", pol_elements(k, 1), pol_elements(k, 2), &
159 0 : "Imag pol.", pol_elements(k, 1), pol_elements(k, 2)
160 : END DO
161 : WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)') &
162 0 : "Real pol.", pol_elements(n_elems, 1), pol_elements(n_elems, 2), &
163 0 : "Imag pol.", pol_elements(n_elems, 1), pol_elements(n_elems, 2)
164 0 : DO i = 1, SIZE(x_eval)
165 0 : IF (pol_unit == bse_unit) THEN
166 : ! Print the stdout preline
167 0 : WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|"
168 : ELSE
169 : ! omega in a.u.
170 0 : WRITE (pol_unit, '(E20.8E3)', advance="no") REAL(x_eval(i), kind=dp)
171 : END IF
172 : ! Common values
173 0 : WRITE (pol_unit, '(E20.8E3)', advance="no") REAL(x_eval(i), kind=dp)*evolt
174 0 : DO k = 1, n_elems - 1
175 : WRITE (pol_unit, '(E20.8E3,E20.8E3)', advance="no") &
176 0 : 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 0 : REAL(polarizability_refined(n_elems, i)), AIMAG(polarizability_refined(n_elems, i))
181 : END DO
182 0 : 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 0 : 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 6 : 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 6 : max_fit = SIZE(x_fit)
223 6 : n_eval = SIZE(x_eval)
224 :
225 : ! Search for the fit start and end indices
226 6 : fit_start = -1
227 6 : 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 6 : IF (fit_e_max < 0) fit_end = max_fit
232 2232 : DO i = 1, max_fit
233 2232 : IF (fit_start == -1 .AND. REAL(x_fit(i)) >= fit_e_min) fit_start = i
234 2232 : IF (fit_end == -1 .AND. REAL(x_fit(i)) > fit_e_max) fit_end = i - 1
235 2232 : IF (fit_start > 0 .AND. fit_end > 0) EXIT
236 : END DO
237 6 : IF (fit_start == -1) fit_start = 1
238 6 : IF (fit_end == -1) fit_end = max_fit
239 6 : n_fit = fit_end - fit_start + 1
240 :
241 6 : n_pade = n_fit/2
242 6 : IF (PRESENT(n_pade_opt)) n_pade = n_pade_opt
243 :
244 : ! Warn about a large number of Padé parameters
245 6 : IF (n_pade > 1000) THEN
246 0 : 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 6 : enforce_symmetry="conjugate")
252 :
253 : ! Check whetner the splice is needed or not
254 6000 : y_eval(1:n_eval) = evaluate_thiele_pade_at(pade_params, x_eval)
255 :
256 6 : 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 6 : END SUBROUTINE greenx_refine_ft
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 204 : 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 204 : regularization=regularization_minimax)
313 : ! Factor 4 is hard-coded in the RPA weights in the internal CP2K minimax routines
314 1484 : wj(:) = wj(:)*4.0_dp
315 204 : IF (ierr == 0) THEN
316 74 : IF (unit_nr > 0) THEN
317 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
318 37 : "GREENX MINIMAX_INFO| Number of integration points:", num_integ_points
319 : WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") &
320 37 : "GREENX MINIMAX_INFO| Gap (Emin):", Emin
321 : WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") &
322 37 : "GREENX MINIMAX_INFO| Maximum eigenvalue difference (Emax):", Emax
323 : WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.4)") &
324 37 : "GREENX MINIMAX_INFO| Energy range (Emax/Emin):", Emax/Emin
325 : WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") &
326 37 : "GREENX MINIMAX_INFO| Frequency grid (scaled):", "Weights", "Abscissas"
327 425 : DO gi = 1, num_integ_points
328 425 : 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 37 : "GREENX MINIMAX_INFO| Time grid (scaled):", "Weights", "Abscissas"
332 425 : DO gi = 1, num_integ_points
333 425 : WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(gi), tau_tj(gi)
334 : END DO
335 37 : CALL m_flush(unit_nr)
336 : END IF
337 : ELSE
338 130 : IF (unit_nr > 0) THEN
339 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75)") &
340 65 : "GREENX MINIMAX_INFO| Grid not available, use internal CP2K grids"
341 65 : CALL m_flush(unit_nr)
342 : END IF
343 130 : IF (ALLOCATED(tau_tj)) THEN
344 130 : DEALLOCATE (tau_tj)
345 : END IF
346 130 : IF (ALLOCATED(tau_wj)) THEN
347 130 : DEALLOCATE (tau_wj)
348 : END IF
349 130 : IF (ALLOCATED(tj)) THEN
350 130 : DEALLOCATE (tj)
351 : END IF
352 130 : IF (ALLOCATED(wj)) THEN
353 130 : DEALLOCATE (wj)
354 : END IF
355 130 : IF (ALLOCATED(weights_cos_tf_t_to_w)) THEN
356 0 : DEALLOCATE (weights_cos_tf_t_to_w)
357 : END IF
358 130 : IF (ALLOCATED(weights_cos_tf_w_to_t)) THEN
359 0 : DEALLOCATE (weights_cos_tf_w_to_t)
360 : END IF
361 130 : IF (ALLOCATED(weights_sin_tf_t_to_w)) THEN
362 0 : 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 204 : END SUBROUTINE greenx_get_minimax_grid
382 :
383 : END MODULE greenx_interface
|
