ddx 1.0.2
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
26type ddx_switching_type
27 integer :: ncav
30 real(dp), allocatable :: u_ni(:,:)
33 real(dp), allocatable :: u_i_cav(:)
35 real(dp), allocatable :: d_ni(:,:)
37 real(dp), allocatable :: f_ni(:,:)
40 real(dp), allocatable :: grad_u_ni(:,:,:)
41 real(dp), allocatable :: dr_u_ni(:,:)
42end type ddx_switching_type
43
45type ddx_constants_type
47 integer :: nbasis
50 integer :: n
57 integer :: dmax
59 integer :: nscales
62 real(dp), allocatable :: vscales(:)
65 real(dp), allocatable :: v4pi2lp1(:)
70 real(dp), allocatable :: vscales_rel(:)
74 integer :: nfact
76 real(dp), allocatable :: vfact(:)
80 real(dp), allocatable :: vcnk(:)
84 real(dp), allocatable :: m2l_ztranslate_coef(:, :, :)
88 real(dp), allocatable :: m2l_ztranslate_adj_coef(:, :, :)
90 real(dp), allocatable :: cgrid(:, :)
92 real(dp), allocatable :: wgrid(:)
98 integer :: vgrid_dmax
100 integer :: vgrid_nbasis
103 real(dp), allocatable :: vgrid(:, :)
107 real(dp), allocatable :: vwgrid(:, :)
109 real(dp), allocatable :: vgrid2(:, :)
111 integer :: lmax0
113 integer :: nbasis0
116 real(dp), allocatable :: pchi(:, :, :)
119 real(dp), allocatable :: c_ik(:, :)
121 real(dp), allocatable :: si_ri(:, :)
124 real(dp), allocatable :: di_ri(:, :)
126 real(dp), allocatable :: sk_ri(:, :)
129 real(dp), allocatable :: dk_ri(:, :)
131 real(dp), allocatable :: termimat(:, :)
133 real(dp), allocatable :: b(:,:,:)
135 real(dp), allocatable :: l(:,:,:)
138 integer :: nngmax
140 integer, allocatable :: inl(:)
143 integer, allocatable :: nl(:)
145 integer, allocatable :: itrnl(:)
148 real(dp), allocatable :: fi(:, :)
151 real(dp), allocatable :: ui(:, :)
154 real(dp), allocatable :: ui_cav(:)
157 real(dp), allocatable :: zi(:, :, :)
160 real(dp), allocatable :: zi_dr(:, :)
162 integer :: ncav
164 integer, allocatable :: ncav_sph(:)
166 real(dp), allocatable :: ccav(:, :)
168 integer, allocatable :: icav_ia(:)
170 integer, allocatable :: icav_ja(:)
173 real(dp), allocatable :: rx_prc(:, :, :)
174 !! Cluster tree information that is allocated and computed only if fmm=1.
177 integer, allocatable :: order(:)
179 integer :: nclusters
182 integer, allocatable :: cluster(:, :)
185 integer, allocatable :: children(:, :)
188 integer, allocatable :: parent(:)
191 real(dp), allocatable :: cnode(:, :)
194 real(dp), allocatable :: rnode(:)
197 integer, allocatable :: snode(:)
200 real(dp), allocatable :: sk_rnode(:, :)
203 real(dp), allocatable :: si_rnode(:, :)
205 integer :: nnfar
207 integer :: nnnear
210 integer, allocatable :: nfar(:)
213 integer, allocatable :: nnear(:)
216 integer, allocatable :: far(:)
219 integer, allocatable :: near(:)
223 integer, allocatable :: sfar(:)
227 integer, allocatable :: snear(:)
230 integer :: nnear_m2p
233 integer :: m2p_lmax
236 integer :: m2p_nbasis
240 integer :: grad_nbasis
242 real(dp) :: inner_tol
245 logical :: dodiag
247 type(ddx_switching_type) :: switching
248
249end type ddx_constants_type
250
251contains
252
259subroutine constants_init(params, constants, ddx_error)
260 use complex_bessel
261 !! Inputs
262 type(ddx_params_type), intent(in) :: params
263 !! Outputs
264 type(ddx_constants_type), intent(out) :: constants
265 type(ddx_error_type), intent(inout) :: ddx_error
266 !! Local variables
267 integer :: i, alloc_size, l, indl, igrid, isph, l0, &
268 & NZ, ierr, info, tmp_pmax
269 real(dp) :: rho, ctheta, stheta, cphi, sphi, termi, termk, &
270 & s1, s2
271 real(dp), allocatable :: vplm(:), vcos(:), vsin(:), vylm(:), SK_rijn(:), &
272 & DK_rijn(:)
273 complex(dp), allocatable :: bessel_work(:)
274 complex(dp) :: z
275 !! The code
276 ! Clean ddx_error state of constants to proceed with geometry
277 if (ddx_error % flag .ne. 0) then
278 call update_error(ddx_error, "constants_init received input in error " // &
279 & " state, exiting")
280 return
281 end if
282 ! activate inner iterations diagonal in mvp for debugging purposes only.
283 ! could be useful for different linear solvers.
284 constants % dodiag = .false.
285 ! Maximal number of modeling spherical harmonics
286 constants % nbasis = (params % lmax+1) ** 2
287 ! Maximal number of modeling degrees of freedom
288 constants % n = params % nsph * constants % nbasis
289 ! Calculate dmax, vgrid_dmax, m2p_lmax, m2p_nbasis and grad_nbasis
290 if (params % fmm .eq. 0) then
291 if (params % force .eq. 1) then
292 constants % dmax = params % lmax + 1
293 constants % vgrid_dmax = params % lmax + 1
294 else
295 constants % dmax = params % lmax
296 constants % vgrid_dmax = params % lmax
297 end if
298 ! Other constants are not referenced if fmm=0
299 constants % m2p_lmax = -1
300 constants % m2p_nbasis = -1
301 constants % grad_nbasis = -1
302 else
303 ! If forces are required then we need the M2P of a degree lmax+1 for
304 ! the near-field analytical gradients
305 if (params % force .eq. 1) then
306 constants % dmax = max(params % pm+params % pl+1, &
307 & params % lmax+1)
308 constants % m2p_lmax = params % lmax + 1
309 constants % grad_nbasis = (params % lmax+2) ** 2
310 constants % vgrid_dmax = max(params % pl, params % lmax) + 1
311 else
312 constants % dmax = max(params % pm+params % pl, &
313 & params % lmax)
314 constants % m2p_lmax = params % lmax
315 constants % grad_nbasis = -1
316 constants % vgrid_dmax = max(params % pl, params % lmax)
317 end if
318 constants % m2p_nbasis = (constants % m2p_lmax+1) ** 2
319 end if
320 ! Compute sizes of vgrid, vfact and vscales
321 constants % vgrid_nbasis = (constants % vgrid_dmax+1) ** 2
322 constants % nfact = 2*constants % dmax+1
323 constants % nscales = (constants % dmax+1) ** 2
324 ! Allocate space for scaling factors of spherical harmonics
325 allocate(constants % vscales(constants % nscales), stat=info)
326 if (info .ne. 0) then
327 call update_error(ddx_error, "constants_init: `vscales` allocation " &
328 & // "failed")
329 return
330 end if
331 allocate(constants % v4pi2lp1(constants % dmax+1), stat=info)
332 if (info .ne. 0) then
333 call update_error(ddx_error, "constants_init: `v4pi2lp1` allocation " &
334 & // "failed")
335 return
336 end if
337 allocate(constants % vscales_rel(constants % nscales), stat=info)
338 if (info .ne. 0) then
339 call update_error(ddx_error, "constants_init: `vscales_rel` " // &
340 & "allocation failed")
341 return
342 end if
343 ! Compute scaling factors of spherical harmonics
344 call ylmscale(constants % dmax, constants % vscales, &
345 & constants % v4pi2lp1, constants % vscales_rel)
346 ! Allocate square roots of factorials
347 allocate(constants % vfact(constants % nfact), stat=info)
348 if (info .ne. 0) then
349 call update_error(ddx_error, "constants_init: `vfact` allocation failed")
350 return
351 end if
352 ! Compute square roots of factorials
353 constants % vfact(1) = 1
354 do i = 2, constants % nfact
355 constants % vfact(i) = constants % vfact(i-1) * sqrt(dble(i-1))
356 end do
357 ! Allocate square roots of combinatorial numbers C_n^k and M2L OZ
358 ! translation coefficients
359 if (params % fmm .eq. 1 .or. &
360 & (params % model .eq. 3 .and. params % force .eq. 1)) then
361 alloc_size = 2*constants % dmax + 1
362 alloc_size = alloc_size * (constants % dmax+1)
363 ! Allocate C_n^k
364 allocate(constants % vcnk(alloc_size), stat=info)
365 if (info .ne. 0) then
366 call update_error(ddx_error, "constants_init: `vcnk` allocation " &
367 & // "failed")
368 return
369 end if
370 ! Allocate M2L OZ coefficients
371 allocate(constants % m2l_ztranslate_coef(&
372 & (params % pm+1), (params % pl+1), (params % pl+1)), &
373 & stat=info)
374 if (info .ne. 0) then
375 call update_error(ddx_error, "constants_init: " // &
376 & "`m2l_ztranslate_coef` allocation failed")
377 return
378 end if
379 ! Allocate adjoint M2L OZ coefficients
380 allocate(constants % m2l_ztranslate_adj_coef(&
381 & (params % pl+1), (params % pl+1), (params % pm+1)), &
382 & stat=info)
383 if (info .ne. 0) then
384 call update_error(ddx_error, "constants_init: " // &
385 & "`m2l_ztranslate_adj_coef` allocation failed")
386 return
387 end if
388 ! Compute combinatorial numbers C_n^k and M2L OZ translate coefficients
389 call fmm_constants(constants % dmax, params % pm, &
390 & params % pl, constants % vcnk, &
391 & constants % m2l_ztranslate_coef, &
392 & constants % m2l_ztranslate_adj_coef)
393 end if
394 ! Allocate space for Lebedev grid coordinates and weights
395 allocate(constants % cgrid(3, params % ngrid), &
396 & constants % wgrid(params % ngrid), stat=info)
397 if (info .ne. 0) then
398 call update_error(ddx_error, "constants_init: `cgrid` and `wgrid` " // &
399 & "allocations failed")
400 return
401 end if
402 ! Get weights and coordinates of Lebedev grid points
403 call llgrid(params % ngrid, constants % wgrid, constants % cgrid)
404 ! Allocate space for values of non-weighted and weighted spherical
405 ! harmonics and L2P at Lebedev grid points
406 allocate(constants % vgrid(constants % vgrid_nbasis, params % ngrid), &
407 & constants % vwgrid(constants % vgrid_nbasis, params % ngrid), &
408 & stat=info)
409 if (info .ne. 0) then
410 call update_error(ddx_error, "constants_init: `vgrid`, `wgrid` and" // &
411 & " allocations failed")
412 return
413 end if
414 allocate(vplm(constants % vgrid_nbasis), vcos(constants % vgrid_dmax+1), &
415 & vsin(constants % vgrid_dmax+1), stat=info)
416 if (info .ne. 0) then
417 call update_error(ddx_error, "constants_init: `vplm`, `vcos` and " // &
418 & "`vsin` allocations failed")
419 return
420 end if
421 ! Compute non-weighted and weighted spherical harmonics and the single
422 ! layer operator at Lebedev grid points
423 do igrid = 1, params % ngrid
424 call ylmbas(constants % cgrid(:, igrid), rho, ctheta, stheta, cphi, &
425 & sphi, constants % vgrid_dmax, constants % vscales, &
426 & constants % vgrid(:, igrid), vplm, vcos, vsin)
427 constants % vwgrid(:, igrid) = constants % wgrid(igrid) * &
428 & constants % vgrid(:, igrid)
429 end do
430 if (params % fmm .eq. 1) then
431 allocate(constants % vgrid2( &
432 & constants % vgrid_nbasis, params % ngrid), &
433 & stat=info)
434 if (info .ne. 0) then
435 call update_error(ddx_error, "constants_init: `vgrid2` " // &
436 & "allocation failed")
437 return
438 end if
439 do igrid = 1, params % ngrid
440 do l = 0, constants % vgrid_dmax
441 indl = l*l + l + 1
442 constants % vgrid2(indl-l:indl+l, igrid) = &
443 & constants % vgrid(indl-l:indl+l, igrid) / &
444 & constants % vscales(indl)**2
445 end do
446 end do
447 end if
448
449 ! Generate geometry-related constants (required by the LPB code)
450 call constants_geometry_init(params, constants, ddx_error)
451 if (ddx_error % flag .ne. 0) then
452 call update_error(ddx_error, "constants_geometry_init returned an " // &
453 & "error, exiting")
454 return
455 end if
456
457 ! Precompute LPB-related constants
458 if (params % model .eq. 3) then
459 constants % lmax0 = params % lmax
460 constants % nbasis0 = constants % nbasis
461 allocate(vylm(constants % vgrid_nbasis), &
462 & sk_rijn(0:constants % lmax0), dk_rijn(0:constants % lmax0), &
463 & bessel_work(constants % dmax+2), stat=info)
464 if (info .ne. 0) then
465 call update_error(ddx_error, "constants_init: `vylm`, `SK_rijn` " // &
466 & "and `DK_rijn` allocations failed")
467 return
468 end if
469 allocate(constants % SI_ri(0:constants % dmax+1, params % nsph), &
470 & constants % DI_ri(0:constants % dmax+1, params % nsph), &
471 & constants % SK_ri(0:params % lmax+1, params % nsph), &
472 & constants % DK_ri(0:params % lmax+1, params % nsph), &
473 & constants % Pchi(constants % nbasis, constants % nbasis0, params % nsph), &
474 & constants % C_ik(0:params % lmax, params % nsph), &
475 & constants % termimat(0:params % lmax, params % nsph), &
476 & stat=info)
477 if (info .ne. 0) then
478 call update_error(ddx_error, "constants_init: allocation failed")
479 return
480 end if
481 sk_rijn = zero
482 dk_rijn = zero
483 do isph = 1, params % nsph
484 ! We compute Bessel functions of degrees 0..lmax+1 because the
485 ! largest degree is required for forces
486 call modified_spherical_bessel_first_kind(constants % dmax+1, &
487 & params % rsph(isph)*params % kappa,&
488 & constants % SI_ri(:, isph), constants % DI_ri(:, isph), &
489 & bessel_work)
490 call modified_spherical_bessel_second_kind(params % lmax+1, &
491 & params % rsph(isph)*params % kappa, &
492 & constants % SK_ri(:, isph), constants % DK_ri(:, isph), &
493 & bessel_work)
494 ! Compute matrix PU_i^e(x_in)
495 ! Previous implementation in update_rhs. Made it in allocate_model, so as to use
496 ! it in Forces as well.
497 call mkpmat(params, constants, isph, constants % Pchi(:, :, isph))
498 ! Compute i'_l(r_i)/i_l(r_i)
499 do l = 0, params % lmax
500 constants % termimat(l, isph) = constants % DI_ri(l, isph) / &
501 & constants % SI_ri(l, isph) * params % kappa
502 end do
503 ! Compute (i'_l0/i_l0 - k'_l0/k_l0)^(-1) is computed in Eq.(97)
504 do l0 = 0, constants % lmax0
505 termi = constants % DI_ri(l0, isph) / &
506 & constants % SI_ri(l0, isph) * params % kappa
507 termk = constants % DK_ri(l0, isph)/ &
508 & constants % SK_ri(l0, isph) * params % kappa
509 constants % C_ik(l0, isph) = one / (termi-termk)
510 end do
511 end do
512 if (params % fmm .eq. 1) then
513 tmp_pmax = max(params % pm, params % pl)
514 allocate(constants % SI_rnode(tmp_pmax+1, constants % nclusters), &
515 & constants % SK_rnode(tmp_pmax+1, constants % nclusters), &
516 & stat=info)
517 if (info.ne.0) then
518 call update_error(ddx_error, "constants_init: `si_rnode, " // &
519 & "sk_rnode allocation failed")
520 return
521 end if
522 do i = 1, constants % nclusters
523 z = constants % rnode(i) * params % kappa
524 s1 = sqrt(two/(pi*real(z)))
525 s2 = sqrt(pi/(two*real(z)))
526 call cbesk(z, pt5, 1, 1, bessel_work(1), nz, ierr)
527 constants % SK_rnode(1, i) = s1 * real(bessel_work(1))
528 call cbesi(z, pt5, 1, 1, bessel_work(1), nz, ierr)
529 constants % SI_rnode(1, i) = s2 * real(bessel_work(1))
530 if (tmp_pmax .gt. 0) then
531 call cbesk(z, 1.5d0, 1, tmp_pmax, bessel_work(2:), nz, ierr)
532 constants % SK_rnode(2:, i) = s1 * real(bessel_work(2:tmp_pmax+1))
533 call cbesi(z, 1.5d0, 1, tmp_pmax, bessel_work(2:), nz, ierr)
534 constants % SI_rnode(2:, i) = s2 * real(bessel_work(2:tmp_pmax+1))
535 end if
536 end do
537 end if
538 deallocate(vylm, sk_rijn, dk_rijn, bessel_work, stat=info)
539 if (info .ne. 0) then
540 call update_error(ddx_error, "constants_init: `vylm`, `SK_rijn` " // &
541 & "and `DK_rijn` deallocations failed")
542 return
543 end if
544 ! if doing incore build nonzero blocks of B
545 if (params % matvecmem .eq. 1) then
546 call build_b(constants, params, ddx_error)
547 if (ddx_error % flag .ne. 0) then
548 call update_error(ddx_error, "build_b returned an " // &
549 & "error, exiting")
550 return
551 end if
552 end if
553 end if
554 deallocate(vplm, vcos, vsin, stat=info)
555 if (info .ne. 0) then
556 call update_error(ddx_error, "constants_init: `vplm`, `vcos` and " // &
557 & "`vsin` deallocations failed")
558 return
559 end if
560 ! if doing incore build nonzero blocks of L
561 if (params % matvecmem .eq. 1) then
562 call build_itrnl(constants, params, ddx_error)
563 if (ddx_error % flag .ne. 0) then
564 call update_error(ddx_error, "build_itrnl returned an " // &
565 & "error, exiting")
566 return
567 end if
568 call build_l(constants, params, ddx_error)
569 if (ddx_error % flag .ne. 0) then
570 call update_error(ddx_error, "build_l returned an " // &
571 & "error, exiting")
572 return
573 end if
574 end if
575
576end subroutine constants_init
577
579subroutine build_itrnl(constants, params, ddx_error)
580 implicit none
581 type(ddx_params_type), intent(in) :: params
582 type(ddx_constants_type), intent(inout) :: constants
583 type(ddx_error_type), intent(inout) :: ddx_error
584 integer :: isph, ij, jsph, ji, istat
585
586 allocate(constants % itrnl(constants % inl(params % nsph + 1)), stat=istat)
587 if (istat.ne.0) then
588 call update_error(ddx_error, 'Allocation failed in build_itrnl')
589 return
590 end if
591
592 do isph = 1, params % nsph
593 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
594 jsph = constants % nl(ij)
595 do ji = constants % inl(jsph), constants % inl(jsph + 1) - 1
596 if (constants % nl(ji) .eq. isph) exit
597 end do
598 constants % itrnl(ij) = ji
599 end do
600 end do
601end subroutine build_itrnl
602
604subroutine build_l(constants, params, ddx_error)
605 implicit none
606 type(ddx_params_type), intent(in) :: params
607 type(ddx_constants_type), intent(inout) :: constants
608 type(ddx_error_type), intent(inout) :: ddx_error
609 integer :: isph, ij, jsph, igrid, l, m, ind, info
610 real(dp), dimension(3) :: vij, sij
611 real(dp) :: vvij, tij, xij, oij, rho, ctheta, stheta, cphi, sphi, &
612 & fac, tt, thigh, d_i, f_i
613 real(dp), dimension(constants % nbasis) :: vylm, vplm
614 real(dp), dimension(params % lmax + 1) :: vcos, vsin
615 real(dp), dimension(constants % nbasis, params % ngrid) :: scratch
616
617 allocate(constants % l(constants % nbasis, constants % nbasis, &
618 & constants % inl(params % nsph + 1)), stat=info)
619 if (info .ne. 0) then
620 call update_error(ddx_error, 'Allocation failed in build_l')
621 return
622 end if
623
624 thigh = one + pt5*(params % se + one)*params % eta
625
626 if (params%switching.eq.0) then
627 !$omp parallel do default(none) shared(params,constants,thigh) &
628 !$omp private(isph,ij,jsph,scratch,igrid,vij,vvij,tij,sij,xij,oij, &
629 !$omp rho,ctheta,stheta,cphi,sphi,vylm,vplm,vcos,vsin,l,fac,ind,m,tt)
630 do isph = 1, params % nsph
631 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
632 jsph = constants % nl(ij)
633 scratch = zero
634 do igrid = 1, params % ngrid
635 if (constants % ui(igrid, isph).eq.one) cycle
636 vij = params % csph(:, isph) &
637 & + params % rsph(isph)*constants % cgrid(:,igrid) &
638 & - params % csph(:, jsph)
639 vvij = sqrt(dot_product(vij, vij))
640 tij = vvij/params % rsph(jsph)
641 if (tij.lt.thigh) then
642 if (tij.ne.zero) then
643 sij = vij/vvij
644 else
645 sij = one
646 end if
647 xij = fsw(tij, params % se, params % eta)
648 if (constants % fi(igrid, isph).gt.one) then
649 oij = xij/constants % fi(igrid, isph)
650 else
651 oij = xij
652 end if
653 call ylmbas(sij, rho, ctheta, stheta, cphi, sphi, &
654 & params % lmax, constants % vscales, vylm, vplm, &
655 & vcos, vsin)
656 tt = oij
657 do l = 0, params % lmax
658 ind = l*l + l + 1
659 fac = - tt/(constants % vscales(ind)**2)
660 do m = -l, l
661 scratch(ind + m, igrid) = fac*vylm(ind + m)
662 end do
663 tt = tt*tij
664 end do
665 end if
666 end do
667 call dgemm('n', 't', constants % nbasis, constants % nbasis, params % ngrid, &
668 & one, constants % vwgrid, constants % vgrid_nbasis, scratch, &
669 & constants % nbasis, zero, constants % l(:,:,ij), constants % nbasis)
670 end do
671 end do
672 else
673 !$omp parallel do default(none) shared(params,constants,thigh) &
674 !$omp private(isph,ij,jsph,scratch,igrid,vij,vvij,tij,sij,xij, &
675 !$omp oij,rho,ctheta,stheta,cphi,sphi,vylm,vplm,vcos,vsin,l, &
676 !$omp d_i,f_i,fac,ind,m,tt)
677 do isph = 1, params % nsph
678 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
679 jsph = constants % nl(ij)
680 scratch = zero
681 do igrid = 1, params % ngrid
682 d_i = constants%switching%d_ni(igrid,isph)
683 f_i = constants%switching%f_ni(igrid,isph)
684 if (constants % ui(igrid, isph).eq.one) cycle
685 vij = params % csph(:, isph) &
686 & + params % rsph(isph)*constants % cgrid(:,igrid) &
687 & - params % csph(:, jsph)
688 vvij = sqrt(dot_product(vij, vij))
689 tij = vvij/params % rsph(jsph)
690 if (tij.lt.thigh) then
691 if (tij.ne.zero) then
692 sij = vij/vvij
693 else
694 sij = one
695 end if
696 xij = fsw(tij, params % se, params % eta)
697 oij = xij/(d_i + f_i)
698 call ylmbas(sij, rho, ctheta, stheta, cphi, sphi, &
699 & params % lmax, constants % vscales, vylm, vplm, &
700 & vcos, vsin)
701 tt = oij
702 do l = 0, params % lmax
703 ind = l*l + l + 1
704 fac = - tt/(constants % vscales(ind)**2)
705 do m = -l, l
706 scratch(ind + m, igrid) = fac*vylm(ind + m)
707 end do
708 tt = tt*tij
709 end do
710 end if
711 end do
712 call dgemm('n', 't', constants % nbasis, constants % nbasis, params % ngrid, &
713 & one, constants % vwgrid, constants % vgrid_nbasis, scratch, &
714 & constants % nbasis, zero, constants % l(:,:,ij), constants % nbasis)
715 end do
716 end do
717 end if
718end subroutine build_l
719
721subroutine build_b(constants, params, ddx_error)
722 implicit none
723 type(ddx_params_type), intent(in) :: params
724 type(ddx_constants_type), intent(inout) :: constants
725 type(ddx_error_type), intent(inout) :: ddx_error
726 integer :: isph, ij, jsph, igrid, l, m, ind
727 real(dp), dimension(3) :: vij, sij
728 real(dp) :: vvij, tij, xij, oij, rho, ctheta, stheta, cphi, sphi, &
729 & fac, vvtij, thigh
730 real(dp), dimension(constants % nbasis) :: vylm, vplm
731 real(dp), dimension(params % lmax + 1) :: vcos, vsin
732 complex(dp), dimension(max(2, params % lmax + 1)) :: bessel_work
733 real(dp), dimension(0:params % lmax) :: SI_rijn, DI_rijn
734 real(dp), dimension(constants % nbasis, params % ngrid) :: scratch
735 integer :: info
736 real(dp) :: f_i, d_i
737
738 allocate(constants % b(constants % nbasis, constants % nbasis, &
739 & constants % inl(params % nsph + 1)), stat=info)
740 if (info.ne.0) then
741 call update_error(ddx_error, "Allocation failed in build_b")
742 return
743 end if
744
745 thigh = one + pt5*(params % se + one)*params % eta
746
747 if (params%switching.eq.0) then
748 !$omp parallel do default(none) shared(params,constants,thigh) &
749 !$omp private(isph,ij,jsph,scratch,igrid,vij,vvij,tij,sij,xij,oij, &
750 !$omp rho,ctheta,stheta,cphi,sphi,vylm,vplm,vcos,vsin,si_rijn,di_rijn, &
751 !$omp vvtij,l,fac,ind,m,bessel_work)
752 do isph = 1, params % nsph
753 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
754 jsph = constants % nl(ij)
755 scratch = zero
756 do igrid = 1, params % ngrid
757 if (constants % ui(igrid, isph).eq.one) cycle
758 vij = params % csph(:, isph) &
759 & + params % rsph(isph)*constants % cgrid(:,igrid) &
760 & - params % csph(:, jsph)
761 vvij = sqrt(dot_product(vij, vij))
762 tij = vvij/params % rsph(jsph)
763 if (tij.lt.thigh) then
764 if (tij.ne.zero) then
765 sij = vij/vvij
766 else
767 sij = one
768 end if
769 xij = fsw(tij, params % se, params % eta)
770 if (constants % fi(igrid, isph).gt.one) then
771 oij = xij/constants % fi(igrid, isph)
772 else
773 oij = xij
774 end if
775 call ylmbas(sij, rho, ctheta, stheta, cphi, sphi, &
776 & params % lmax, constants % vscales, vylm, vplm, &
777 & vcos, vsin)
778 si_rijn = 0
779 di_rijn = 0
780 vvtij = vvij*params % kappa
781 call modified_spherical_bessel_first_kind(params % lmax, &
782 & vvtij, si_rijn, di_rijn, bessel_work)
783 do l = 0, params % lmax
784 fac = - oij*si_rijn(l)/constants % SI_ri(l, jsph)
785 ind = l*l + l + 1
786 do m = -l, l
787 scratch(ind + m, igrid) = fac*vylm(ind + m)
788 end do
789 end do
790 end if
791 end do
792 call dgemm('n', 't', constants % nbasis, constants % nbasis, params % ngrid, &
793 & one, constants % vwgrid, constants % vgrid_nbasis, scratch, &
794 & constants % nbasis, zero, constants % b(:,:,ij), constants % nbasis)
795 end do
796 end do
797 else
798 !$omp parallel do default(none) shared(params,constants,thigh) &
799 !$omp private(isph,ij,jsph,scratch,igrid,vij,vvij,tij,sij,xij,oij, &
800 !$omp rho,ctheta,stheta,cphi,sphi,vylm,vplm,vcos,vsin,si_rijn,di_rijn, &
801 !$omp vvtij,l,fac,ind,m,bessel_work,d_i,f_i)
802 do isph = 1, params % nsph
803 do ij = constants % inl(isph), constants % inl(isph + 1) - 1
804 jsph = constants % nl(ij)
805 scratch = zero
806 do igrid = 1, params % ngrid
807 d_i = constants%switching%d_ni(igrid,isph)
808 f_i = constants%switching%f_ni(igrid,isph)
809 if (constants % ui(igrid, isph).eq.one) cycle
810 vij = params % csph(:, isph) &
811 & + params % rsph(isph)*constants % cgrid(:,igrid) &
812 & - params % csph(:, jsph)
813 vvij = sqrt(dot_product(vij, vij))
814 tij = vvij/params % rsph(jsph)
815 if (tij.lt.thigh) then
816 if (tij.ne.zero) then
817 sij = vij/vvij
818 else
819 sij = one
820 end if
821 xij = fsw(tij, params % se, params % eta)
822 oij = xij/(d_i + f_i)
823 call ylmbas(sij, rho, ctheta, stheta, cphi, sphi, &
824 & params % lmax, constants % vscales, vylm, vplm, &
825 & vcos, vsin)
826 si_rijn = 0
827 di_rijn = 0
828 vvtij = vvij*params % kappa
829 call modified_spherical_bessel_first_kind(params % lmax, &
830 & vvtij, si_rijn, di_rijn, bessel_work)
831 do l = 0, params % lmax
832 fac = - oij*si_rijn(l)/constants % SI_ri(l, jsph)
833 ind = l*l + l + 1
834 do m = -l, l
835 scratch(ind + m, igrid) = fac*vylm(ind + m)
836 end do
837 end do
838 end if
839 end do
840 call dgemm('n', 't', constants % nbasis, constants % nbasis, params % ngrid, &
841 & one, constants % vwgrid, constants % vgrid_nbasis, scratch, &
842 & constants % nbasis, zero, constants % b(:,:,ij), constants % nbasis)
843 end do
844 end do
845 end if
846end subroutine build_b
847
849subroutine mkpmat(params, constants, isph, pmat)
850 type(ddx_params_type), intent(in) :: params
851 type(ddx_constants_type), intent(in) :: constants
852 integer, intent(in) :: isph
853 real(dp), dimension(constants % nbasis, constants % nbasis0), intent(inout) :: pmat
854 integer :: l, m, ind, l0, m0, ind0, its
855 real(dp) :: f, f0
856 pmat(:,:) = zero
857 do its = 1, params % ngrid
858 if (constants % ui(its,isph).ne.0) then
859 do l = 0, params % lmax
860 ind = l*l + l + 1
861 do m = -l,l
862 f = constants % wgrid(its) * constants % vgrid(ind+m,its) * constants % ui(its,isph)
863 do l0 = 0, constants % lmax0
864 ind0 = l0*l0 + l0 + 1
865 do m0 = -l0, l0
866 f0 = constants % vgrid(ind0+m0,its)
867 pmat(ind+m,ind0+m0) = pmat(ind+m,ind0+m0) + f * f0
868 end do
869 end do
870 end do
871 end do
872 end if
873 end do
874end subroutine mkpmat
875
885subroutine constants_geometry_init(params, constants, ddx_error)
886 !! Inputs
887 type(ddx_params_type), intent(in) :: params
888 !! Outputs
889 type(ddx_constants_type), intent(inout) :: constants
890 type(ddx_error_type), intent(inout) :: ddx_error
891 !! Local variables
892 real(dp) :: swthr, v(3), maxv, ssqv, vv, t
893 integer :: i, isph, jsph, inear, igrid, iwork, jwork, lwork, &
894 & old_lwork, icav, info
895 integer, allocatable :: work(:, :), tmp_work(:, :)
896 real(dp) :: start_time
897 !! The code
898 ! Prepare FMM structures if needed
899 start_time = omp_get_wtime()
900 if (params % fmm .eq. 1) then
901 ! Allocate space for a cluster tree
902 allocate(constants % order(params % nsph), stat=info)
903 if (info .ne. 0) then
904 call update_error(ddx_error, "constants_geometry_init: `order` " &
905 & // "allocation failed")
906 return
907 end if
908 constants % nclusters = 2*params % nsph - 1
909 allocate(constants % cluster(2, constants % nclusters), stat=info)
910 if (info .ne. 0) then
911 call update_error(ddx_error, "constants_geometry_init: `cluster` " &
912 & // "allocation failed")
913 return
914 end if
915 allocate(constants % children(2, constants % nclusters), stat=info)
916 if (info .ne. 0) then
917 call update_error(ddx_error, "constants_geometry_init: " // &
918 & "`children` allocation failed")
919 return
920 endif
921 allocate(constants % parent(constants % nclusters), stat=info)
922 if (info .ne. 0) then
923 call update_error(ddx_error, "constants_geometry_init: `parent` " &
924 & // "allocation failed")
925 return
926 endif
927 allocate(constants % cnode(3, constants % nclusters), stat=info)
928 if (info .ne. 0) then
929 call update_error(ddx_error, "constants_geometry_init: `cnode` " &
930 & // "allocation failed")
931 return
932 endif
933 allocate(constants % rnode(constants % nclusters), stat=info)
934 if (info .ne. 0) then
935 call update_error(ddx_error, "constants_geometry_init: `rnode` " &
936 & // "allocation failed")
937 return
938 endif
939 allocate(constants % snode(params % nsph), stat=info)
940 if (info .ne. 0) then
941 call update_error(ddx_error, "constants_geometry_init: `snode` " &
942 & // "allocation failed")
943 return
944 endif
945 allocate(constants % nfar(constants % nclusters), stat=info)
946 if (info .ne. 0) then
947 call update_error(ddx_error, "constants_geometry_init: `nfar` " &
948 & // "allocation failed")
949 return
950 endif
951 allocate(constants % nnear(constants % nclusters), stat=info)
952 if (info .ne. 0) then
953 call update_error(ddx_error, "constants_geometry_init: `nnear` " &
954 & // "allocation failed")
955 return
956 endif
957 ! Get the tree
958 call tree_rib_build(params % nsph, params % csph, params % rsph, &
959 & constants % order, constants % cluster, constants % children, &
960 & constants % parent, constants % cnode, constants % rnode, &
961 & constants % snode, ddx_error)
962 if (ddx_error % flag .ne. 0) then
963 call update_error(ddx_error, "tree_rib_build returned an " // &
964 & "ddx_error, exiting")
965 return
966 end if
967 ! Get number of far and near admissible pairs
968 iwork = 0
969 jwork = 1
970 lwork = 1
971 allocate(work(3, lwork), stat=info)
972 if (info .ne. 0) then
973 call update_error(ddx_error, "constants_geometry_init: `work` " &
974 & // "allocation failed")
975 return
976 end if
977 do while (iwork .le. jwork)
978 allocate(tmp_work(3, lwork), stat=info)
979 if (info .ne. 0) then
980 call update_error(ddx_error, "constants_geometry_init: " // &
981 & "`tmp_work` allocation failed")
982 return
983 end if
984 tmp_work = work
985 deallocate(work, stat=info)
986 if (info .ne. 0) then
987 call update_error(ddx_error, "constants_geometry_init: " // &
988 & "`work` deallocation failed")
989 return
990 end if
991 old_lwork = lwork
992 lwork = old_lwork + 1000*params % nsph
993 allocate(work(3, lwork), stat=info)
994 if (info .ne. 0) then
995 call update_error(ddx_error, "constants_geometry_init: " // &
996 & "`work` allocation failed")
997 return
998 end if
999 work(:, 1:old_lwork) = tmp_work
1000 deallocate(tmp_work, stat=info)
1001 if (info .ne. 0) then
1002 call update_error(ddx_error, "constants_geometry_init: " // &
1003 & "`tmp_work` deallocation failed")
1004 return
1005 end if
1006 call tree_get_farnear_work(constants % nclusters, &
1007 & constants % children, constants % cnode, &
1008 & constants % rnode, lwork, iwork, jwork, work, &
1009 & constants % nnfar, constants % nfar, constants % nnnear, &
1010 & constants % nnear)
1011 end do
1012 allocate(constants % sfar(constants % nclusters+1), &
1013 & constants % snear(constants % nclusters+1), stat=info)
1014 if (info .ne. 0) then
1015 call update_error(ddx_error, "constants_geometry_init: `sfar` " // &
1016 & "and `snear` allocations failed")
1017 return
1018 end if
1019 allocate(constants % far(constants % nnfar), stat=info)
1020 if (info .ne. 0) then
1021 call update_error(ddx_error, "constants_geometry_init: `far` " // &
1022 & "allocation failed")
1023 return
1024 end if
1025 allocate(constants % near(constants % nnnear), stat=info)
1026 if (info .ne. 0) then
1027 call update_error(ddx_error, "constants_geometry_init: `near` " // &
1028 & "allocation failed")
1029 return
1030 end if
1031 call tree_get_farnear(jwork, lwork, work, constants % nclusters, &
1032 & constants % nnfar, constants % nfar, constants % sfar, &
1033 & constants % far, constants % nnnear, constants % nnear, &
1034 & constants % snear, constants % near)
1035 deallocate(work, stat=info)
1036 if (info .ne. 0) then
1037 call update_error(ddx_error, "constants_geometry_init: `work` " // &
1038 & "deallocation failed")
1039 return
1040 end if
1041 end if
1042 ! Upper bound of switch region. Defines intersection criterion for spheres
1043 swthr = one + (params % se+one)*params % eta/two
1044 ! Assemble neighbor list
1045 if (params % fmm .eq. 1) then
1046 call neighbor_list_init_fmm(params, constants, ddx_error)
1047 else
1048 call neighbor_list_init(params, constants, ddx_error)
1049 end if
1050 if (ddx_error % flag .ne. 0) then
1051 call update_error(ddx_error, "Neighbor list construction failed")
1052 return
1053 end if
1054 ! Allocate space for characteristic functions fi and ui
1055 allocate(constants % fi(params % ngrid, params % nsph), &
1056 & constants % ui(params % ngrid, params % nsph), stat=info)
1057 if (info .ne. 0) then
1058 call update_error(ddx_error, "`fi` and `ui` allocations failed")
1059 return
1060 end if
1061 constants % fi = zero
1062 constants % ui = zero
1063 ! Allocate space for force-related arrays
1064 if (params % force .eq. 1) then
1065 allocate(constants % zi(3, params % ngrid, params % nsph), stat=info)
1066 if (info .ne. 0) then
1067 call update_error(ddx_error, "`zi` allocation failed")
1068 return
1069 endif
1070 constants % zi = zero
1071
1072 allocate(constants % zi_dr(params % ngrid, params % nsph), stat=info)
1073 if (info .ne. 0) then
1074 call update_error(ddx_error, "`zi_dr` allocation failed")
1075 return
1076 endif
1077 constants % zi_dr = zero
1078 end if
1079
1080 ! Build arrays fi, ui, zi
1081 !$omp parallel do default(none) shared(params,constants,swthr) &
1082 !$omp private(isph,igrid,jsph,v,maxv,ssqv,vv,t) schedule(dynamic)
1083 do isph = 1, params % nsph
1084 do igrid = 1, params % ngrid
1085 ! Loop over neighbours of i-th sphere
1086 do inear = constants % inl(isph), constants % inl(isph+1)-1
1087 ! Neighbour's index
1088 jsph = constants % nl(inear)
1089 ! Compute t_n^ij
1090 v = params % csph(:, isph) - params % csph(:, jsph) + &
1091 & params % rsph(isph)*constants % cgrid(:, igrid)
1092 maxv = max(abs(v(1)), abs(v(2)), abs(v(3)))
1093 if (maxv .gt. zero) then
1094 ssqv = (v(1)/maxv)**2 + (v(2)/maxv)**2 + (v(3)/maxv)**2
1095 else
1096 ssqv = zero
1097 end if
1098 vv = maxv * sqrt(ssqv)
1099 t = vv / params % rsph(jsph)
1100 ! Accumulate characteristic function \chi(t_n^ij)
1101 constants % fi(igrid, isph) = constants % fi(igrid, isph) + &
1102 & fsw(t, params % se, params % eta)
1103 ! Check if gradients are required
1104 if (params % force .eq. 1) then
1105 ! Check if t_n^ij belongs to switch region
1106 if ((t .lt. swthr) .and. (t .gt. swthr-params % eta)) then
1107 ! Accumulate gradient of characteristic function \chi
1108 vv = dfsw(t, params % se, params % eta) / &
1109 & params % rsph(jsph) / vv
1110 constants % zi(:, igrid, isph) = &
1111 & constants % zi(:, igrid, isph) + vv*v
1112 constants % zi_dr(igrid, isph) = &
1113 & constants % zi_dr(igrid, isph) &
1114 & + vv*dot_product(v, constants % cgrid(:, igrid))
1115 end if
1116 end if
1117 enddo
1118 ! Compute characteristic function of a molecular surface ui
1119 if (constants % fi(igrid, isph) .le. one) then
1120 constants % ui(igrid, isph) = one - constants % fi(igrid, isph)
1121 end if
1122 enddo
1123 enddo
1124
1125 if (params%switching.eq.1) then
1126 call switching_init(params, constants, constants%switching, ddx_error)
1127 constants % ui = constants%switching%u_ni
1128 if (params%force.eq.1) then
1129 constants % zi = constants%switching%grad_u_ni
1130 constants % zi_dr = constants%switching%dr_u_ni
1131 end if
1132 end if
1133 if (ddx_error % flag .ne. 0) then
1134 call update_error(ddx_error, "switching_init returned an " // &
1135 & "error, exiting")
1136 return
1137 end if
1138
1139 ! Build cavity array. At first get total count for each sphere
1140 allocate(constants % ncav_sph(params % nsph), stat=info)
1141 if (info .ne. 0) then
1142 call update_error(ddx_error, "`ncav_sph` allocation failed")
1143 return
1144 endif
1145 !$omp parallel do default(none) shared(params,constants) &
1146 !$omp private(isph,i) schedule(dynamic)
1147 do isph = 1, params % nsph
1148 constants % ncav_sph(isph) = 0
1149 ! Loop over integration points
1150 do i = 1, params % ngrid
1151 ! Positive contribution from integration point
1152 if (constants % ui(i, isph) .gt. zero) then
1153 constants % ncav_sph(isph) = constants % ncav_sph(isph) + 1
1154 end if
1155 end do
1156 end do
1157 constants % ncav = sum(constants % ncav_sph)
1158 ! Allocate cavity array and CSR format for indexes of cavities
1159 allocate(constants % ccav(3, constants % ncav), &
1160 & constants % icav_ia(params % nsph+1), &
1161 & constants % icav_ja(constants % ncav), stat=info)
1162 if (info .ne. 0) then
1163 call update_error(ddx_error, "`ccav`, `icav_ia` and " // &
1164 & "`icav_ja` allocations failed")
1165 return
1166 endif
1167 ! Allocate space for characteristic functions ui at cavity points
1168 allocate(constants % ui_cav(constants % ncav), stat=info)
1169 if (info .ne. 0) then
1170 call update_error(ddx_error, "`ui_cav` allocations failed")
1171 return
1172 end if
1173 ! Get actual cavity coordinates and indexes in CSR format and fill in
1174 ! ui_cav aray
1175 constants % icav_ia(1) = 1
1176 icav = 1
1177 do isph = 1, params % nsph
1178 constants % icav_ia(isph+1) = constants % icav_ia(isph) + &
1179 & constants % ncav_sph(isph)
1180 ! Loop over integration points
1181 do igrid = 1, params % ngrid
1182 ! Ignore zero contribution
1183 if (constants % ui(igrid, isph) .eq. zero) cycle
1184 ! Store coordinates
1185 constants % ccav(:, icav) = params % csph(:, isph) + &
1186 & params % rsph(isph)* &
1187 & constants % cgrid(:, igrid)
1188 ! Store index
1189 constants % icav_ja(icav) = igrid
1190 ! Store characteristic function
1191 constants % ui_cav(icav) = constants % ui(igrid, isph)
1192 ! Advance cavity array index
1193 icav = icav + 1
1194 end do
1195 end do
1196 ! Compute diagonal preconditioner for PCM equations
1197 if (params % model .eq. 2) then
1198 allocate(constants % rx_prc( &
1199 & constants % nbasis, constants % nbasis, params % nsph), &
1200 & stat=info)
1201 if (info .ne. 0) then
1202 call update_error(ddx_error, "constants_geometry_init: `rx_prc` " &
1203 & // "allocation failed")
1204 return
1205 endif
1206 call mkprec(params % lmax, constants % nbasis, params % nsph, &
1207 & params % ngrid, params % eps, constants % ui, &
1208 & constants % wgrid, constants % vgrid, constants % vgrid_nbasis, &
1209 & constants % rx_prc, ddx_error)
1210 if (info .ne. 0) then
1211 call update_error(ddx_error, "mkprec returned an ddx_error, exiting")
1212 return
1213 end if
1214 end if
1215end subroutine constants_geometry_init
1216
1218subroutine neighbor_list_init(params, constants, ddx_error)
1219 type(ddx_params_type), intent(in) :: params
1220 type(ddx_constants_type), intent(inout) :: constants
1221 type(ddx_error_type), intent(inout) :: ddx_error
1222 real(dp) :: swthr, v(3), maxv, ssqv, vv, r
1223 integer :: nngmax, i, lnl, isph, jsph, info
1224 integer, allocatable :: tmp_nl(:)
1225 ! Upper bound of switch region. Defines intersection criterion for spheres
1226 swthr = (params % se+one)*params % eta/two
1227 ! Build list of neighbours in CSR format
1228 nngmax = 1
1229 allocate(constants % inl(params % nsph+1), &
1230 & constants % nl(params % nsph*nngmax), stat=info)
1231 if (info .ne. 0) then
1232 call update_error(ddx_error, "`inl` and `nl` allocations failed")
1233 return
1234 end if
1235 i = 1
1236 lnl = 0
1237 do isph = 1, params % nsph
1238 constants % inl(isph) = lnl + 1
1239 do jsph = 1, params % nsph
1240 if (isph .ne. jsph) then
1241 v = params % csph(:, isph)-params % csph(:, jsph)
1242 maxv = max(abs(v(1)), abs(v(2)), abs(v(3)))
1243 if (maxv .gt. zero) then
1244 ssqv = (v(1)/maxv)**2 + (v(2)/maxv)**2 + (v(3)/maxv)**2
1245 else
1246 ssqv = zero
1247 end if
1248 vv = maxv * sqrt(ssqv)
1249 ! Take regularization parameter into account with respect to
1250 ! shift se. It is described properly by the upper bound of a
1251 ! switch region `swthr`. Note that we take the largest
1252 ! switching region as we need a symmetric neighbor list to
1253 ! use it also for the adjoint products.
1254 r = params % rsph(isph) + params % rsph(jsph) &
1255 & + swthr * max(params % rsph(isph), params % rsph(jsph))
1256 if (vv .le. r) then
1257 constants % nl(i) = jsph
1258 i = i + 1
1259 lnl = lnl + 1
1260 ! Extend ddx_data % nl if needed
1261 if (i .gt. params % nsph*nngmax) then
1262 allocate(tmp_nl(params % nsph*nngmax), stat=info)
1263 if (info .ne. 0) then
1264 call update_error(ddx_error, "`tmp_nl` " // &
1265 & "allocation failed")
1266 return
1267 end if
1268 tmp_nl(1:params % nsph*nngmax) = &
1269 & constants % nl(1:params % nsph*nngmax)
1270 deallocate(constants % nl, stat=info)
1271 if (info .ne. 0) then
1272 call update_error(ddx_error, "`nl` " // &
1273 & "deallocation failed")
1274 return
1275 end if
1276 nngmax = nngmax + 10
1277 allocate(constants % nl(params % nsph*nngmax), &
1278 & stat=info)
1279 if (info .ne. 0) then
1280 call update_error(ddx_error, "`nl` " // &
1281 & "allocation failed")
1282 return
1283 end if
1284 constants % nl(1:params % nsph*(nngmax-10)) = &
1285 & tmp_nl(1:params % nsph*(nngmax-10))
1286 deallocate(tmp_nl, stat=info)
1287 if (info .ne. 0) then
1288 call update_error(ddx_error, "`tmp_nl` " // &
1289 & "deallocation failed")
1290 return
1291 end if
1292 end if
1293 end if
1294 end if
1295 end do
1296 end do
1297 constants % inl(params % nsph+1) = lnl+1
1298 constants % nngmax = nngmax
1299end subroutine neighbor_list_init
1300
1303subroutine neighbor_list_init_fmm(params, constants, ddx_error)
1304 type(ddx_params_type), intent(in) :: params
1305 type(ddx_constants_type), intent(inout) :: constants
1306 type(ddx_error_type), intent(inout) :: ddx_error
1307 real(dp) :: swthr, v(3), maxv, ssqv, vv, r
1308 integer :: nngmax, i, lnl, isph, jsph, inode, jnode, j, k, info
1309 integer, allocatable :: tmp_nl(:)
1310 ! Upper bound of switch region. Defines intersection criterion for spheres
1311 swthr = (params % se+one)*params % eta/two
1312 ! Build list of neighbours in CSR format
1313 nngmax = 1
1314 allocate(constants % inl(params % nsph+1), &
1315 & constants % nl(params % nsph*nngmax), stat=info)
1316 if (info .ne. 0) then
1317 call update_error(ddx_error, "`inl` and `nl` allocations failed")
1318 return
1319 end if
1320 i = 1
1321 lnl = 0
1322 ! loop over leaf clusters
1323 do isph = 1, params % nsph
1324 constants % inl(isph) = lnl + 1
1325 inode = constants % snode(isph)
1326 ! loop over near clusters
1327 do j = constants % snear(inode), constants % snear(inode+1) - 1
1328 jnode = constants % near(j)
1329 do k = constants % cluster(1, jnode), constants % cluster(2, jnode)
1330 jsph = constants % order(k)
1331 if (isph .ne. jsph) then
1332 v = params % csph(:, isph)-params % csph(:, jsph)
1333 maxv = max(abs(v(1)), abs(v(2)), abs(v(3)))
1334 if (maxv .gt. zero) then
1335 ssqv = (v(1)/maxv)**2 + (v(2)/maxv)**2 + (v(3)/maxv)**2
1336 else
1337 ssqv = zero
1338 end if
1339 vv = maxv * sqrt(ssqv)
1340 ! Take regularization parameter into account with respect to
1341 ! shift se. It is described properly by the upper bound of a
1342 ! switch region `swthr`. Note that we take the largest
1343 ! switching region as we need a symmetric neighbor list to
1344 ! use it also for the adjoint products.
1345 r = params % rsph(isph) + params % rsph(jsph) &
1346 & + swthr * max(params % rsph(isph), params % rsph(jsph))
1347 if (vv .le. r) then
1348 constants % nl(i) = jsph
1349 i = i + 1
1350 lnl = lnl + 1
1351 ! Extend ddx_data % nl if needed
1352 if (i .gt. params % nsph*nngmax) then
1353 allocate(tmp_nl(params % nsph*nngmax), stat=info)
1354 if (info .ne. 0) then
1355 call update_error(ddx_error, "`tmp_nl` " // &
1356 & "allocation failed")
1357 return
1358 end if
1359 tmp_nl(1:params % nsph*nngmax) = &
1360 & constants % nl(1:params % nsph*nngmax)
1361 deallocate(constants % nl, stat=info)
1362 if (info .ne. 0) then
1363 call update_error(ddx_error, "`nl` " // &
1364 & "deallocation failed")
1365 return
1366 end if
1367 nngmax = nngmax + 10
1368 allocate(constants % nl(params % nsph*nngmax), &
1369 & stat=info)
1370 if (info .ne. 0) then
1371 call update_error(ddx_error, "`nl` " // &
1372 & "allocation failed")
1373 return
1374 end if
1375 constants % nl(1:params % nsph*(nngmax-10)) = &
1376 & tmp_nl(1:params % nsph*(nngmax-10))
1377 deallocate(tmp_nl, stat=info)
1378 if (info .ne. 0) then
1379 call update_error(ddx_error, "`tmp_nl` " // &
1380 & "deallocation failed")
1381 return
1382 end if
1383 end if
1384 end if
1385 end if
1386 end do
1387 end do
1388 end do
1389 constants % inl(params % nsph+1) = lnl+1
1390 constants % nngmax = nngmax
1391end subroutine neighbor_list_init_fmm
1392
1420real(dp) function fsw(t, se, eta)
1421 ! Inputs
1422 real(dp), intent(in) :: t, se, eta
1423 ! Local variables
1424 real(dp) :: a, b, c, x
1425 ! Apply shift:
1426 ! se = 0 => t - eta/2 [ CENTERED ]
1427 ! se = 1 => t - eta [ EXTERIOR ]
1428 ! se = -1 => t [ INTERIOR ]
1429 x = t - (se*pt5 + pt5)*eta
1430 ! a must be in range (0,eta) for the switch region
1431 a = one - x
1432 ! Define switch function chi
1433 ! a <= 0 corresponds to the outside region
1434 if (a .le. zero) then
1435 fsw = zero
1436 ! a >= eta corresponds to the inside region
1437 else if (a .ge. eta) then
1438 fsw = one
1439 ! Switch region wher x is in range (1-eta,1)
1440 else
1441 ! Normalize a to region (0,1)
1442 a = a / eta
1443 b = 6d0*a - 15
1444 c = b*a + 10
1445 fsw = a * a * a * c
1446 end if
1447end function fsw
1448
1475real(dp) function dfsw(t, se, eta)
1476 ! Inputs
1477 real(dp), intent(in) :: t, se, eta
1478 ! Local variables
1479 real(dp) :: flow, x
1480 real(dp), parameter :: f30=30.0d0
1481 ! Apply shift:
1482 ! s = 0 => t - eta/2 [ CENTERED ]
1483 ! s = 1 => t - eta [ EXTERIOR ]
1484 ! s = -1 => t [ INTERIOR ]
1485 x = t - (se + 1.d0)*eta / 2.d0
1486 ! Lower bound of switch region
1487 flow = one - eta
1488 ! Define derivative of switch function chi
1489 if (x .ge. one) then
1490 dfsw = zero
1491 else if (x .le. flow) then
1492 dfsw = zero
1493 else
1494 dfsw = -f30 * (( (one-x)*(x-one+eta) )**2) / (eta**5)
1495 endif
1496end function dfsw
1497
1506real(dp) function ratio_d(d, chi)
1507 real(dp), intent(in) :: d, chi
1508 real(dp) :: den
1509 den = one - chi
1510 if (d .le. zero) then
1511 ratio_d = zero
1512 else if (den .le. d) then
1513 ratio_d = one
1514 else
1515 ratio_d = d/den
1516 end if
1517end function ratio_d
1518
1522subroutine mkprec(lmax, nbasis, nsph, ngrid, eps, ui, wgrid, vgrid, &
1523 & vgrid_nbasis, rx_prc, ddx_error)
1524 !! Inputs
1525 integer, intent(in) :: lmax, nbasis, nsph, ngrid, vgrid_nbasis
1526 real(dp), intent(in) :: eps, ui(ngrid, nsph), wgrid(ngrid), &
1527 & vgrid(vgrid_nbasis, ngrid)
1528 !! Output
1529 real(dp), intent(out) :: rx_prc(nbasis, nbasis, nsph)
1530 type(ddx_error_type), intent(inout) :: ddx_error
1531 !! Local variables
1532 integer :: info, isph, lm, l1, m1, ind1, igrid
1533 real(dp) :: f, f1
1534 integer, allocatable :: ipiv(:)
1535 real(dp), allocatable :: work(:)
1536 external :: dgetrf, dgetri
1537 ! Allocation of temporaries
1538 allocate(ipiv(nbasis), work(nbasis**2), stat=info)
1539 if (info .ne. 0) then
1540 call update_error(ddx_error, "mkprec: `ipiv` and `work` allocation failed")
1541 return
1542 endif
1543 ! Init
1544 rx_prc = zero
1545 ! Off-diagonal part
1546 do isph = 1, nsph
1547 do igrid = 1, ngrid
1548 f = twopi * ui(igrid, isph) * wgrid(igrid)
1549 do l1 = 0, lmax
1550 ind1 = l1*l1 + l1 + 1
1551 do m1 = -l1, l1
1552 f1 = f * vgrid(ind1+m1, igrid) / (two*dble(l1) + one)
1553 do lm = 1, nbasis
1554 rx_prc(lm, ind1+m1, isph) = &
1555 & rx_prc(lm, ind1+m1, isph) + f1*vgrid(lm, igrid)
1556 end do
1557 end do
1558 end do
1559 end do
1560 end do
1561 ! add diagonal
1562 f = twopi * (eps+one) / (eps-one)
1563 do isph = 1, nsph
1564 do lm = 1, nbasis
1565 rx_prc(lm, lm, isph) = rx_prc(lm, lm, isph) + f
1566 end do
1567 end do
1568 ! invert the blocks
1569 do isph = 1, nsph
1570 call dgetrf(nbasis, nbasis, rx_prc(:, :, isph), nbasis, ipiv, info)
1571 if (info .ne. 0) then
1572 call update_error(ddx_error, "mkprec: dgetrf failed")
1573 return
1574 end if
1575 call dgetri(nbasis, rx_prc(:, :, isph), nbasis, ipiv, work, &
1576 & nbasis**2, info)
1577 if (info .ne. 0) then
1578 call update_error(ddx_error, "mkprec: dgetri failed")
1579 return
1580 end if
1581 end do
1582 ! Cleanup temporaries
1583 deallocate(work, ipiv, stat=info)
1584 if (info .ne. 0) then
1585 call update_error(ddx_error, "mkprec: `ipiv` and `work` " // &
1586 & "deallocation failed")
1587 return
1588 end if
1589end subroutine mkprec
1590
1613subroutine tree_rib_build(nsph, csph, rsph, order, cluster, children, parent, &
1614 & cnode, rnode, snode, ddx_error)
1615 ! Inputs
1616 integer, intent(in) :: nsph
1617 real(dp), intent(in) :: csph(3, nsph), rsph(nsph)
1618 ! Outputs
1619 integer, intent(out) :: order(nsph), cluster(2, 2*nsph-1), &
1620 & children(2, 2*nsph-1), parent(2*nsph-1), snode(nsph)
1621 real(dp), intent(out) :: cnode(3, 2*nsph-1), rnode(2*nsph-1)
1622 type(ddx_error_type), intent(inout) :: ddx_error
1623 ! Local variables
1624 integer :: nclusters, i, j, n, s, e, div
1625 real(dp) :: r, r1, r2, c(3), c1(3), c2(3), d, maxc
1626 !! At first construct the tree
1627 nclusters = 2*nsph - 1
1628 ! Init the root node
1629 cluster(1, 1) = 1
1630 cluster(2, 1) = nsph
1631 parent(1) = 0
1632 ! Index of the first unassigned node
1633 j = 2
1634 ! Init spheres ordering
1635 do i = 1, nsph
1636 order(i) = i
1637 end do
1638 ! Divide nodes until a single spheres
1639 do i = 1, nclusters
1640 ! The first sphere in i-th node
1641 s = cluster(1, i)
1642 ! The last sphere in i-th node
1643 e = cluster(2, i)
1644 ! Number of spheres in i-th node
1645 n = e - s + 1
1646 ! Divide only if there are 2 or more spheres
1647 if (n .gt. 1) then
1648 ! Use inertial bisection to reorder spheres and cut into 2 halfs
1649 call tree_rib_node_bisect(nsph, csph, n, order(s:e), div, &
1650 & ddx_error)
1651 if (ddx_error % flag .ne. 0) then
1652 call update_error(ddx_error, "tree_rib_node_bisect returned " // &
1653 & "an ddx_error, exiting")
1654 return
1655 end if
1656 ! Assign the first half to the j-th node
1657 cluster(1, j) = s
1658 cluster(2, j) = s + div - 1
1659 ! Assign the second half to the (j+1)-th node
1660 cluster(1, j+1) = s + div
1661 cluster(2, j+1) = e
1662 ! Update list of children of i-th node
1663 children(1, i) = j
1664 children(2, i) = j + 1
1665 ! Set parents of new j-th and (j+1)-th node
1666 parent(j) = i
1667 parent(j+1) = i
1668 ! Shift index of the first unassigned node
1669 j = j + 2
1670 ! Set information for a leaf node
1671 else
1672 ! No children nodes
1673 children(:, i) = 0
1674 ! i-th node contains sphere ind(s)
1675 snode(order(s)) = i
1676 end if
1677 end do
1678 !! Compute bounding spheres of each node of the tree
1679 ! Bottom-to-top pass over all nodes of the tree
1680 do i = nclusters, 1, -1
1681 ! In case of a leaf node just use corresponding input sphere as a
1682 ! bounding sphere of the node
1683 if (children(1, i) .eq. 0) then
1684 ! Get correct index of the corresponding input sphere
1685 j = order(cluster(1, i))
1686 ! Get corresponding center and radius
1687 cnode(:, i) = csph(:, j)
1688 rnode(i) = rsph(j)
1689 ! In case of a non-leaf node get minimal sphere that contains bounding
1690 ! spheres of children nodes
1691 else
1692 ! The first child
1693 j = children(1, i)
1694 c1 = cnode(:, j)
1695 r1 = rnode(j)
1696 ! The second child
1697 j = children(2, i)
1698 c2 = cnode(:, j)
1699 r2 = rnode(j)
1700 ! Distance between centers of bounding spheres of children nodes
1701 c = c1 - c2
1702 ! Compute distance by scale and sum of scaled squares
1703 maxc = max(abs(c(1)), abs(c(2)), abs(c(3)))
1704 if (maxc .eq. zero) then
1705 d = zero
1706 else
1707 d = (c(1)/maxc)**2 + (c(2)/maxc)**2 + (c(3)/maxc)**2
1708 d = maxc * sqrt(d)
1709 end if
1710 ! If sphere #2 is completely inside sphere #1 use the first sphere
1711 if (r1 .ge. (r2+d)) then
1712 c = c1
1713 r = r1
1714 ! If sphere #1 is completely inside sphere #2 use the second sphere
1715 else if(r2 .ge. (r1+d)) then
1716 c = c2
1717 r = r2
1718 ! Otherwise use special formula to find a minimal sphere
1719 else
1720 r = (r1+r2+d) / 2
1721 c = c2 + c/d*(r-r2)
1722 end if
1723 ! Assign computed bounding sphere
1724 cnode(:, i) = c
1725 rnode(i) = r
1726 end if
1727 end do
1728end subroutine tree_rib_build
1729
1740subroutine tree_rib_node_bisect(nsph, csph, n, order, div, ddx_error)
1741 ! Inputs
1742 integer, intent(in) :: nsph, n
1743 real(dp), intent(in) :: csph(3, nsph)
1744 ! Outputs
1745 integer, intent(inout) :: order(n)
1746 integer, intent(out) :: div
1747 type(ddx_error_type), intent(inout) :: ddx_error
1748 ! Local variables
1749 real(dp) :: c(3), s(3)
1750 real(dp), allocatable :: tmp_csph(:, :), work(:)
1751 external :: dgesvd
1752 integer :: i, l, r, lwork, info, istat
1753 integer, allocatable :: tmp_order(:)
1754
1755 allocate(tmp_csph(3, n), tmp_order(n), stat=istat)
1756 if (istat.ne.0) then
1757 call update_error(ddx_error, 'Allocation failed in tree_node_bisect')
1758 return
1759 end if
1760
1761 ! Get average coordinate
1762 c = zero
1763 do i = 1, n
1764 c = c + csph(:, order(i))
1765 end do
1766 c = c / n
1767 ! Get coordinates minus average in a contiguous array
1768 do i = 1, n
1769 tmp_csph(:, i) = csph(:, order(i)) - c
1770 end do
1771 !! Find right singular vectors
1772 ! Get proper size of temporary workspace
1773 lwork = -1
1774 call dgesvd('N', 'O', 3, n, tmp_csph, 3, s, tmp_csph, 3, tmp_csph, 3, &
1775 & s, lwork, info)
1776 lwork = int(s(1))
1777 allocate(work(lwork), stat=istat)
1778 if (istat.ne.0) then
1779 call update_error(ddx_error, 'Allocation failed in tree_node_bisect')
1780 return
1781 end if
1782 ! Get right singular vectors
1783 call dgesvd('N', 'O', 3, n, tmp_csph, 3, s, tmp_csph, 3, tmp_csph, 3, &
1784 & work, lwork, info)
1785 if (info.ne.0) then
1786 call update_error(ddx_error, 'DGESVD failed in tree_node_bisect')
1787 return
1788 end if
1789 deallocate(work, stat=istat)
1790 if (istat.ne.0) then
1791 call update_error(ddx_error, 'Deallocation failed in tree_node_bisect')
1792 return
1793 end if
1794 !! Sort spheres by sign of the above scalar product, which is equal to
1795 !! the leading right singular vector scaled by the leading singular value.
1796 !! However, we only care about sign, so we take into account only the
1797 !! leading right singular vector.
1798 ! First empty index from the left
1799 l = 1
1800 ! First empty index from the right
1801 r = n
1802 ! Cycle over values of the singular vector
1803 do i = 1, n
1804 ! Positive scalar products are moved to the beginning of `order`
1805 if (tmp_csph(1, i) .ge. zero) then
1806 tmp_order(l) = order(i)
1807 l = l + 1
1808 ! Negative scalar products are moved to the end of `order`
1809 else
1810 tmp_order(r) = order(i)
1811 r = r - 1
1812 end if
1813 end do
1814 ! Set divider and update order
1815 div = r
1816 order = tmp_order
1817 deallocate(tmp_csph, tmp_order, stat=istat)
1818 if (istat.ne.0) then
1819 call update_error(ddx_error, 'Deallocation failed in tree_node_bisect')
1820 return
1821 end if
1822end subroutine tree_rib_node_bisect
1823
1848subroutine tree_get_farnear_work(n, children, cnode, rnode, lwork, iwork, &
1849 & jwork, work, nnfar, nfar, nnnear, nnear)
1850 ! Inputs
1851 integer, intent(in) :: n, children(2, n), lwork
1852 real(dp), intent(in) :: cnode(3, n), rnode(n)
1853 ! Outputs
1854 integer, intent(inout) :: iwork, jwork, work(3, lwork)
1855 integer, intent(out) :: nnfar, nfar(n), nnnear, nnear(n)
1856 ! Local variables
1857 integer :: j(2), npairs, k1, k2
1858 real(dp) :: c(3), r, d, dmax, dssq
1859 ! iwork is current temporary item in work array to process
1860 if (iwork .eq. 0) then
1861 work(1, 1) = 1
1862 work(2, 1) = 1
1863 iwork = 1
1864 jwork = 1
1865 end if
1866 ! jwork is total amount of temporary items in work array
1867 do while (iwork .le. jwork)
1868 j = work(1:2, iwork)
1869 c = cnode(:, j(1)) - cnode(:, j(2))
1870 r = rnode(j(1)) + rnode(j(2)) + max(rnode(j(1)), rnode(j(2)))
1871 !r = rnode(j(1)) + rnode(j(2))
1872 dmax = max(abs(c(1)), abs(c(2)), abs(c(3)))
1873 if (dmax .gt. zero) then
1874 dssq = (c(1)/dmax)**2 + (c(2)/dmax)**2 + (c(3)/dmax)**2
1875 else
1876 dssq = zero
1877 end if
1878 d = dmax * sqrt(dssq)
1879 !d = sqrt(c(1)**2 + c(2)**2 + c(3)**2)
1880 ! If node has no children then assume itself for purpose of finding
1881 ! far-field and near-filed interactions with children nodes of another
1882 ! node
1883 npairs = max(1, children(2, j(1))-children(1, j(1))+1) * &
1884 & max(1, children(2, j(2))-children(1, j(2))+1)
1885 if (d .ge. r) then
1886 ! Mark as far admissible pair
1887 work(3, iwork) = 1
1888 else if (npairs .eq. 1) then
1889 ! Mark as near admissible pair if both nodes are leaves
1890 work(3, iwork) = 2
1891 else if (jwork+npairs .gt. lwork) then
1892 ! Exit procedure, since work array was too small
1893 return
1894 else
1895 ! Mark as non-admissible pair and check all pairs of children nodes
1896 ! or pairs of one node (if it is a leaf node) with children of
1897 ! another node
1898 work(3, iwork) = 0
1899 if (children(1, j(1)) .eq. 0) then
1900 k1 = j(1)
1901 do k2 = children(1, j(2)), children(2, j(2))
1902 jwork = jwork + 1
1903 work(1, jwork) = k1
1904 work(2, jwork) = k2
1905 end do
1906 else if(children(1, j(2)) .eq. 0) then
1907 k2 = j(2)
1908 do k1 = children(1, j(1)), children(2, j(1))
1909 jwork = jwork + 1
1910 work(1, jwork) = k1
1911 work(2, jwork) = k2
1912 end do
1913 else
1914 do k1 = children(1, j(1)), children(2, j(1))
1915 do k2 = children(1, j(2)), children(2, j(2))
1916 jwork = jwork + 1
1917 work(1, jwork) = k1
1918 work(2, jwork) = k2
1919 end do
1920 end do
1921 end if
1922 end if
1923 iwork = iwork + 1
1924 end do
1925 nfar = 0
1926 nnear = 0
1927 do iwork = 1, jwork
1928 if (work(3, iwork) .eq. 1) then
1929 nfar(work(1, iwork)) = nfar(work(1, iwork)) + 1
1930 else if (work(3, iwork) .eq. 2) then
1931 nnear(work(1, iwork)) = nnear(work(1, iwork)) + 1
1932 end if
1933 end do
1934 iwork = jwork + 1
1935 nnfar = sum(nfar)
1936 nnnear = sum(nnear)
1937end subroutine tree_get_farnear_work
1938
1941subroutine tree_get_farnear(jwork, lwork, work, n, nnfar, nfar, sfar, far, &
1942 & nnnear, nnear, snear, near)
1943! Parameters:
1944! jwork: Total number of checked pairs in work array
1945! lwork: Total length of work array
1946! work: Work array itself
1947! n: Number of nodes
1948! nnfar: Total number of all far-field interactions
1949! nfar: Number of far-field interactions of each node
1950! sfar: Index in far array of first far-field node for each node
1951! far: Indexes of far-field nodes
1952! nnnear: Total number of all near-field interactions
1953! nnear: Number of near-field interactions of each node
1954! snear: Index in near array of first near-field node for each node
1955! near: Indexes of near-field nodes
1956 integer, intent(in) :: jwork, lwork, work(3, lwork), n, nnfar, nnnear
1957 integer, intent(in) :: nfar(n), nnear(n)
1958 integer, intent(out) :: sfar(n+1), far(nnfar), snear(n+1), near(nnnear)
1959 integer :: i, j
1960 integer, allocatable :: cfar(:), cnear(:)
1961
1962 allocate(cfar(n+1), cnear(n+1))
1963 sfar(1) = 1
1964 snear(1) = 1
1965 do i = 2, n+1
1966 sfar(i) = sfar(i-1) + nfar(i-1)
1967 snear(i) = snear(i-1) + nnear(i-1)
1968 end do
1969 cfar = sfar
1970 cnear = snear
1971 do i = 1, jwork
1972 if (work(3, i) .eq. 1) then
1973 ! Far
1974 j = work(1, i)
1975 if ((j .gt. n) .or. (j .le. 0)) then
1976 write(*,*) "ALARM", j
1977 end if
1978 far(cfar(j)) = work(2, i)
1979 cfar(j) = cfar(j) + 1
1980 else if (work(3, i) .eq. 2) then
1981 ! Near
1982 j = work(1, i)
1983 if ((j .gt. n) .or. (j .le. 0)) then
1984 write(*,*) "ALARM", j
1985 end if
1986 near(cnear(j)) = work(2, i)
1987 cnear(j) = cnear(j) + 1
1988 end if
1989 end do
1990
1991 deallocate(cfar, cnear)
1992
1993end subroutine tree_get_farnear
1994
2000subroutine constants_free(constants, ddx_error)
2001 implicit none
2002 type(ddx_constants_type), intent(out) :: constants
2003 type(ddx_error_type), intent(inout) :: ddx_error
2004 integer :: istat
2005
2006 istat = 0
2007
2008 if (allocated(constants % vscales)) then
2009 deallocate(constants % vscales, stat=istat)
2010 if (istat .ne. 0) then
2011 call update_error(ddx_error, "`vscales` deallocation failed!")
2012 end if
2013 end if
2014 if (allocated(constants % v4pi2lp1)) then
2015 deallocate(constants % v4pi2lp1, stat=istat)
2016 if (istat .ne. 0) then
2017 call update_error(ddx_error, "`v4pi2lp1` deallocation failed!")
2018 end if
2019 end if
2020 if (allocated(constants % vscales_rel)) then
2021 deallocate(constants % vscales_rel, stat=istat)
2022 if (istat .ne. 0) then
2023 call update_error(ddx_error, "`vscales_rel` deallocation failed!")
2024 end if
2025 end if
2026 if (allocated(constants % vfact)) then
2027 deallocate(constants % vfact, stat=istat)
2028 if (istat .ne. 0) then
2029 call update_error(ddx_error, "`vfact` deallocation failed!")
2030 end if
2031 end if
2032 if (allocated(constants % vcnk)) then
2033 deallocate(constants % vcnk, stat=istat)
2034 if (istat .ne. 0) then
2035 call update_error(ddx_error, "`vcnk` deallocation failed!")
2036 end if
2037 end if
2038 if (allocated(constants % m2l_ztranslate_coef)) then
2039 deallocate(constants % m2l_ztranslate_coef, stat=istat)
2040 if (istat .ne. 0) then
2041 call update_error(ddx_error, &
2042 & "`m2l_ztranslate_coef` deallocation failed!")
2043 end if
2044 end if
2045 if (allocated(constants % m2l_ztranslate_adj_coef)) then
2046 deallocate(constants % m2l_ztranslate_adj_coef, stat=istat)
2047 if (istat .ne. 0) then
2048 call update_error(ddx_error, &
2049 & "`m2l_ztranslate_adj_coef` deallocation failed!")
2050 end if
2051 end if
2052 if (allocated(constants % cgrid)) then
2053 deallocate(constants % cgrid, stat=istat)
2054 if (istat .ne. 0) then
2055 call update_error(ddx_error, "`cgrid` deallocation failed!")
2056 end if
2057 end if
2058 if (allocated(constants % wgrid)) then
2059 deallocate(constants % wgrid, stat=istat)
2060 if (istat .ne. 0) then
2061 call update_error(ddx_error, "`wgrid` deallocation failed!")
2062 end if
2063 end if
2064 if (allocated(constants % vgrid)) then
2065 deallocate(constants % vgrid, stat=istat)
2066 if (istat .ne. 0) then
2067 call update_error(ddx_error, "`vgrid` deallocation failed!")
2068 end if
2069 end if
2070 if (allocated(constants % vwgrid)) then
2071 deallocate(constants % vwgrid, stat=istat)
2072 if (istat .ne. 0) then
2073 call update_error(ddx_error, "`vwgrid` deallocation failed!")
2074 end if
2075 end if
2076 if (allocated(constants % vgrid2)) then
2077 deallocate(constants % vgrid2, stat=istat)
2078 if (istat .ne. 0) then
2079 call update_error(ddx_error, "`vgrid2` deallocation failed!")
2080 end if
2081 end if
2082 if (allocated(constants % pchi)) then
2083 deallocate(constants % pchi, stat=istat)
2084 if (istat .ne. 0) then
2085 call update_error(ddx_error, "`pchi` deallocation failed!")
2086 end if
2087 end if
2088 if (allocated(constants % c_ik)) then
2089 deallocate(constants % c_ik, stat=istat)
2090 if (istat .ne. 0) then
2091 call update_error(ddx_error, "`c_ik` deallocation failed!")
2092 end if
2093 end if
2094 if (allocated(constants % si_ri)) then
2095 deallocate(constants % si_ri, stat=istat)
2096 if (istat .ne. 0) then
2097 call update_error(ddx_error, "`si_ri` deallocation failed!")
2098 end if
2099 end if
2100 if (allocated(constants % di_ri)) then
2101 deallocate(constants % di_ri, stat=istat)
2102 if (istat .ne. 0) then
2103 call update_error(ddx_error, "`di_ri` deallocation failed!")
2104 end if
2105 end if
2106 if (allocated(constants % sk_ri)) then
2107 deallocate(constants % sk_ri, stat=istat)
2108 if (istat .ne. 0) then
2109 call update_error(ddx_error, "`sk_ri` deallocation failed!")
2110 end if
2111 end if
2112 if (allocated(constants % dk_ri)) then
2113 deallocate(constants % dk_ri, stat=istat)
2114 if (istat .ne. 0) then
2115 call update_error(ddx_error, "`dk_ri` deallocation failed!")
2116 end if
2117 end if
2118 if (allocated(constants % termimat)) then
2119 deallocate(constants % termimat, stat=istat)
2120 if (istat .ne. 0) then
2121 call update_error(ddx_error, "`termimat` deallocation failed!")
2122 end if
2123 end if
2124 if (allocated(constants % b)) then
2125 deallocate(constants % b, stat=istat)
2126 if (istat .ne. 0) then
2127 call update_error(ddx_error, "`b` deallocation failed!")
2128 end if
2129 end if
2130 if (allocated(constants % l)) then
2131 deallocate(constants % l, stat=istat)
2132 if (istat .ne. 0) then
2133 call update_error(ddx_error, "`l` deallocation failed!")
2134 end if
2135 end if
2136 if (allocated(constants % inl)) then
2137 deallocate(constants % inl, stat=istat)
2138 if (istat .ne. 0) then
2139 call update_error(ddx_error, "`inl` deallocation failed!")
2140 end if
2141 end if
2142 if (allocated(constants % nl)) then
2143 deallocate(constants % nl, stat=istat)
2144 if (istat .ne. 0) then
2145 call update_error(ddx_error, "`nl` deallocation failed!")
2146 end if
2147 end if
2148 if (allocated(constants % itrnl)) then
2149 deallocate(constants % itrnl, stat=istat)
2150 if (istat .ne. 0) then
2151 call update_error(ddx_error, "`itrnl` deallocation failed!")
2152 end if
2153 end if
2154 if (allocated(constants % fi)) then
2155 deallocate(constants % fi, stat=istat)
2156 if (istat .ne. 0) then
2157 call update_error(ddx_error, "`fi` deallocation failed!")
2158 end if
2159 end if
2160 if (allocated(constants % ui)) then
2161 deallocate(constants % ui, stat=istat)
2162 if (istat .ne. 0) then
2163 call update_error(ddx_error, "`ui` deallocation failed!")
2164 end if
2165 end if
2166 if (allocated(constants % ui_cav)) then
2167 deallocate(constants % ui_cav, stat=istat)
2168 if (istat .ne. 0) then
2169 call update_error(ddx_error, "`ui_cav` deallocation failed!")
2170 end if
2171 end if
2172 if (allocated(constants % zi)) then
2173 deallocate(constants % zi, stat=istat)
2174 if (istat .ne. 0) then
2175 call update_error(ddx_error, "`zi` deallocation failed!")
2176 end if
2177 end if
2178 if (allocated(constants % zi_dr)) then
2179 deallocate(constants % zi_dr, stat=istat)
2180 if (istat .ne. 0) then
2181 call update_error(ddx_error, "`zi_dr` deallocation failed!")
2182 end if
2183 end if
2184 if (allocated(constants % ncav_sph)) then
2185 deallocate(constants % ncav_sph, stat=istat)
2186 if (istat .ne. 0) then
2187 call update_error(ddx_error, "`ncav_sph` deallocation failed!")
2188 end if
2189 end if
2190 if (allocated(constants % ccav)) then
2191 deallocate(constants % ccav, stat=istat)
2192 if (istat .ne. 0) then
2193 call update_error(ddx_error, "`ccav` deallocation failed!")
2194 end if
2195 end if
2196 if (allocated(constants % icav_ia)) then
2197 deallocate(constants % icav_ia, stat=istat)
2198 if (istat .ne. 0) then
2199 call update_error(ddx_error, "`icav_ia` deallocation failed!")
2200 end if
2201 end if
2202 if (allocated(constants % icav_ja)) then
2203 deallocate(constants % icav_ja, stat=istat)
2204 if (istat .ne. 0) then
2205 call update_error(ddx_error, "`icav_ja` deallocation failed!")
2206 end if
2207 end if
2208 if (allocated(constants % rx_prc)) then
2209 deallocate(constants % rx_prc, stat=istat)
2210 if (istat .ne. 0) then
2211 call update_error(ddx_error, "`rx_prc` deallocation failed!")
2212 end if
2213 end if
2214 if (allocated(constants % order)) then
2215 deallocate(constants % order, stat=istat)
2216 if (istat .ne. 0) then
2217 call update_error(ddx_error, "`order` deallocation failed!")
2218 end if
2219 end if
2220 if (allocated(constants % cluster)) then
2221 deallocate(constants % cluster, stat=istat)
2222 if (istat .ne. 0) then
2223 call update_error(ddx_error, "`cluster` deallocation failed!")
2224 end if
2225 end if
2226 if (allocated(constants % children)) then
2227 deallocate(constants % children, stat=istat)
2228 if (istat .ne. 0) then
2229 call update_error(ddx_error, "`children` deallocation failed!")
2230 end if
2231 end if
2232 if (allocated(constants % parent)) then
2233 deallocate(constants % parent, stat=istat)
2234 if (istat .ne. 0) then
2235 call update_error(ddx_error, "`parent` deallocation failed!")
2236 end if
2237 end if
2238 if (allocated(constants % cnode)) then
2239 deallocate(constants % cnode, stat=istat)
2240 if (istat .ne. 0) then
2241 call update_error(ddx_error, "`cnode` deallocation failed!")
2242 end if
2243 end if
2244 if (allocated(constants % rnode)) then
2245 deallocate(constants % rnode, stat=istat)
2246 if (istat .ne. 0) then
2247 call update_error(ddx_error, "`rnode` deallocation failed!")
2248 end if
2249 end if
2250 if (allocated(constants % snode)) then
2251 deallocate(constants % snode, stat=istat)
2252 if (istat .ne. 0) then
2253 call update_error(ddx_error, "`snode` deallocation failed!")
2254 end if
2255 end if
2256 if (allocated(constants % sk_rnode)) then
2257 deallocate(constants % sk_rnode, stat=istat)
2258 if (istat .ne. 0) then
2259 call update_error(ddx_error, "`sk_rnode` deallocation failed!")
2260 end if
2261 end if
2262 if (allocated(constants % si_rnode)) then
2263 deallocate(constants % si_rnode, stat=istat)
2264 if (istat .ne. 0) then
2265 call update_error(ddx_error, "`si_rnode` deallocation failed!")
2266 end if
2267 end if
2268 if (allocated(constants % nfar)) then
2269 deallocate(constants % nfar, stat=istat)
2270 if (istat .ne. 0) then
2271 call update_error(ddx_error, "`nfar` deallocation failed!")
2272 end if
2273 end if
2274 if (allocated(constants % nnear)) then
2275 deallocate(constants % nnear, stat=istat)
2276 if (istat .ne. 0) then
2277 call update_error(ddx_error, "`nnear` deallocation failed!")
2278 end if
2279 end if
2280 if (allocated(constants % far)) then
2281 deallocate(constants % far, stat=istat)
2282 if (istat .ne. 0) then
2283 call update_error(ddx_error, "`far` deallocation failed!")
2284 end if
2285 end if
2286 if (allocated(constants % near)) then
2287 deallocate(constants % near, stat=istat)
2288 if (istat .ne. 0) then
2289 call update_error(ddx_error, "`near` deallocation failed!")
2290 end if
2291 end if
2292 if (allocated(constants % sfar)) then
2293 deallocate(constants % sfar, stat=istat)
2294 if (istat .ne. 0) then
2295 call update_error(ddx_error, "`sfar` deallocation failed!")
2296 end if
2297 end if
2298 if (allocated(constants % snear)) then
2299 deallocate(constants % snear, stat=istat)
2300 if (istat .ne. 0) then
2301 call update_error(ddx_error, "`snear` deallocation failed!")
2302 end if
2303 end if
2304 if (allocated(constants % m2l_ztranslate_coef)) then
2305 deallocate(constants % m2l_ztranslate_coef, stat=istat)
2306 if (istat .ne. 0) then
2307 call update_error(ddx_error, "`m2l_ztranslate_coef` " // &
2308 & "deallocation failed!")
2309 end if
2310 end if
2311 if (allocated(constants % m2l_ztranslate_adj_coef)) then
2312 deallocate(constants % m2l_ztranslate_adj_coef, stat=istat)
2313 if (istat .ne. 0) then
2314 call update_error(ddx_error, "`m2l_ztranslate_adj_coef` " // &
2315 & "deallocation failed!")
2316 end if
2317 end if
2318
2319 call switching_free(constants%switching, ddx_error)
2320
2321end subroutine constants_free
2322
2323subroutine switching_init(params, constants, switching, ddx_error)
2324 type(ddx_params_type), intent(in) :: params
2325 type(ddx_constants_type), intent(in) :: constants
2326 type(ddx_switching_type), intent(inout) :: switching
2327 type(ddx_error_type), intent(inout) :: ddx_error
2328 integer :: info, isph, n, icav, ncav
2329 real(dp) :: d_ni, f_ni, u_ni
2330
2331 allocate(switching%u_ni(params%ngrid, params%nsph), &
2332 & switching%d_ni(params%ngrid, params%nsph), &
2333 & switching%f_ni(params%ngrid, params%nsph), &
2334 & stat=info)
2335 if (info.ne.0) then
2336 call update_error(ddx_error, &
2337 & "Allocation failed in switching_init")
2338 return
2339 end if
2340
2341 do isph = 1, params%nsph
2342 do n = 1, params%ngrid
2343 call compute_d_and_f(params, constants, isph, n, d_ni, f_ni)
2344 switching%d_ni(n, isph) = d_ni
2345 switching%f_ni(n, isph) = f_ni
2346 end do
2347 end do
2348
2349 ncav = 0
2350 do isph = 1, params%nsph
2351 do n = 1, params%ngrid
2352 u_ni = compute_u(params, constants, isph, n)
2353 if (u_ni.ne.zero) ncav = ncav + 1
2354 end do
2355 end do
2356 switching%ncav = ncav
2357
2358 allocate(switching%u_i_cav(switching%ncav), stat=info)
2359 if (info.ne.0) then
2360 call update_error(ddx_error, &
2361 & "Allocation failed in switching_init")
2362 return
2363 end if
2364
2365 icav = 0
2366 do isph = 1, params%nsph
2367 do n = 1, params%ngrid
2368 u_ni = compute_u(params, constants, isph, n)
2369 switching%u_ni(n, isph) = u_ni
2370 if (u_ni.ne.zero) then
2371 icav = icav + 1
2372 switching%u_i_cav(icav) = u_ni
2373 end if
2374 end do
2375 end do
2376
2377 if (params%force.eq.1) then
2378 allocate(switching%grad_u_ni(3, params%ngrid, params%nsph), &
2379 & switching%dr_u_ni(params%ngrid, params%nsph), &
2380 & stat=info)
2381 if (info.ne.0) then
2382 call update_error(ddx_error, &
2383 & "Allocation failed in switching_init")
2384 return
2385 end if
2386
2387 do isph = 1, params%nsph
2388 do n = 1, params%ngrid
2389 call compute_grad_i_u_i(params, constants, isph, n, &
2390 & switching%grad_u_ni(:, n, isph), &
2391 & switching%dr_u_ni(n, isph))
2392 end do
2393 end do
2394 end if
2395
2396end subroutine switching_init
2397
2398subroutine switching_free(switching, ddx_error)
2399 type(ddx_switching_type), intent(inout) :: switching
2400 type(ddx_error_type), intent(inout) :: ddx_error
2401 integer :: info
2402 if (allocated(switching%u_ni)) then
2403 deallocate(switching%u_ni, stat=info)
2404 if (info.ne.0) then
2405 call update_error(ddx_error, &
2406 & "Deallocation failed in switching_free")
2407 return
2408 end if
2409 end if
2410 if (allocated(switching%u_i_cav)) then
2411 deallocate(switching%u_i_cav, stat=info)
2412 if (info.ne.0) then
2413 call update_error(ddx_error, &
2414 & "Deallocation failed in switching_free")
2415 return
2416 end if
2417 end if
2418 if (allocated(switching%d_ni)) then
2419 deallocate(switching%d_ni, stat=info)
2420 if (info.ne.0) then
2421 call update_error(ddx_error, &
2422 & "Deallocation failed in switching_free")
2423 return
2424 end if
2425 end if
2426 if (allocated(switching%f_ni)) then
2427 deallocate(switching%f_ni, stat=info)
2428 if (info.ne.0) then
2429 call update_error(ddx_error, &
2430 & "Deallocation failed in switching_free")
2431 return
2432 end if
2433 end if
2434 if (allocated(switching%grad_u_ni)) then
2435 deallocate(switching%grad_u_ni, stat=info)
2436 if (info.ne.0) then
2437 call update_error(ddx_error, &
2438 & "Deallocation failed in switching_free")
2439 return
2440 end if
2441 end if
2442 if (allocated(switching%dr_u_ni)) then
2443 deallocate(switching%dr_u_ni, stat=info)
2444 if (info.ne.0) then
2445 call update_error(ddx_error, &
2446 & "Deallocation failed in switching_free")
2447 return
2448 end if
2449 end if
2450end subroutine switching_free
2451
2452real(dp) function compute_t(params, constants, isph, jsph, n)
2453 type(ddx_params_type), intent(in) :: params
2454 type(ddx_constants_type), intent(in) :: constants
2455 integer, intent(in) :: isph, jsph, n
2456 associate(cgrid => constants%cgrid, rsph => params%rsph, &
2457 & csph => params%csph)
2458 compute_t = sqrt((cgrid(1,n)*rsph(isph) + csph(1,isph) &
2459 & - csph(1,jsph))**2 + (cgrid(2,n)*rsph(isph) &
2460 & + csph(2,isph) - csph(2,jsph))**2 &
2461 & + (cgrid(3,n)*rsph(isph) + csph(3,isph) &
2462 & - csph(3,jsph))**2)/rsph(jsph)
2463 end associate
2464end function compute_t
2465
2466real(dp) function compute_chi(params, constants, isph, jsph, n)
2467 type(ddx_params_type), intent(in) :: params
2468 type(ddx_constants_type), intent(in) :: constants
2469 integer, intent(in) :: isph, jsph, n
2470 real(dp) :: t_ijn
2471 t_ijn = compute_t(params, constants, isph, jsph, n)
2472 compute_chi = fsw(t_ijn, params%se, params%eta)
2473end function compute_chi
2474
2475subroutine compute_d_and_f(params, constants, isph, n, d, f)
2476 type(ddx_params_type), intent(in) :: params
2477 type(ddx_constants_type), intent(in) :: constants
2478 integer, intent(in) :: isph, n
2479 real(dp), intent(out) :: d, f
2480 integer :: ij, jsph
2481 real(dp) :: chi_ijn
2482 d = one
2483 f = zero
2484 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2485 jsph = constants%nl(ij)
2486 chi_ijn = compute_chi(params, constants, isph, jsph, n)
2487 d = d*(one - chi_ijn)
2488 f = f + chi_ijn
2489 end do
2490end subroutine compute_d_and_f
2491
2492real(dp) function compute_d(params, constants, isph, n)
2493 type(ddx_params_type), intent(in) :: params
2494 type(ddx_constants_type), intent(in) :: constants
2495 integer, intent(in) :: isph, n
2496 integer :: ij, jsph
2497 real(dp) :: chi_ijn
2498 compute_d = one
2499 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2500 jsph = constants%nl(ij)
2501 chi_ijn = compute_chi(params, constants, isph, jsph, n)
2502 compute_d = compute_d*(one - chi_ijn)
2503 end do
2504end function
2505
2506real(dp) function compute_f(params, constants, isph, n)
2507 type(ddx_params_type), intent(in) :: params
2508 type(ddx_constants_type), intent(in) :: constants
2509 integer, intent(in) :: isph, n
2510 integer :: ij, jsph
2511
2512 compute_f = zero
2513 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2514 jsph = constants%nl(ij)
2515 compute_f = compute_f + compute_chi(params, constants, isph, jsph, n)
2516 end do
2517end function compute_f
2518
2519real(dp) function compute_omega(params, constants, isph, jsph, n)
2520 type(ddx_params_type), intent(in) :: params
2521 type(ddx_constants_type), intent(in) :: constants
2522 integer, intent(in) :: isph, jsph, n
2523 real(dp) :: chi_nij, d_ni, f_ni
2524
2525 chi_nij = compute_chi(params, constants, isph, jsph, n)
2526 d_ni = constants%switching%d_ni(n, isph)
2527 f_ni = constants%switching%f_ni(n, isph)
2528 compute_omega = chi_nij/(d_ni + f_ni)
2529end function compute_omega
2530
2531real(dp) function compute_u(params, constants, isph, n)
2532 type(ddx_params_type), intent(in) :: params
2533 type(ddx_constants_type), intent(in) :: constants
2534 integer, intent(in) :: isph, n
2535 integer :: ij, jsph
2536 compute_u = one
2537 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2538 jsph = constants%nl(ij)
2539 compute_u = compute_u - compute_omega(params, constants, isph, jsph, n)
2540 end do
2541 ! Remove roundoff so that the `u.eq.zero` checks are reliable. The exact u
2542 ! is in [0, 1], so a negative u is removed too: the cavity points are
2543 ! counted with `u.gt.zero` but stored for `u.ne.zero`.
2544 if (compute_u .le. 1d-15) compute_u = zero
2545end function compute_u
2546
2547subroutine compute_grad_i_u_i(params, constants, isph, n, grad, dr)
2548 type(ddx_params_type), intent(in) :: params
2549 type(ddx_constants_type), intent(in) :: constants
2550 real(dp), intent(out) :: grad(3), dr
2551 integer, intent(in) :: isph, n
2552 integer :: ij, jsph
2553 real(dp) :: vij(3), vvij, tij, sij(3), grad_tij(3), dr_tij, &
2554 & d_i, f_i, chi_ij, a, c, grad_p
2555
2556 grad(:) = zero
2557 dr = zero
2558 d_i = constants%switching%d_ni(n, isph)
2559 f_i = constants%switching%f_ni(n, isph)
2560 a = one/(d_i + f_i)**2
2561 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2562 jsph = constants%nl(ij)
2563 vij = params % csph(:,isph) + &
2564 & params % rsph(isph)*constants % cgrid(:,n) - &
2565 & params % csph(:,jsph)
2566 vvij = sqrt(vij(1)*vij(1) + vij(2)*vij(2) + vij(3)*vij(3))
2567 tij = vvij/params % rsph(jsph)
2568 sij = vij/vvij
2569
2570 grad_tij(:) = sij(:)/params%rsph(jsph)
2571 dr_tij = dot_product(vij, constants%cgrid(:,n)) &
2572 & /(vvij*params%rsph(jsph))
2573
2574 grad_p = dfsw(tij, params % se, params % eta)
2575 chi_ij = fsw(tij, params % se, params % eta)
2576
2577 c = a*(d_i + f_i*ratio_d(d_i, chi_ij))*grad_p
2578
2579 grad(:) = grad(:) + c*grad_tij(:)
2580 dr = dr + c*dr_tij
2581 end do
2582
2583end subroutine compute_grad_i_u_i
2584
2585end module ddx_constants
2586
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:580
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:722
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:1850
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:1219
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:1524
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:605
ddx_constants::ratio_d
real(dp) function ratio_d(d, chi)
Bounded evaluation of d/(1 - chi_j) with .
Definition: ddx_constants.f90:1507
ddx_constants::dfsw
real(dp) function dfsw(t, se, eta)
Derivative of a switching function.
Definition: ddx_constants.f90:1476
ddx_constants::constants_init
subroutine constants_init(params, constants, ddx_error)
Compute all necessary constants.
Definition: ddx_constants.f90:260
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:1304
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:886
ddx_constants::constants_free
subroutine constants_free(constants, ddx_error)
Deallocate the constants.
Definition: ddx_constants.f90:2001
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:1615
ddx_constants::fsw
real(dp) function fsw(t, se, eta)
Switching function.
Definition: ddx_constants.f90:1421
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:1741
ddx_constants::mkpmat
subroutine mkpmat(params, constants, isph, pmat)
Computation of P_chi.
Definition: ddx_constants.f90:850
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:1943
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:45
ddx_constants::ddx_switching_type
Definition: ddx_constants.f90:26
ddx_parameters::ddx_params_type
Type to check and store user input parameters.
Definition: ddx_parameters.f90:25