ddx 0.6.8
Libary for domain-decomposition methods for polarizable continuum models
ddx_constants.f90
Go to the documentation of this file.
1
15
17module ddx_constants
18! Get ddx_params_type and all compile-time definitions
19use ddx_parameters
20! Get harmonics-related functions
21use ddx_harmonics
22use omp_lib, only : omp_get_wtime
23
24implicit none
25
26
28type ddx_constants_type
30 integer :: nbasis
33 integer :: n
40 integer :: dmax
42 integer :: nscales
45 real(dp), allocatable :: vscales(:)
48 real(dp), allocatable :: v4pi2lp1(:)
53 real(dp), allocatable :: vscales_rel(:)
57 integer :: nfact
59 real(dp), allocatable :: vfact(:)
63 real(dp), allocatable :: vcnk(:)
67 real(dp), allocatable :: m2l_ztranslate_coef(:, :, :)
71 real(dp), allocatable :: m2l_ztranslate_adj_coef(:, :, :)
73 real(dp), allocatable :: cgrid(:, :)
75 real(dp), allocatable :: wgrid(:)
81 integer :: vgrid_dmax
83 integer :: vgrid_nbasis
86 real(dp), allocatable :: vgrid(:, :)
90 real(dp), allocatable :: vwgrid(:, :)
92 real(dp), allocatable :: vgrid2(:, :)
94 integer :: lmax0
96 integer :: nbasis0
99 real(dp), allocatable :: pchi(:, :, :)
102 real(dp), allocatable :: c_ik(:, :)
104 real(dp), allocatable :: si_ri(:, :)
107 real(dp), allocatable :: di_ri(:, :)
109 real(dp), allocatable :: sk_ri(:, :)
112 real(dp), allocatable :: dk_ri(:, :)
114 real(dp), allocatable :: termimat(:, :)
116 real(dp), allocatable :: b(:,:,:)
118 real(dp), allocatable :: l(:,:,:)
121 integer :: nngmax
123 integer, allocatable :: inl(:)
126 integer, allocatable :: nl(:)
128 integer, allocatable :: itrnl(:)
131 real(dp), allocatable :: fi(:, :)
134 real(dp), allocatable :: ui(:, :)
137 real(dp), allocatable :: ui_cav(:)
140 real(dp), allocatable :: zi(:, :, :)
143 real(dp), allocatable :: zi_dr(:, :)
145 integer :: ncav
147 integer, allocatable :: ncav_sph(:)
149 real(dp), allocatable :: ccav(:, :)
151 integer, allocatable :: icav_ia(:)
153 integer, allocatable :: icav_ja(:)
156 real(dp), allocatable :: rx_prc(:, :, :)
157 !! Cluster tree information that is allocated and computed only if fmm=1.
160 integer, allocatable :: order(:)
162 integer :: nclusters
165 integer, allocatable :: cluster(:, :)
168 integer, allocatable :: children(:, :)
171 integer, allocatable :: parent(:)
174 real(dp), allocatable :: cnode(:, :)
177 real(dp), allocatable :: rnode(:)
180 integer, allocatable :: snode(:)
183 real(dp), allocatable :: sk_rnode(:, :)
186 real(dp), allocatable :: si_rnode(:, :)
188 integer :: nnfar
190 integer :: nnnear
193 integer, allocatable :: nfar(:)
196 integer, allocatable :: nnear(:)
199 integer, allocatable :: far(:)
202 integer, allocatable :: near(:)
206 integer, allocatable :: sfar(:)
210 integer, allocatable :: snear(:)
213 integer :: nnear_m2p
216 integer :: m2p_lmax
219 integer :: m2p_nbasis
223 integer :: grad_nbasis
225 real(dp) :: inner_tol
228 logical :: dodiag
229end type ddx_constants_type
230
231contains
232
239subroutine constants_init(params, constants, ddx_error)
240 use complex_bessel
241 !! Inputs
242 type(ddx_params_type), intent(in) :: params
243 !! Outputs
244 type(ddx_constants_type), intent(out) :: constants
245 type(ddx_error_type), intent(inout) :: ddx_error
246 !! Local variables
247 integer :: i, alloc_size, l, indl, igrid, isph, l0, &
248 & NZ, ierr, info, tmp_pmax
249 real(dp) :: rho, ctheta, stheta, cphi, sphi, termi, termk, &
250 & s1, s2
251 real(dp), allocatable :: vplm(:), vcos(:), vsin(:), vylm(:), SK_rijn(:), &
252 & DK_rijn(:)
253 complex(dp), allocatable :: bessel_work(:)
254 complex(dp) :: z
255 !! The code
256 ! Clean ddx_error state of constants to proceed with geometry
257 if (ddx_error % flag .ne. 0) then
258 call update_error(ddx_error, "constants_init received input in error " // &
259 & " state, exiting")
260 return
261 end if
262 ! activate inner iterations diagonal in mvp for debugging purposes only.
263 ! could be useful for different linear solvers.
264 constants % dodiag = .false.
265 ! Maximal number of modeling spherical harmonics
266 constants % nbasis = (params % lmax+1) ** 2
267 ! Maximal number of modeling degrees of freedom
268 constants % n = params % nsph * constants % nbasis
269 ! Calculate dmax, vgrid_dmax, m2p_lmax, m2p_nbasis and grad_nbasis
270 if (params % fmm .eq. 0) then
271 if (params % force .eq. 1) then
272 constants % dmax = params % lmax + 1
273 constants % vgrid_dmax = params % lmax + 1
274 else
275 constants % dmax = params % lmax
276 constants % vgrid_dmax = params % lmax
277 end if
278 ! Other constants are not referenced if fmm=0
279 constants % m2p_lmax = -1
280 constants % m2p_nbasis = -1
281 constants % grad_nbasis = -1
282 else
283 ! If forces are required then we need the M2P of a degree lmax+1 for
284 ! the near-field analytical gradients
285 if (params % force .eq. 1) then
286 constants % dmax = max(params % pm+params % pl+1, &
287 & params % lmax+1)
288 constants % m2p_lmax = params % lmax + 1
289 constants % grad_nbasis = (params % lmax+2) ** 2
290 constants % vgrid_dmax = max(params % pl, params % lmax) + 1
291 else
292 constants % dmax = max(params % pm+params % pl, &
293 & params % lmax)
294 constants % m2p_lmax = params % lmax
295 constants % grad_nbasis = -1
296 constants % vgrid_dmax = max(params % pl, params % lmax)
297 end if
298 constants % m2p_nbasis = (constants % m2p_lmax+1) ** 2
299 end if
300 ! Compute sizes of vgrid, vfact and vscales
301 constants % vgrid_nbasis = (constants % vgrid_dmax+1) ** 2
302 constants % nfact = 2*constants % dmax+1
303 constants % nscales = (constants % dmax+1) ** 2
304 ! Allocate space for scaling factors of spherical harmonics
305 allocate(constants % vscales(constants % nscales), stat=info)
306 if (info .ne. 0) then
307 call update_error(ddx_error, "constants_init: `vscales` allocation " &
308 & // "failed")
309 return
310 end if
311 allocate(constants % v4pi2lp1(constants % dmax+1), stat=info)
312 if (info .ne. 0) then
313 call update_error(ddx_error, "constants_init: `v4pi2lp1` allocation " &
314 & // "failed")
315 return
316 end if
317 allocate(constants % vscales_rel(constants % nscales), stat=info)
318 if (info .ne. 0) then
319 call update_error(ddx_error, "constants_init: `vscales_rel` " // &
320 & "allocation failed")
321 return
322 end if
323 ! Compute scaling factors of spherical harmonics
324 call ylmscale(constants % dmax, constants % vscales, &
325 & constants % v4pi2lp1, constants % vscales_rel)
326 ! Allocate square roots of factorials
327 allocate(constants % vfact(constants % nfact), stat=info)
328 if (info .ne. 0) then
329 call update_error(ddx_error, "constants_init: `vfact` allocation failed")
330 return
331 end if
332 ! Compute square roots of factorials
333 constants % vfact(1) = 1
334 do i = 2, constants % nfact
335 constants % vfact(i) = constants % vfact(i-1) * sqrt(dble(i-1))
336 end do
337 ! Allocate square roots of combinatorial numbers C_n^k and M2L OZ
338 ! translation coefficients
339 if (params % fmm .eq. 1 .or. &
340 & (params % model .eq. 3 .and. params % force .eq. 1)) then
341 alloc_size = 2*constants % dmax + 1
342 alloc_size = alloc_size * (constants % dmax+1)
343 ! Allocate C_n^k
344 allocate(constants % vcnk(alloc_size), stat=info)
345 if (info .ne. 0) then
346 call update_error(ddx_error, "constants_init: `vcnk` allocation " &
347 & // "failed")
348 return
349 end if
350 ! Allocate M2L OZ coefficients
351 allocate(constants % m2l_ztranslate_coef(&
352 & (params % pm+1), (params % pl+1), (params % pl+1)), &
353 & stat=info)
354 if (info .ne. 0) then
355 call update_error(ddx_error, "constants_init: " // &
356 & "`m2l_ztranslate_coef` allocation failed")
357 return
358 end if
359 ! Allocate adjoint M2L OZ coefficients
360 allocate(constants % m2l_ztranslate_adj_coef(&
361 & (params % pl+1), (params % pl+1), (params % pm+1)), &
362 & stat=info)
363 if (info .ne. 0) then
364 call update_error(ddx_error, "constants_init: " // &
365 & "`m2l_ztranslate_adj_coef` allocation failed")
366 return
367 end if
368 ! Compute combinatorial numbers C_n^k and M2L OZ translate coefficients
369 call fmm_constants(constants % dmax, params % pm, &
370 & params % pl, constants % vcnk, &
371 & constants % m2l_ztranslate_coef, &
372 & constants % m2l_ztranslate_adj_coef)
373 end if
374 ! Allocate space for Lebedev grid coordinates and weights
375 allocate(constants % cgrid(3, params % ngrid), &
376 & constants % wgrid(params % ngrid), stat=info)
377 if (info .ne. 0) then
378 call update_error(ddx_error, "constants_init: `cgrid` and `wgrid` " // &
379 & "allocations failed")
380 return
381 end if
382 ! Get weights and coordinates of Lebedev grid points
383 call llgrid(params % ngrid, constants % wgrid, constants % cgrid)
384 ! Allocate space for values of non-weighted and weighted spherical
385 ! harmonics and L2P at Lebedev grid points
386 allocate(constants % vgrid(constants % vgrid_nbasis, params % ngrid), &
387 & constants % vwgrid(constants % vgrid_nbasis, params % ngrid), &
388 & stat=info)
389 if (info .ne. 0) then
390 call update_error(ddx_error, "constants_init: `vgrid`, `wgrid` and" // &
391 & " allocations failed")
392 return
393 end if
394 allocate(vplm(constants % vgrid_nbasis), vcos(constants % vgrid_dmax+1), &
395 & vsin(constants % vgrid_dmax+1), stat=info)
396 if (info .ne. 0) then
397 call update_error(ddx_error, "constants_init: `vplm`, `vcos` and " // &
398 & "`vsin` allocations failed")
399 return
400 end if
401 ! Compute non-weighted and weighted spherical harmonics and the single
402 ! layer operator at Lebedev grid points
403 do igrid = 1, params % ngrid
404 call ylmbas(constants % cgrid(:, igrid), rho, ctheta, stheta, cphi, &
405 & sphi, constants % vgrid_dmax, constants % vscales, &
406 & constants % vgrid(:, igrid), vplm, vcos, vsin)
407 constants % vwgrid(:, igrid) = constants % wgrid(igrid) * &
408 & constants % vgrid(:, igrid)
409 end do
410 if (params % fmm .eq. 1) then
411 allocate(constants % vgrid2( &
412 & constants % vgrid_nbasis, params % ngrid), &
413 & stat=info)
414 if (info .ne. 0) then
415 call update_error(ddx_error, "constants_init: `vgrid2` " // &
416 & "allocation failed")
417 return
418 end if
419 do igrid = 1, params % ngrid
420 do l = 0, constants % vgrid_dmax
421 indl = l*l + l + 1
422 constants % vgrid2(indl-l:indl+l, igrid) = &
423 & constants % vgrid(indl-l:indl+l, igrid) / &
424 & constants % vscales(indl)**2
425 end do
426 end do
427 end if
428
429 ! Generate geometry-related constants (required by the LPB code)
430 call constants_geometry_init(params, constants, ddx_error)
431 if (ddx_error % flag .ne. 0) then
432 call update_error(ddx_error, "constants_geometry_init returned an " // &
433 & "error, exiting")
434 return
435 end if
436
437 ! Precompute LPB-related constants
438 if (params % model .eq. 3) then
439 constants % lmax0 = params % lmax
440 constants % nbasis0 = constants % nbasis
441 allocate(vylm(constants % vgrid_nbasis), &
442 & sk_rijn(0:constants % lmax0), dk_rijn(0:constants % lmax0), &
443 & bessel_work(constants % dmax+2), stat=info)
444 if (info .ne. 0) then
445 call update_error(ddx_error, "constants_init: `vylm`, `SK_rijn` " // &
446 & "and `DK_rijn` allocations failed")
447 return
448 end if
449 allocate(constants % SI_ri(0:constants % dmax+1, params % nsph), &
450 & constants % DI_ri(0:constants % dmax+1, params % nsph), &
451 & constants % SK_ri(0:params % lmax+1, params % nsph), &
452 & constants % DK_ri(0:params % lmax+1, params % nsph), &
453 & constants % Pchi(constants % nbasis, constants % nbasis0, params % nsph), &
454 & constants % C_ik(0:params % lmax, params % nsph), &
455 & constants % termimat(0:params % lmax, params % nsph), &
456 & stat=info)
457 if (info .ne. 0) then
458 call update_error(ddx_error, "constants_init: allocation failed")
459 return
460 end if
461 sk_rijn = zero
462 dk_rijn = zero
463 do isph = 1, params % nsph
464 ! We compute Bessel functions of degrees 0..lmax+1 because the
465 ! largest degree is required for forces
466 call modified_spherical_bessel_first_kind(constants % dmax+1, &
467 & params % rsph(isph)*params % kappa,&
468 & constants % SI_ri(:, isph), constants % DI_ri(:, isph), &
469 & bessel_work)
470 call modified_spherical_bessel_second_kind(params % lmax+1, &
471 & params % rsph(isph)*params % kappa, &
472 & constants % SK_ri(:, isph), constants % DK_ri(:, isph), &
473 & bessel_work)
474 ! Compute matrix PU_i^e(x_in)
475 ! Previous implementation in update_rhs. Made it in allocate_model, so as to use
476 ! it in Forces as well.
477 call mkpmat(params, constants, isph, constants % Pchi(:, :, isph))
478 ! Compute i'_l(r_i)/i_l(r_i)
479 do l = 0, params % lmax
480 constants % termimat(l, isph) = constants % DI_ri(l, isph) / &
481 & constants % SI_ri(l, isph) * params % kappa
482 end do
483 ! Compute (i'_l0/i_l0 - k'_l0/k_l0)^(-1) is computed in Eq.(97)
484 do l0 = 0, constants % lmax0
485 termi = constants % DI_ri(l0, isph) / &
486 & constants % SI_ri(l0, isph) * params % kappa
487 termk = constants % DK_ri(l0, isph)/ &
488 & constants % SK_ri(l0, isph) * params % kappa
489 constants % C_ik(l0, isph) = one / (termi-termk)
490 end do
491 end do
492 if (params % fmm .eq. 1) then
493 tmp_pmax = max(params % pm, params % pl)
494 allocate(constants % SI_rnode(tmp_pmax+1, constants % nclusters), &
495 & constants % SK_rnode(tmp_pmax+1, constants % nclusters), &
496 & stat=info)
497 if (info.ne.0) then
498 call update_error(ddx_error, "constants_init: `si_rnode, " // &
499 & "sk_rnode allocation failed")
500 return
501 end if
502 do i = 1, constants % nclusters
503 z = constants % rnode(i) * params % kappa
504 s1 = sqrt(two/(pi*real(z)))
505 s2 = sqrt(pi/(two*real(z)))
506 call cbesk(z, pt5, 1, 1, bessel_work(1), nz, ierr)
507 constants % SK_rnode(1, i) = s1 * real(bessel_work(1))
508 call cbesi(z, pt5, 1, 1, bessel_work(1), nz, ierr)
509 constants % SI_rnode(1, i) = s2 * real(bessel_work(1))
510 if (tmp_pmax .gt. 0) then
511 call cbesk(z, 1.5d0, 1, tmp_pmax, bessel_work(2:), nz, ierr)
512 constants % SK_rnode(2:, i) = s1 * real(bessel_work(2:tmp_pmax+1))
513 call cbesi(z, 1.5d0, 1, tmp_pmax, bessel_work(2:), nz, ierr)
514 constants % SI_rnode(2:, i) = s2 * real(bessel_work(2:tmp_pmax+1))
515 end if
516 end do
517 end if
518 deallocate(vylm, sk_rijn, dk_rijn, bessel_work, stat=info)
519 if (info .ne. 0) then
520 call update_error(ddx_error, "constants_init: `vylm`, `SK_rijn` " // &
521 & "and `DK_rijn` deallocations failed")
522 return
523 end if
524 ! if doing incore build nonzero blocks of B
525 if (params % matvecmem .eq. 1) then
526 call build_b(constants, params, ddx_error)
527 if (ddx_error % flag .ne. 0) then
528 call update_error(ddx_error, "build_b returned an " // &
529 & "error, exiting")
530 return
531 end if
532 end if
533 end if
534 deallocate(vplm, vcos, vsin, stat=info)
535 if (info .ne. 0) then
536 call update_error(ddx_error, "constants_init: `vplm`, `vcos` and " // &
537 & "`vsin` deallocations failed")
538 return
539 end if
540 ! if doing incore build nonzero blocks of L
541 if (params % matvecmem .eq. 1) then
542 call build_itrnl(constants, params, ddx_error)
543 if (ddx_error % flag .ne. 0) then
544 call update_error(ddx_error, "build_itrnl returned an " // &
545 & "error, exiting")
546 return
547 end if
548 call build_l(constants, params, ddx_error)
549 if (ddx_error % flag .ne. 0) then
550 call update_error(ddx_error, "build_l returned an " // &
551 & "error, exiting")
552 return
553 end if
554 end if
555end subroutine constants_init
556
558subroutine build_itrnl(constants, params, ddx_error)
559 implicit none
560 type(ddx_params_type), intent(in) :: params
561 type(ddx_constants_type), intent(inout) :: constants
562 type(ddx_error_type), intent(inout) :: ddx_error
563 integer :: isph, ij, jsph, ji, istat
564
565 allocate(constants % itrnl(constants % inl(params % nsph + 1)), stat=istat)
566 if (istat.ne.0) then
567 call update_error(ddx_error, 'Allocation failed in build_itrnl')
568 return
569 end if
570
571 do isph = 1, params % nsph
572 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
573 jsph = constants % nl(ij)
574 do ji = constants % inl(jsph), constants % inl(jsph + 1) - 1
575 if (constants % nl(ji) .eq. isph) exit
576 end do
577 constants % itrnl(ij) = ji
578 end do
579 end do
580end subroutine build_itrnl
581
583subroutine build_l(constants, params, ddx_error)
584 implicit none
585 type(ddx_params_type), intent(in) :: params
586 type(ddx_constants_type), intent(inout) :: constants
587 type(ddx_error_type), intent(inout) :: ddx_error
588 integer :: isph, ij, jsph, igrid, l, m, ind, info
589 real(dp), dimension(3) :: vij, sij
590 real(dp) :: vvij, tij, xij, oij, rho, ctheta, stheta, cphi, sphi, &
591 & fac, tt, thigh
592 real(dp), dimension(constants % nbasis) :: vylm, vplm
593 real(dp), dimension(params % lmax + 1) :: vcos, vsin
594 real(dp), dimension(constants % nbasis, params % ngrid) :: scratch
595 real(dp) :: t
596
597 allocate(constants % l(constants % nbasis, constants % nbasis, &
598 & constants % inl(params % nsph + 1)), stat=info)
599 if (info .ne. 0) then
600 call update_error(ddx_error, 'Allocation failed in build_l')
601 return
602 end if
603
604 thigh = one + pt5*(params % se + one)*params % eta
605
606 t = omp_get_wtime()
607 !$omp parallel do default(none) shared(params,constants,thigh) &
608 !$omp private(isph,ij,jsph,scratch,igrid,vij,vvij,tij,sij,xij,oij, &
609 !$omp rho,ctheta,stheta,cphi,sphi,vylm,vplm,vcos,vsin,l,fac,ind,m,tt)
610 do isph = 1, params % nsph
611 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
612 jsph = constants % nl(ij)
613 scratch = zero
614 do igrid = 1, params % ngrid
615 if (constants % ui(igrid, isph).eq.one) cycle
616 vij = params % csph(:, isph) &
617 & + params % rsph(isph)*constants % cgrid(:,igrid) &
618 & - params % csph(:, jsph)
619 vvij = sqrt(dot_product(vij, vij))
620 tij = vvij/params % rsph(jsph)
621 if (tij.lt.thigh) then
622 if (tij.ne.zero) then
623 sij = vij/vvij
624 else
625 sij = one
626 end if
627 xij = fsw(tij, params % se, params % eta)
628 if (constants % fi(igrid, isph).gt.one) then
629 oij = xij/constants % fi(igrid, isph)
630 else
631 oij = xij
632 end if
633 call ylmbas(sij, rho, ctheta, stheta, cphi, sphi, &
634 & params % lmax, constants % vscales, vylm, vplm, &
635 & vcos, vsin)
636 tt = oij
637 do l = 0, params % lmax
638 ind = l*l + l + 1
639 fac = - tt/(constants % vscales(ind)**2)
640 do m = -l, l
641 scratch(ind + m, igrid) = fac*vylm(ind + m)
642 end do
643 tt = tt*tij
644 end do
645 end if
646 end do
647 call dgemm('n', 't', constants % nbasis, constants % nbasis, params % ngrid, &
648 & one, constants % vwgrid, constants % vgrid_nbasis, scratch, &
649 & constants % nbasis, zero, constants % l(:,:,ij), constants % nbasis)
650 end do
651 end do
652end subroutine build_l
653
655subroutine build_b(constants, params, ddx_error)
656 implicit none
657 type(ddx_params_type), intent(in) :: params
658 type(ddx_constants_type), intent(inout) :: constants
659 type(ddx_error_type), intent(inout) :: ddx_error
660 integer :: isph, ij, jsph, igrid, l, m, ind
661 real(dp), dimension(3) :: vij, sij
662 real(dp) :: vvij, tij, xij, oij, rho, ctheta, stheta, cphi, sphi, &
663 & fac, vvtij, thigh
664 real(dp), dimension(constants % nbasis) :: vylm, vplm
665 real(dp), dimension(params % lmax + 1) :: vcos, vsin
666 complex(dp), dimension(max(2, params % lmax + 1)) :: bessel_work
667 real(dp), dimension(0:params % lmax) :: SI_rijn, DI_rijn
668 real(dp), dimension(constants % nbasis, params % ngrid) :: scratch
669 real(dp) :: t
670 integer :: info
671
672 allocate(constants % b(constants % nbasis, constants % nbasis, &
673 & constants % inl(params % nsph + 1)), stat=info)
674 if (info.ne.0) then
675 call update_error(ddx_error, "Allocation failed in build_b")
676 return
677 end if
678
679 thigh = one + pt5*(params % se + one)*params % eta
680
681 t = omp_get_wtime()
682 !$omp parallel do default(none) shared(params,constants,thigh) &
683 !$omp private(isph,ij,jsph,scratch,igrid,vij,vvij,tij,sij,xij,oij, &
684 !$omp rho,ctheta,stheta,cphi,sphi,vylm,vplm,vcos,vsin,si_rijn,di_rijn, &
685 !$omp vvtij,l,fac,ind,m,bessel_work)
686 do isph = 1, params % nsph
687 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
688 jsph = constants % nl(ij)
689 scratch = zero
690 do igrid = 1, params % ngrid
691 if (constants % ui(igrid, isph).eq.one) cycle
692 vij = params % csph(:, isph) &
693 & + params % rsph(isph)*constants % cgrid(:,igrid) &
694 & - params % csph(:, jsph)
695 vvij = sqrt(dot_product(vij, vij))
696 tij = vvij/params % rsph(jsph)
697 if (tij.lt.thigh) then
698 if (tij.ne.zero) then
699 sij = vij/vvij
700 else
701 sij = one
702 end if
703 xij = fsw(tij, params % se, params % eta)
704 if (constants % fi(igrid, isph).gt.one) then
705 oij = xij/constants % fi(igrid, isph)
706 else
707 oij = xij
708 end if
709 call ylmbas(sij, rho, ctheta, stheta, cphi, sphi, &
710 & params % lmax, constants % vscales, vylm, vplm, &
711 & vcos, vsin)
712 si_rijn = 0
713 di_rijn = 0
714 vvtij = vvij*params % kappa
715 call modified_spherical_bessel_first_kind(params % lmax, &
716 & vvtij, si_rijn, di_rijn, bessel_work)
717 do l = 0, params % lmax
718 fac = - oij*si_rijn(l)/constants % SI_ri(l, jsph)
719 ind = l*l + l + 1
720 do m = -l, l
721 scratch(ind + m, igrid) = fac*vylm(ind + m)
722 end do
723 end do
724 end if
725 end do
726 call dgemm('n', 't', constants % nbasis, constants % nbasis, params % ngrid, &
727 & one, constants % vwgrid, constants % vgrid_nbasis, scratch, &
728 & constants % nbasis, zero, constants % b(:,:,ij), constants % nbasis)
729 end do
730 end do
731end subroutine build_b
732
734subroutine mkpmat(params, constants, isph, pmat)
735 type(ddx_params_type), intent(in) :: params
736 type(ddx_constants_type), intent(in) :: constants
737 integer, intent(in) :: isph
738 real(dp), dimension(constants % nbasis, constants % nbasis0), intent(inout) :: pmat
739 integer :: l, m, ind, l0, m0, ind0, its
740 real(dp) :: f, f0
741 pmat(:,:) = zero
742 do its = 1, params % ngrid
743 if (constants % ui(its,isph).ne.0) then
744 do l = 0, params % lmax
745 ind = l*l + l + 1
746 do m = -l,l
747 f = constants % wgrid(its) * constants % vgrid(ind+m,its) * constants % ui(its,isph)
748 do l0 = 0, constants % lmax0
749 ind0 = l0*l0 + l0 + 1
750 do m0 = -l0, l0
751 f0 = constants % vgrid(ind0+m0,its)
752 pmat(ind+m,ind0+m0) = pmat(ind+m,ind0+m0) + f * f0
753 end do
754 end do
755 end do
756 end do
757 end if
758 end do
759end subroutine mkpmat
760
770subroutine constants_geometry_init(params, constants, ddx_error)
771 !! Inputs
772 type(ddx_params_type), intent(in) :: params
773 !! Outputs
774 type(ddx_constants_type), intent(inout) :: constants
775 type(ddx_error_type), intent(inout) :: ddx_error
776 !! Local variables
777 real(dp) :: swthr, v(3), maxv, ssqv, vv, t
778 integer :: i, isph, jsph, inear, igrid, iwork, jwork, lwork, &
779 & old_lwork, icav, info
780 integer, allocatable :: work(:, :), tmp_work(:, :)
781 real(dp) :: start_time
782 !! The code
783 ! Prepare FMM structures if needed
784 start_time = omp_get_wtime()
785 if (params % fmm .eq. 1) then
786 ! Allocate space for a cluster tree
787 allocate(constants % order(params % nsph), stat=info)
788 if (info .ne. 0) then
789 call update_error(ddx_error, "constants_geometry_init: `order` " &
790 & // "allocation failed")
791 return
792 end if
793 constants % nclusters = 2*params % nsph - 1
794 allocate(constants % cluster(2, constants % nclusters), stat=info)
795 if (info .ne. 0) then
796 call update_error(ddx_error, "constants_geometry_init: `cluster` " &
797 & // "allocation failed")
798 return
799 end if
800 allocate(constants % children(2, constants % nclusters), stat=info)
801 if (info .ne. 0) then
802 call update_error(ddx_error, "constants_geometry_init: " // &
803 & "`children` allocation failed")
804 return
805 endif
806 allocate(constants % parent(constants % nclusters), stat=info)
807 if (info .ne. 0) then
808 call update_error(ddx_error, "constants_geometry_init: `parent` " &
809 & // "allocation failed")
810 return
811 endif
812 allocate(constants % cnode(3, constants % nclusters), stat=info)
813 if (info .ne. 0) then
814 call update_error(ddx_error, "constants_geometry_init: `cnode` " &
815 & // "allocation failed")
816 return
817 endif
818 allocate(constants % rnode(constants % nclusters), stat=info)
819 if (info .ne. 0) then
820 call update_error(ddx_error, "constants_geometry_init: `rnode` " &
821 & // "allocation failed")
822 return
823 endif
824 allocate(constants % snode(params % nsph), stat=info)
825 if (info .ne. 0) then
826 call update_error(ddx_error, "constants_geometry_init: `snode` " &
827 & // "allocation failed")
828 return
829 endif
830 allocate(constants % nfar(constants % nclusters), stat=info)
831 if (info .ne. 0) then
832 call update_error(ddx_error, "constants_geometry_init: `nfar` " &
833 & // "allocation failed")
834 return
835 endif
836 allocate(constants % nnear(constants % nclusters), stat=info)
837 if (info .ne. 0) then
838 call update_error(ddx_error, "constants_geometry_init: `nnear` " &
839 & // "allocation failed")
840 return
841 endif
842 ! Get the tree
843 call tree_rib_build(params % nsph, params % csph, params % rsph, &
844 & constants % order, constants % cluster, constants % children, &
845 & constants % parent, constants % cnode, constants % rnode, &
846 & constants % snode, ddx_error)
847 if (ddx_error % flag .ne. 0) then
848 call update_error(ddx_error, "tree_rib_build returned an " // &
849 & "ddx_error, exiting")
850 return
851 end if
852 ! Get number of far and near admissible pairs
853 iwork = 0
854 jwork = 1
855 lwork = 1
856 allocate(work(3, lwork), stat=info)
857 if (info .ne. 0) then
858 call update_error(ddx_error, "constants_geometry_init: `work` " &
859 & // "allocation failed")
860 return
861 end if
862 do while (iwork .le. jwork)
863 allocate(tmp_work(3, lwork), stat=info)
864 if (info .ne. 0) then
865 call update_error(ddx_error, "constants_geometry_init: " // &
866 & "`tmp_work` allocation failed")
867 return
868 end if
869 tmp_work = work
870 deallocate(work, stat=info)
871 if (info .ne. 0) then
872 call update_error(ddx_error, "constants_geometry_init: " // &
873 & "`work` deallocation failed")
874 return
875 end if
876 old_lwork = lwork
877 lwork = old_lwork + 1000*params % nsph
878 allocate(work(3, lwork), stat=info)
879 if (info .ne. 0) then
880 call update_error(ddx_error, "constants_geometry_init: " // &
881 & "`work` allocation failed")
882 return
883 end if
884 work(:, 1:old_lwork) = tmp_work
885 deallocate(tmp_work, stat=info)
886 if (info .ne. 0) then
887 call update_error(ddx_error, "constants_geometry_init: " // &
888 & "`tmp_work` deallocation failed")
889 return
890 end if
891 call tree_get_farnear_work(constants % nclusters, &
892 & constants % children, constants % cnode, &
893 & constants % rnode, lwork, iwork, jwork, work, &
894 & constants % nnfar, constants % nfar, constants % nnnear, &
895 & constants % nnear)
896 end do
897 allocate(constants % sfar(constants % nclusters+1), &
898 & constants % snear(constants % nclusters+1), stat=info)
899 if (info .ne. 0) then
900 call update_error(ddx_error, "constants_geometry_init: `sfar` " // &
901 & "and `snear` allocations failed")
902 return
903 end if
904 allocate(constants % far(constants % nnfar), stat=info)
905 if (info .ne. 0) then
906 call update_error(ddx_error, "constants_geometry_init: `far` " // &
907 & "allocation failed")
908 return
909 end if
910 allocate(constants % near(constants % nnnear), stat=info)
911 if (info .ne. 0) then
912 call update_error(ddx_error, "constants_geometry_init: `near` " // &
913 & "allocation failed")
914 return
915 end if
916 call tree_get_farnear(jwork, lwork, work, constants % nclusters, &
917 & constants % nnfar, constants % nfar, constants % sfar, &
918 & constants % far, constants % nnnear, constants % nnear, &
919 & constants % snear, constants % near)
920 deallocate(work, stat=info)
921 if (info .ne. 0) then
922 call update_error(ddx_error, "constants_geometry_init: `work` " // &
923 & "deallocation failed")
924 return
925 end if
926 end if
927 ! Upper bound of switch region. Defines intersection criterion for spheres
928 swthr = one + (params % se+one)*params % eta/two
929 ! Assemble neighbor list
930 if (params % fmm .eq. 1) then
931 call neighbor_list_init_fmm(params, constants, ddx_error)
932 else
933 call neighbor_list_init(params, constants, ddx_error)
934 end if
935 if (ddx_error % flag .ne. 0) then
936 call update_error(ddx_error, "Neighbor list construction failed")
937 return
938 end if
939 ! Allocate space for characteristic functions fi and ui
940 allocate(constants % fi(params % ngrid, params % nsph), &
941 & constants % ui(params % ngrid, params % nsph), stat=info)
942 if (info .ne. 0) then
943 call update_error(ddx_error, "`fi` and `ui` allocations failed")
944 return
945 end if
946 constants % fi = zero
947 constants % ui = zero
948 ! Allocate space for force-related arrays
949 if (params % force .eq. 1) then
950 allocate(constants % zi(3, params % ngrid, params % nsph), stat=info)
951 if (info .ne. 0) then
952 call update_error(ddx_error, "`zi` allocation failed")
953 return
954 endif
955 constants % zi = zero
956
957 allocate(constants % zi_dr(params % ngrid, params % nsph), stat=info)
958 if (info .ne. 0) then
959 call update_error(ddx_error, "`zi_dr` allocation failed")
960 return
961 endif
962 constants % zi_dr = zero
963 end if
964 ! Build arrays fi, ui, zi
965 !$omp parallel do default(none) shared(params,constants,swthr) &
966 !$omp private(isph,igrid,jsph,v,maxv,ssqv,vv,t) schedule(dynamic)
967 do isph = 1, params % nsph
968 do igrid = 1, params % ngrid
969 ! Loop over neighbours of i-th sphere
970 do inear = constants % inl(isph), constants % inl(isph+1)-1
971 ! Neighbour's index
972 jsph = constants % nl(inear)
973 ! Compute t_n^ij
974 v = params % csph(:, isph) - params % csph(:, jsph) + &
975 & params % rsph(isph)*constants % cgrid(:, igrid)
976 maxv = max(abs(v(1)), abs(v(2)), abs(v(3)))
977 if (maxv .gt. zero) then
978 ssqv = (v(1)/maxv)**2 + (v(2)/maxv)**2 + (v(3)/maxv)**2
979 else
980 ssqv = zero
981 end if
982 vv = maxv * sqrt(ssqv)
983 t = vv / params % rsph(jsph)
984 ! Accumulate characteristic function \chi(t_n^ij)
985 constants % fi(igrid, isph) = constants % fi(igrid, isph) + &
986 & fsw(t, params % se, params % eta)
987 ! Check if gradients are required
988 if (params % force .eq. 1) then
989 ! Check if t_n^ij belongs to switch region
990 if ((t .lt. swthr) .and. (t .gt. swthr-params % eta)) then
991 ! Accumulate gradient of characteristic function \chi
992 vv = dfsw(t, params % se, params % eta) / &
993 & params % rsph(jsph) / vv
994 constants % zi(:, igrid, isph) = &
995 & constants % zi(:, igrid, isph) + vv*v
996 constants % zi_dr(igrid, isph) = &
997 & constants % zi_dr(igrid, isph) &
998 & + vv*dot_product(v, constants % cgrid(:, igrid))
999 end if
1000 end if
1001 enddo
1002 ! Compute characteristic function of a molecular surface ui
1003 if (constants % fi(igrid, isph) .le. one) then
1004 constants % ui(igrid, isph) = one - constants % fi(igrid, isph)
1005 end if
1006 enddo
1007 enddo
1008 ! Build cavity array. At first get total count for each sphere
1009 allocate(constants % ncav_sph(params % nsph), stat=info)
1010 if (info .ne. 0) then
1011 call update_error(ddx_error, "`ncav_sph` allocation failed")
1012 return
1013 endif
1014 !$omp parallel do default(none) shared(params,constants) &
1015 !$omp private(isph,i) schedule(dynamic)
1016 do isph = 1, params % nsph
1017 constants % ncav_sph(isph) = 0
1018 ! Loop over integration points
1019 do i = 1, params % ngrid
1020 ! Positive contribution from integration point
1021 if (constants % ui(i, isph) .gt. zero) then
1022 constants % ncav_sph(isph) = constants % ncav_sph(isph) + 1
1023 end if
1024 end do
1025 end do
1026 constants % ncav = sum(constants % ncav_sph)
1027 ! Allocate cavity array and CSR format for indexes of cavities
1028 allocate(constants % ccav(3, constants % ncav), &
1029 & constants % icav_ia(params % nsph+1), &
1030 & constants % icav_ja(constants % ncav), stat=info)
1031 if (info .ne. 0) then
1032 call update_error(ddx_error, "`ccav`, `icav_ia` and " // &
1033 & "`icav_ja` allocations failed")
1034 return
1035 endif
1036 ! Allocate space for characteristic functions ui at cavity points
1037 allocate(constants % ui_cav(constants % ncav), stat=info)
1038 if (info .ne. 0) then
1039 call update_error(ddx_error, "`ui_cav` allocations failed")
1040 return
1041 end if
1042 ! Get actual cavity coordinates and indexes in CSR format and fill in
1043 ! ui_cav aray
1044 constants % icav_ia(1) = 1
1045 icav = 1
1046 do isph = 1, params % nsph
1047 constants % icav_ia(isph+1) = constants % icav_ia(isph) + &
1048 & constants % ncav_sph(isph)
1049 ! Loop over integration points
1050 do igrid = 1, params % ngrid
1051 ! Ignore zero contribution
1052 if (constants % ui(igrid, isph) .eq. zero) cycle
1053 ! Store coordinates
1054 constants % ccav(:, icav) = params % csph(:, isph) + &
1055 & params % rsph(isph)* &
1056 & constants % cgrid(:, igrid)
1057 ! Store index
1058 constants % icav_ja(icav) = igrid
1059 ! Store characteristic function
1060 constants % ui_cav(icav) = constants % ui(igrid, isph)
1061 ! Advance cavity array index
1062 icav = icav + 1
1063 end do
1064 end do
1065 ! Compute diagonal preconditioner for PCM equations
1066 if (params % model .eq. 2) then
1067 allocate(constants % rx_prc( &
1068 & constants % nbasis, constants % nbasis, params % nsph), &
1069 & stat=info)
1070 if (info .ne. 0) then
1071 call update_error(ddx_error, "constants_geometry_init: `rx_prc` " &
1072 & // "allocation failed")
1073 return
1074 endif
1075 call mkprec(params % lmax, constants % nbasis, params % nsph, &
1076 & params % ngrid, params % eps, constants % ui, &
1077 & constants % wgrid, constants % vgrid, constants % vgrid_nbasis, &
1078 & constants % rx_prc, ddx_error)
1079 if (info .ne. 0) then
1080 call update_error(ddx_error, "mkprec returned an ddx_error, exiting")
1081 return
1082 end if
1083 end if
1084end subroutine constants_geometry_init
1085
1087subroutine neighbor_list_init(params, constants, ddx_error)
1088 type(ddx_params_type), intent(in) :: params
1089 type(ddx_constants_type), intent(inout) :: constants
1090 type(ddx_error_type), intent(inout) :: ddx_error
1091 real(dp) :: swthr, v(3), maxv, ssqv, vv, r
1092 integer :: nngmax, i, lnl, isph, jsph, info
1093 integer, allocatable :: tmp_nl(:)
1094 ! Upper bound of switch region. Defines intersection criterion for spheres
1095 swthr = (params % se+one)*params % eta/two
1096 ! Build list of neighbours in CSR format
1097 nngmax = 1
1098 allocate(constants % inl(params % nsph+1), &
1099 & constants % nl(params % nsph*nngmax), stat=info)
1100 if (info .ne. 0) then
1101 call update_error(ddx_error, "`inl` and `nl` allocations failed")
1102 return
1103 end if
1104 i = 1
1105 lnl = 0
1106 do isph = 1, params % nsph
1107 constants % inl(isph) = lnl + 1
1108 do jsph = 1, params % nsph
1109 if (isph .ne. jsph) then
1110 v = params % csph(:, isph)-params % csph(:, jsph)
1111 maxv = max(abs(v(1)), abs(v(2)), abs(v(3)))
1112 if (maxv .gt. zero) then
1113 ssqv = (v(1)/maxv)**2 + (v(2)/maxv)**2 + (v(3)/maxv)**2
1114 else
1115 ssqv = zero
1116 end if
1117 vv = maxv * sqrt(ssqv)
1118 ! Take regularization parameter into account with respect to
1119 ! shift se. It is described properly by the upper bound of a
1120 ! switch region `swthr`. Note that we take the largest
1121 ! switching region as we need a symmetric neighbor list to
1122 ! use it also for the adjoint products.
1123 r = params % rsph(isph) + params % rsph(jsph) &
1124 & + swthr * max(params % rsph(isph), params % rsph(jsph))
1125 if (vv .le. r) then
1126 constants % nl(i) = jsph
1127 i = i + 1
1128 lnl = lnl + 1
1129 ! Extend ddx_data % nl if needed
1130 if (i .gt. params % nsph*nngmax) then
1131 allocate(tmp_nl(params % nsph*nngmax), stat=info)
1132 if (info .ne. 0) then
1133 call update_error(ddx_error, "`tmp_nl` " // &
1134 & "allocation failed")
1135 return
1136 end if
1137 tmp_nl(1:params % nsph*nngmax) = &
1138 & constants % nl(1:params % nsph*nngmax)
1139 deallocate(constants % nl, stat=info)
1140 if (info .ne. 0) then
1141 call update_error(ddx_error, "`nl` " // &
1142 & "deallocation failed")
1143 return
1144 end if
1145 nngmax = nngmax + 10
1146 allocate(constants % nl(params % nsph*nngmax), &
1147 & stat=info)
1148 if (info .ne. 0) then
1149 call update_error(ddx_error, "`nl` " // &
1150 & "allocation failed")
1151 return
1152 end if
1153 constants % nl(1:params % nsph*(nngmax-10)) = &
1154 & tmp_nl(1:params % nsph*(nngmax-10))
1155 deallocate(tmp_nl, stat=info)
1156 if (info .ne. 0) then
1157 call update_error(ddx_error, "`tmp_nl` " // &
1158 & "deallocation failed")
1159 return
1160 end if
1161 end if
1162 end if
1163 end if
1164 end do
1165 end do
1166 constants % inl(params % nsph+1) = lnl+1
1167 constants % nngmax = nngmax
1168end subroutine neighbor_list_init
1169
1172subroutine neighbor_list_init_fmm(params, constants, ddx_error)
1173 type(ddx_params_type), intent(in) :: params
1174 type(ddx_constants_type), intent(inout) :: constants
1175 type(ddx_error_type), intent(inout) :: ddx_error
1176 real(dp) :: swthr, v(3), maxv, ssqv, vv, r
1177 integer :: nngmax, i, lnl, isph, jsph, inode, jnode, j, k, info
1178 integer, allocatable :: tmp_nl(:)
1179 ! Upper bound of switch region. Defines intersection criterion for spheres
1180 swthr = (params % se+one)*params % eta/two
1181 ! Build list of neighbours in CSR format
1182 nngmax = 1
1183 allocate(constants % inl(params % nsph+1), &
1184 & constants % nl(params % nsph*nngmax), stat=info)
1185 if (info .ne. 0) then
1186 call update_error(ddx_error, "`inl` and `nl` allocations failed")
1187 return
1188 end if
1189 i = 1
1190 lnl = 0
1191 ! loop over leaf clusters
1192 do isph = 1, params % nsph
1193 constants % inl(isph) = lnl + 1
1194 inode = constants % snode(isph)
1195 ! loop over near clusters
1196 do j = constants % snear(inode), constants % snear(inode+1) - 1
1197 jnode = constants % near(j)
1198 do k = constants % cluster(1, jnode), constants % cluster(2, jnode)
1199 jsph = constants % order(k)
1200 if (isph .ne. jsph) then
1201 v = params % csph(:, isph)-params % csph(:, jsph)
1202 maxv = max(abs(v(1)), abs(v(2)), abs(v(3)))
1203 if (maxv .gt. zero) then
1204 ssqv = (v(1)/maxv)**2 + (v(2)/maxv)**2 + (v(3)/maxv)**2
1205 else
1206 ssqv = zero
1207 end if
1208 vv = maxv * sqrt(ssqv)
1209 ! Take regularization parameter into account with respect to
1210 ! shift se. It is described properly by the upper bound of a
1211 ! switch region `swthr`. Note that we take the largest
1212 ! switching region as we need a symmetric neighbor list to
1213 ! use it also for the adjoint products.
1214 r = params % rsph(isph) + params % rsph(jsph) &
1215 & + swthr * max(params % rsph(isph), params % rsph(jsph))
1216 if (vv .le. r) then
1217 constants % nl(i) = jsph
1218 i = i + 1
1219 lnl = lnl + 1
1220 ! Extend ddx_data % nl if needed
1221 if (i .gt. params % nsph*nngmax) then
1222 allocate(tmp_nl(params % nsph*nngmax), stat=info)
1223 if (info .ne. 0) then
1224 call update_error(ddx_error, "`tmp_nl` " // &
1225 & "allocation failed")
1226 return
1227 end if
1228 tmp_nl(1:params % nsph*nngmax) = &
1229 & constants % nl(1:params % nsph*nngmax)
1230 deallocate(constants % nl, stat=info)
1231 if (info .ne. 0) then
1232 call update_error(ddx_error, "`nl` " // &
1233 & "deallocation failed")
1234 return
1235 end if
1236 nngmax = nngmax + 10
1237 allocate(constants % nl(params % nsph*nngmax), &
1238 & stat=info)
1239 if (info .ne. 0) then
1240 call update_error(ddx_error, "`nl` " // &
1241 & "allocation failed")
1242 return
1243 end if
1244 constants % nl(1:params % nsph*(nngmax-10)) = &
1245 & tmp_nl(1:params % nsph*(nngmax-10))
1246 deallocate(tmp_nl, stat=info)
1247 if (info .ne. 0) then
1248 call update_error(ddx_error, "`tmp_nl` " // &
1249 & "deallocation failed")
1250 return
1251 end if
1252 end if
1253 end if
1254 end if
1255 end do
1256 end do
1257 end do
1258 constants % inl(params % nsph+1) = lnl+1
1259 constants % nngmax = nngmax
1260end subroutine neighbor_list_init_fmm
1261
1289real(dp) function fsw(t, se, eta)
1290 ! Inputs
1291 real(dp), intent(in) :: t, se, eta
1292 ! Local variables
1293 real(dp) :: a, b, c, x
1294 ! Apply shift:
1295 ! se = 0 => t - eta/2 [ CENTERED ]
1296 ! se = 1 => t - eta [ EXTERIOR ]
1297 ! se = -1 => t [ INTERIOR ]
1298 x = t - (se*pt5 + pt5)*eta
1299 ! a must be in range (0,eta) for the switch region
1300 a = one - x
1301 ! Define switch function chi
1302 ! a <= 0 corresponds to the outside region
1303 if (a .le. zero) then
1304 fsw = zero
1305 ! a >= eta corresponds to the inside region
1306 else if (a .ge. eta) then
1307 fsw = one
1308 ! Switch region wher x is in range (1-eta,1)
1309 else
1310 ! Normalize a to region (0,1)
1311 a = a / eta
1312 b = 6d0*a - 15
1313 c = b*a + 10
1314 fsw = a * a * a * c
1315 end if
1316end function fsw
1317
1344real(dp) function dfsw(t, se, eta)
1345 ! Inputs
1346 real(dp), intent(in) :: t, se, eta
1347 ! Local variables
1348 real(dp) :: flow, x
1349 real(dp), parameter :: f30=30.0d0
1350 ! Apply shift:
1351 ! s = 0 => t - eta/2 [ CENTERED ]
1352 ! s = 1 => t - eta [ EXTERIOR ]
1353 ! s = -1 => t [ INTERIOR ]
1354 x = t - (se + 1.d0)*eta / 2.d0
1355 ! Lower bound of switch region
1356 flow = one - eta
1357 ! Define derivative of switch function chi
1358 if (x .ge. one) then
1359 dfsw = zero
1360 else if (x .le. flow) then
1361 dfsw = zero
1362 else
1363 dfsw = -f30 * (( (one-x)*(x-one+eta) )**2) / (eta**5)
1364 endif
1365end function dfsw
1366
1370subroutine mkprec(lmax, nbasis, nsph, ngrid, eps, ui, wgrid, vgrid, &
1371 & vgrid_nbasis, rx_prc, ddx_error)
1372 !! Inputs
1373 integer, intent(in) :: lmax, nbasis, nsph, ngrid, vgrid_nbasis
1374 real(dp), intent(in) :: eps, ui(ngrid, nsph), wgrid(ngrid), &
1375 & vgrid(vgrid_nbasis, ngrid)
1376 !! Output
1377 real(dp), intent(out) :: rx_prc(nbasis, nbasis, nsph)
1378 type(ddx_error_type), intent(inout) :: ddx_error
1379 !! Local variables
1380 integer :: info, isph, lm, l1, m1, ind1, igrid
1381 real(dp) :: f, f1
1382 integer, allocatable :: ipiv(:)
1383 real(dp), allocatable :: work(:)
1384 external :: dgetrf, dgetri
1385 ! Allocation of temporaries
1386 allocate(ipiv(nbasis), work(nbasis**2), stat=info)
1387 if (info .ne. 0) then
1388 call update_error(ddx_error, "mkprec: `ipiv` and `work` allocation failed")
1389 return
1390 endif
1391 ! Init
1392 rx_prc = zero
1393 ! Off-diagonal part
1394 do isph = 1, nsph
1395 do igrid = 1, ngrid
1396 f = twopi * ui(igrid, isph) * wgrid(igrid)
1397 do l1 = 0, lmax
1398 ind1 = l1*l1 + l1 + 1
1399 do m1 = -l1, l1
1400 f1 = f * vgrid(ind1+m1, igrid) / (two*dble(l1) + one)
1401 do lm = 1, nbasis
1402 rx_prc(lm, ind1+m1, isph) = &
1403 & rx_prc(lm, ind1+m1, isph) + f1*vgrid(lm, igrid)
1404 end do
1405 end do
1406 end do
1407 end do
1408 end do
1409 ! add diagonal
1410 f = twopi * (eps+one) / (eps-one)
1411 do isph = 1, nsph
1412 do lm = 1, nbasis
1413 rx_prc(lm, lm, isph) = rx_prc(lm, lm, isph) + f
1414 end do
1415 end do
1416 ! invert the blocks
1417 do isph = 1, nsph
1418 call dgetrf(nbasis, nbasis, rx_prc(:, :, isph), nbasis, ipiv, info)
1419 if (info .ne. 0) then
1420 call update_error(ddx_error, "mkprec: dgetrf failed")
1421 return
1422 end if
1423 call dgetri(nbasis, rx_prc(:, :, isph), nbasis, ipiv, work, &
1424 & nbasis**2, info)
1425 if (info .ne. 0) then
1426 call update_error(ddx_error, "mkprec: dgetri failed")
1427 return
1428 end if
1429 end do
1430 ! Cleanup temporaries
1431 deallocate(work, ipiv, stat=info)
1432 if (info .ne. 0) then
1433 call update_error(ddx_error, "mkprec: `ipiv` and `work` " // &
1434 & "deallocation failed")
1435 return
1436 end if
1437end subroutine mkprec
1438
1461subroutine tree_rib_build(nsph, csph, rsph, order, cluster, children, parent, &
1462 & cnode, rnode, snode, ddx_error)
1463 ! Inputs
1464 integer, intent(in) :: nsph
1465 real(dp), intent(in) :: csph(3, nsph), rsph(nsph)
1466 ! Outputs
1467 integer, intent(out) :: order(nsph), cluster(2, 2*nsph-1), &
1468 & children(2, 2*nsph-1), parent(2*nsph-1), snode(nsph)
1469 real(dp), intent(out) :: cnode(3, 2*nsph-1), rnode(2*nsph-1)
1470 type(ddx_error_type), intent(inout) :: ddx_error
1471 ! Local variables
1472 integer :: nclusters, i, j, n, s, e, div
1473 real(dp) :: r, r1, r2, c(3), c1(3), c2(3), d, maxc
1474 !! At first construct the tree
1475 nclusters = 2*nsph - 1
1476 ! Init the root node
1477 cluster(1, 1) = 1
1478 cluster(2, 1) = nsph
1479 parent(1) = 0
1480 ! Index of the first unassigned node
1481 j = 2
1482 ! Init spheres ordering
1483 do i = 1, nsph
1484 order(i) = i
1485 end do
1486 ! Divide nodes until a single spheres
1487 do i = 1, nclusters
1488 ! The first sphere in i-th node
1489 s = cluster(1, i)
1490 ! The last sphere in i-th node
1491 e = cluster(2, i)
1492 ! Number of spheres in i-th node
1493 n = e - s + 1
1494 ! Divide only if there are 2 or more spheres
1495 if (n .gt. 1) then
1496 ! Use inertial bisection to reorder spheres and cut into 2 halfs
1497 call tree_rib_node_bisect(nsph, csph, n, order(s:e), div, &
1498 & ddx_error)
1499 if (ddx_error % flag .ne. 0) then
1500 call update_error(ddx_error, "tree_rib_node_bisect returned " // &
1501 & "an ddx_error, exiting")
1502 return
1503 end if
1504 ! Assign the first half to the j-th node
1505 cluster(1, j) = s
1506 cluster(2, j) = s + div - 1
1507 ! Assign the second half to the (j+1)-th node
1508 cluster(1, j+1) = s + div
1509 cluster(2, j+1) = e
1510 ! Update list of children of i-th node
1511 children(1, i) = j
1512 children(2, i) = j + 1
1513 ! Set parents of new j-th and (j+1)-th node
1514 parent(j) = i
1515 parent(j+1) = i
1516 ! Shift index of the first unassigned node
1517 j = j + 2
1518 ! Set information for a leaf node
1519 else
1520 ! No children nodes
1521 children(:, i) = 0
1522 ! i-th node contains sphere ind(s)
1523 snode(order(s)) = i
1524 end if
1525 end do
1526 !! Compute bounding spheres of each node of the tree
1527 ! Bottom-to-top pass over all nodes of the tree
1528 do i = nclusters, 1, -1
1529 ! In case of a leaf node just use corresponding input sphere as a
1530 ! bounding sphere of the node
1531 if (children(1, i) .eq. 0) then
1532 ! Get correct index of the corresponding input sphere
1533 j = order(cluster(1, i))
1534 ! Get corresponding center and radius
1535 cnode(:, i) = csph(:, j)
1536 rnode(i) = rsph(j)
1537 ! In case of a non-leaf node get minimal sphere that contains bounding
1538 ! spheres of children nodes
1539 else
1540 ! The first child
1541 j = children(1, i)
1542 c1 = cnode(:, j)
1543 r1 = rnode(j)
1544 ! The second child
1545 j = children(2, i)
1546 c2 = cnode(:, j)
1547 r2 = rnode(j)
1548 ! Distance between centers of bounding spheres of children nodes
1549 c = c1 - c2
1550 ! Compute distance by scale and sum of scaled squares
1551 maxc = max(abs(c(1)), abs(c(2)), abs(c(3)))
1552 if (maxc .eq. zero) then
1553 d = zero
1554 else
1555 d = (c(1)/maxc)**2 + (c(2)/maxc)**2 + (c(3)/maxc)**2
1556 d = maxc * sqrt(d)
1557 end if
1558 ! If sphere #2 is completely inside sphere #1 use the first sphere
1559 if (r1 .ge. (r2+d)) then
1560 c = c1
1561 r = r1
1562 ! If sphere #1 is completely inside sphere #2 use the second sphere
1563 else if(r2 .ge. (r1+d)) then
1564 c = c2
1565 r = r2
1566 ! Otherwise use special formula to find a minimal sphere
1567 else
1568 r = (r1+r2+d) / 2
1569 c = c2 + c/d*(r-r2)
1570 end if
1571 ! Assign computed bounding sphere
1572 cnode(:, i) = c
1573 rnode(i) = r
1574 end if
1575 end do
1576end subroutine tree_rib_build
1577
1588subroutine tree_rib_node_bisect(nsph, csph, n, order, div, ddx_error)
1589 ! Inputs
1590 integer, intent(in) :: nsph, n
1591 real(dp), intent(in) :: csph(3, nsph)
1592 ! Outputs
1593 integer, intent(inout) :: order(n)
1594 integer, intent(out) :: div
1595 type(ddx_error_type), intent(inout) :: ddx_error
1596 ! Local variables
1597 real(dp) :: c(3), s(3)
1598 real(dp), allocatable :: tmp_csph(:, :), work(:)
1599 external :: dgesvd
1600 integer :: i, l, r, lwork, info, istat
1601 integer, allocatable :: tmp_order(:)
1602
1603 allocate(tmp_csph(3, n), tmp_order(n), stat=istat)
1604 if (istat.ne.0) then
1605 call update_error(ddx_error, 'Allocation failed in tree_node_bisect')
1606 return
1607 end if
1608
1609 ! Get average coordinate
1610 c = zero
1611 do i = 1, n
1612 c = c + csph(:, order(i))
1613 end do
1614 c = c / n
1615 ! Get coordinates minus average in a contiguous array
1616 do i = 1, n
1617 tmp_csph(:, i) = csph(:, order(i)) - c
1618 end do
1619 !! Find right singular vectors
1620 ! Get proper size of temporary workspace
1621 lwork = -1
1622 call dgesvd('N', 'O', 3, n, tmp_csph, 3, s, tmp_csph, 3, tmp_csph, 3, &
1623 & s, lwork, info)
1624 lwork = int(s(1))
1625 allocate(work(lwork), stat=istat)
1626 if (istat.ne.0) then
1627 call update_error(ddx_error, 'Allocation failed in tree_node_bisect')
1628 return
1629 end if
1630 ! Get right singular vectors
1631 call dgesvd('N', 'O', 3, n, tmp_csph, 3, s, tmp_csph, 3, tmp_csph, 3, &
1632 & work, lwork, info)
1633 if (info.ne.0) then
1634 call update_error(ddx_error, 'DGESVD failed in tree_node_bisect')
1635 return
1636 end if
1637 deallocate(work, stat=istat)
1638 if (istat.ne.0) then
1639 call update_error(ddx_error, 'Deallocation failed in tree_node_bisect')
1640 return
1641 end if
1642 !! Sort spheres by sign of the above scalar product, which is equal to
1643 !! the leading right singular vector scaled by the leading singular value.
1644 !! However, we only care about sign, so we take into account only the
1645 !! leading right singular vector.
1646 ! First empty index from the left
1647 l = 1
1648 ! First empty index from the right
1649 r = n
1650 ! Cycle over values of the singular vector
1651 do i = 1, n
1652 ! Positive scalar products are moved to the beginning of `order`
1653 if (tmp_csph(1, i) .ge. zero) then
1654 tmp_order(l) = order(i)
1655 l = l + 1
1656 ! Negative scalar products are moved to the end of `order`
1657 else
1658 tmp_order(r) = order(i)
1659 r = r - 1
1660 end if
1661 end do
1662 ! Set divider and update order
1663 div = r
1664 order = tmp_order
1665 deallocate(tmp_csph, tmp_order, stat=istat)
1666 if (istat.ne.0) then
1667 call update_error(ddx_error, 'Deallocation failed in tree_node_bisect')
1668 return
1669 end if
1670end subroutine tree_rib_node_bisect
1671
1696subroutine tree_get_farnear_work(n, children, cnode, rnode, lwork, iwork, &
1697 & jwork, work, nnfar, nfar, nnnear, nnear)
1698 ! Inputs
1699 integer, intent(in) :: n, children(2, n), lwork
1700 real(dp), intent(in) :: cnode(3, n), rnode(n)
1701 ! Outputs
1702 integer, intent(inout) :: iwork, jwork, work(3, lwork)
1703 integer, intent(out) :: nnfar, nfar(n), nnnear, nnear(n)
1704 ! Local variables
1705 integer :: j(2), npairs, k1, k2
1706 real(dp) :: c(3), r, d, dmax, dssq
1707 ! iwork is current temporary item in work array to process
1708 if (iwork .eq. 0) then
1709 work(1, 1) = 1
1710 work(2, 1) = 1
1711 iwork = 1
1712 jwork = 1
1713 end if
1714 ! jwork is total amount of temporary items in work array
1715 do while (iwork .le. jwork)
1716 j = work(1:2, iwork)
1717 c = cnode(:, j(1)) - cnode(:, j(2))
1718 r = rnode(j(1)) + rnode(j(2)) + max(rnode(j(1)), rnode(j(2)))
1719 !r = rnode(j(1)) + rnode(j(2))
1720 dmax = max(abs(c(1)), abs(c(2)), abs(c(3)))
1721 if (dmax .gt. zero) then
1722 dssq = (c(1)/dmax)**2 + (c(2)/dmax)**2 + (c(3)/dmax)**2
1723 else
1724 dssq = zero
1725 end if
1726 d = dmax * sqrt(dssq)
1727 !d = sqrt(c(1)**2 + c(2)**2 + c(3)**2)
1728 ! If node has no children then assume itself for purpose of finding
1729 ! far-field and near-filed interactions with children nodes of another
1730 ! node
1731 npairs = max(1, children(2, j(1))-children(1, j(1))+1) * &
1732 & max(1, children(2, j(2))-children(1, j(2))+1)
1733 if (d .ge. r) then
1734 ! Mark as far admissible pair
1735 work(3, iwork) = 1
1736 else if (npairs .eq. 1) then
1737 ! Mark as near admissible pair if both nodes are leaves
1738 work(3, iwork) = 2
1739 else if (jwork+npairs .gt. lwork) then
1740 ! Exit procedure, since work array was too small
1741 return
1742 else
1743 ! Mark as non-admissible pair and check all pairs of children nodes
1744 ! or pairs of one node (if it is a leaf node) with children of
1745 ! another node
1746 work(3, iwork) = 0
1747 if (children(1, j(1)) .eq. 0) then
1748 k1 = j(1)
1749 do k2 = children(1, j(2)), children(2, j(2))
1750 jwork = jwork + 1
1751 work(1, jwork) = k1
1752 work(2, jwork) = k2
1753 end do
1754 else if(children(1, j(2)) .eq. 0) then
1755 k2 = j(2)
1756 do k1 = children(1, j(1)), children(2, j(1))
1757 jwork = jwork + 1
1758 work(1, jwork) = k1
1759 work(2, jwork) = k2
1760 end do
1761 else
1762 do k1 = children(1, j(1)), children(2, j(1))
1763 do k2 = children(1, j(2)), children(2, j(2))
1764 jwork = jwork + 1
1765 work(1, jwork) = k1
1766 work(2, jwork) = k2
1767 end do
1768 end do
1769 end if
1770 end if
1771 iwork = iwork + 1
1772 end do
1773 nfar = 0
1774 nnear = 0
1775 do iwork = 1, jwork
1776 if (work(3, iwork) .eq. 1) then
1777 nfar(work(1, iwork)) = nfar(work(1, iwork)) + 1
1778 else if (work(3, iwork) .eq. 2) then
1779 nnear(work(1, iwork)) = nnear(work(1, iwork)) + 1
1780 end if
1781 end do
1782 iwork = jwork + 1
1783 nnfar = sum(nfar)
1784 nnnear = sum(nnear)
1785end subroutine tree_get_farnear_work
1786
1789subroutine tree_get_farnear(jwork, lwork, work, n, nnfar, nfar, sfar, far, &
1790 & nnnear, nnear, snear, near)
1791! Parameters:
1792! jwork: Total number of checked pairs in work array
1793! lwork: Total length of work array
1794! work: Work array itself
1795! n: Number of nodes
1796! nnfar: Total number of all far-field interactions
1797! nfar: Number of far-field interactions of each node
1798! sfar: Index in far array of first far-field node for each node
1799! far: Indexes of far-field nodes
1800! nnnear: Total number of all near-field interactions
1801! nnear: Number of near-field interactions of each node
1802! snear: Index in near array of first near-field node for each node
1803! near: Indexes of near-field nodes
1804 integer, intent(in) :: jwork, lwork, work(3, lwork), n, nnfar, nnnear
1805 integer, intent(in) :: nfar(n), nnear(n)
1806 integer, intent(out) :: sfar(n+1), far(nnfar), snear(n+1), near(nnnear)
1807 integer :: i, j
1808 integer, allocatable :: cfar(:), cnear(:)
1809
1810 allocate(cfar(n+1), cnear(n+1))
1811 sfar(1) = 1
1812 snear(1) = 1
1813 do i = 2, n+1
1814 sfar(i) = sfar(i-1) + nfar(i-1)
1815 snear(i) = snear(i-1) + nnear(i-1)
1816 end do
1817 cfar = sfar
1818 cnear = snear
1819 do i = 1, jwork
1820 if (work(3, i) .eq. 1) then
1821 ! Far
1822 j = work(1, i)
1823 if ((j .gt. n) .or. (j .le. 0)) then
1824 write(*,*) "ALARM", j
1825 end if
1826 far(cfar(j)) = work(2, i)
1827 cfar(j) = cfar(j) + 1
1828 else if (work(3, i) .eq. 2) then
1829 ! Near
1830 j = work(1, i)
1831 if ((j .gt. n) .or. (j .le. 0)) then
1832 write(*,*) "ALARM", j
1833 end if
1834 near(cnear(j)) = work(2, i)
1835 cnear(j) = cnear(j) + 1
1836 end if
1837 end do
1838
1839 deallocate(cfar, cnear)
1840
1841end subroutine tree_get_farnear
1842
1848subroutine constants_free(constants, ddx_error)
1849 implicit none
1850 type(ddx_constants_type), intent(out) :: constants
1851 type(ddx_error_type), intent(inout) :: ddx_error
1852 integer :: istat
1853
1854 istat = 0
1855
1856 if (allocated(constants % vscales)) then
1857 deallocate(constants % vscales, stat=istat)
1858 if (istat .ne. 0) then
1859 call update_error(ddx_error, "`vscales` deallocation failed!")
1860 end if
1861 end if
1862 if (allocated(constants % v4pi2lp1)) then
1863 deallocate(constants % v4pi2lp1, stat=istat)
1864 if (istat .ne. 0) then
1865 call update_error(ddx_error, "`v4pi2lp1` deallocation failed!")
1866 end if
1867 end if
1868 if (allocated(constants % vscales_rel)) then
1869 deallocate(constants % vscales_rel, stat=istat)
1870 if (istat .ne. 0) then
1871 call update_error(ddx_error, "`vscales_rel` deallocation failed!")
1872 end if
1873 end if
1874 if (allocated(constants % vfact)) then
1875 deallocate(constants % vfact, stat=istat)
1876 if (istat .ne. 0) then
1877 call update_error(ddx_error, "`vfact` deallocation failed!")
1878 end if
1879 end if
1880 if (allocated(constants % vcnk)) then
1881 deallocate(constants % vcnk, stat=istat)
1882 if (istat .ne. 0) then
1883 call update_error(ddx_error, "`vcnk` deallocation failed!")
1884 end if
1885 end if
1886 if (allocated(constants % m2l_ztranslate_coef)) then
1887 deallocate(constants % m2l_ztranslate_coef, stat=istat)
1888 if (istat .ne. 0) then
1889 call update_error(ddx_error, &
1890 & "`m2l_ztranslate_coef` deallocation failed!")
1891 end if
1892 end if
1893 if (allocated(constants % m2l_ztranslate_adj_coef)) then
1894 deallocate(constants % m2l_ztranslate_adj_coef, stat=istat)
1895 if (istat .ne. 0) then
1896 call update_error(ddx_error, &
1897 & "`m2l_ztranslate_adj_coef` deallocation failed!")
1898 end if
1899 end if
1900 if (allocated(constants % cgrid)) then
1901 deallocate(constants % cgrid, stat=istat)
1902 if (istat .ne. 0) then
1903 call update_error(ddx_error, "`cgrid` deallocation failed!")
1904 end if
1905 end if
1906 if (allocated(constants % wgrid)) then
1907 deallocate(constants % wgrid, stat=istat)
1908 if (istat .ne. 0) then
1909 call update_error(ddx_error, "`wgrid` deallocation failed!")
1910 end if
1911 end if
1912 if (allocated(constants % vgrid)) then
1913 deallocate(constants % vgrid, stat=istat)
1914 if (istat .ne. 0) then
1915 call update_error(ddx_error, "`vgrid` deallocation failed!")
1916 end if
1917 end if
1918 if (allocated(constants % vwgrid)) then
1919 deallocate(constants % vwgrid, stat=istat)
1920 if (istat .ne. 0) then
1921 call update_error(ddx_error, "`vwgrid` deallocation failed!")
1922 end if
1923 end if
1924 if (allocated(constants % vgrid2)) then
1925 deallocate(constants % vgrid2, stat=istat)
1926 if (istat .ne. 0) then
1927 call update_error(ddx_error, "`vgrid2` deallocation failed!")
1928 end if
1929 end if
1930 if (allocated(constants % pchi)) then
1931 deallocate(constants % pchi, stat=istat)
1932 if (istat .ne. 0) then
1933 call update_error(ddx_error, "`pchi` deallocation failed!")
1934 end if
1935 end if
1936 if (allocated(constants % c_ik)) then
1937 deallocate(constants % c_ik, stat=istat)
1938 if (istat .ne. 0) then
1939 call update_error(ddx_error, "`c_ik` deallocation failed!")
1940 end if
1941 end if
1942 if (allocated(constants % si_ri)) then
1943 deallocate(constants % si_ri, stat=istat)
1944 if (istat .ne. 0) then
1945 call update_error(ddx_error, "`si_ri` deallocation failed!")
1946 end if
1947 end if
1948 if (allocated(constants % di_ri)) then
1949 deallocate(constants % di_ri, stat=istat)
1950 if (istat .ne. 0) then
1951 call update_error(ddx_error, "`di_ri` deallocation failed!")
1952 end if
1953 end if
1954 if (allocated(constants % sk_ri)) then
1955 deallocate(constants % sk_ri, stat=istat)
1956 if (istat .ne. 0) then
1957 call update_error(ddx_error, "`sk_ri` deallocation failed!")
1958 end if
1959 end if
1960 if (allocated(constants % dk_ri)) then
1961 deallocate(constants % dk_ri, stat=istat)
1962 if (istat .ne. 0) then
1963 call update_error(ddx_error, "`dk_ri` deallocation failed!")
1964 end if
1965 end if
1966 if (allocated(constants % termimat)) then
1967 deallocate(constants % termimat, stat=istat)
1968 if (istat .ne. 0) then
1969 call update_error(ddx_error, "`termimat` deallocation failed!")
1970 end if
1971 end if
1972 if (allocated(constants % b)) then
1973 deallocate(constants % b, stat=istat)
1974 if (istat .ne. 0) then
1975 call update_error(ddx_error, "`b` deallocation failed!")
1976 end if
1977 end if
1978 if (allocated(constants % l)) then
1979 deallocate(constants % l, stat=istat)
1980 if (istat .ne. 0) then
1981 call update_error(ddx_error, "`l` deallocation failed!")
1982 end if
1983 end if
1984 if (allocated(constants % inl)) then
1985 deallocate(constants % inl, stat=istat)
1986 if (istat .ne. 0) then
1987 call update_error(ddx_error, "`inl` deallocation failed!")
1988 end if
1989 end if
1990 if (allocated(constants % nl)) then
1991 deallocate(constants % nl, stat=istat)
1992 if (istat .ne. 0) then
1993 call update_error(ddx_error, "`nl` deallocation failed!")
1994 end if
1995 end if
1996 if (allocated(constants % itrnl)) then
1997 deallocate(constants % itrnl, stat=istat)
1998 if (istat .ne. 0) then
1999 call update_error(ddx_error, "`itrnl` deallocation failed!")
2000 end if
2001 end if
2002 if (allocated(constants % fi)) then
2003 deallocate(constants % fi, stat=istat)
2004 if (istat .ne. 0) then
2005 call update_error(ddx_error, "`fi` deallocation failed!")
2006 end if
2007 end if
2008 if (allocated(constants % ui)) then
2009 deallocate(constants % ui, stat=istat)
2010 if (istat .ne. 0) then
2011 call update_error(ddx_error, "`ui` deallocation failed!")
2012 end if
2013 end if
2014 if (allocated(constants % ui_cav)) then
2015 deallocate(constants % ui_cav, stat=istat)
2016 if (istat .ne. 0) then
2017 call update_error(ddx_error, "`ui_cav` deallocation failed!")
2018 end if
2019 end if
2020 if (allocated(constants % zi)) then
2021 deallocate(constants % zi, stat=istat)
2022 if (istat .ne. 0) then
2023 call update_error(ddx_error, "`zi` deallocation failed!")
2024 end if
2025 end if
2026 if (allocated(constants % zi_dr)) then
2027 deallocate(constants % zi_dr, stat=istat)
2028 if (istat .ne. 0) then
2029 call update_error(ddx_error, "`zi_dr` deallocation failed!")
2030 end if
2031 end if
2032 if (allocated(constants % ncav_sph)) then
2033 deallocate(constants % ncav_sph, stat=istat)
2034 if (istat .ne. 0) then
2035 call update_error(ddx_error, "`ncav_sph` deallocation failed!")
2036 end if
2037 end if
2038 if (allocated(constants % ccav)) then
2039 deallocate(constants % ccav, stat=istat)
2040 if (istat .ne. 0) then
2041 call update_error(ddx_error, "`ccav` deallocation failed!")
2042 end if
2043 end if
2044 if (allocated(constants % icav_ia)) then
2045 deallocate(constants % icav_ia, stat=istat)
2046 if (istat .ne. 0) then
2047 call update_error(ddx_error, "`icav_ia` deallocation failed!")
2048 end if
2049 end if
2050 if (allocated(constants % icav_ja)) then
2051 deallocate(constants % icav_ja, stat=istat)
2052 if (istat .ne. 0) then
2053 call update_error(ddx_error, "`icav_ja` deallocation failed!")
2054 end if
2055 end if
2056 if (allocated(constants % rx_prc)) then
2057 deallocate(constants % rx_prc, stat=istat)
2058 if (istat .ne. 0) then
2059 call update_error(ddx_error, "`rx_prc` deallocation failed!")
2060 end if
2061 end if
2062 if (allocated(constants % order)) then
2063 deallocate(constants % order, stat=istat)
2064 if (istat .ne. 0) then
2065 call update_error(ddx_error, "`order` deallocation failed!")
2066 end if
2067 end if
2068 if (allocated(constants % cluster)) then
2069 deallocate(constants % cluster, stat=istat)
2070 if (istat .ne. 0) then
2071 call update_error(ddx_error, "`cluster` deallocation failed!")
2072 end if
2073 end if
2074 if (allocated(constants % children)) then
2075 deallocate(constants % children, stat=istat)
2076 if (istat .ne. 0) then
2077 call update_error(ddx_error, "`children` deallocation failed!")
2078 end if
2079 end if
2080 if (allocated(constants % parent)) then
2081 deallocate(constants % parent, stat=istat)
2082 if (istat .ne. 0) then
2083 call update_error(ddx_error, "`parent` deallocation failed!")
2084 end if
2085 end if
2086 if (allocated(constants % cnode)) then
2087 deallocate(constants % cnode, stat=istat)
2088 if (istat .ne. 0) then
2089 call update_error(ddx_error, "`cnode` deallocation failed!")
2090 end if
2091 end if
2092 if (allocated(constants % rnode)) then
2093 deallocate(constants % rnode, stat=istat)
2094 if (istat .ne. 0) then
2095 call update_error(ddx_error, "`rnode` deallocation failed!")
2096 end if
2097 end if
2098 if (allocated(constants % snode)) then
2099 deallocate(constants % snode, stat=istat)
2100 if (istat .ne. 0) then
2101 call update_error(ddx_error, "`snode` deallocation failed!")
2102 end if
2103 end if
2104 if (allocated(constants % sk_rnode)) then
2105 deallocate(constants % sk_rnode, stat=istat)
2106 if (istat .ne. 0) then
2107 call update_error(ddx_error, "`sk_rnode` deallocation failed!")
2108 end if
2109 end if
2110 if (allocated(constants % si_rnode)) then
2111 deallocate(constants % si_rnode, stat=istat)
2112 if (istat .ne. 0) then
2113 call update_error(ddx_error, "`si_rnode` deallocation failed!")
2114 end if
2115 end if
2116 if (allocated(constants % nfar)) then
2117 deallocate(constants % nfar, stat=istat)
2118 if (istat .ne. 0) then
2119 call update_error(ddx_error, "`nfar` deallocation failed!")
2120 end if
2121 end if
2122 if (allocated(constants % nnear)) then
2123 deallocate(constants % nnear, stat=istat)
2124 if (istat .ne. 0) then
2125 call update_error(ddx_error, "`nnear` deallocation failed!")
2126 end if
2127 end if
2128 if (allocated(constants % far)) then
2129 deallocate(constants % far, stat=istat)
2130 if (istat .ne. 0) then
2131 call update_error(ddx_error, "`far` deallocation failed!")
2132 end if
2133 end if
2134 if (allocated(constants % near)) then
2135 deallocate(constants % near, stat=istat)
2136 if (istat .ne. 0) then
2137 call update_error(ddx_error, "`near` deallocation failed!")
2138 end if
2139 end if
2140 if (allocated(constants % sfar)) then
2141 deallocate(constants % sfar, stat=istat)
2142 if (istat .ne. 0) then
2143 call update_error(ddx_error, "`sfar` deallocation failed!")
2144 end if
2145 end if
2146 if (allocated(constants % snear)) then
2147 deallocate(constants % snear, stat=istat)
2148 if (istat .ne. 0) then
2149 call update_error(ddx_error, "`snear` deallocation failed!")
2150 end if
2151 end if
2152 if (allocated(constants % m2l_ztranslate_coef)) then
2153 deallocate(constants % m2l_ztranslate_coef, stat=istat)
2154 if (istat .ne. 0) then
2155 call update_error(ddx_error, "`m2l_ztranslate_coef` " // &
2156 & "deallocation failed!")
2157 end if
2158 end if
2159 if (allocated(constants % m2l_ztranslate_adj_coef)) then
2160 deallocate(constants % m2l_ztranslate_adj_coef, stat=istat)
2161 if (istat .ne. 0) then
2162 call update_error(ddx_error, "`m2l_ztranslate_adj_coef` " // &
2163 & "deallocation failed!")
2164 end if
2165 end if
2166end subroutine constants_free
2167
2168end module ddx_constants
2169
ddx_constants
Run-time constants.
Definition: ddx_constants.f90:17
ddx_constants::build_itrnl
subroutine build_itrnl(constants, params, ddx_error)
Build the transposed neighbor list.
Definition: ddx_constants.f90:559
ddx_constants::build_b
subroutine build_b(constants, params, ddx_error)
Allocate and build the HSP sparse matrix, only if incore is set.
Definition: ddx_constants.f90:656
ddx_constants::tree_get_farnear_work
subroutine tree_get_farnear_work(n, children, cnode, rnode, lwork, iwork, jwork, work, nnfar, nfar, nnnear, nnear)
Find near and far admissible pairs of tree nodes and store it in work array.
Definition: ddx_constants.f90:1698
ddx_constants::neighbor_list_init
subroutine neighbor_list_init(params, constants, ddx_error)
Build the neighbor list using a N^2 code.
Definition: ddx_constants.f90:1088
ddx_constants::mkprec
subroutine mkprec(lmax, nbasis, nsph, ngrid, eps, ui, wgrid, vgrid, vgrid_nbasis, rx_prc, ddx_error)
Compute preconditioner assemble the diagonal blocks of the reps matrix then invert them to build the ...
Definition: ddx_constants.f90:1372
ddx_constants::build_l
subroutine build_l(constants, params, ddx_error)
Allocate and build the ddCOSMO sparse matrix, only if incore is set.
Definition: ddx_constants.f90:584
ddx_constants::dfsw
real(dp) function dfsw(t, se, eta)
Derivative of a switching function.
Definition: ddx_constants.f90:1345
ddx_constants::constants_init
subroutine constants_init(params, constants, ddx_error)
Compute all necessary constants.
Definition: ddx_constants.f90:240
ddx_constants::neighbor_list_init_fmm
subroutine neighbor_list_init_fmm(params, constants, ddx_error)
Build the neighbor list using a linear scaling code (only if FMMs are enabled)
Definition: ddx_constants.f90:1173
ddx_constants::constants_geometry_init
subroutine constants_geometry_init(params, constants, ddx_error)
Initialize geometry-related constants like list of neighbouring spheres.
Definition: ddx_constants.f90:771
ddx_constants::constants_free
subroutine constants_free(constants, ddx_error)
Deallocate the constants.
Definition: ddx_constants.f90:1849
ddx_constants::tree_rib_build
subroutine tree_rib_build(nsph, csph, rsph, order, cluster, children, parent, cnode, rnode, snode, ddx_error)
Build a recursive inertial binary tree.
Definition: ddx_constants.f90:1463
ddx_constants::fsw
real(dp) function fsw(t, se, eta)
Switching function.
Definition: ddx_constants.f90:1290
ddx_constants::tree_rib_node_bisect
subroutine tree_rib_node_bisect(nsph, csph, n, order, div, ddx_error)
Divide given cluster of spheres into two subclusters by inertial bisection.
Definition: ddx_constants.f90:1589
ddx_constants::mkpmat
subroutine mkpmat(params, constants, isph, pmat)
Computation of P_chi.
Definition: ddx_constants.f90:735
ddx_constants::tree_get_farnear
subroutine tree_get_farnear(jwork, lwork, work, n, nnfar, nfar, sfar, far, nnnear, nnear, snear, near)
Get near and far admissible pairs from work array of tree_get_farnear_work Works only for binary tree...
Definition: ddx_constants.f90:1791
ddx_harmonics
Harmonics-related core routines.
Definition: ddx_harmonics.f90:16
ddx_harmonics::ylmscale
subroutine ylmscale(p, vscales, v4pi2lp1, vscales_rel)
Compute scaling factors of real normalized spherical harmonics.
Definition: ddx_harmonics.f90:38
ddx_harmonics::ylmbas
subroutine ylmbas(x, rho, ctheta, stheta, cphi, sphi, p, vscales, vylm, vplm, vcos, vsin)
Compute all spherical harmonics up to a given degree at a given point.
Definition: ddx_harmonics.f90:216
ddx_harmonics::fmm_constants
subroutine fmm_constants(dmax, pm, pl, vcnk, m2l_ztranslate_coef, m2l_ztranslate_adj_coef)
Compute FMM-related constants.
Definition: ddx_harmonics.f90:82
ddx_parameters
Module to treat properly user input parameters.
Definition: ddx_parameters.f90:15
ddx_constants::ddx_constants_type
Container for precomputed constants.
Definition: ddx_constants.f90:28
ddx_parameters::ddx_params_type
Type to check and store user input parameters.
Definition: ddx_parameters.f90:25