Branch data Line data Source code
1 : : ! **************************************************************************************************
2 : : ! Copyright (C) 2020-2023 Green-X library
3 : : ! This file is distributed under the terms of the APACHE2 License.
4 : : !
5 : : ! **************************************************************************************************
6 : : !> \brief Routines to calculate frequency and time grids (integration points and weights)
7 : : !> for correlation methods as well as weights for the inhomogeneous cosine/sine transform.
8 : : !>
9 : : !> NB: When dealing with runtime exceptions, we set ierr to a non-zero value and return immediately
10 : : ! to the caller so we don't need to goto to a cleanup section at the end of the procedure.
11 : : ! Assume -std=f2008: i.e. allocatable arrays are automatically deallocated when going out of scope.
12 : : !> reference: [https://doi.org/10.1021/ct5001268](https://doi.org/10.1021/ct5001268)
13 : : !> reference: [https://doi.org/10.1103/PhysRevB.94.165109](https://doi.org/10.1103/PhysRevB.94.165109)
14 : : ! **************************************************************************************************
15 : :
16 : : module minimax_grids
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : : #include "abi_common.h"
21 : :
22 : : #include "gx_common.h"
23 : : use defs_basis, only: dp, pi
24 : : use m_errors
25 : : !use kinds, only: dp
26 : : !use error_handling, only: register_exc
27 : : !use constants, only: pi
28 : : use minimax_tau, only: get_points_weights_tau
29 : : use minimax_omega, only: get_points_weights_omega
30 : : use minimax_utils, only: cosine_wt, cosine_tw, sine_tw
31 : : !use lapack_interfaces, only: dgemm, dgesdd
32 : :
33 : : implicit none
34 : :
35 : : private
36 : :
37 : : !> Main entry point for client code.
38 : : public :: gx_minimax_grid, gx_minimax_grid_frequency
39 : :
40 : : contains
41 : :
42 : : !> \brief Compute minimax grid for GW calculation on imaginary time/frequency domain.
43 : : !! @param[in] num_points: Number of mesh points.
44 : : !! @param[in] e_min: Minimum transition energy. Arbitrary units as we only need e_max/e_min
45 : : !! @param[in] e_max: Maximum transition energy.
46 : : !! @param[out] tau_points: imaginary time grid points
47 : : !! @param[out] tau_weights: weights for imaginary time grid points weights
48 : : !! @param[out] omega_points: imaginary frequency grid points
49 : : !! @param[out] omega_weights: weights for imaginary frequency grid points
50 : : !! @param[out] cosft_wt: weights for tau -> omega cosine transform. cos(w*t) factor is included.
51 : : !! @param[out] cosft_tw: weights for omega -> tau cosine transform. cos(w*t) factor is included.
52 : : !! @param[out] sinft_wt: weights for tau -> omega sine transform. sin(w*t) factor is included.
53 : : !! @param[out] max_errors: Max error for the three kind of transforms (same order as previous args)
54 : : !! @param[out] cosft_duality_error. Max_{ij} |AB - I| where A and B are the cosft_wt and cosft_tw matrices.
55 : : !! @param[out] ierr: Exit status
56 : : !! @param[in] bare_cos_sin_weights: whether the cosine and sine weights are multiplied by cos and sin term, optional
57 : 2 : subroutine gx_minimax_grid(num_points, e_min, e_max, &
58 : : tau_points, tau_weights, omega_points, omega_weights, &
59 : : cosft_wt, cosft_tw, sinft_wt, &
60 : : max_errors, cosft_duality_error, ierr, bare_cos_sin_weights, regterm)
61 : :
62 : : integer, intent(in) :: num_points
63 : : real(kind=dp), intent(in) :: e_min, e_max
64 : : real(kind=dp), allocatable, dimension(:), &
65 : : intent(out) :: tau_points, tau_weights
66 : : real(kind=dp), allocatable, dimension(:), &
67 : : intent(out) :: omega_points, omega_weights
68 : : real(kind=dp), allocatable, dimension(:, :), &
69 : : intent(out) :: cosft_wt, cosft_tw, sinft_wt
70 : : real(kind=dp), intent(out) :: max_errors(3), cosft_duality_error
71 : : logical, intent(in), optional :: bare_cos_sin_weights
72 : : integer, intent(out) :: ierr
73 : : real(kind=dp),optional,intent(in) :: regterm
74 : :
75 : : ! Internal variables
76 : : logical :: my_bare_cos_sin_weights
77 : : integer, parameter :: cos_t_to_cos_w = 1
78 : : integer, parameter :: cos_w_to_cos_t = 2
79 : : integer, parameter :: sin_t_to_sin_w = 3
80 : : integer :: i_point, j_point
81 : : real(kind=dp) :: e_range, scaling, regterm__
82 : 2 : real(kind=dp), dimension(:), allocatable :: x_tw
83 : 2 : real(kind=dp), dimension(:, :), allocatable :: mat
84 : 2 : real(kind=dp), dimension(:, :), allocatable :: tmp_cosft_wt, tmp_cosft_tw
85 : :
86 : 2 : my_bare_cos_sin_weights = .false.
87 [ - + ]: 2 : if (present(bare_cos_sin_weights)) then
88 : 0 : my_bare_cos_sin_weights = bare_cos_sin_weights
89 : : endif
90 : :
91 [ - + ]: 2 : regterm__ = 0.0_dp; if (present(regterm)) regterm__ = regterm
92 : :
93 : : ! Begin work
94 : 2 : e_range = e_max/e_min
95 : 2 : ierr = 0
96 : :
97 : : ! Allocations
98 [ + - - + : 6 : allocate (x_tw(2*num_points))
+ - + - +
- - + -
+ ]
99 [ + - ]: 2 : if (.not. allocated(omega_points)) then
100 [ - + + - : 6 : allocate (omega_points(num_points))
+ - + - -
+ - + -
+ ]
101 : : end if
102 [ + - ]: 2 : if (.not. allocated(omega_weights)) then
103 [ - + + - : 6 : allocate (omega_weights(num_points))
+ - + - -
+ - + -
+ ]
104 : : end if
105 [ + - ]: 2 : if (.not. allocated(tau_points)) then
106 [ - + + - : 6 : allocate (tau_points(num_points))
+ - + - -
+ - + -
+ ]
107 : : end if
108 [ + - ]: 2 : if (.not. allocated(tau_weights)) then
109 [ - + + - : 6 : allocate (tau_weights(num_points))
+ - + - -
+ - + -
+ ]
110 : : end if
111 : :
112 : : ! Get the frequency grid points and weights
113 : 2 : call get_points_weights_omega(num_points, e_range, x_tw, ierr)
114 [ + - ]: 2 : if (ierr /= 0) return
115 : :
116 : : ! Scale the frequency grid points and weights from [1,R] to [e_min,e_max]
117 : : ! Note: the frequency grid points and weights include a factor of two
118 : 2 : scaling = e_min
119 [ + + ]: 32 : omega_points(:) = x_tw(1: num_points) *scaling
120 [ + + ]: 32 : omega_weights(:) = x_tw(num_points+1: 2* num_points) *scaling
121 : :
122 : : ! Get the time grid points and weights
123 : 2 : call get_points_weights_tau(num_points, e_range, x_tw, ierr)
124 [ + - ]: 2 : if (ierr /= 0) return
125 : :
126 : : ! Scale the time grid points and weights from [1,R] to [e_min,e_max]
127 : 2 : scaling = 2.0_dp *e_min
128 [ + + ]: 32 : tau_points(:) = x_tw(1:num_points) / scaling
129 [ + + ]: 32 : tau_weights(:) = x_tw(num_points+1:2*num_points) / scaling
130 : :
131 [ - + + - : 8 : allocate (cosft_wt(num_points, num_points))
- + + - +
- + - - +
- + - + ]
132 [ - + + - : 6 : allocate (cosft_tw(num_points, num_points))
- + + - -
+ - + -
+ ]
133 [ - + + - : 6 : allocate (sinft_wt(num_points, num_points))
- + + - -
+ - + -
+ ]
134 [ - + + - : 6 : allocate (tmp_cosft_wt(num_points, num_points))
- + + - -
+ - + ]
135 [ - + + - : 6 : allocate (tmp_cosft_tw(num_points, num_points))
- + + - -
+ - + ]
136 : :
137 : : ! get the weights for the cosine transform W^c(it) -> W^c(iw)
138 : : call get_transformation_weights(num_points, tau_points, omega_points, cosft_wt, e_min, e_max, &
139 : 2 : max_errors(1), cos_t_to_cos_w, regterm__, ierr)
140 [ + - ]: 2 : if (ierr /= 0) return
141 : :
142 : : ! get the weights for the cosine transform W^c(iw) -> W^c(it)
143 : : call get_transformation_weights(num_points, tau_points, omega_points, cosft_tw, e_min, e_max, &
144 : 2 : max_errors(2), cos_w_to_cos_t, regterm__, ierr)
145 [ + - ]: 2 : if (ierr /= 0) return
146 : :
147 : : ! get the weights for the sine transform Sigma^sin(it) -> Sigma^sin(iw) (PRB 94, 165109 (2016), Eq. 71)
148 : : call get_transformation_weights(num_points, tau_points, omega_points, sinft_wt, e_min, e_max, &
149 : 2 : max_errors(3), sin_t_to_sin_w, regterm__, ierr)
150 [ + - ]: 2 : if (ierr /= 0) return
151 : :
152 : : ! Compute the actual weights used for the inhomogeneous cosine/ FT and check whether
153 : : ! the two matrices for the forward/backward transform are the inverse of each other.
154 [ + - ]: 2 : if(.not.my_bare_cos_sin_weights) then
155 [ + + ]: 32 : do j_point = 1, num_points
156 [ + + ]: 532 : do i_point = 1, num_points
157 : 500 : cosft_wt(j_point, i_point) = cosft_wt(j_point, i_point)*cos(tau_points(i_point)*omega_points(j_point))
158 : 500 : cosft_tw(i_point, j_point) = cosft_tw(i_point, j_point)*cos(tau_points(i_point)*omega_points(j_point))
159 : 530 : sinft_wt(j_point, i_point) = sinft_wt(j_point, i_point)*sin(tau_points(i_point)*omega_points(j_point))
160 : : end do
161 : : end do
162 : : else
163 [ # # ]: 0 : do j_point = 1, num_points
164 [ # # ]: 0 : do i_point = 1, num_points
165 : 0 : tmp_cosft_wt(j_point, i_point) = cosft_wt(j_point, i_point)*cos(tau_points(i_point)*omega_points(j_point))
166 : 0 : tmp_cosft_tw(i_point, j_point) = cosft_tw(i_point, j_point)*cos(tau_points(i_point)*omega_points(j_point))
167 : : end do
168 : : end do
169 : : end if
170 : :
171 [ - + + - : 6 : allocate (mat(num_points, num_points))
- + + - -
+ - + ]
172 [ + - ]: 2 : if(.not.my_bare_cos_sin_weights) then
173 [ + - + + : 37062 : mat(:,:) = matmul(cosft_wt, cosft_tw)
+ + + + +
+ + + + -
+ - + - +
- ]
174 : : else
175 [ # # # # : 0 : mat(:,:) = matmul(tmp_cosft_wt, tmp_cosft_tw)
# # # # #
# # # ]
176 : : endif
177 : :
178 [ + + ]: 32 : do i_point = 1, num_points
179 : 32 : mat(i_point, i_point) = mat(i_point, i_point) - 1.0_dp
180 : : end do
181 [ + + + + : 532 : cosft_duality_error = maxval(abs(mat))
+ + + - -
+ - - ]
182 : :
183 [ - + ]: 2 : deallocate (mat)
184 [ - + ]: 2 : deallocate (x_tw)
185 [ - + - + ]: 2 : deallocate (tmp_cosft_wt,tmp_cosft_tw)
186 : :
187 [ - - - - : 4 : end subroutine gx_minimax_grid
- - ]
188 : :
189 : : !> \brief Retrieves the frequency grid for a canonical GW/RPA calculation
190 : : !! @param[in] num_points: Number of mesh points.
191 : : !! @param[in] e_min: Minimum transition energy. Arbitrary units as we only need e_max/e_min
192 : : !! @param[in] e_max: Maximum transition energy.
193 : : !! @param[out] omega_points: imaginary frequency grid points
194 : : !! @param[out] omega_weights: weights for imaginary frequency grid points
195 : : !! @param[out] ierr: Exit status
196 : 0 : subroutine gx_minimax_grid_frequency (num_points, e_min, e_max, omega_points, omega_weights, ierr)
197 : : integer, intent(in) :: num_points
198 : : real(kind=dp), intent(in) :: e_min, e_max
199 : : real(kind=dp), allocatable, dimension(:), &
200 : : intent(out) :: omega_points(:), omega_weights(:)
201 : : integer, intent(out) :: ierr
202 : :
203 : : ! Internal variables
204 : : real(kind=dp) :: e_range, scaling
205 : 0 : real(kind=dp), dimension(:), allocatable :: x_tw
206 : :
207 : : ! Begin work
208 : 0 : e_range = e_max/e_min
209 : 0 : ierr = 0
210 : :
211 : : ! Allocations
212 [ # # # # : 0 : allocate (x_tw(2*num_points))
# # # # #
# # # #
# ]
213 [ # # ]: 0 : if (.not. allocated(omega_points)) then
214 [ # # # # : 0 : allocate (omega_points(num_points))
# # # # #
# # # #
# ]
215 : : end if
216 [ # # ]: 0 : if (.not. allocated(omega_weights)) then
217 [ # # # # : 0 : allocate (omega_weights(num_points))
# # # # #
# # # #
# ]
218 : : end if
219 : :
220 : : ! Get the frequency grid points and weights
221 : 0 : call get_points_weights_omega(num_points, e_range, x_tw, ierr)
222 [ # # ]: 0 : if (ierr /= 0) return
223 : :
224 : : ! Scale the frequency grid points and weights from [1,R] to [e_min,e_max]
225 : : ! Note: the frequency grid points and weights include a factor of two
226 : 0 : scaling = e_min
227 [ # # ]: 0 : omega_points(:) = x_tw(1: num_points) *scaling
228 [ # # ]: 0 : omega_weights(:) = x_tw(num_points+1: 2* num_points) *scaling
229 : :
230 [ # # ]: 0 : deallocate (x_tw)
231 : :
232 [ # # ]: 0 : end subroutine gx_minimax_grid_frequency
233 : :
234 : :
235 : : !> \brief Get the weights eiter for the cosine/sin transformation for tau to omega or viceversa
236 : : !! @param[in] num_points: Number of mesh points.
237 : : !! @param[in] tau_points: imaginary time grid points
238 : : !! @param[in] omega_points: imaginary frequency grid points
239 : : !! @param[inout] weights: corresponding tranformation weights
240 : : !! @param[in] e_min: Minimum transition energy.
241 : : !! @param[in] e_max: Maximum transition energy.
242 : : !! @param[inout] max_error: Max error for the transform
243 : : !! @param[in] transformation type : 1 the cosine transform cos(it) -> cos(iw)
244 : : !! : 2 the cosine transform cos(iw) -> cos(it)
245 : : !! : 3 the sine transform sin(it) -> sin(iw)
246 : : !! @param[in] ierr: exit status
247 : 6 : subroutine get_transformation_weights(num_points, tau_points, omega_points, weights, e_min, e_max, &
248 : : max_error, transformation_type, regterm, ierr)
249 : :
250 : : integer, intent(in) :: num_points
251 : : real(kind=dp), allocatable, dimension(:), &
252 : : intent(in) :: tau_points
253 : : real(kind=dp), allocatable, dimension(:), &
254 : : intent(in) :: omega_points
255 : : real(kind=dp), allocatable, dimension(:, :), &
256 : : intent(inout) :: weights
257 : : real(kind=dp), intent(in) :: e_min, e_max
258 : : real(kind=dp), intent(inout) :: max_error
259 : : integer, intent(in) :: transformation_type
260 : : integer, intent(out) :: ierr
261 : : real(kind=dp),intent(in) :: regterm
262 : :
263 : : ! Internal variables
264 : : integer :: i_node, i_point, j_point, k_point, &
265 : : num_x_nodes
266 : : integer, parameter :: nodes_factor = 200
267 : : real(kind=dp) :: current_point, x_factor
268 : : real(kind=dp), allocatable, dimension(:) :: weights_work, x_mu, psi
269 : : real(kind=dp), allocatable, dimension(:, :) :: mat_A
270 : :
271 : : integer :: lwork
272 : 6 : integer, allocatable, dimension(:) :: iwork
273 : 6 : real(kind=dp), allocatable, dimension(:) :: vec_S, vec_UT_psi, work
274 : 6 : real(kind=dp), allocatable, dimension(:, :) :: mat_U, mat_VT, mat_VT_s
275 : :
276 : : ! Begin work
277 : 6 : ierr = 0
278 : :
279 [ - + + - : 108 : allocate (weights_work(num_points), source=0.0_dp)
+ - + - -
+ - + +
+ ]
280 : :
281 : : ! compute the number of x nodes per magnitude points per 10-interval
282 : 6 : num_x_nodes = (int(log10(e_max/e_min)) + 1)*nodes_factor
283 : :
284 : : ! make sure that the number of x nodes are at least as many integration points
285 : 6 : num_x_nodes = max(num_x_nodes, num_points)
286 : :
287 [ - + + - : 3018 : allocate (x_mu(num_x_nodes), source=0.0_dp)
+ - + - -
+ - + +
+ ]
288 [ - + + - : 3012 : allocate (psi(num_x_nodes), source=0.0_dp)
- + - + +
+ ]
289 : :
290 : : ! Allocations for the BLAS routines
291 : : ! double the value nessary for 'A' to achieve good performance
292 : 6 : lwork = 8*num_points*num_points + 12*num_points + 2*num_x_nodes
293 [ + - - + : 738 : allocate (iwork(8*num_points), source=0)
+ - + - +
- - + - +
+ + ]
294 [ - + + - : 19098 : allocate (work(lwork), source=0.0_dp)
+ - + - -
+ - + +
+ ]
295 : :
296 [ - + + - : 42114 : allocate (mat_A(num_x_nodes, num_points), source=0.0_dp)
- + + - +
- + - - +
- + + + +
+ ]
297 [ - + + - : 1563024 : allocate (mat_U(num_x_nodes, num_x_nodes), source=0.0_dp)
- + + - +
- + - - +
- + + + +
+ ]
298 [ - + + - : 42108 : allocate (mat_VT(num_x_nodes, num_points), source=0.0_dp)
- + + - -
+ - + + +
+ + ]
299 [ - + + - : 45018 : allocate (mat_VT_s(num_points, num_x_nodes), source=0.0_dp)
- + + - -
+ - + + +
+ + ]
300 [ - + + - : 102 : allocate (vec_S(num_points), source=0.0_dp)
- + - + +
+ ]
301 [ - + + - : 3012 : allocate (vec_UT_psi(num_x_nodes), source=0.0_dp)
- + - + +
+ ]
302 : :
303 : : ! set the x-mu logarithmically in the interval [e_min,e_max]
304 : 6 : x_factor = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=dp) - 1.0_dp))
305 [ + + ]: 3006 : do i_node = 1, num_x_nodes
306 : 3006 : x_mu(i_node) = e_min*x_factor**(i_node - 1)
307 : : end do
308 : :
309 : 6 : current_point = 0.0_dp
310 : 6 : max_error = 0.0_dp
311 : :
312 : : ! loop over all grid points
313 [ + + ]: 96 : do i_point = 1, num_points
314 : : ! calculate psi and mat_A
315 : : call calculate_psi_and_mat_A(num_points, tau_points, omega_points, num_x_nodes, x_mu, psi, &
316 : 90 : mat_A, i_point, current_point, transformation_type)
317 : :
318 : : ! Singular value decomposition of mat_A = U*Sigma*V^T
319 : : call dgesdd('A', num_x_nodes, num_points, mat_A, num_x_nodes, vec_S, mat_U, num_x_nodes, &
320 : 90 : mat_VT, num_x_nodes, work, lwork, iwork, ierr)
321 : :
322 [ - + ]: 90 : if (ierr /= 0) then
323 : 0 : _REGISTER_EXC("DGESDD returned ierr != 0")
324 : 0 : return
325 : : end if
326 : :
327 : : ! integration weights = (V Sigma^-1 U^T)*psi
328 : :
329 : : ! 1) V * Sigma^-1
330 : : !regterm = 0.01_dp
331 [ + + ]: 1590 : do j_point = 1, num_points
332 [ + + ]: 28590 : do k_point = 1, num_points
333 [ - + ]: 28500 : if (regterm > tiny(regterm)) then
334 : : mat_VT_s(k_point, j_point) = mat_VT(j_point, k_point)*vec_S(j_point) / &
335 : 0 : (vec_S(j_point)**2+regterm**2)
336 : : else
337 : 27000 : mat_VT_s(k_point, j_point) = mat_VT(j_point, k_point) / vec_S(j_point)
338 : : end if
339 : : end do ! k_point
340 : : end do ! j_point
341 : :
342 : : ! 2) (U^T)*psi
343 : : call dgemm('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, psi, num_x_nodes, &
344 : 90 : 0.0_dp, vec_UT_psi, num_x_nodes)
345 : :
346 : : ! 3) (V*Sigma^-1) * (U^T*psi)
347 : : call dgemm('N', 'N', num_points, 1, num_x_nodes, 1.0_dp, mat_VT_s, num_points, vec_UT_psi, &
348 : 90 : num_x_nodes, 0.0_dp, weights_work, num_points)
349 : :
350 [ + + ]: 1590 : weights(i_point, :) = weights_work(:)
351 : :
352 : : ! calculate the maximum error of the fitting
353 : : call calculate_max_error(num_points, tau_points, omega_points, weights_work, num_x_nodes, x_mu, &
354 : 96 : psi, current_point, max_error, transformation_type)
355 : : end do ! i_point
356 : :
357 [ - + - + : 6 : deallocate (x_mu, psi, mat_A, weights_work, vec_S, mat_U, mat_VT, work, iwork, mat_VT_s, vec_UT_psi)
- + - + -
+ - + - +
- + - + -
+ - + ]
358 : :
359 [ - + - + : 6 : end subroutine get_transformation_weights
- + - + -
+ - + - +
- + - + -
+ - + ]
360 : :
361 : : !> \brief Calculate the auxiliary matrix for cosine/sin transformation for tau to omega or viceversa
362 : : !! @param[in] num_points: Number of mesh points.
363 : : !! @param[in] tau_points: imaginary time grid points
364 : : !! @param[in] omega_points: imaginary frequency grid points
365 : : !! @param[in] num_x_nodes: Number of node in the interval [e_min,e_max]
366 : : !! @param[in] x_mu : Transition energy (nodes in the interval [e_min,e_max])
367 : : !! @param[inout] psi: corresponding auxiliary function (see transformation type definition)
368 : : !! @param[inout] mat_A: auxiliary matrix (see transformation type definition)
369 : : !! @param[in] i_point: pointer for the current grid point
370 : : !! @param[inout] current_point: current grid point ether omega(i_point) or tau_(i_point)
371 : : !! @param[in] transformation type :
372 : : !! (1) the cosine transform cos(it) -> cos(iw): psi(omega,x), mat_A = cos(omega*tau)*psi(tau,x)
373 : : !! (2) the cosine transform cos(iw) -> cos(it): psi(tau,x) , mat_A = cos(omega*tau)*psi(omega,x)
374 : : !! (3) the sine transform sin(it) -> sin(iw): psi(omega,x), mat_A = sin(omega*tau)*psi(tau,x)
375 : 90 : subroutine calculate_psi_and_mat_A(num_points, tau_points, omega_points, num_x_nodes, x_mu, psi, &
376 : : mat_A, i_point, current_point, transformation_type)
377 : :
378 : : integer, intent(in) :: num_points, num_x_nodes, i_point
379 : : real(kind=dp), allocatable, dimension(:), &
380 : : intent(in) :: tau_points, omega_points, x_mu
381 : : real(kind=dp), allocatable, dimension(:), &
382 : : intent(inout) :: psi
383 : : real(kind=dp), allocatable, dimension(:, :), &
384 : : intent(inout) :: mat_A
385 : : real(kind=dp), intent(inout) :: current_point
386 : : integer, intent(in) :: transformation_type
387 : :
388 : : ! Internal variables
389 : : integer :: i_node, j_point
390 : : real(kind=dp) :: tau, omega
391 : :
392 : : ! Begin work
393 : :
394 : : ! the cosine transform cos(it) -> cos(iw)
395 [ + + ]: 90 : if (transformation_type == cosine_tw) then
396 : 30 : omega = omega_points(i_point)
397 : 30 : current_point = omega
398 : :
399 : : ! psi(omega_k,x) = 2x/(x^2+omega_k^2)
400 [ + + ]: 14030 : do i_node = 1, num_x_nodes
401 : 14030 : psi(i_node) = 2.0_dp*x_mu(i_node)/((x_mu(i_node))**2 + omega**2)
402 : : end do
403 : :
404 : : ! mat_A = cos(omega_k * tau) psi(tau,x)
405 [ + + ]: 530 : do j_point = 1, num_points
406 : 500 : tau = tau_points(j_point)
407 [ + + ]: 220530 : do i_node = 1, num_x_nodes
408 : 220500 : mat_A(i_node, j_point) = cos(omega*tau)*exp(-x_mu(i_node)*tau)
409 : : end do
410 : : end do
411 : :
412 : : ! the cosine transform cos(iw) -> cos(it)
413 [ + + ]: 60 : else if (transformation_type == cosine_wt) then
414 : 30 : tau = tau_points(i_point)
415 : 30 : current_point = tau
416 : :
417 : : ! psi(tau_k,x) = =exp(-x*|tau_k|)
418 [ + + ]: 14030 : do i_node = 1, num_x_nodes
419 : 14030 : psi(i_node) = exp(-x_mu(i_node)*tau)
420 : : end do
421 : :
422 : : ! mat_A = cos(tau_k,omega) psi(omega,x)
423 [ + + ]: 530 : do j_point = 1, num_points
424 : 500 : omega = omega_points(j_point)
425 [ + + ]: 220530 : do i_node = 1, num_x_nodes
426 : 220500 : mat_A(i_node, j_point) = cos(tau*omega)*2.0_dp*x_mu(i_node)/(x_mu(i_node)**2 + omega**2)
427 : : end do
428 : : end do
429 : :
430 : : ! the sine transform sin(it) -> sin(iw)
431 [ + - ]: 30 : else if (transformation_type == sine_tw) then
432 : 30 : omega = omega_points(i_point)
433 : 30 : current_point = omega
434 : :
435 : : ! psi(omega_k,x) = 2*omega_k/(x^2+omega_k^2)
436 [ + + ]: 14030 : do i_node = 1, num_x_nodes
437 : 14030 : psi(i_node) = 2.0_dp*omega/(x_mu(i_node)**2 + omega**2)
438 : : end do
439 : :
440 : : ! mat_A = sin(omega_k,tau)*phi(tau,x)
441 [ + + ]: 530 : do j_point = 1, num_points
442 : 500 : tau = tau_points(j_point)
443 [ + + ]: 220530 : do i_node = 1, num_x_nodes
444 : 220500 : mat_A(i_node, j_point) = sin(omega*tau)*exp(-x_mu(i_node)*tau)
445 : : end do
446 : : end do
447 : : end if
448 : :
449 : 90 : end subroutine calculate_psi_and_mat_A
450 : :
451 : : !> /brief calculate the error of the fit function for the cosine/sin transformation for tau to omega or viceversa
452 : : !! @param[in] num_points: Number of mesh points.
453 : : !! @param[in] tau_points: imaginary time grid points
454 : : !! @param[in] omega_points: imaginary frequency grid points
455 : : !! @param[in] weights_work: work vector for the transformation weights
456 : : !! @param[in] num_x_nodes: Number of node in the interval [e_min,e_max]
457 : : !! @param[in] x_mu : Transition energy (nodes in the interval [e_min,e_max])
458 : : !! @param[in] psi: corresponding auxiliary function (see transformation type definition)
459 : : !! @param[in] current_point: current grid point ether omega(i_point) or tau_(i_point)
460 : : !! @param[out] max_error: Max error for the transform.
461 : : !! @param[in] transformation type : 1 fit function for the cosine transform cos(it) -> cos(iw), psi(omeaga,x)
462 : : !! : 2 fit function for the cosine transform cos(iw) -> cos(it), psi(tau,x)
463 : : !! : 3 fit function for the sine transform sin(it) -> sin(iw), psi(omega,x)
464 : 90 : subroutine calculate_max_error(num_points, tau_points, omega_points, weights_work, num_x_nodes, x_mu, &
465 : : psi, current_point, max_error, transformation_type)
466 : :
467 : : real(kind=dp), intent(out) :: max_error
468 : : real(kind=dp), intent(in) :: current_point
469 : : real(kind=dp), allocatable, dimension(:), &
470 : : intent(in) :: tau_points, omega_points, x_mu, psi, &
471 : : weights_work
472 : : integer, intent(in) :: num_points, num_x_nodes
473 : : integer, intent(in) :: transformation_type
474 : :
475 : : ! Internal variables
476 : : integer :: i_node,i_point
477 : : real(kind=dp) :: func_val, func_val_temp, max_error_tmp, &
478 : : tau, omega, x_val
479 : :
480 : : ! Begin work
481 : 90 : max_error_tmp = 0.0_dp
482 : :
483 : : ! the cosine transform cos(it) -> cos(iw)
484 [ + + ]: 90 : if (transformation_type == cosine_tw) then
485 : 30 : omega=current_point
486 : :
487 [ + + ]: 14030 : do i_node = 1, num_x_nodes
488 : 14000 : func_val = 0.0_dp
489 : : ! calculate value of the fit function f(x) = f(x) + weights(omega)cos(omega*tau)psi(tau.x)
490 [ + + ]: 234000 : do i_point = 1, num_points
491 : 220000 : tau = tau_points(i_point)
492 : 234000 : func_val = func_val + weights_work(i_point)*cos(omega*tau)*exp(-x_mu(i_node)*tau)
493 : : end do
494 : :
495 : 14030 : if (abs(psi(i_node) - func_val) > max_error_tmp) then
496 : : max_error_tmp = abs(psi(i_node) - func_val)
497 : : func_val_temp = func_val
498 : : end if
499 : : end do
500 : :
501 : : ! the cosine transform cos(iw) -> cos(it)
502 [ + + ]: 60 : else if (transformation_type == cosine_wt) then
503 : 30 : tau = current_point
504 : :
505 [ + + ]: 14030 : do i_node = 1, num_x_nodes
506 : 14000 : func_val = 0.0_dp
507 : 14000 : x_val=x_mu(i_node)
508 : : ! calculate value of the fit function f(x) = f(x) + weights(tau)cos(omega*tau)psi(omega.x)
509 [ + + ]: 234000 : do i_point = 1, num_points
510 : 220000 : omega = omega_points(i_point)
511 : 234000 : func_val = func_val + weights_work(i_point)*cos(tau*omega)*2.0_dp*x_val/(x_val**2 + omega**2)
512 : : end do
513 : :
514 : 14030 : if (abs(psi(i_node) - func_val) > max_error_tmp) then
515 : : max_error_tmp = abs(psi(i_node) - func_val)
516 : : func_val_temp = func_val
517 : : end if
518 : : end do
519 : :
520 : : ! the sine transform sin(it) -> sin(iw)
521 [ + - ]: 30 : else if (transformation_type == sine_tw) then
522 : 30 : omega = current_point
523 : :
524 [ + + ]: 14030 : do i_node = 1, num_x_nodes
525 : 14000 : func_val = 0.0_dp
526 : : ! calculate value of the fit function f(x) = f(x) + weights(omega)sin(omega*tau)psi(tau.x)
527 [ + + ]: 234000 : do i_point = 1, num_points
528 : 220000 : tau = tau_points(i_point)
529 : 234000 : func_val = func_val + weights_work(i_point)*sin(omega*tau)*exp(-x_mu(i_node)*tau)
530 : : end do
531 : :
532 : 14030 : if (abs(psi(i_node) - func_val) > max_error_tmp) then
533 : : max_error_tmp = abs(psi(i_node) - func_val)
534 : : func_val_temp = func_val
535 : : end if
536 : : end do
537 : : end if
538 : :
539 [ + + ]: 90 : if (max_error_tmp > max_error) then
540 : 13 : max_error = max_error_tmp
541 : : end if
542 : :
543 : 90 : end subroutine calculate_max_error
544 : :
545 : : end module minimax_grids
|
