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