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