Line data Source code
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 50 : 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 30 : POINTER :: w_data, w_data_my
136 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
137 :
138 30 : CALL timeset(routineN, handle)
139 :
140 30 : CPASSERT(nnodes > 4)
141 :
142 : ! ensure that MOD(nnodes-1, 4) == 0
143 30 : nnodes = 4*((nnodes - 1)/4) + 1
144 :
145 30 : sr_env%shape_id = shape_id
146 30 : sr_env%a = a
147 30 : sr_env%b = b
148 30 : sr_env%conv = conv
149 30 : sr_env%error = HUGE(0.0_dp)
150 30 : sr_env%error_conv = 0.0_dp
151 :
152 30 : NULLIFY (sr_env%error_fm, sr_env%weights)
153 30 : CALL cp_fm_get_info(weights, local_data=w_data, nrow_local=nrows, ncol_local=ncols, matrix_struct=fm_struct)
154 30 : ALLOCATE (sr_env%error_fm, sr_env%weights)
155 30 : CALL cp_fm_create(sr_env%error_fm, fm_struct)
156 30 : CALL cp_fm_create(sr_env%weights, fm_struct)
157 30 : 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 704 : DO icol = 1, ncols
161 8769 : DO irow = 1, nrows
162 8739 : w_data_my(irow, icol) = ABS(w_data(irow, icol))/15.0_dp
163 : END DO
164 : END DO
165 :
166 30 : NULLIFY (sr_env%integral, sr_env%integral_conv)
167 30 : NULLIFY (sr_env%integral_abc, sr_env%integral_cde, sr_env%integral_ace)
168 :
169 90 : ALLOCATE (sr_env%tnodes(nnodes))
170 :
171 30 : IF (PRESENT(tnodes_restart)) THEN
172 2880 : sr_env%tnodes(1:nnodes) = tnodes_restart(1:nnodes)
173 : ELSE
174 10 : CALL equidistant_nodes_a_b(-1.0_dp, 1.0_dp, nnodes, sr_env%tnodes)
175 : END IF
176 30 : CALL rescale_normalised_nodes(nnodes, sr_env%tnodes, a, b, shape_id, xnodes)
177 :
178 30 : CALL timestop(handle)
179 110 : 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 30 : 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 30 : CALL timeset(routineN, handle)
195 30 : IF (ALLOCATED(sr_env%subintervals)) THEN
196 0 : DO interval = SIZE(sr_env%subintervals), 1, -1
197 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fa)
198 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fb)
199 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fc)
200 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fd)
201 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fe)
202 : END DO
203 :
204 0 : DEALLOCATE (sr_env%subintervals)
205 : END IF
206 :
207 30 : IF (ASSOCIATED(sr_env%integral)) THEN
208 30 : CALL cp_cfm_release(sr_env%integral)
209 30 : DEALLOCATE (sr_env%integral)
210 : NULLIFY (sr_env%integral)
211 : END IF
212 30 : IF (ASSOCIATED(sr_env%integral_conv)) THEN
213 30 : CALL cp_cfm_release(sr_env%integral_conv)
214 30 : DEALLOCATE (sr_env%integral_conv)
215 : NULLIFY (sr_env%integral_conv)
216 : END IF
217 30 : IF (ASSOCIATED(sr_env%integral_abc)) THEN
218 30 : CALL cp_cfm_release(sr_env%integral_abc)
219 30 : DEALLOCATE (sr_env%integral_abc)
220 : NULLIFY (sr_env%integral_abc)
221 : END IF
222 30 : IF (ASSOCIATED(sr_env%integral_cde)) THEN
223 30 : CALL cp_cfm_release(sr_env%integral_cde)
224 30 : DEALLOCATE (sr_env%integral_cde)
225 : NULLIFY (sr_env%integral_cde)
226 : END IF
227 30 : IF (ASSOCIATED(sr_env%integral_ace)) THEN
228 30 : CALL cp_cfm_release(sr_env%integral_ace)
229 30 : DEALLOCATE (sr_env%integral_ace)
230 : NULLIFY (sr_env%integral_ace)
231 : END IF
232 30 : IF (ASSOCIATED(sr_env%error_fm)) THEN
233 30 : CALL cp_fm_release(sr_env%error_fm)
234 30 : DEALLOCATE (sr_env%error_fm)
235 : NULLIFY (sr_env%error_fm)
236 : END IF
237 30 : IF (ASSOCIATED(sr_env%weights)) THEN
238 30 : CALL cp_fm_release(sr_env%weights)
239 30 : DEALLOCATE (sr_env%weights)
240 : NULLIFY (sr_env%weights)
241 : END IF
242 :
243 30 : IF (ALLOCATED(sr_env%tnodes)) DEALLOCATE (sr_env%tnodes)
244 :
245 30 : CALL timestop(handle)
246 30 : 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 68 : 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 68 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes, tnodes_old
267 :
268 68 : CALL timeset(routineN, handle)
269 204 : ALLOCATE (tnodes(nnodes))
270 :
271 68 : CALL simpsonrule_get_next_nodes_real(sr_env, tnodes, nnodes)
272 68 : IF (nnodes > 0) THEN
273 68 : CALL MOVE_ALLOC(sr_env%tnodes, tnodes_old)
274 68 : nnodes_old = SIZE(tnodes_old)
275 :
276 204 : ALLOCATE (sr_env%tnodes(nnodes_old + nnodes))
277 6040 : sr_env%tnodes(1:nnodes_old) = tnodes_old(1:nnodes_old)
278 1108 : sr_env%tnodes(nnodes_old + 1:nnodes_old + nnodes) = tnodes(1:nnodes)
279 68 : DEALLOCATE (tnodes_old)
280 :
281 68 : CALL rescale_normalised_nodes(nnodes, tnodes, sr_env%a, sr_env%b, sr_env%shape_id, xnodes_next)
282 : END IF
283 :
284 68 : DEALLOCATE (tnodes)
285 68 : CALL timestop(handle)
286 136 : 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 68 : 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 68 : CALL timeset(routineN, handle)
306 :
307 68 : IF (ALLOCATED(sr_env%subintervals)) THEN
308 68 : nintervals = SIZE(sr_env%subintervals)
309 : ELSE
310 : nintervals = 0
311 : END IF
312 :
313 68 : IF (nintervals > 0) THEN
314 68 : IF (SIZE(xnodes_unity) < 4*nintervals) &
315 52 : nintervals = SIZE(xnodes_unity)/4
316 :
317 328 : DO interval = 1, nintervals
318 : xnodes_unity(4*interval - 3) = 0.125_dp* &
319 260 : (7.0_dp*sr_env%subintervals(interval)%lb + sr_env%subintervals(interval)%ub)
320 : xnodes_unity(4*interval - 2) = 0.125_dp* &
321 260 : (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 260 : (3.0_dp*sr_env%subintervals(interval)%lb + 5.0_dp*sr_env%subintervals(interval)%ub)
324 328 : 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 68 : nnodes = 4*nintervals
329 68 : CALL timestop(handle)
330 68 : 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 98 : 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 98 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: zscale
347 : COMPLEX(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
348 98 : POINTER :: error_zdata
349 : INTEGER :: handle, interval, ipoint, jpoint, &
350 : nintervals, nintervals_exist, npoints
351 98 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
352 : LOGICAL :: interval_converged, interval_exists
353 : REAL(kind=dp) :: my_bound, rscale
354 98 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: errors
355 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
356 98 : POINTER :: error_rdata
357 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
358 : TYPE(simpsonrule_subinterval_type), ALLOCATABLE, &
359 98 : DIMENSION(:) :: subintervals
360 :
361 98 : CALL timeset(routineN, handle)
362 :
363 98 : npoints = SIZE(zdata_next)
364 98 : 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 68 : 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 30 : 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 98 : nintervals_exist = SIZE(sr_env%tnodes)
379 98 : CPASSERT(nintervals_exist >= npoints)
380 294 : ALLOCATE (zscale(npoints))
381 :
382 : CALL rescale_normalised_nodes(npoints, sr_env%tnodes(nintervals_exist - npoints + 1:nintervals_exist), &
383 98 : sr_env%a, sr_env%b, sr_env%shape_id, weights=zscale)
384 :
385 : ! rescale integrand values
386 4128 : DO ipoint = 1, npoints
387 4128 : CALL cp_cfm_scale(zscale(ipoint), zdata_next(ipoint))
388 : END DO
389 :
390 98 : DEALLOCATE (zscale)
391 :
392 : ! insert new points
393 98 : nintervals = npoints/4
394 98 : IF (ASSOCIATED(sr_env%integral)) THEN
395 : ! subdivide existing intervals
396 68 : nintervals_exist = SIZE(sr_env%subintervals)
397 68 : CPASSERT(nintervals <= nintervals_exist)
398 :
399 1624 : ALLOCATE (subintervals(nintervals_exist + nintervals))
400 :
401 328 : DO interval = 1, nintervals
402 260 : subintervals(2*interval - 1)%lb = sr_env%subintervals(interval)%lb
403 260 : subintervals(2*interval - 1)%ub = 0.5_dp*(sr_env%subintervals(interval)%lb + sr_env%subintervals(interval)%ub)
404 260 : subintervals(2*interval - 1)%conv = 0.5_dp*sr_env%subintervals(interval)%conv
405 260 : subintervals(2*interval - 1)%fa = sr_env%subintervals(interval)%fa
406 260 : subintervals(2*interval - 1)%fb = zdata_next(4*interval - 3)
407 260 : subintervals(2*interval - 1)%fc = sr_env%subintervals(interval)%fb
408 260 : subintervals(2*interval - 1)%fd = zdata_next(4*interval - 2)
409 260 : subintervals(2*interval - 1)%fe = sr_env%subintervals(interval)%fc
410 :
411 260 : subintervals(2*interval)%lb = subintervals(2*interval - 1)%ub
412 260 : subintervals(2*interval)%ub = sr_env%subintervals(interval)%ub
413 260 : subintervals(2*interval)%conv = subintervals(2*interval - 1)%conv
414 260 : subintervals(2*interval)%fa = sr_env%subintervals(interval)%fc
415 260 : subintervals(2*interval)%fb = zdata_next(4*interval - 1)
416 260 : subintervals(2*interval)%fc = sr_env%subintervals(interval)%fd
417 260 : subintervals(2*interval)%fd = zdata_next(4*interval)
418 260 : subintervals(2*interval)%fe = sr_env%subintervals(interval)%fe
419 :
420 1368 : zdata_next(4*interval - 3:4*interval) = cfm_null
421 : END DO
422 :
423 968 : DO interval = nintervals + 1, nintervals_exist
424 968 : subintervals(interval + nintervals) = sr_env%subintervals(interval)
425 : END DO
426 68 : DEALLOCATE (sr_env%subintervals)
427 : ELSE
428 : ! first time -- allocate matrices and create a new set of intervals
429 30 : CALL cp_cfm_get_info(zdata_next(1), matrix_struct=fm_struct)
430 : ALLOCATE (sr_env%integral, sr_env%integral_conv, &
431 30 : sr_env%integral_abc, sr_env%integral_cde, sr_env%integral_ace)
432 30 : CALL cp_cfm_create(sr_env%integral, fm_struct)
433 30 : CALL cp_cfm_create(sr_env%integral_conv, fm_struct)
434 30 : CALL cp_cfm_create(sr_env%integral_abc, fm_struct)
435 30 : CALL cp_cfm_create(sr_env%integral_cde, fm_struct)
436 30 : CALL cp_cfm_create(sr_env%integral_ace, fm_struct)
437 :
438 30 : CALL cp_cfm_set_all(sr_env%integral_conv, z_zero)
439 :
440 830 : ALLOCATE (subintervals(nintervals))
441 :
442 30 : rscale = 1.0_dp/REAL(nintervals, kind=dp)
443 :
444 770 : DO interval = 1, nintervals
445 : ! lower bound: point with indices 1, 5, 9, ..., 4*nintervals+1
446 740 : subintervals(interval)%lb = sr_env%tnodes(4*interval - 3)
447 740 : subintervals(interval)%ub = sr_env%tnodes(4*interval + 1)
448 740 : subintervals(interval)%conv = rscale*sr_env%conv
449 :
450 740 : subintervals(interval)%fa = zdata_next(4*interval - 3)
451 740 : subintervals(interval)%fb = zdata_next(4*interval - 2)
452 740 : subintervals(interval)%fc = zdata_next(4*interval - 1)
453 740 : subintervals(interval)%fd = zdata_next(4*interval)
454 770 : 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 4128 : zdata_next(1:npoints) = cfm_null
461 :
462 98 : CALL cp_fm_get_info(sr_env%error_fm, local_data=error_rdata)
463 98 : CALL cp_cfm_get_info(sr_env%integral_ace, local_data=error_zdata)
464 :
465 : ! do actual integration
466 98 : CALL cp_cfm_to_cfm(sr_env%integral_conv, sr_env%integral)
467 98 : sr_env%error = sr_env%error_conv
468 98 : nintervals_exist = SIZE(subintervals)
469 :
470 2258 : DO interval = 1, nintervals_exist
471 2160 : 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 2160 : -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 2160 : 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 2160 : 0.25_dp*rscale)
487 :
488 2160 : CALL cp_cfm_scale_and_add(z_one, sr_env%integral_abc, z_one, sr_env%integral_cde)
489 2160 : 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 2160 : 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 674220 : error_rdata(:, :) = ABS(error_zdata(:, :))
505 2160 : CALL cp_fm_trace(sr_env%error_fm, sr_env%weights, subintervals(interval)%error)
506 :
507 2160 : sr_env%error = sr_env%error + subintervals(interval)%error
508 :
509 : ! add contributions from converged subintervals, so we could drop them afterward
510 2258 : IF (subintervals(interval)%error <= subintervals(interval)%conv) THEN
511 766 : CALL cp_cfm_scale_and_add(z_one, sr_env%integral_conv, z_one, sr_env%integral_abc)
512 766 : sr_env%error_conv = sr_env%error_conv + subintervals(interval)%error
513 : END IF
514 : END DO
515 :
516 98 : 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 894 : DO interval = nintervals_exist, 1, -1
528 864 : interval_exists = .FALSE.
529 864 : my_bound = subintervals(interval)%lb
530 18824 : DO jpoint = 1, nintervals_exist
531 18824 : IF (subintervals(jpoint)%ub == my_bound) THEN
532 : interval_exists = .TRUE.
533 : EXIT
534 : END IF
535 : END DO
536 864 : IF (.NOT. interval_exists) THEN
537 : ! interval does not exist anymore, so it is safe to release the matrix
538 58 : 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 864 : CALL cp_cfm_release(subintervals(interval)%fb)
543 864 : CALL cp_cfm_release(subintervals(interval)%fc)
544 864 : CALL cp_cfm_release(subintervals(interval)%fd)
545 894 : CALL cp_cfm_release(subintervals(interval)%fe)
546 : END DO
547 : ELSE
548 : ! sort subinterval according to their convergence, and drop convergent ones
549 340 : ALLOCATE (errors(nintervals_exist), inds(nintervals_exist))
550 :
551 1364 : nintervals = 0
552 1364 : DO interval = 1, nintervals_exist
553 1296 : errors(interval) = subintervals(interval)%error
554 :
555 1296 : IF (subintervals(interval)%error > subintervals(interval)%conv) &
556 1228 : nintervals = nintervals + 1
557 : END DO
558 :
559 68 : CALL sort(errors, nintervals_exist, inds)
560 :
561 68 : IF (nintervals > 0) &
562 1364 : ALLOCATE (sr_env%subintervals(nintervals))
563 :
564 : nintervals = 0
565 1364 : DO ipoint = nintervals_exist, 1, -1
566 1296 : interval = inds(ipoint)
567 :
568 1364 : IF (subintervals(interval)%error > subintervals(interval)%conv) THEN
569 1160 : nintervals = nintervals + 1
570 :
571 1160 : 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 136 : interval_exists = .FALSE.
576 136 : my_bound = subintervals(interval)%lb
577 1538 : DO jpoint = 1, nintervals_exist
578 1538 : IF (subintervals(jpoint)%ub == my_bound) THEN
579 : interval_exists = .TRUE.
580 : EXIT
581 : END IF
582 : END DO
583 136 : IF (.NOT. interval_exists) THEN
584 : ! interval does not exist anymore, so it is safe to release the matrix
585 24 : 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 136 : CALL cp_cfm_release(subintervals(interval)%fb)
590 136 : CALL cp_cfm_release(subintervals(interval)%fc)
591 136 : 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 136 : interval_exists = .FALSE.
596 136 : interval_converged = .FALSE.
597 136 : my_bound = subintervals(interval)%ub
598 1670 : DO jpoint = 1, nintervals_exist
599 1670 : IF (subintervals(jpoint)%lb == my_bound) THEN
600 112 : interval_exists = .TRUE.
601 112 : IF (subintervals(jpoint)%error <= subintervals(jpoint)%conv) interval_converged = .TRUE.
602 : EXIT
603 : END IF
604 : END DO
605 136 : IF (.NOT. interval_exists .OR. interval_converged) THEN
606 84 : CALL cp_cfm_release(subintervals(interval)%fe)
607 : END IF
608 : END IF
609 : END DO
610 :
611 68 : DEALLOCATE (errors, inds)
612 : END IF
613 :
614 98 : DEALLOCATE (subintervals)
615 :
616 98 : CALL timestop(handle)
617 98 : 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 6480 : 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 6480 : CALL cp_cfm_to_cfm(fa, integral)
634 6480 : CALL cp_cfm_scale_and_add(z_one, integral, z_four, fb)
635 6480 : CALL cp_cfm_scale_and_add(z_one, integral, z_one, fc)
636 6480 : CALL cp_cfm_scale(length/6.0_dp, integral)
637 6480 : 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 0 : 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 0 : rscale = length/90.0_dp
661 :
662 0 : CALL cp_cfm_to_cfm(fc, integral)
663 0 : CALL cp_cfm_scale(12.0_dp*rscale, integral)
664 :
665 0 : CALL cp_cfm_to_cfm(fa, work)
666 0 : CALL cp_cfm_scale_and_add(z_one, work, z_one, fe)
667 0 : CALL cp_cfm_scale(7.0_dp*rscale, work)
668 0 : CALL cp_cfm_scale_and_add(z_one, integral, z_one, work)
669 :
670 0 : CALL cp_cfm_to_cfm(fb, work)
671 0 : CALL cp_cfm_scale_and_add(z_one, work, z_one, fd)
672 0 : CALL cp_cfm_scale(32.0_dp*rscale, work)
673 0 : CALL cp_cfm_scale_and_add(z_one, integral, z_one, work)
674 0 : END SUBROUTINE do_boole_rule
675 0 : END MODULE negf_integr_simpson
|
