1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Adaptive Clenshaw-Curtis quadrature algorithm to integrate a complex-valued function in 8!> a complex plane 9!> \par History 10!> * 05.2017 created [Sergey Chulkov] 11! ************************************************************************************************** 12MODULE negf_integr_cc 13 USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale,& 14 cp_cfm_scale_and_add 15 USE cp_cfm_types, ONLY: cp_cfm_create,& 16 cp_cfm_get_info,& 17 cp_cfm_p_type,& 18 cp_cfm_release,& 19 cp_cfm_type 20 USE cp_fm_basic_linalg, ONLY: cp_fm_trace 21 USE cp_fm_struct, ONLY: cp_fm_struct_equivalent,& 22 cp_fm_struct_type 23 USE cp_fm_types, ONLY: cp_fm_create,& 24 cp_fm_get_info,& 25 cp_fm_release,& 26 cp_fm_type 27 USE fft_tools, ONLY: fft_fw1d 28 USE kahan_sum, ONLY: accurate_sum 29 USE kinds, ONLY: dp,& 30 int_8 31 USE mathconstants, ONLY: z_one 32 USE negf_integr_utils, ONLY: contour_shape_arc,& 33 contour_shape_linear,& 34 equidistant_nodes_a_b,& 35 rescale_nodes_cos,& 36 rescale_normalised_nodes 37#include "./base/base_uses.f90" 38 39 IMPLICIT NONE 40 PRIVATE 41 42 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_integr_cc' 43 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE. 44 45 INTEGER, PARAMETER, PUBLIC :: cc_interval_full = 0, & 46 cc_interval_half = 1 47 48 INTEGER, PARAMETER, PUBLIC :: cc_shape_linear = contour_shape_linear, & 49 cc_shape_arc = contour_shape_arc 50 51 PUBLIC :: ccquad_type 52 53 PUBLIC :: ccquad_init, & 54 ccquad_release, & 55 ccquad_double_number_of_points, & 56 ccquad_reduce_and_append_zdata, & 57 ccquad_refine_integral 58 59! ************************************************************************************************** 60!> \brief Adaptive Clenshaw-Curtis environment. 61! ************************************************************************************************** 62 TYPE ccquad_type 63 !> integration lower and upper bounds 64 COMPLEX(kind=dp) :: a, b 65 !> integration interval: 66 !> cc_interval_full -- [a .. b], 67 !> grid density: 'a' .. . . . . . .. 'b'; 68 !> cc_interval_half -- [a .. 2b-a], assuming int_{b}^{2b-a} f(x) dx = 0, 69 !> grid density: 'a' .. . . . 'b' 70 INTEGER :: interval_id 71 !> integration shape 72 INTEGER :: shape_id 73 !> estimated error 74 REAL(kind=dp) :: error 75 !> approximate integral value 76 TYPE(cp_cfm_type), POINTER :: integral 77 !> error estimate for every element of the 'integral' matrix 78 TYPE(cp_fm_type), POINTER :: error_fm 79 !> weights associated with matrix elements; the 'error' variable contains the value Trace(error_fm * weights) 80 TYPE(cp_fm_type), POINTER :: weights 81 !> integrand value at grid points. Due to symmetry of Clenshaw-Curtis quadratures, 82 !> we only need to keep the left half-interval 83 TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:) :: zdata_cache 84 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes 85 END TYPE ccquad_type 86 87CONTAINS 88 89! ************************************************************************************************** 90!> \brief Initialise a Clenshaw-Curtis quadrature environment variable. 91!> \param cc_env environment variable to initialise 92!> \param xnodes points at which an integrand needs to be computed (initialised on exit) 93!> \param nnodes initial number of points to compute (initialised on exit) 94!> \param a integral lower bound 95!> \param b integral upper bound 96!> \param interval_id full [-1 .. 1] or half [-1 .. 0] interval 97!> \param shape_id shape of a curve along which the integral will be evaluated 98!> \param weights weights associated with matrix elements; used to compute cumulative error 99!> \param tnodes_restart list of nodes over the interval [-1 .. 1] from a previous integral evaluation. 100!> If present, the same set of 'xnodes' will be used to compute this integral. 101!> \par History 102!> * 05.2017 created [Sergey Chulkov] 103!> \note Clenshaw-Curtis quadratures are defined on the interval [-1 .. 1] and have non-uniforms node 104!> distribution which is symmetric and much sparse about 0. When the half-interval [-1 .. 0] 105!> is requested, the integrand value on another subinterval (0 .. 1] is assumed to be zero. 106!> Half interval mode is typically useful for rapidly decaying integrands (e.g. multiplied by 107!> Fermi function), so we do not actually need a fine grid spacing on this tail. 108! ************************************************************************************************** 109 SUBROUTINE ccquad_init(cc_env, xnodes, nnodes, a, b, interval_id, shape_id, weights, tnodes_restart) 110 TYPE(ccquad_type), INTENT(out) :: cc_env 111 INTEGER, INTENT(inout) :: nnodes 112 COMPLEX(kind=dp), DIMENSION(nnodes), INTENT(out) :: xnodes 113 COMPLEX(kind=dp), INTENT(in) :: a, b 114 INTEGER, INTENT(in) :: interval_id, shape_id 115 TYPE(cp_fm_type), POINTER :: weights 116 REAL(kind=dp), DIMENSION(nnodes), INTENT(in), & 117 OPTIONAL :: tnodes_restart 118 119 CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_init', routineP = moduleN//':'//routineN 120 121 INTEGER :: handle, icol, ipoint, irow, ncols, & 122 nnodes_half, nrows 123 REAL(kind=dp), DIMENSION(:, :), POINTER :: w_data, w_data_my 124 TYPE(cp_fm_struct_type), POINTER :: fm_struct 125 126 CALL timeset(routineN, handle) 127 128 CPASSERT(nnodes > 2) 129 CPASSERT(ASSOCIATED(weights)) 130 131 ! ensure that MOD(nnodes-1, 2) == 0 132 nnodes = 2*((nnodes - 1)/2) + 1 133 134 cc_env%interval_id = interval_id 135 cc_env%shape_id = shape_id 136 cc_env%a = a 137 cc_env%b = b 138 cc_env%error = HUGE(0.0_dp) 139 140 NULLIFY (cc_env%integral, cc_env%error_fm, cc_env%weights) 141 CALL cp_fm_get_info(weights, local_data=w_data, nrow_local=nrows, ncol_local=ncols, matrix_struct=fm_struct) 142 CALL cp_fm_create(cc_env%weights, fm_struct) 143 CALL cp_fm_get_info(cc_env%weights, local_data=w_data_my) 144 145 ! use the explicit loop to avoid temporary arrays 146 DO icol = 1, ncols 147 DO irow = 1, nrows 148 w_data_my(irow, icol) = ABS(w_data(irow, icol)) 149 END DO 150 END DO 151 152 SELECT CASE (interval_id) 153 CASE (cc_interval_full) 154 nnodes_half = nnodes/2 + 1 155 CASE (cc_interval_half) 156 nnodes_half = nnodes 157 CASE DEFAULT 158 CPABORT("Unimplemented interval type") 159 END SELECT 160 161 ALLOCATE (cc_env%tnodes(nnodes)) 162 163 IF (PRESENT(tnodes_restart)) THEN 164 cc_env%tnodes(1:nnodes) = tnodes_restart(1:nnodes) 165 ELSE 166 CALL equidistant_nodes_a_b(-1.0_dp, 0.0_dp, nnodes_half, cc_env%tnodes) 167 168 ! rescale all but the end-points, as they are transformed into themselves (-1.0 -> -1.0; 0.0 -> 0.0). 169 ! Moreover, by applying this rescaling transformation to the end-points we cannot guarantee the exact 170 ! result due to rounding errors in evaluation of COS function. 171 IF (nnodes_half > 2) & 172 CALL rescale_nodes_cos(nnodes_half - 2, cc_env%tnodes(2:)) 173 174 SELECT CASE (interval_id) 175 CASE (cc_interval_full) 176 ! reflect symmetric nodes 177 DO ipoint = nnodes_half - 1, 1, -1 178 cc_env%tnodes(nnodes_half + ipoint) = -cc_env%tnodes(nnodes_half - ipoint) 179 END DO 180 CASE (cc_interval_half) 181 ! rescale half-interval : [-1 .. 0] -> [-1 .. 1] 182 cc_env%tnodes(1:nnodes_half) = 2.0_dp*cc_env%tnodes(1:nnodes_half) + 1.0_dp 183 END SELECT 184 END IF 185 186 CALL rescale_normalised_nodes(nnodes, cc_env%tnodes, a, b, shape_id, xnodes) 187 188 CALL timestop(handle) 189 END SUBROUTINE ccquad_init 190 191! ************************************************************************************************** 192!> \brief Release a Clenshaw-Curtis quadrature environment variable. 193!> \param cc_env environment variable to release (modified on exit) 194!> \par History 195!> * 05.2017 created [Sergey Chulkov] 196! ************************************************************************************************** 197 SUBROUTINE ccquad_release(cc_env) 198 TYPE(ccquad_type), INTENT(inout) :: cc_env 199 200 CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_release', routineP = moduleN//':'//routineN 201 202 INTEGER :: handle, ipoint 203 204 CALL timeset(routineN, handle) 205 206 IF (ASSOCIATED(cc_env%error_fm)) & 207 CALL cp_fm_release(cc_env%error_fm) 208 209 IF (ASSOCIATED(cc_env%weights)) & 210 CALL cp_fm_release(cc_env%weights) 211 212 IF (ASSOCIATED(cc_env%integral)) & 213 CALL cp_cfm_release(cc_env%integral) 214 215 IF (ALLOCATED(cc_env%zdata_cache)) THEN 216 DO ipoint = SIZE(cc_env%zdata_cache), 1, -1 217 IF (ASSOCIATED(cc_env%zdata_cache(ipoint)%matrix)) & 218 CALL cp_cfm_release(cc_env%zdata_cache(ipoint)%matrix) 219 END DO 220 221 DEALLOCATE (cc_env%zdata_cache) 222 END IF 223 224 IF (ALLOCATED(cc_env%tnodes)) DEALLOCATE (cc_env%tnodes) 225 226 CALL timestop(handle) 227 END SUBROUTINE ccquad_release 228 229! ************************************************************************************************** 230!> \brief Get the next set of points at which the integrand needs to be computed. These points are 231!> then can be used to refine the integral approximation. 232!> \param cc_env environment variable (modified on exit) 233!> \param xnodes_next set of additional points (allocated and initialised on exit) 234!> \par History 235!> * 05.2017 created [Sergey Chulkov] 236! ************************************************************************************************** 237 SUBROUTINE ccquad_double_number_of_points(cc_env, xnodes_next) 238 TYPE(ccquad_type), INTENT(inout) :: cc_env 239 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:), & 240 INTENT(inout) :: xnodes_next 241 242 CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_double_number_of_points', & 243 routineP = moduleN//':'//routineN 244 245 INTEGER :: handle, ipoint, nnodes_exist, & 246 nnodes_half, nnodes_next 247 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes, tnodes_old 248 249 CALL timeset(routineN, handle) 250 251 CPASSERT(.NOT. ALLOCATED(xnodes_next)) 252 CPASSERT(ASSOCIATED(cc_env%integral)) 253 CPASSERT(ASSOCIATED(cc_env%error_fm)) 254 CPASSERT(ALLOCATED(cc_env%zdata_cache)) 255 256 ! due to symmetry of Clenshaw-Curtis quadratures, we only need to keep the left half-interval [-1 .. 0] 257 nnodes_exist = SIZE(cc_env%zdata_cache) 258 ! new nodes will be placed between the existed ones, so the number of nodes 259 ! on the left half-interval [-1 .. 0] is equal to nnodes_exist - 1 260 nnodes_half = nnodes_exist - 1 261 262 SELECT CASE (cc_env%interval_id) 263 CASE (cc_interval_full) 264 ! double number of nodes as we have 2 half-intervals [-1 .. 0] and [0 .. 1] 265 nnodes_next = 2*nnodes_half 266 CASE (cc_interval_half) 267 nnodes_next = nnodes_half 268 CASE DEFAULT 269 CPABORT("Unimplemented interval type") 270 END SELECT 271 272 ALLOCATE (xnodes_next(nnodes_next)) 273 ALLOCATE (tnodes(nnodes_next)) 274 275 CALL equidistant_nodes_a_b(0.5_dp/REAL(nnodes_half, kind=dp) - 1.0_dp, & 276 -0.5_dp/REAL(nnodes_half, kind=dp), & 277 nnodes_half, tnodes) 278 279 CALL rescale_nodes_cos(nnodes_half, tnodes) 280 281 SELECT CASE (cc_env%interval_id) 282 CASE (cc_interval_full) 283 ! reflect symmetric nodes 284 DO ipoint = 1, nnodes_half 285 tnodes(nnodes_half + ipoint) = -tnodes(nnodes_half - ipoint + 1) 286 END DO 287 CASE (cc_interval_half) 288 ! rescale half-interval : [-1 .. 0] -> [-1 .. 1] 289 tnodes(1:nnodes_half) = 2.0_dp*tnodes(1:nnodes_half) + 1.0_dp 290 END SELECT 291 292 ! append new tnodes to the cache 293 CALL MOVE_ALLOC(cc_env%tnodes, tnodes_old) 294 nnodes_exist = SIZE(tnodes_old) 295 296 ALLOCATE (cc_env%tnodes(nnodes_exist + nnodes_next)) 297 cc_env%tnodes(1:nnodes_exist) = tnodes_old(1:nnodes_exist) 298 cc_env%tnodes(nnodes_exist + 1:nnodes_exist + nnodes_next) = tnodes(1:nnodes_next) 299 DEALLOCATE (tnodes_old) 300 301 ! rescale nodes [-1 .. 1] -> [a .. b] according to the shape 302 CALL rescale_normalised_nodes(nnodes_next, tnodes, cc_env%a, cc_env%b, cc_env%shape_id, xnodes_next) 303 304 DEALLOCATE (tnodes) 305 CALL timestop(handle) 306 END SUBROUTINE ccquad_double_number_of_points 307 308! ************************************************************************************************** 309!> \brief Prepare Clenshaw-Curtis environment for the subsequent refinement of the integral. 310!> \param cc_env environment variable (modified on exit) 311!> \param zdata_next additional integrand value at additional points (modified on exit) 312!> \par History 313!> * 05.2017 created [Sergey Chulkov] 314!> \note Due to symmetry of Clenshaw-Curtis quadratures (weight(x) == weight(-x)), we do not need to 315!> keep all the matrices from 'zdata_next', only 'zdata_next(x) + zdata_next(-x)' is needed. 316!> In order to reduce the number of matrix allocations, we move some of the matrices from the 317!> end of the 'zdata_new' array to the 'cc_env%zdata_cache' array, and nullify the corresponding 318!> pointers at 'zdata_next' array. So the calling subroutine need to release the remained 319!> matrices or reuse them but taking into account the missed ones. 320! ************************************************************************************************** 321 SUBROUTINE ccquad_reduce_and_append_zdata(cc_env, zdata_next) 322 TYPE(ccquad_type), INTENT(inout) :: cc_env 323 TYPE(cp_cfm_p_type), DIMENSION(:), INTENT(inout) :: zdata_next 324 325 CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_reduce_and_append_zdata', & 326 routineP = moduleN//':'//routineN 327 328 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: zscale 329 INTEGER :: handle, ipoint, nnodes_exist, & 330 nnodes_half, nnodes_next 331 TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:) :: zdata_tmp 332 333 CALL timeset(routineN, handle) 334 335 nnodes_next = SIZE(zdata_next) 336 CPASSERT(nnodes_next > 0) 337 338 ! compute weights of new points on a complex contour according to their values of the 't' parameter 339 nnodes_exist = SIZE(cc_env%tnodes) 340 CPASSERT(nnodes_exist >= nnodes_next) 341 342 ALLOCATE (zscale(nnodes_next)) 343 CALL rescale_normalised_nodes(nnodes_next, cc_env%tnodes(nnodes_exist - nnodes_next + 1:nnodes_exist), & 344 cc_env%a, cc_env%b, cc_env%shape_id, weights=zscale) 345 346 IF (cc_env%interval_id == cc_interval_half) zscale(:) = 2.0_dp*zscale(:) 347 348 ! rescale integrand values 349 DO ipoint = 1, nnodes_next 350 CALL cp_cfm_scale(zscale(ipoint), zdata_next(ipoint)%matrix) 351 END DO 352 DEALLOCATE (zscale) 353 354 ! squash points with the same clenshaw-curtis weights together 355 IF (ALLOCATED(cc_env%zdata_cache)) THEN 356 nnodes_exist = SIZE(cc_env%zdata_cache) 357 ELSE 358 nnodes_exist = 0 359 END IF 360 361 SELECT CASE (cc_env%interval_id) 362 CASE (cc_interval_full) 363 IF (ALLOCATED(cc_env%zdata_cache)) THEN 364 CPASSERT(nnodes_exist == nnodes_next/2 + 1) 365 nnodes_half = nnodes_exist - 1 366 ELSE 367 CPASSERT(MOD(nnodes_next, 2) == 1) 368 nnodes_half = nnodes_next/2 + 1 369 END IF 370 CASE (cc_interval_half) 371 IF (ALLOCATED(cc_env%zdata_cache)) THEN 372 CPASSERT(nnodes_exist == nnodes_next + 1) 373 END IF 374 375 nnodes_half = nnodes_next 376 END SELECT 377 378 IF (cc_env%interval_id == cc_interval_full) THEN 379 DO ipoint = nnodes_next/2, 1, -1 380 CALL cp_cfm_scale_and_add(z_one, zdata_next(ipoint)%matrix, z_one, zdata_next(nnodes_next - ipoint + 1)%matrix) 381 END DO 382 END IF 383 384 IF (ALLOCATED(cc_env%zdata_cache)) THEN 385 ! note that nnodes_half+1 == nnodes_exist for both half- and full-intervals 386 ALLOCATE (zdata_tmp(nnodes_half + nnodes_exist)) 387 388 DO ipoint = 1, nnodes_half 389 zdata_tmp(2*ipoint - 1)%matrix => cc_env%zdata_cache(ipoint)%matrix 390 zdata_tmp(2*ipoint)%matrix => zdata_next(ipoint)%matrix 391 NULLIFY (zdata_next(ipoint)%matrix) 392 END DO 393 zdata_tmp(nnodes_half + nnodes_exist)%matrix => cc_env%zdata_cache(nnodes_exist)%matrix 394 395 DEALLOCATE (cc_env%zdata_cache) 396 CALL MOVE_ALLOC(zdata_tmp, cc_env%zdata_cache) 397 ELSE 398 CALL cp_cfm_scale(2.0_dp, zdata_next(nnodes_half)%matrix) 399 400 ALLOCATE (cc_env%zdata_cache(nnodes_half)) 401 402 DO ipoint = 1, nnodes_half 403 cc_env%zdata_cache(ipoint)%matrix => zdata_next(ipoint)%matrix 404 NULLIFY (zdata_next(ipoint)%matrix) 405 END DO 406 END IF 407 408 CALL timestop(handle) 409 END SUBROUTINE ccquad_reduce_and_append_zdata 410 411! ************************************************************************************************** 412!> \brief Refine approximated integral. 413!> \param cc_env environment variable (modified on exit) 414!> \par History 415!> * 05.2017 created [Sergey Chulkov] 416! ************************************************************************************************** 417 SUBROUTINE ccquad_refine_integral(cc_env) 418 TYPE(ccquad_type), INTENT(inout) :: cc_env 419 420 CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_refine_integral', & 421 routineP = moduleN//':'//routineN 422 423 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ztmp, ztmp_dct 424 INTEGER :: handle, icol, ipoint, irow, ncols_local, nintervals, nintervals_half, & 425 nintervals_half_plus_1, nintervals_half_plus_2, nintervals_plus_2, nrows_local, stat 426 LOGICAL :: equiv 427 REAL(kind=dp) :: rscale 428 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights 429 TYPE(cp_fm_struct_type), POINTER :: fm_struct 430 431! TYPE(fft_plan_type) :: fft_plan 432! INTEGER(kind=int_8) :: plan 433 434 CALL timeset(routineN, handle) 435 436 CPASSERT(ALLOCATED(cc_env%zdata_cache)) 437 438 nintervals_half_plus_1 = SIZE(cc_env%zdata_cache) 439 nintervals_half = nintervals_half_plus_1 - 1 440 nintervals_half_plus_2 = nintervals_half_plus_1 + 1 441 nintervals = 2*nintervals_half 442 nintervals_plus_2 = nintervals + 2 443 CPASSERT(nintervals_half > 1) 444 445 IF (.NOT. ASSOCIATED(cc_env%integral)) THEN 446 CALL cp_cfm_get_info(cc_env%zdata_cache(1)%matrix, matrix_struct=fm_struct) 447 equiv = cp_fm_struct_equivalent(fm_struct, cc_env%weights%matrix_struct) 448 CPASSERT(equiv) 449 450 CALL cp_cfm_create(cc_env%integral, fm_struct) 451 CALL cp_fm_create(cc_env%error_fm, fm_struct) 452 END IF 453 454 IF (debug_this_module) THEN 455 DO ipoint = 1, nintervals_half_plus_1 456 equiv = cp_fm_struct_equivalent(cc_env%zdata_cache(ipoint)%matrix%matrix_struct, cc_env%integral%matrix_struct) 457 CPASSERT(equiv) 458 END DO 459 END IF 460 461 CALL cp_cfm_get_info(cc_env%integral, nrow_local=nrows_local, ncol_local=ncols_local) 462 463 ALLOCATE (weights(nintervals_half)) 464 465 ! omit the trivial weights(1) = 0.5 466 DO ipoint = 2, nintervals_half 467 rscale = REAL(2*(ipoint - 1), kind=dp) 468 weights(ipoint) = 1.0_dp/(1.0_dp - rscale*rscale) 469 END DO 470 ! weights(1) <- weights(intervals_half + 1) 471 rscale = REAL(nintervals, kind=dp) 472 weights(1) = 1.0_dp/(1.0_dp - rscale*rscale) 473 474 ! 1.0 / nintervals 475 rscale = 1.0_dp/rscale 476 477 ALLOCATE (ztmp(nintervals, nrows_local, ncols_local)) 478 ALLOCATE (ztmp_dct(nintervals, nrows_local, ncols_local)) 479 480!$OMP PARALLEL DO DEFAULT(NONE), PRIVATE(icol, ipoint, irow), & 481!$OMP SHARED(cc_env, ncols_local, nintervals_half, nintervals_half_plus_1, nintervals_half_plus_2, nrows_local, ztmp) 482 DO icol = 1, ncols_local 483 DO irow = 1, nrows_local 484 DO ipoint = 1, nintervals_half_plus_1 485 ztmp(ipoint, irow, icol) = cc_env%zdata_cache(ipoint)%matrix%local_data(irow, icol) 486 END DO 487 488 DO ipoint = 2, nintervals_half 489 ztmp(nintervals_half + ipoint, irow, icol) = ztmp(nintervals_half_plus_2 - ipoint, irow, icol) 490 END DO 491 END DO 492 END DO 493!$OMP END PARALLEL DO 494 495 CALL fft_fw1d(nintervals, nrows_local*ncols_local, .FALSE., ztmp, ztmp_dct, 1.0_dp, stat) 496 IF (stat /= 0) THEN 497 CALL cp_abort(__LOCATION__, & 498 "An FFT library is required for Clenshaw-Curtis quadrature. "// & 499 "You can use an alternative integration method instead.") 500 END IF 501 502!$OMP PARALLEL DO DEFAULT(NONE), PRIVATE(icol, ipoint, irow), & 503!$OMP SHARED(cc_env, rscale, ncols_local, nintervals_half, nintervals_half_plus_1, nintervals_plus_2), & 504!$OMP SHARED(nrows_local, weights, ztmp_dct) 505 DO icol = 1, ncols_local 506 DO irow = 1, nrows_local 507 ztmp_dct(1, irow, icol) = 0.5_dp*ztmp_dct(1, irow, icol) 508 DO ipoint = 2, nintervals_half 509 ztmp_dct(ipoint, irow, icol) = 0.5_dp*weights(ipoint)*(ztmp_dct(ipoint, irow, icol) + & 510 ztmp_dct(nintervals_plus_2 - ipoint, irow, icol)) 511 END DO 512 ztmp_dct(nintervals_half_plus_1, irow, icol) = weights(1)*ztmp_dct(nintervals_half_plus_1, irow, icol) 513 514 cc_env%integral%local_data(irow, icol) = rscale*accurate_sum(ztmp_dct(1:nintervals_half_plus_1, irow, icol)) 515 cc_env%error_fm%local_data(irow, icol) = rscale*ABS(ztmp_dct(nintervals_half_plus_1, irow, icol)) 516 END DO 517 END DO 518!$OMP END PARALLEL DO 519 520 DEALLOCATE (ztmp, ztmp_dct) 521 522 CALL cp_fm_trace(cc_env%error_fm, cc_env%weights, cc_env%error) 523 524 DEALLOCATE (weights) 525 CALL timestop(handle) 526 END SUBROUTINE ccquad_refine_integral 527 528END MODULE negf_integr_cc 529