(git:d2a9ebd)
Loading...
Searching...
No Matches
negf_integr_simpson.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 Adaptive Simpson's rule algorithm to integrate a complex-valued function in a complex plane
10! **************************************************************************************************
11MODULE negf_integr_simpson
12 USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale,&
13 cp_cfm_scale_and_add
14 USE cp_cfm_types, ONLY: cp_cfm_create,&
15 cp_cfm_get_info,&
16 cp_cfm_release,&
17 cp_cfm_set_all,&
18 cp_cfm_to_cfm,&
19 cp_cfm_type
20 USE cp_fm_basic_linalg, ONLY: cp_fm_trace
21 USE cp_fm_struct, ONLY: cp_fm_struct_type
22 USE cp_fm_types, ONLY: cp_fm_create,&
23 cp_fm_get_info,&
24 cp_fm_release,&
25 cp_fm_type
26 USE kinds, ONLY: dp
27 USE mathconstants, ONLY: pi,&
28 z_one,&
29 z_zero
30 USE negf_integr_utils, ONLY: contour_shape_arc,&
31 contour_shape_linear,&
32 equidistant_nodes_a_b,&
33 rescale_normalised_nodes
34 USE util, ONLY: sort
35#include "./base/base_uses.f90"
36
37 IMPLICIT NONE
38 PRIVATE
39
40 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_integr_simpson'
41 ! adaptive Simpson method requires 5 points per subinterval for the error estimate.
42 ! So, in principle, at the end we can compute the value of the integral using
43 ! Boole's rule and possibly improve the actual accuracy by up to one order of magnitude.
44 LOGICAL, PARAMETER, PRIVATE :: is_boole = .false.
45
46 INTEGER, PARAMETER, PUBLIC :: sr_shape_linear = contour_shape_linear, &
47 sr_shape_arc = contour_shape_arc
48
49 PUBLIC :: simpsonrule_type
50 PUBLIC :: simpsonrule_init, simpsonrule_release, simpsonrule_get_next_nodes, simpsonrule_refine_integral
51
52! **************************************************************************************************
53!> \brief A structure to store data for non-converged sub-interval.
54! **************************************************************************************************
55 TYPE simpsonrule_subinterval_type
56 !> unscaled lower and upper boundaries within the interval [-1 .. 1]
57 REAL(kind=dp) :: lb = -1.0_dp, ub = -1.0_dp
58 !> target accuracy for this sub-interval
59 REAL(kind=dp) :: conv = -1.0_dp
60 !> estimated error value on this sub-interval
61 REAL(kind=dp) :: error = -1.0_dp
62 !> integrand values at equally spaced points [a, b, c, d, e] located on the curve shape([lb .. ub])
63 TYPE(cp_cfm_type) :: fa = cp_cfm_type(), fb = cp_cfm_type(), fc = cp_cfm_type(), &
64 fd = cp_cfm_type(), fe = cp_cfm_type()
65 END TYPE simpsonrule_subinterval_type
66
67! **************************************************************************************************
68!> \brief A structure to store data needed for adaptive Simpson's rule algorithm.
69! **************************************************************************************************
70 TYPE simpsonrule_type
71 !> lower and upper boundaries of the curve on the complex plane
72 COMPLEX(kind=dp) :: a = z_zero, b = z_zero
73 !> ID number which determines the shape of a curve along which the integral will be evaluated
74 INTEGER :: shape_id = -1
75 !> target accuracy
76 REAL(kind=dp) :: conv = -1.0_dp
77 !> estimated error value on the entire integration interval,
78 !> as well as on converged sub-intervals only
79 REAL(kind=dp) :: error = -1.0_dp, error_conv = -1.0_dp
80 !> the estimated value of the integral on the entire interval
81 TYPE(cp_cfm_type), POINTER :: integral => null()
82 !> work matrix to store the contribution to the integral on converged sub-intervals
83 TYPE(cp_cfm_type), POINTER :: integral_conv => null()
84 !> work matrices which stores approximated integral computed by using a/b/c, c/d/e, and a/c/e points respectively
85 TYPE(cp_cfm_type), POINTER :: integral_abc => null(), integral_cde => null(), integral_ace => null()
86 !> work matrix to temporarily store error estimate of the integral on a sub-interval for every matrix element
87 TYPE(cp_fm_type), POINTER :: error_fm => null()
88 !> weights associated with matrix elements; the final error is computed as Trace(error_fm * weights)
89 TYPE(cp_fm_type), POINTER :: weights => null()
90 ! non-converged sub-intervals
91 TYPE(simpsonrule_subinterval_type), ALLOCATABLE, &
92 DIMENSION(:) :: subintervals
93 !> complete list of nodes over the normalised interval [-1 .. 1] needed to restart
94 !> Useful when a series of similar integrals need to be computed at an identical set
95 !> of points, so intermediate quantities can be saved and reused.
96 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes
97 END TYPE simpsonrule_type
98
99 COMPLEX(kind=dp), PARAMETER, PRIVATE :: z_four = 4.0_dp*z_one
100
101CONTAINS
102
103! **************************************************************************************************
104!> \brief Initialise a Simpson's rule environment variable.
105!> \param sr_env Simpson's rule environment (initialised on exit)
106!> \param xnodes points at which an integrand needs to be computed (initialised on exit)
107!> \param nnodes initial number of points to compute (initialised on exit)
108!> \param a integral lower boundary
109!> \param b integral upper boundary
110!> \param shape_id shape of a curve along which the integral will be evaluated
111!> \param conv convergence threshold
112!> \param weights weights associated with matrix elements; used to compute cumulative error
113!> \param tnodes_restart list of nodes over the interval [-1 .. 1] from a previous integral evaluation.
114!> If present, the same set of 'xnodes' will be used to compute this integral.
115!> \par History
116!> * 05.2017 created [Sergey Chulkov]
117!> \note When we integrate the retarded Green's function times the Fermi function over the energy
118!> domain and pass the overlap matrix (S) as the 'weights' matrix, the convergence threshold
119!> ('conv') becomes the maximum error in the total number of electrons multiplied by pi.
120! **************************************************************************************************
121 SUBROUTINE simpsonrule_init(sr_env, xnodes, nnodes, a, b, shape_id, conv, weights, tnodes_restart)
122 TYPE(simpsonrule_type), INTENT(out) :: sr_env
123 INTEGER, INTENT(inout) :: nnodes
124 COMPLEX(kind=dp), DIMENSION(nnodes), INTENT(out) :: xnodes
125 COMPLEX(kind=dp), INTENT(in) :: a, b
126 INTEGER, INTENT(in) :: shape_id
127 REAL(kind=dp), INTENT(in) :: conv
128 TYPE(cp_fm_type), INTENT(IN) :: weights
129 REAL(kind=dp), DIMENSION(nnodes), INTENT(in), &
130 OPTIONAL :: tnodes_restart
131
132 CHARACTER(len=*), PARAMETER :: routinen = 'simpsonrule_init'
133
134 INTEGER :: handle, icol, irow, ncols, nrows
135 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
136 POINTER :: w_data, w_data_my
137 TYPE(cp_fm_struct_type), POINTER :: fm_struct
138
139 CALL timeset(routinen, handle)
140
141 cpassert(nnodes > 4)
142
143 ! ensure that MOD(nnodes-1, 4) == 0
144 nnodes = 4*((nnodes - 1)/4) + 1
145
146 sr_env%shape_id = shape_id
147 sr_env%a = a
148 sr_env%b = b
149 sr_env%conv = conv
150 sr_env%error = huge(0.0_dp)
151 sr_env%error_conv = 0.0_dp
152
153 NULLIFY (sr_env%error_fm, sr_env%weights)
154 CALL cp_fm_get_info(weights, local_data=w_data, nrow_local=nrows, ncol_local=ncols, matrix_struct=fm_struct)
155 ALLOCATE (sr_env%error_fm, sr_env%weights)
156 CALL cp_fm_create(sr_env%error_fm, fm_struct)
157 CALL cp_fm_create(sr_env%weights, fm_struct)
158 CALL cp_fm_get_info(sr_env%weights, local_data=w_data_my)
159
160 ! use the explicit loop to avoid temporary arrays. The magic constant 15.0 is due to Simpson's rule error analysis.
161 DO icol = 1, ncols
162 DO irow = 1, nrows
163 w_data_my(irow, icol) = abs(w_data(irow, icol))/15.0_dp
164 END DO
165 END DO
166
167 NULLIFY (sr_env%integral, sr_env%integral_conv)
168 NULLIFY (sr_env%integral_abc, sr_env%integral_cde, sr_env%integral_ace)
169
170 ALLOCATE (sr_env%tnodes(nnodes))
171
172 IF (PRESENT(tnodes_restart)) THEN
173 sr_env%tnodes(1:nnodes) = tnodes_restart(1:nnodes)
174 ELSE
175 CALL equidistant_nodes_a_b(-1.0_dp, 1.0_dp, nnodes, sr_env%tnodes)
176 END IF
177 CALL rescale_normalised_nodes(nnodes, sr_env%tnodes, a, b, shape_id, xnodes)
178
179 CALL timestop(handle)
180 END SUBROUTINE simpsonrule_init
181
182! **************************************************************************************************
183!> \brief Release a Simpson's rule environment variable.
184!> \param sr_env Simpson's rule environment (modified on exit)
185!> \par History
186!> * 05.2017 created [Sergey Chulkov]
187! **************************************************************************************************
188 SUBROUTINE simpsonrule_release(sr_env)
189 TYPE(simpsonrule_type), INTENT(inout) :: sr_env
190
191 CHARACTER(len=*), PARAMETER :: routinen = 'simpsonrule_release'
192
193 INTEGER :: handle, interval
194
195 CALL timeset(routinen, handle)
196 IF (ALLOCATED(sr_env%subintervals)) THEN
197 DO interval = SIZE(sr_env%subintervals), 1, -1
198 CALL cp_cfm_release(sr_env%subintervals(interval)%fa)
199 CALL cp_cfm_release(sr_env%subintervals(interval)%fb)
200 CALL cp_cfm_release(sr_env%subintervals(interval)%fc)
201 CALL cp_cfm_release(sr_env%subintervals(interval)%fd)
202 CALL cp_cfm_release(sr_env%subintervals(interval)%fe)
203 END DO
204
205 DEALLOCATE (sr_env%subintervals)
206 END IF
207
208 IF (ASSOCIATED(sr_env%integral)) THEN
209 CALL cp_cfm_release(sr_env%integral)
210 DEALLOCATE (sr_env%integral)
211 NULLIFY (sr_env%integral)
212 END IF
213 IF (ASSOCIATED(sr_env%integral_conv)) THEN
214 CALL cp_cfm_release(sr_env%integral_conv)
215 DEALLOCATE (sr_env%integral_conv)
216 NULLIFY (sr_env%integral_conv)
217 END IF
218 IF (ASSOCIATED(sr_env%integral_abc)) THEN
219 CALL cp_cfm_release(sr_env%integral_abc)
220 DEALLOCATE (sr_env%integral_abc)
221 NULLIFY (sr_env%integral_abc)
222 END IF
223 IF (ASSOCIATED(sr_env%integral_cde)) THEN
224 CALL cp_cfm_release(sr_env%integral_cde)
225 DEALLOCATE (sr_env%integral_cde)
226 NULLIFY (sr_env%integral_cde)
227 END IF
228 IF (ASSOCIATED(sr_env%integral_ace)) THEN
229 CALL cp_cfm_release(sr_env%integral_ace)
230 DEALLOCATE (sr_env%integral_ace)
231 NULLIFY (sr_env%integral_ace)
232 END IF
233 IF (ASSOCIATED(sr_env%error_fm)) THEN
234 CALL cp_fm_release(sr_env%error_fm)
235 DEALLOCATE (sr_env%error_fm)
236 NULLIFY (sr_env%error_fm)
237 END IF
238 IF (ASSOCIATED(sr_env%weights)) THEN
239 CALL cp_fm_release(sr_env%weights)
240 DEALLOCATE (sr_env%weights)
241 NULLIFY (sr_env%weights)
242 END IF
243
244 IF (ALLOCATED(sr_env%tnodes)) DEALLOCATE (sr_env%tnodes)
245
246 CALL timestop(handle)
247 END SUBROUTINE simpsonrule_release
248
249! **************************************************************************************************
250!> \brief Get the next set of nodes where to compute integrand.
251!> \param sr_env Simpson's rule environment (modified on exit)
252!> \param xnodes_next list of additional points (initialised on exit)
253!> \param nnodes actual number of points to compute (modified on exit)
254!> \par History
255!> * 05.2017 created [Sergey Chulkov]
256!> \note The number of nodes returned is limited by the initial value of the nnodes variable;
257!> un exit nnodes == 0 means that the target accuracy has been achieved.
258! **************************************************************************************************
259 SUBROUTINE simpsonrule_get_next_nodes(sr_env, xnodes_next, nnodes)
260 TYPE(simpsonrule_type), INTENT(inout) :: sr_env
261 INTEGER, INTENT(inout) :: nnodes
262 COMPLEX(kind=dp), DIMENSION(nnodes), INTENT(out) :: xnodes_next
263
264 CHARACTER(len=*), PARAMETER :: routinen = 'simpsonrule_get_next_nodes'
265
266 INTEGER :: handle, nnodes_old
267 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes, tnodes_old
268
269 CALL timeset(routinen, handle)
270 ALLOCATE (tnodes(nnodes))
271
272 CALL simpsonrule_get_next_nodes_real(sr_env, tnodes, nnodes)
273 IF (nnodes > 0) THEN
274 CALL move_alloc(sr_env%tnodes, tnodes_old)
275 nnodes_old = SIZE(tnodes_old)
276
277 ALLOCATE (sr_env%tnodes(nnodes_old + nnodes))
278 sr_env%tnodes(1:nnodes_old) = tnodes_old(1:nnodes_old)
279 sr_env%tnodes(nnodes_old + 1:nnodes_old + nnodes) = tnodes(1:nnodes)
280 DEALLOCATE (tnodes_old)
281
282 CALL rescale_normalised_nodes(nnodes, tnodes, sr_env%a, sr_env%b, sr_env%shape_id, xnodes_next)
283 END IF
284
285 DEALLOCATE (tnodes)
286 CALL timestop(handle)
287 END SUBROUTINE simpsonrule_get_next_nodes
288
289! **************************************************************************************************
290!> \brief Low level routine that returns unscaled nodes on interval [-1 .. 1].
291!> \param sr_env Simpson's rule environment
292!> \param xnodes_unity list of additional unscaled nodes (initialised on exit)
293!> \param nnodes actual number of points to compute (initialised on exit)
294!> \par History
295!> * 05.2017 created [Sergey Chulkov]
296! **************************************************************************************************
297 SUBROUTINE simpsonrule_get_next_nodes_real(sr_env, xnodes_unity, nnodes)
298 TYPE(simpsonrule_type), INTENT(in) :: sr_env
299 REAL(kind=dp), DIMENSION(:), INTENT(out) :: xnodes_unity
300 INTEGER, INTENT(out) :: nnodes
301
302 CHARACTER(len=*), PARAMETER :: routinen = 'simpsonrule_get_next_nodes_real'
303
304 INTEGER :: handle, interval, nintervals
305
306 CALL timeset(routinen, handle)
307
308 IF (ALLOCATED(sr_env%subintervals)) THEN
309 nintervals = SIZE(sr_env%subintervals)
310 ELSE
311 nintervals = 0
312 END IF
313
314 IF (nintervals > 0) THEN
315 IF (SIZE(xnodes_unity) < 4*nintervals) THEN
316 nintervals = SIZE(xnodes_unity)/4
317 END IF
318
319 DO interval = 1, nintervals
320 xnodes_unity(4*interval - 3) = 0.125_dp* &
321 (7.0_dp*sr_env%subintervals(interval)%lb + sr_env%subintervals(interval)%ub)
322 xnodes_unity(4*interval - 2) = 0.125_dp* &
323 (5.0_dp*sr_env%subintervals(interval)%lb + 3.0_dp*sr_env%subintervals(interval)%ub)
324 xnodes_unity(4*interval - 1) = 0.125_dp* &
325 (3.0_dp*sr_env%subintervals(interval)%lb + 5.0_dp*sr_env%subintervals(interval)%ub)
326 xnodes_unity(4*interval) = 0.125_dp*(sr_env%subintervals(interval)%lb + 7.0_dp*sr_env%subintervals(interval)%ub)
327 END DO
328 END IF
329
330 nnodes = 4*nintervals
331 CALL timestop(handle)
332 END SUBROUTINE simpsonrule_get_next_nodes_real
333
334! **************************************************************************************************
335!> \brief Compute integral using the simpson's rules.
336!> \param sr_env Simpson's rule environment
337!> \param zdata_next precomputed integrand values at points xnodes_next (nullified on exit)
338!> \par History
339!> * 05.2017 created [Sergey Chulkov]
340! **************************************************************************************************
341 SUBROUTINE simpsonrule_refine_integral(sr_env, zdata_next)
342 TYPE(simpsonrule_type), INTENT(inout) :: sr_env
343 TYPE(cp_cfm_type), DIMENSION(:), INTENT(inout) :: zdata_next
344
345 CHARACTER(len=*), PARAMETER :: routinen = 'simpsonrule_refine_integral'
346 TYPE(cp_cfm_type), PARAMETER :: cfm_null = cp_cfm_type()
347
348 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: zscale
349 COMPLEX(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
350 POINTER :: error_zdata
351 INTEGER :: handle, interval, ipoint, jpoint, &
352 nintervals, nintervals_exist, npoints
353 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
354 LOGICAL :: interval_converged, interval_exists
355 REAL(kind=dp) :: my_bound, rscale
356 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: errors
357 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
358 POINTER :: error_rdata
359 TYPE(cp_fm_struct_type), POINTER :: fm_struct
360 TYPE(simpsonrule_subinterval_type), ALLOCATABLE, &
361 DIMENSION(:) :: subintervals
362
363 CALL timeset(routinen, handle)
364
365 npoints = SIZE(zdata_next)
366 IF (ASSOCIATED(sr_env%integral)) THEN
367 ! we need 4 new points per subinterval (p, q, r, s)
368 ! p q r s
369 ! a . b . c . d . e
370 cpassert(npoints > 0 .AND. mod(npoints, 4) == 0)
371 ELSE
372 ! first call: need 4*n+1 points
373 ! a1 b1 c1 d1 e1
374 ! a2 b2 c2 d2 e2
375 ! a3 b3 c3 d3 e3
376 cpassert(npoints > 1 .AND. mod(npoints, 4) == 1)
377 END IF
378
379 ! compute weights of new points on a complex contour according to their values of the 't' parameter
380 nintervals_exist = SIZE(sr_env%tnodes)
381 cpassert(nintervals_exist >= npoints)
382 ALLOCATE (zscale(npoints))
383
384 CALL rescale_normalised_nodes(npoints, sr_env%tnodes(nintervals_exist - npoints + 1:nintervals_exist), &
385 sr_env%a, sr_env%b, sr_env%shape_id, weights=zscale)
386
387 ! rescale integrand values
388 DO ipoint = 1, npoints
389 CALL cp_cfm_scale(zscale(ipoint), zdata_next(ipoint))
390 END DO
391
392 DEALLOCATE (zscale)
393
394 ! insert new points
395 nintervals = npoints/4
396 IF (ASSOCIATED(sr_env%integral)) THEN
397 ! subdivide existing intervals
398 nintervals_exist = SIZE(sr_env%subintervals)
399 cpassert(nintervals <= nintervals_exist)
400
401 ALLOCATE (subintervals(nintervals_exist + nintervals))
402
403 DO interval = 1, nintervals
404 subintervals(2*interval - 1)%lb = sr_env%subintervals(interval)%lb
405 subintervals(2*interval - 1)%ub = 0.5_dp*(sr_env%subintervals(interval)%lb + sr_env%subintervals(interval)%ub)
406 subintervals(2*interval - 1)%conv = 0.5_dp*sr_env%subintervals(interval)%conv
407 subintervals(2*interval - 1)%fa = sr_env%subintervals(interval)%fa
408 subintervals(2*interval - 1)%fb = zdata_next(4*interval - 3)
409 subintervals(2*interval - 1)%fc = sr_env%subintervals(interval)%fb
410 subintervals(2*interval - 1)%fd = zdata_next(4*interval - 2)
411 subintervals(2*interval - 1)%fe = sr_env%subintervals(interval)%fc
412
413 subintervals(2*interval)%lb = subintervals(2*interval - 1)%ub
414 subintervals(2*interval)%ub = sr_env%subintervals(interval)%ub
415 subintervals(2*interval)%conv = subintervals(2*interval - 1)%conv
416 subintervals(2*interval)%fa = sr_env%subintervals(interval)%fc
417 subintervals(2*interval)%fb = zdata_next(4*interval - 1)
418 subintervals(2*interval)%fc = sr_env%subintervals(interval)%fd
419 subintervals(2*interval)%fd = zdata_next(4*interval)
420 subintervals(2*interval)%fe = sr_env%subintervals(interval)%fe
421
422 zdata_next(4*interval - 3:4*interval) = cfm_null
423 END DO
424
425 DO interval = nintervals + 1, nintervals_exist
426 subintervals(interval + nintervals) = sr_env%subintervals(interval)
427 END DO
428 DEALLOCATE (sr_env%subintervals)
429 ELSE
430 ! first time -- allocate matrices and create a new set of intervals
431 CALL cp_cfm_get_info(zdata_next(1), matrix_struct=fm_struct)
432 ALLOCATE (sr_env%integral, sr_env%integral_conv, &
433 sr_env%integral_abc, sr_env%integral_cde, sr_env%integral_ace)
434 CALL cp_cfm_create(sr_env%integral, fm_struct)
435 CALL cp_cfm_create(sr_env%integral_conv, fm_struct)
436 CALL cp_cfm_create(sr_env%integral_abc, fm_struct)
437 CALL cp_cfm_create(sr_env%integral_cde, fm_struct)
438 CALL cp_cfm_create(sr_env%integral_ace, fm_struct)
439
440 CALL cp_cfm_set_all(sr_env%integral_conv, z_zero)
441
442 ALLOCATE (subintervals(nintervals))
443
444 rscale = 1.0_dp/real(nintervals, kind=dp)
445
446 DO interval = 1, nintervals
447 ! lower bound: point with indices 1, 5, 9, ..., 4*nintervals+1
448 subintervals(interval)%lb = sr_env%tnodes(4*interval - 3)
449 subintervals(interval)%ub = sr_env%tnodes(4*interval + 1)
450 subintervals(interval)%conv = rscale*sr_env%conv
451
452 subintervals(interval)%fa = zdata_next(4*interval - 3)
453 subintervals(interval)%fb = zdata_next(4*interval - 2)
454 subintervals(interval)%fc = zdata_next(4*interval - 1)
455 subintervals(interval)%fd = zdata_next(4*interval)
456 subintervals(interval)%fe = zdata_next(4*interval + 1)
457 END DO
458 END IF
459
460 ! we kept the originals matrices for internal use, so set the matrix to null
461 ! to prevent alteration of the matrices from the outside
462 zdata_next(1:npoints) = cfm_null
463
464 CALL cp_fm_get_info(sr_env%error_fm, local_data=error_rdata)
465 CALL cp_cfm_get_info(sr_env%integral_ace, local_data=error_zdata)
466
467 ! do actual integration
468 CALL cp_cfm_to_cfm(sr_env%integral_conv, sr_env%integral)
469 sr_env%error = sr_env%error_conv
470 nintervals_exist = SIZE(subintervals)
471
472 DO interval = 1, nintervals_exist
473 rscale = subintervals(interval)%ub - subintervals(interval)%lb
474 CALL do_simpson_rule(sr_env%integral_ace, &
475 subintervals(interval)%fa, &
476 subintervals(interval)%fc, &
477 subintervals(interval)%fe, &
478 -0.5_dp*rscale)
479 CALL do_simpson_rule(sr_env%integral_abc, &
480 subintervals(interval)%fa, &
481 subintervals(interval)%fb, &
482 subintervals(interval)%fc, &
483 0.25_dp*rscale)
484 CALL do_simpson_rule(sr_env%integral_cde, &
485 subintervals(interval)%fc, &
486 subintervals(interval)%fd, &
487 subintervals(interval)%fe, &
488 0.25_dp*rscale)
489
490 CALL cp_cfm_scale_and_add(z_one, sr_env%integral_abc, z_one, sr_env%integral_cde)
491 CALL cp_cfm_scale_and_add(z_one, sr_env%integral_ace, z_one, sr_env%integral_abc)
492
493 IF (is_boole) THEN
494 CALL do_boole_rule(sr_env%integral_abc, &
495 subintervals(interval)%fa, &
496 subintervals(interval)%fb, &
497 subintervals(interval)%fc, &
498 subintervals(interval)%fd, &
499 subintervals(interval)%fe, &
500 0.5_dp*rscale, sr_env%integral_cde)
501 END IF
502
503 CALL cp_cfm_scale_and_add(z_one, sr_env%integral, z_one, sr_env%integral_abc)
504
505 ! sr_env%error_fm = ABS(sr_env%integral_ace); no temporary arrays as pointers have different types
506 error_rdata(:, :) = abs(error_zdata(:, :))
507 CALL cp_fm_trace(sr_env%error_fm, sr_env%weights, subintervals(interval)%error)
508
509 sr_env%error = sr_env%error + subintervals(interval)%error
510
511 ! add contributions from converged subintervals, so we could drop them afterward
512 IF (subintervals(interval)%error <= subintervals(interval)%conv) THEN
513 CALL cp_cfm_scale_and_add(z_one, sr_env%integral_conv, z_one, sr_env%integral_abc)
514 sr_env%error_conv = sr_env%error_conv + subintervals(interval)%error
515 END IF
516 END DO
517
518 IF (sr_env%error <= sr_env%conv) THEN
519 ! We have already reached the target accuracy, so we can drop all subintervals
520 ! (even those where local convergence has not been achieved). From now on environment
521 ! components 'sr_env%error' and 'sr_env%integral_conv' hold incorrect values,
522 ! but they should not been accessed from the outside anyway
523 ! (uncomment the following two lines if they are actually need)
524
525 ! sr_env%error_conv = sr_env%error
526 ! CALL cp_cfm_to_cfm(sr_env%integral, sr_env%integral_conv)
527
528 ! Only deallocate the fa component explicitly if there is no interval to the left from it
529 DO interval = nintervals_exist, 1, -1
530 interval_exists = .false.
531 my_bound = subintervals(interval)%lb
532 DO jpoint = 1, nintervals_exist
533 IF (subintervals(jpoint)%ub == my_bound) THEN
534 interval_exists = .true.
535 EXIT
536 END IF
537 END DO
538 IF (.NOT. interval_exists) THEN
539 ! interval does not exist anymore, so it is safe to release the matrix
540 CALL cp_cfm_release(subintervals(interval)%fa)
541 ELSE IF (interval_converged) THEN
542 ! the interval still exists and will be released with fe
543 END IF
544 CALL cp_cfm_release(subintervals(interval)%fb)
545 CALL cp_cfm_release(subintervals(interval)%fc)
546 CALL cp_cfm_release(subintervals(interval)%fd)
547 CALL cp_cfm_release(subintervals(interval)%fe)
548 END DO
549 ELSE
550 ! sort subinterval according to their convergence, and drop convergent ones
551 ALLOCATE (errors(nintervals_exist), inds(nintervals_exist))
552
553 nintervals = 0
554 DO interval = 1, nintervals_exist
555 errors(interval) = subintervals(interval)%error
556
557 IF (subintervals(interval)%error > subintervals(interval)%conv) THEN
558 nintervals = nintervals + 1
559 END IF
560 END DO
561
562 CALL sort(errors, nintervals_exist, inds)
563
564 IF (nintervals > 0) THEN
565 ALLOCATE (sr_env%subintervals(nintervals))
566 END IF
567
568 nintervals = 0
569 DO ipoint = nintervals_exist, 1, -1
570 interval = inds(ipoint)
571
572 IF (subintervals(interval)%error > subintervals(interval)%conv) THEN
573 nintervals = nintervals + 1
574
575 sr_env%subintervals(nintervals) = subintervals(interval)
576 ELSE
577 ! Release matrices of converged intervals. Special cases: left and right boundary
578 ! Check whether the neighboring interval still exists and if it does, check for its convergence
579 interval_exists = .false.
580 my_bound = subintervals(interval)%lb
581 DO jpoint = 1, nintervals_exist
582 IF (subintervals(jpoint)%ub == my_bound) THEN
583 interval_exists = .true.
584 EXIT
585 END IF
586 END DO
587 IF (.NOT. interval_exists) THEN
588 ! interval does not exist anymore, so it is safe to release the matrix
589 CALL cp_cfm_release(subintervals(interval)%fa)
590 ELSE IF (interval_converged) THEN
591 ! the interval still exists and will be released with fe
592 END IF
593 CALL cp_cfm_release(subintervals(interval)%fb)
594 CALL cp_cfm_release(subintervals(interval)%fc)
595 CALL cp_cfm_release(subintervals(interval)%fd)
596
597 ! Right border: Check for the existence and the convergence of the interval
598 ! If the right interval does not exist or has converged, release the matrix
599 interval_exists = .false.
600 interval_converged = .false.
601 my_bound = subintervals(interval)%ub
602 DO jpoint = 1, nintervals_exist
603 IF (subintervals(jpoint)%lb == my_bound) THEN
604 interval_exists = .true.
605 IF (subintervals(jpoint)%error <= subintervals(jpoint)%conv) interval_converged = .true.
606 EXIT
607 END IF
608 END DO
609 IF (.NOT. interval_exists .OR. interval_converged) THEN
610 CALL cp_cfm_release(subintervals(interval)%fe)
611 END IF
612 END IF
613 END DO
614
615 DEALLOCATE (errors, inds)
616 END IF
617
618 DEALLOCATE (subintervals)
619
620 CALL timestop(handle)
621 END SUBROUTINE simpsonrule_refine_integral
622
623! **************************************************************************************************
624!> \brief Approximate value of the integral on subinterval [a .. c] using the Simpson's rule.
625!> \param integral approximated integral = length / 6 * (fa + 4*fb + fc) (initialised on exit)
626!> \param fa integrand value at point a
627!> \param fb integrand value at point b = (a + c) / 2
628!> \param fc integrand value at point c
629!> \param length distance between points a and c [ABS(c-a)]
630!> \par History
631!> * 05.2017 created [Sergey Chulkov]
632! **************************************************************************************************
633 SUBROUTINE do_simpson_rule(integral, fa, fb, fc, length)
634 TYPE(cp_cfm_type), INTENT(IN) :: integral, fa, fb, fc
635 REAL(kind=dp), INTENT(in) :: length
636
637 CALL cp_cfm_to_cfm(fa, integral)
638 CALL cp_cfm_scale_and_add(z_one, integral, z_four, fb)
639 CALL cp_cfm_scale_and_add(z_one, integral, z_one, fc)
640 CALL cp_cfm_scale(length/6.0_dp, integral)
641 END SUBROUTINE do_simpson_rule
642
643! **************************************************************************************************
644!> \brief Approximate value of the integral on subinterval [a .. e] using the Boole's rule.
645!> \param integral approximated integral = length / 90 * (7*fa + 32*fb + 12*fc + 32*fd + 7*fe)
646!> (initialised on exit)
647!> \param fa integrand value at point a
648!> \param fb integrand value at point b = a + (e-a)/4
649!> \param fc integrand value at point c = a + (e-a)/2
650!> \param fd integrand value at point d = a + 3*(e-a)/4
651!> \param fe integrand value at point e
652!> \param length distance between points a and e [ABS(e-a)]
653!> \param work work matrix
654!> \par History
655!> * 05.2017 created [Sergey Chulkov]
656! **************************************************************************************************
657 SUBROUTINE do_boole_rule(integral, fa, fb, fc, fd, fe, length, work)
658 TYPE(cp_cfm_type), INTENT(IN) :: integral, fa, fb, fc, fd, fe
659 REAL(kind=dp), INTENT(in) :: length
660 TYPE(cp_cfm_type), INTENT(IN) :: work
661
662 REAL(kind=dp) :: rscale
663
664 rscale = length/90.0_dp
665
666 CALL cp_cfm_to_cfm(fc, integral)
667 CALL cp_cfm_scale(12.0_dp*rscale, integral)
668
669 CALL cp_cfm_to_cfm(fa, work)
670 CALL cp_cfm_scale_and_add(z_one, work, z_one, fe)
671 CALL cp_cfm_scale(7.0_dp*rscale, work)
672 CALL cp_cfm_scale_and_add(z_one, integral, z_one, work)
673
674 CALL cp_cfm_to_cfm(fb, work)
675 CALL cp_cfm_scale_and_add(z_one, work, z_one, fd)
676 CALL cp_cfm_scale(32.0_dp*rscale, work)
677 CALL cp_cfm_scale_and_add(z_one, integral, z_one, work)
678 END SUBROUTINE do_boole_rule
679END MODULE negf_integr_simpson
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:60
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:56
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:78
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
negf_integr_utils::equidistant_nodes_a_b
Definition negf_integr_utils.F:30
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:240
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:631
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:202
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
negf_integr_simpson
Adaptive Simpson's rule algorithm to integrate a complex-valued function in a complex plane.
Definition negf_integr_simpson.F:11
negf_integr_simpson::sr_shape_arc
integer, parameter, public sr_shape_arc
Definition negf_integr_simpson.F:46
negf_integr_simpson::simpsonrule_refine_integral
subroutine, public simpsonrule_refine_integral(sr_env, zdata_next)
Compute integral using the simpson's rules.
Definition negf_integr_simpson.F:342
negf_integr_simpson::simpsonrule_init
subroutine, public simpsonrule_init(sr_env, xnodes, nnodes, a, b, shape_id, conv, weights, tnodes_restart)
Initialise a Simpson's rule environment variable.
Definition negf_integr_simpson.F:122
negf_integr_simpson::simpsonrule_release
subroutine, public simpsonrule_release(sr_env)
Release a Simpson's rule environment variable.
Definition negf_integr_simpson.F:189
negf_integr_simpson::sr_shape_linear
integer, parameter, public sr_shape_linear
Definition negf_integr_simpson.F:46
negf_integr_simpson::simpsonrule_get_next_nodes
subroutine, public simpsonrule_get_next_nodes(sr_env, xnodes_next, nnodes)
Get the next set of nodes where to compute integrand.
Definition negf_integr_simpson.F:260
negf_integr_utils
Helper functions for integration routines.
Definition negf_integr_utils.F:13
negf_integr_utils::rescale_normalised_nodes
subroutine, public rescale_normalised_nodes(nnodes, tnodes, a, b, shape_id, xnodes, weights)
Definition negf_integr_utils.F:87
negf_integr_utils::contour_shape_arc
integer, parameter, public contour_shape_arc
Definition negf_integr_utils.F:27
negf_integr_utils::contour_shape_linear
integer, parameter, public contour_shape_linear
Definition negf_integr_utils.F:27
util
All kind of helpful little routines.
Definition util.F:14
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
negf_integr_simpson::simpsonrule_type
A structure to store data needed for adaptive Simpson's rule algorithm.
Definition negf_integr_simpson.F:70