Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 = -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 :
101 : CONTAINS
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 50 : 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 30 : POINTER :: w_data, w_data_my
137 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
138 :
139 30 : CALL timeset(routineN, handle)
140 :
141 30 : CPASSERT(nnodes > 4)
142 :
143 : ! ensure that MOD(nnodes-1, 4) == 0
144 30 : nnodes = 4*((nnodes - 1)/4) + 1
145 :
146 30 : sr_env%shape_id = shape_id
147 30 : sr_env%a = a
148 30 : sr_env%b = b
149 30 : sr_env%conv = conv
150 30 : sr_env%error = HUGE(0.0_dp)
151 30 : sr_env%error_conv = 0.0_dp
152 :
153 30 : NULLIFY (sr_env%error_fm, sr_env%weights)
154 30 : CALL cp_fm_get_info(weights, local_data=w_data, nrow_local=nrows, ncol_local=ncols, matrix_struct=fm_struct)
155 30 : ALLOCATE (sr_env%error_fm, sr_env%weights)
156 30 : CALL cp_fm_create(sr_env%error_fm, fm_struct)
157 30 : CALL cp_fm_create(sr_env%weights, fm_struct)
158 30 : 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 704 : DO icol = 1, ncols
162 8769 : DO irow = 1, nrows
163 8739 : w_data_my(irow, icol) = ABS(w_data(irow, icol))/15.0_dp
164 : END DO
165 : END DO
166 :
167 30 : NULLIFY (sr_env%integral, sr_env%integral_conv)
168 30 : NULLIFY (sr_env%integral_abc, sr_env%integral_cde, sr_env%integral_ace)
169 :
170 90 : ALLOCATE (sr_env%tnodes(nnodes))
171 :
172 30 : IF (PRESENT(tnodes_restart)) THEN
173 2880 : sr_env%tnodes(1:nnodes) = tnodes_restart(1:nnodes)
174 : ELSE
175 10 : CALL equidistant_nodes_a_b(-1.0_dp, 1.0_dp, nnodes, sr_env%tnodes)
176 : END IF
177 30 : CALL rescale_normalised_nodes(nnodes, sr_env%tnodes, a, b, shape_id, xnodes)
178 :
179 30 : CALL timestop(handle)
180 110 : 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 30 : 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 30 : CALL timeset(routineN, handle)
196 30 : IF (ALLOCATED(sr_env%subintervals)) THEN
197 0 : DO interval = SIZE(sr_env%subintervals), 1, -1
198 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fa)
199 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fb)
200 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fc)
201 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fd)
202 0 : CALL cp_cfm_release(sr_env%subintervals(interval)%fe)
203 : END DO
204 :
205 0 : DEALLOCATE (sr_env%subintervals)
206 : END IF
207 :
208 30 : IF (ASSOCIATED(sr_env%integral)) THEN
209 30 : CALL cp_cfm_release(sr_env%integral)
210 30 : DEALLOCATE (sr_env%integral)
211 : NULLIFY (sr_env%integral)
212 : END IF
213 30 : IF (ASSOCIATED(sr_env%integral_conv)) THEN
214 30 : CALL cp_cfm_release(sr_env%integral_conv)
215 30 : DEALLOCATE (sr_env%integral_conv)
216 : NULLIFY (sr_env%integral_conv)
217 : END IF
218 30 : IF (ASSOCIATED(sr_env%integral_abc)) THEN
219 30 : CALL cp_cfm_release(sr_env%integral_abc)
220 30 : DEALLOCATE (sr_env%integral_abc)
221 : NULLIFY (sr_env%integral_abc)
222 : END IF
223 30 : IF (ASSOCIATED(sr_env%integral_cde)) THEN
224 30 : CALL cp_cfm_release(sr_env%integral_cde)
225 30 : DEALLOCATE (sr_env%integral_cde)
226 : NULLIFY (sr_env%integral_cde)
227 : END IF
228 30 : IF (ASSOCIATED(sr_env%integral_ace)) THEN
229 30 : CALL cp_cfm_release(sr_env%integral_ace)
230 30 : DEALLOCATE (sr_env%integral_ace)
231 : NULLIFY (sr_env%integral_ace)
232 : END IF
233 30 : IF (ASSOCIATED(sr_env%error_fm)) THEN
234 30 : CALL cp_fm_release(sr_env%error_fm)
235 30 : DEALLOCATE (sr_env%error_fm)
236 : NULLIFY (sr_env%error_fm)
237 : END IF
238 30 : IF (ASSOCIATED(sr_env%weights)) THEN
239 30 : CALL cp_fm_release(sr_env%weights)
240 30 : DEALLOCATE (sr_env%weights)
241 : NULLIFY (sr_env%weights)
242 : END IF
243 :
244 30 : IF (ALLOCATED(sr_env%tnodes)) DEALLOCATE (sr_env%tnodes)
245 :
246 30 : CALL timestop(handle)
247 30 : 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 68 : 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 68 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes, tnodes_old
268 :
269 68 : CALL timeset(routineN, handle)
270 204 : ALLOCATE (tnodes(nnodes))
271 :
272 68 : CALL simpsonrule_get_next_nodes_real(sr_env, tnodes, nnodes)
273 68 : IF (nnodes > 0) THEN
274 68 : CALL MOVE_ALLOC(sr_env%tnodes, tnodes_old)
275 68 : nnodes_old = SIZE(tnodes_old)
276 :
277 204 : ALLOCATE (sr_env%tnodes(nnodes_old + nnodes))
278 6040 : sr_env%tnodes(1:nnodes_old) = tnodes_old(1:nnodes_old)
279 1108 : sr_env%tnodes(nnodes_old + 1:nnodes_old + nnodes) = tnodes(1:nnodes)
280 68 : DEALLOCATE (tnodes_old)
281 :
282 68 : CALL rescale_normalised_nodes(nnodes, tnodes, sr_env%a, sr_env%b, sr_env%shape_id, xnodes_next)
283 : END IF
284 :
285 68 : DEALLOCATE (tnodes)
286 68 : CALL timestop(handle)
287 136 : 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 68 : 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 68 : CALL timeset(routineN, handle)
307 :
308 68 : IF (ALLOCATED(sr_env%subintervals)) THEN
309 68 : nintervals = SIZE(sr_env%subintervals)
310 : ELSE
311 : nintervals = 0
312 : END IF
313 :
314 68 : IF (nintervals > 0) THEN
315 68 : IF (SIZE(xnodes_unity) < 4*nintervals) &
316 52 : nintervals = SIZE(xnodes_unity)/4
317 :
318 328 : DO interval = 1, nintervals
319 : xnodes_unity(4*interval - 3) = 0.125_dp* &
320 260 : (7.0_dp*sr_env%subintervals(interval)%lb + sr_env%subintervals(interval)%ub)
321 : xnodes_unity(4*interval - 2) = 0.125_dp* &
322 260 : (5.0_dp*sr_env%subintervals(interval)%lb + 3.0_dp*sr_env%subintervals(interval)%ub)
323 : xnodes_unity(4*interval - 1) = 0.125_dp* &
324 260 : (3.0_dp*sr_env%subintervals(interval)%lb + 5.0_dp*sr_env%subintervals(interval)%ub)
325 328 : xnodes_unity(4*interval) = 0.125_dp*(sr_env%subintervals(interval)%lb + 7.0_dp*sr_env%subintervals(interval)%ub)
326 : END DO
327 : END IF
328 :
329 68 : nnodes = 4*nintervals
330 68 : CALL timestop(handle)
331 68 : END SUBROUTINE simpsonrule_get_next_nodes_real
332 :
333 : ! **************************************************************************************************
334 : !> \brief Compute integral using the simpson's rules.
335 : !> \param sr_env Simpson's rule environment
336 : !> \param zdata_next precomputed integrand values at points xnodes_next (nullified on exit)
337 : !> \par History
338 : !> * 05.2017 created [Sergey Chulkov]
339 : ! **************************************************************************************************
340 98 : SUBROUTINE simpsonrule_refine_integral(sr_env, zdata_next)
341 : TYPE(simpsonrule_type), INTENT(inout) :: sr_env
342 : TYPE(cp_cfm_type), DIMENSION(:), INTENT(inout) :: zdata_next
343 :
344 : CHARACTER(len=*), PARAMETER :: routineN = 'simpsonrule_refine_integral'
345 : TYPE(cp_cfm_type), PARAMETER :: cfm_null = cp_cfm_type()
346 :
347 98 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: zscale
348 : COMPLEX(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
349 98 : POINTER :: error_zdata
350 : INTEGER :: handle, interval, ipoint, jpoint, &
351 : nintervals, nintervals_exist, npoints
352 98 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
353 : LOGICAL :: interval_converged, interval_exists
354 : REAL(kind=dp) :: my_bound, rscale
355 98 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: errors
356 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
357 98 : POINTER :: error_rdata
358 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
359 : TYPE(simpsonrule_subinterval_type), ALLOCATABLE, &
360 98 : DIMENSION(:) :: subintervals
361 :
362 98 : CALL timeset(routineN, handle)
363 :
364 98 : npoints = SIZE(zdata_next)
365 98 : IF (ASSOCIATED(sr_env%integral)) THEN
366 : ! we need 4 new points per subinterval (p, q, r, s)
367 : ! p q r s
368 : ! a . b . c . d . e
369 68 : CPASSERT(npoints > 0 .AND. MOD(npoints, 4) == 0)
370 : ELSE
371 : ! first call: need 4*n+1 points
372 : ! a1 b1 c1 d1 e1
373 : ! a2 b2 c2 d2 e2
374 : ! a3 b3 c3 d3 e3
375 30 : CPASSERT(npoints > 1 .AND. MOD(npoints, 4) == 1)
376 : END IF
377 :
378 : ! compute weights of new points on a complex contour according to their values of the 't' parameter
379 98 : nintervals_exist = SIZE(sr_env%tnodes)
380 98 : CPASSERT(nintervals_exist >= npoints)
381 294 : ALLOCATE (zscale(npoints))
382 :
383 : CALL rescale_normalised_nodes(npoints, sr_env%tnodes(nintervals_exist - npoints + 1:nintervals_exist), &
384 98 : sr_env%a, sr_env%b, sr_env%shape_id, weights=zscale)
385 :
386 : ! rescale integrand values
387 4128 : DO ipoint = 1, npoints
388 4128 : CALL cp_cfm_scale(zscale(ipoint), zdata_next(ipoint))
389 : END DO
390 :
391 98 : DEALLOCATE (zscale)
392 :
393 : ! insert new points
394 98 : nintervals = npoints/4
395 98 : IF (ASSOCIATED(sr_env%integral)) THEN
396 : ! subdivide existing intervals
397 68 : nintervals_exist = SIZE(sr_env%subintervals)
398 68 : CPASSERT(nintervals <= nintervals_exist)
399 :
400 1624 : ALLOCATE (subintervals(nintervals_exist + nintervals))
401 :
402 328 : DO interval = 1, nintervals
403 260 : subintervals(2*interval - 1)%lb = sr_env%subintervals(interval)%lb
404 260 : subintervals(2*interval - 1)%ub = 0.5_dp*(sr_env%subintervals(interval)%lb + sr_env%subintervals(interval)%ub)
405 260 : subintervals(2*interval - 1)%conv = 0.5_dp*sr_env%subintervals(interval)%conv
406 260 : subintervals(2*interval - 1)%fa = sr_env%subintervals(interval)%fa
407 260 : subintervals(2*interval - 1)%fb = zdata_next(4*interval - 3)
408 260 : subintervals(2*interval - 1)%fc = sr_env%subintervals(interval)%fb
409 260 : subintervals(2*interval - 1)%fd = zdata_next(4*interval - 2)
410 260 : subintervals(2*interval - 1)%fe = sr_env%subintervals(interval)%fc
411 :
412 260 : subintervals(2*interval)%lb = subintervals(2*interval - 1)%ub
413 260 : subintervals(2*interval)%ub = sr_env%subintervals(interval)%ub
414 260 : subintervals(2*interval)%conv = subintervals(2*interval - 1)%conv
415 260 : subintervals(2*interval)%fa = sr_env%subintervals(interval)%fc
416 260 : subintervals(2*interval)%fb = zdata_next(4*interval - 1)
417 260 : subintervals(2*interval)%fc = sr_env%subintervals(interval)%fd
418 260 : subintervals(2*interval)%fd = zdata_next(4*interval)
419 260 : subintervals(2*interval)%fe = sr_env%subintervals(interval)%fe
420 :
421 1368 : zdata_next(4*interval - 3:4*interval) = cfm_null
422 : END DO
423 :
424 968 : DO interval = nintervals + 1, nintervals_exist
425 968 : subintervals(interval + nintervals) = sr_env%subintervals(interval)
426 : END DO
427 68 : DEALLOCATE (sr_env%subintervals)
428 : ELSE
429 : ! first time -- allocate matrices and create a new set of intervals
430 30 : CALL cp_cfm_get_info(zdata_next(1), matrix_struct=fm_struct)
431 : ALLOCATE (sr_env%integral, sr_env%integral_conv, &
432 30 : sr_env%integral_abc, sr_env%integral_cde, sr_env%integral_ace)
433 30 : CALL cp_cfm_create(sr_env%integral, fm_struct)
434 30 : CALL cp_cfm_create(sr_env%integral_conv, fm_struct)
435 30 : CALL cp_cfm_create(sr_env%integral_abc, fm_struct)
436 30 : CALL cp_cfm_create(sr_env%integral_cde, fm_struct)
437 30 : CALL cp_cfm_create(sr_env%integral_ace, fm_struct)
438 :
439 30 : CALL cp_cfm_set_all(sr_env%integral_conv, z_zero)
440 :
441 830 : ALLOCATE (subintervals(nintervals))
442 :
443 30 : rscale = 1.0_dp/REAL(nintervals, kind=dp)
444 :
445 770 : DO interval = 1, nintervals
446 : ! lower bound: point with indices 1, 5, 9, ..., 4*nintervals+1
447 740 : subintervals(interval)%lb = sr_env%tnodes(4*interval - 3)
448 740 : subintervals(interval)%ub = sr_env%tnodes(4*interval + 1)
449 740 : subintervals(interval)%conv = rscale*sr_env%conv
450 :
451 740 : subintervals(interval)%fa = zdata_next(4*interval - 3)
452 740 : subintervals(interval)%fb = zdata_next(4*interval - 2)
453 740 : subintervals(interval)%fc = zdata_next(4*interval - 1)
454 740 : subintervals(interval)%fd = zdata_next(4*interval)
455 770 : subintervals(interval)%fe = zdata_next(4*interval + 1)
456 : END DO
457 : END IF
458 :
459 : ! we kept the originals matrices for internal use, so set the matrix to null
460 : ! to prevent alteration of the matrices from the outside
461 4128 : zdata_next(1:npoints) = cfm_null
462 :
463 98 : CALL cp_fm_get_info(sr_env%error_fm, local_data=error_rdata)
464 98 : CALL cp_cfm_get_info(sr_env%integral_ace, local_data=error_zdata)
465 :
466 : ! do actual integration
467 98 : CALL cp_cfm_to_cfm(sr_env%integral_conv, sr_env%integral)
468 98 : sr_env%error = sr_env%error_conv
469 98 : nintervals_exist = SIZE(subintervals)
470 :
471 2258 : DO interval = 1, nintervals_exist
472 2160 : rscale = subintervals(interval)%ub - subintervals(interval)%lb
473 : CALL do_simpson_rule(sr_env%integral_ace, &
474 : subintervals(interval)%fa, &
475 : subintervals(interval)%fc, &
476 : subintervals(interval)%fe, &
477 2160 : -0.5_dp*rscale)
478 : CALL do_simpson_rule(sr_env%integral_abc, &
479 : subintervals(interval)%fa, &
480 : subintervals(interval)%fb, &
481 : subintervals(interval)%fc, &
482 2160 : 0.25_dp*rscale)
483 : CALL do_simpson_rule(sr_env%integral_cde, &
484 : subintervals(interval)%fc, &
485 : subintervals(interval)%fd, &
486 : subintervals(interval)%fe, &
487 2160 : 0.25_dp*rscale)
488 :
489 2160 : CALL cp_cfm_scale_and_add(z_one, sr_env%integral_abc, z_one, sr_env%integral_cde)
490 2160 : CALL cp_cfm_scale_and_add(z_one, sr_env%integral_ace, z_one, sr_env%integral_abc)
491 :
492 : IF (is_boole) THEN
493 : CALL do_boole_rule(sr_env%integral_abc, &
494 : subintervals(interval)%fa, &
495 : subintervals(interval)%fb, &
496 : subintervals(interval)%fc, &
497 : subintervals(interval)%fd, &
498 : subintervals(interval)%fe, &
499 : 0.5_dp*rscale, sr_env%integral_cde)
500 : END IF
501 :
502 2160 : CALL cp_cfm_scale_and_add(z_one, sr_env%integral, z_one, sr_env%integral_abc)
503 :
504 : ! sr_env%error_fm = ABS(sr_env%integral_ace); no temporary arrays as pointers have different types
505 674220 : error_rdata(:, :) = ABS(error_zdata(:, :))
506 2160 : CALL cp_fm_trace(sr_env%error_fm, sr_env%weights, subintervals(interval)%error)
507 :
508 2160 : sr_env%error = sr_env%error + subintervals(interval)%error
509 :
510 : ! add contributions from converged subintervals, so we could drop them afterward
511 2258 : IF (subintervals(interval)%error <= subintervals(interval)%conv) THEN
512 766 : CALL cp_cfm_scale_and_add(z_one, sr_env%integral_conv, z_one, sr_env%integral_abc)
513 766 : sr_env%error_conv = sr_env%error_conv + subintervals(interval)%error
514 : END IF
515 : END DO
516 :
517 98 : IF (sr_env%error <= sr_env%conv) THEN
518 : ! We have already reached the target accuracy, so we can drop all subintervals
519 : ! (even those where local convergence has not been achieved). From now on environment
520 : ! components 'sr_env%error' and 'sr_env%integral_conv' hold incorrect values,
521 : ! but they should not been accessed from the outside anyway
522 : ! (uncomment the following two lines if they are actually need)
523 :
524 : ! sr_env%error_conv = sr_env%error
525 : ! CALL cp_cfm_to_cfm(sr_env%integral, sr_env%integral_conv)
526 :
527 : ! Only deallocate the fa component explicitly if there is no interval to the left from it
528 894 : DO interval = nintervals_exist, 1, -1
529 864 : interval_exists = .FALSE.
530 864 : my_bound = subintervals(interval)%lb
531 18824 : DO jpoint = 1, nintervals_exist
532 18824 : IF (subintervals(jpoint)%ub == my_bound) THEN
533 : interval_exists = .TRUE.
534 : EXIT
535 : END IF
536 : END DO
537 864 : IF (.NOT. interval_exists) THEN
538 : ! interval does not exist anymore, so it is safe to release the matrix
539 58 : CALL cp_cfm_release(subintervals(interval)%fa)
540 : ELSE IF (interval_converged) THEN
541 : ! the interval still exists and will be released with fe
542 : END IF
543 864 : CALL cp_cfm_release(subintervals(interval)%fb)
544 864 : CALL cp_cfm_release(subintervals(interval)%fc)
545 864 : CALL cp_cfm_release(subintervals(interval)%fd)
546 894 : CALL cp_cfm_release(subintervals(interval)%fe)
547 : END DO
548 : ELSE
549 : ! sort subinterval according to their convergence, and drop convergent ones
550 340 : ALLOCATE (errors(nintervals_exist), inds(nintervals_exist))
551 :
552 1364 : nintervals = 0
553 1364 : DO interval = 1, nintervals_exist
554 1296 : errors(interval) = subintervals(interval)%error
555 :
556 1296 : IF (subintervals(interval)%error > subintervals(interval)%conv) &
557 1228 : nintervals = nintervals + 1
558 : END DO
559 :
560 68 : CALL sort(errors, nintervals_exist, inds)
561 :
562 68 : IF (nintervals > 0) &
563 1364 : ALLOCATE (sr_env%subintervals(nintervals))
564 :
565 : nintervals = 0
566 1364 : DO ipoint = nintervals_exist, 1, -1
567 1296 : interval = inds(ipoint)
568 :
569 1364 : IF (subintervals(interval)%error > subintervals(interval)%conv) THEN
570 1160 : nintervals = nintervals + 1
571 :
572 1160 : sr_env%subintervals(nintervals) = subintervals(interval)
573 : ELSE
574 : ! Release matrices of converged intervals. Special cases: left and right boundary
575 : ! Check whether the neighboring interval still exists and if it does, check for its convergence
576 136 : interval_exists = .FALSE.
577 136 : my_bound = subintervals(interval)%lb
578 1538 : DO jpoint = 1, nintervals_exist
579 1538 : IF (subintervals(jpoint)%ub == my_bound) THEN
580 : interval_exists = .TRUE.
581 : EXIT
582 : END IF
583 : END DO
584 136 : IF (.NOT. interval_exists) THEN
585 : ! interval does not exist anymore, so it is safe to release the matrix
586 24 : CALL cp_cfm_release(subintervals(interval)%fa)
587 : ELSE IF (interval_converged) THEN
588 : ! the interval still exists and will be released with fe
589 : END IF
590 136 : CALL cp_cfm_release(subintervals(interval)%fb)
591 136 : CALL cp_cfm_release(subintervals(interval)%fc)
592 136 : CALL cp_cfm_release(subintervals(interval)%fd)
593 :
594 : ! Right border: Check for the existence and the convergence of the interval
595 : ! If the right interval does not exist or has converged, release the matrix
596 136 : interval_exists = .FALSE.
597 136 : interval_converged = .FALSE.
598 136 : my_bound = subintervals(interval)%ub
599 1670 : DO jpoint = 1, nintervals_exist
600 1670 : IF (subintervals(jpoint)%lb == my_bound) THEN
601 112 : interval_exists = .TRUE.
602 112 : IF (subintervals(jpoint)%error <= subintervals(jpoint)%conv) interval_converged = .TRUE.
603 : EXIT
604 : END IF
605 : END DO
606 136 : IF (.NOT. interval_exists .OR. interval_converged) THEN
607 84 : CALL cp_cfm_release(subintervals(interval)%fe)
608 : END IF
609 : END IF
610 : END DO
611 :
612 68 : DEALLOCATE (errors, inds)
613 : END IF
614 :
615 98 : DEALLOCATE (subintervals)
616 :
617 98 : CALL timestop(handle)
618 98 : END SUBROUTINE simpsonrule_refine_integral
619 :
620 : ! **************************************************************************************************
621 : !> \brief Approximate value of the integral on subinterval [a .. c] using the Simpson's rule.
622 : !> \param integral approximated integral = length / 6 * (fa + 4*fb + fc) (initialised on exit)
623 : !> \param fa integrand value at point a
624 : !> \param fb integrand value at point b = (a + c) / 2
625 : !> \param fc integrand value at point c
626 : !> \param length distance between points a and c [ABS(c-a)]
627 : !> \par History
628 : !> * 05.2017 created [Sergey Chulkov]
629 : ! **************************************************************************************************
630 6480 : SUBROUTINE do_simpson_rule(integral, fa, fb, fc, length)
631 : TYPE(cp_cfm_type), INTENT(IN) :: integral, fa, fb, fc
632 : REAL(kind=dp), INTENT(in) :: length
633 :
634 6480 : CALL cp_cfm_to_cfm(fa, integral)
635 6480 : CALL cp_cfm_scale_and_add(z_one, integral, z_four, fb)
636 6480 : CALL cp_cfm_scale_and_add(z_one, integral, z_one, fc)
637 6480 : CALL cp_cfm_scale(length/6.0_dp, integral)
638 6480 : END SUBROUTINE do_simpson_rule
639 :
640 : ! **************************************************************************************************
641 : !> \brief Approximate value of the integral on subinterval [a .. e] using the Boole's rule.
642 : !> \param integral approximated integral = length / 90 * (7*fa + 32*fb + 12*fc + 32*fd + 7*fe)
643 : !> (initialised on exit)
644 : !> \param fa integrand value at point a
645 : !> \param fb integrand value at point b = a + (e-a)/4
646 : !> \param fc integrand value at point c = a + (e-a)/2
647 : !> \param fd integrand value at point d = a + 3*(e-a)/4
648 : !> \param fe integrand value at point e
649 : !> \param length distance between points a and e [ABS(e-a)]
650 : !> \param work work matrix
651 : !> \par History
652 : !> * 05.2017 created [Sergey Chulkov]
653 : ! **************************************************************************************************
654 0 : SUBROUTINE do_boole_rule(integral, fa, fb, fc, fd, fe, length, work)
655 : TYPE(cp_cfm_type), INTENT(IN) :: integral, fa, fb, fc, fd, fe
656 : REAL(kind=dp), INTENT(in) :: length
657 : TYPE(cp_cfm_type), INTENT(IN) :: work
658 :
659 : REAL(kind=dp) :: rscale
660 :
661 0 : rscale = length/90.0_dp
662 :
663 0 : CALL cp_cfm_to_cfm(fc, integral)
664 0 : CALL cp_cfm_scale(12.0_dp*rscale, integral)
665 :
666 0 : CALL cp_cfm_to_cfm(fa, work)
667 0 : CALL cp_cfm_scale_and_add(z_one, work, z_one, fe)
668 0 : CALL cp_cfm_scale(7.0_dp*rscale, work)
669 0 : CALL cp_cfm_scale_and_add(z_one, integral, z_one, work)
670 :
671 0 : CALL cp_cfm_to_cfm(fb, work)
672 0 : CALL cp_cfm_scale_and_add(z_one, work, z_one, fd)
673 0 : CALL cp_cfm_scale(32.0_dp*rscale, work)
674 0 : CALL cp_cfm_scale_and_add(z_one, integral, z_one, work)
675 0 : END SUBROUTINE do_boole_rule
676 0 : END MODULE negf_integr_simpson
|
