ddx 1.0.0
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
1501subroutine mkprec(lmax, nbasis, nsph, ngrid, eps, ui, wgrid, vgrid, &
1502 & vgrid_nbasis, rx_prc, ddx_error)
1503 !! Inputs
1504 integer, intent(in) :: lmax, nbasis, nsph, ngrid, vgrid_nbasis
1505 real(dp), intent(in) :: eps, ui(ngrid, nsph), wgrid(ngrid), &
1506 & vgrid(vgrid_nbasis, ngrid)
1507 !! Output
1508 real(dp), intent(out) :: rx_prc(nbasis, nbasis, nsph)
1509 type(ddx_error_type), intent(inout) :: ddx_error
1510 !! Local variables
1511 integer :: info, isph, lm, l1, m1, ind1, igrid
1512 real(dp) :: f, f1
1513 integer, allocatable :: ipiv(:)
1514 real(dp), allocatable :: work(:)
1515 external :: dgetrf, dgetri
1516 ! Allocation of temporaries
1517 allocate(ipiv(nbasis), work(nbasis**2), stat=info)
1518 if (info .ne. 0) then
1519 call update_error(ddx_error, "mkprec: `ipiv` and `work` allocation failed")
1520 return
1521 endif
1522 ! Init
1523 rx_prc = zero
1524 ! Off-diagonal part
1525 do isph = 1, nsph
1526 do igrid = 1, ngrid
1527 f = twopi * ui(igrid, isph) * wgrid(igrid)
1528 do l1 = 0, lmax
1529 ind1 = l1*l1 + l1 + 1
1530 do m1 = -l1, l1
1531 f1 = f * vgrid(ind1+m1, igrid) / (two*dble(l1) + one)
1532 do lm = 1, nbasis
1533 rx_prc(lm, ind1+m1, isph) = &
1534 & rx_prc(lm, ind1+m1, isph) + f1*vgrid(lm, igrid)
1535 end do
1536 end do
1537 end do
1538 end do
1539 end do
1540 ! add diagonal
1541 f = twopi * (eps+one) / (eps-one)
1542 do isph = 1, nsph
1543 do lm = 1, nbasis
1544 rx_prc(lm, lm, isph) = rx_prc(lm, lm, isph) + f
1545 end do
1546 end do
1547 ! invert the blocks
1548 do isph = 1, nsph
1549 call dgetrf(nbasis, nbasis, rx_prc(:, :, isph), nbasis, ipiv, info)
1550 if (info .ne. 0) then
1551 call update_error(ddx_error, "mkprec: dgetrf failed")
1552 return
1553 end if
1554 call dgetri(nbasis, rx_prc(:, :, isph), nbasis, ipiv, work, &
1555 & nbasis**2, info)
1556 if (info .ne. 0) then
1557 call update_error(ddx_error, "mkprec: dgetri failed")
1558 return
1559 end if
1560 end do
1561 ! Cleanup temporaries
1562 deallocate(work, ipiv, stat=info)
1563 if (info .ne. 0) then
1564 call update_error(ddx_error, "mkprec: `ipiv` and `work` " // &
1565 & "deallocation failed")
1566 return
1567 end if
1568end subroutine mkprec
1569
1592subroutine tree_rib_build(nsph, csph, rsph, order, cluster, children, parent, &
1593 & cnode, rnode, snode, ddx_error)
1594 ! Inputs
1595 integer, intent(in) :: nsph
1596 real(dp), intent(in) :: csph(3, nsph), rsph(nsph)
1597 ! Outputs
1598 integer, intent(out) :: order(nsph), cluster(2, 2*nsph-1), &
1599 & children(2, 2*nsph-1), parent(2*nsph-1), snode(nsph)
1600 real(dp), intent(out) :: cnode(3, 2*nsph-1), rnode(2*nsph-1)
1601 type(ddx_error_type), intent(inout) :: ddx_error
1602 ! Local variables
1603 integer :: nclusters, i, j, n, s, e, div
1604 real(dp) :: r, r1, r2, c(3), c1(3), c2(3), d, maxc
1605 !! At first construct the tree
1606 nclusters = 2*nsph - 1
1607 ! Init the root node
1608 cluster(1, 1) = 1
1609 cluster(2, 1) = nsph
1610 parent(1) = 0
1611 ! Index of the first unassigned node
1612 j = 2
1613 ! Init spheres ordering
1614 do i = 1, nsph
1615 order(i) = i
1616 end do
1617 ! Divide nodes until a single spheres
1618 do i = 1, nclusters
1619 ! The first sphere in i-th node
1620 s = cluster(1, i)
1621 ! The last sphere in i-th node
1622 e = cluster(2, i)
1623 ! Number of spheres in i-th node
1624 n = e - s + 1
1625 ! Divide only if there are 2 or more spheres
1626 if (n .gt. 1) then
1627 ! Use inertial bisection to reorder spheres and cut into 2 halfs
1628 call tree_rib_node_bisect(nsph, csph, n, order(s:e), div, &
1629 & ddx_error)
1630 if (ddx_error % flag .ne. 0) then
1631 call update_error(ddx_error, "tree_rib_node_bisect returned " // &
1632 & "an ddx_error, exiting")
1633 return
1634 end if
1635 ! Assign the first half to the j-th node
1636 cluster(1, j) = s
1637 cluster(2, j) = s + div - 1
1638 ! Assign the second half to the (j+1)-th node
1639 cluster(1, j+1) = s + div
1640 cluster(2, j+1) = e
1641 ! Update list of children of i-th node
1642 children(1, i) = j
1643 children(2, i) = j + 1
1644 ! Set parents of new j-th and (j+1)-th node
1645 parent(j) = i
1646 parent(j+1) = i
1647 ! Shift index of the first unassigned node
1648 j = j + 2
1649 ! Set information for a leaf node
1650 else
1651 ! No children nodes
1652 children(:, i) = 0
1653 ! i-th node contains sphere ind(s)
1654 snode(order(s)) = i
1655 end if
1656 end do
1657 !! Compute bounding spheres of each node of the tree
1658 ! Bottom-to-top pass over all nodes of the tree
1659 do i = nclusters, 1, -1
1660 ! In case of a leaf node just use corresponding input sphere as a
1661 ! bounding sphere of the node
1662 if (children(1, i) .eq. 0) then
1663 ! Get correct index of the corresponding input sphere
1664 j = order(cluster(1, i))
1665 ! Get corresponding center and radius
1666 cnode(:, i) = csph(:, j)
1667 rnode(i) = rsph(j)
1668 ! In case of a non-leaf node get minimal sphere that contains bounding
1669 ! spheres of children nodes
1670 else
1671 ! The first child
1672 j = children(1, i)
1673 c1 = cnode(:, j)
1674 r1 = rnode(j)
1675 ! The second child
1676 j = children(2, i)
1677 c2 = cnode(:, j)
1678 r2 = rnode(j)
1679 ! Distance between centers of bounding spheres of children nodes
1680 c = c1 - c2
1681 ! Compute distance by scale and sum of scaled squares
1682 maxc = max(abs(c(1)), abs(c(2)), abs(c(3)))
1683 if (maxc .eq. zero) then
1684 d = zero
1685 else
1686 d = (c(1)/maxc)**2 + (c(2)/maxc)**2 + (c(3)/maxc)**2
1687 d = maxc * sqrt(d)
1688 end if
1689 ! If sphere #2 is completely inside sphere #1 use the first sphere
1690 if (r1 .ge. (r2+d)) then
1691 c = c1
1692 r = r1
1693 ! If sphere #1 is completely inside sphere #2 use the second sphere
1694 else if(r2 .ge. (r1+d)) then
1695 c = c2
1696 r = r2
1697 ! Otherwise use special formula to find a minimal sphere
1698 else
1699 r = (r1+r2+d) / 2
1700 c = c2 + c/d*(r-r2)
1701 end if
1702 ! Assign computed bounding sphere
1703 cnode(:, i) = c
1704 rnode(i) = r
1705 end if
1706 end do
1707end subroutine tree_rib_build
1708
1719subroutine tree_rib_node_bisect(nsph, csph, n, order, div, ddx_error)
1720 ! Inputs
1721 integer, intent(in) :: nsph, n
1722 real(dp), intent(in) :: csph(3, nsph)
1723 ! Outputs
1724 integer, intent(inout) :: order(n)
1725 integer, intent(out) :: div
1726 type(ddx_error_type), intent(inout) :: ddx_error
1727 ! Local variables
1728 real(dp) :: c(3), s(3)
1729 real(dp), allocatable :: tmp_csph(:, :), work(:)
1730 external :: dgesvd
1731 integer :: i, l, r, lwork, info, istat
1732 integer, allocatable :: tmp_order(:)
1733
1734 allocate(tmp_csph(3, n), tmp_order(n), stat=istat)
1735 if (istat.ne.0) then
1736 call update_error(ddx_error, 'Allocation failed in tree_node_bisect')
1737 return
1738 end if
1739
1740 ! Get average coordinate
1741 c = zero
1742 do i = 1, n
1743 c = c + csph(:, order(i))
1744 end do
1745 c = c / n
1746 ! Get coordinates minus average in a contiguous array
1747 do i = 1, n
1748 tmp_csph(:, i) = csph(:, order(i)) - c
1749 end do
1750 !! Find right singular vectors
1751 ! Get proper size of temporary workspace
1752 lwork = -1
1753 call dgesvd('N', 'O', 3, n, tmp_csph, 3, s, tmp_csph, 3, tmp_csph, 3, &
1754 & s, lwork, info)
1755 lwork = int(s(1))
1756 allocate(work(lwork), stat=istat)
1757 if (istat.ne.0) then
1758 call update_error(ddx_error, 'Allocation failed in tree_node_bisect')
1759 return
1760 end if
1761 ! Get right singular vectors
1762 call dgesvd('N', 'O', 3, n, tmp_csph, 3, s, tmp_csph, 3, tmp_csph, 3, &
1763 & work, lwork, info)
1764 if (info.ne.0) then
1765 call update_error(ddx_error, 'DGESVD failed in tree_node_bisect')
1766 return
1767 end if
1768 deallocate(work, stat=istat)
1769 if (istat.ne.0) then
1770 call update_error(ddx_error, 'Deallocation failed in tree_node_bisect')
1771 return
1772 end if
1773 !! Sort spheres by sign of the above scalar product, which is equal to
1774 !! the leading right singular vector scaled by the leading singular value.
1775 !! However, we only care about sign, so we take into account only the
1776 !! leading right singular vector.
1777 ! First empty index from the left
1778 l = 1
1779 ! First empty index from the right
1780 r = n
1781 ! Cycle over values of the singular vector
1782 do i = 1, n
1783 ! Positive scalar products are moved to the beginning of `order`
1784 if (tmp_csph(1, i) .ge. zero) then
1785 tmp_order(l) = order(i)
1786 l = l + 1
1787 ! Negative scalar products are moved to the end of `order`
1788 else
1789 tmp_order(r) = order(i)
1790 r = r - 1
1791 end if
1792 end do
1793 ! Set divider and update order
1794 div = r
1795 order = tmp_order
1796 deallocate(tmp_csph, tmp_order, stat=istat)
1797 if (istat.ne.0) then
1798 call update_error(ddx_error, 'Deallocation failed in tree_node_bisect')
1799 return
1800 end if
1801end subroutine tree_rib_node_bisect
1802
1827subroutine tree_get_farnear_work(n, children, cnode, rnode, lwork, iwork, &
1828 & jwork, work, nnfar, nfar, nnnear, nnear)
1829 ! Inputs
1830 integer, intent(in) :: n, children(2, n), lwork
1831 real(dp), intent(in) :: cnode(3, n), rnode(n)
1832 ! Outputs
1833 integer, intent(inout) :: iwork, jwork, work(3, lwork)
1834 integer, intent(out) :: nnfar, nfar(n), nnnear, nnear(n)
1835 ! Local variables
1836 integer :: j(2), npairs, k1, k2
1837 real(dp) :: c(3), r, d, dmax, dssq
1838 ! iwork is current temporary item in work array to process
1839 if (iwork .eq. 0) then
1840 work(1, 1) = 1
1841 work(2, 1) = 1
1842 iwork = 1
1843 jwork = 1
1844 end if
1845 ! jwork is total amount of temporary items in work array
1846 do while (iwork .le. jwork)
1847 j = work(1:2, iwork)
1848 c = cnode(:, j(1)) - cnode(:, j(2))
1849 r = rnode(j(1)) + rnode(j(2)) + max(rnode(j(1)), rnode(j(2)))
1850 !r = rnode(j(1)) + rnode(j(2))
1851 dmax = max(abs(c(1)), abs(c(2)), abs(c(3)))
1852 if (dmax .gt. zero) then
1853 dssq = (c(1)/dmax)**2 + (c(2)/dmax)**2 + (c(3)/dmax)**2
1854 else
1855 dssq = zero
1856 end if
1857 d = dmax * sqrt(dssq)
1858 !d = sqrt(c(1)**2 + c(2)**2 + c(3)**2)
1859 ! If node has no children then assume itself for purpose of finding
1860 ! far-field and near-filed interactions with children nodes of another
1861 ! node
1862 npairs = max(1, children(2, j(1))-children(1, j(1))+1) * &
1863 & max(1, children(2, j(2))-children(1, j(2))+1)
1864 if (d .ge. r) then
1865 ! Mark as far admissible pair
1866 work(3, iwork) = 1
1867 else if (npairs .eq. 1) then
1868 ! Mark as near admissible pair if both nodes are leaves
1869 work(3, iwork) = 2
1870 else if (jwork+npairs .gt. lwork) then
1871 ! Exit procedure, since work array was too small
1872 return
1873 else
1874 ! Mark as non-admissible pair and check all pairs of children nodes
1875 ! or pairs of one node (if it is a leaf node) with children of
1876 ! another node
1877 work(3, iwork) = 0
1878 if (children(1, j(1)) .eq. 0) then
1879 k1 = j(1)
1880 do k2 = children(1, j(2)), children(2, j(2))
1881 jwork = jwork + 1
1882 work(1, jwork) = k1
1883 work(2, jwork) = k2
1884 end do
1885 else if(children(1, j(2)) .eq. 0) then
1886 k2 = j(2)
1887 do k1 = children(1, j(1)), children(2, j(1))
1888 jwork = jwork + 1
1889 work(1, jwork) = k1
1890 work(2, jwork) = k2
1891 end do
1892 else
1893 do k1 = children(1, j(1)), children(2, j(1))
1894 do k2 = children(1, j(2)), children(2, j(2))
1895 jwork = jwork + 1
1896 work(1, jwork) = k1
1897 work(2, jwork) = k2
1898 end do
1899 end do
1900 end if
1901 end if
1902 iwork = iwork + 1
1903 end do
1904 nfar = 0
1905 nnear = 0
1906 do iwork = 1, jwork
1907 if (work(3, iwork) .eq. 1) then
1908 nfar(work(1, iwork)) = nfar(work(1, iwork)) + 1
1909 else if (work(3, iwork) .eq. 2) then
1910 nnear(work(1, iwork)) = nnear(work(1, iwork)) + 1
1911 end if
1912 end do
1913 iwork = jwork + 1
1914 nnfar = sum(nfar)
1915 nnnear = sum(nnear)
1916end subroutine tree_get_farnear_work
1917
1920subroutine tree_get_farnear(jwork, lwork, work, n, nnfar, nfar, sfar, far, &
1921 & nnnear, nnear, snear, near)
1922! Parameters:
1923! jwork: Total number of checked pairs in work array
1924! lwork: Total length of work array
1925! work: Work array itself
1926! n: Number of nodes
1927! nnfar: Total number of all far-field interactions
1928! nfar: Number of far-field interactions of each node
1929! sfar: Index in far array of first far-field node for each node
1930! far: Indexes of far-field nodes
1931! nnnear: Total number of all near-field interactions
1932! nnear: Number of near-field interactions of each node
1933! snear: Index in near array of first near-field node for each node
1934! near: Indexes of near-field nodes
1935 integer, intent(in) :: jwork, lwork, work(3, lwork), n, nnfar, nnnear
1936 integer, intent(in) :: nfar(n), nnear(n)
1937 integer, intent(out) :: sfar(n+1), far(nnfar), snear(n+1), near(nnnear)
1938 integer :: i, j
1939 integer, allocatable :: cfar(:), cnear(:)
1940
1941 allocate(cfar(n+1), cnear(n+1))
1942 sfar(1) = 1
1943 snear(1) = 1
1944 do i = 2, n+1
1945 sfar(i) = sfar(i-1) + nfar(i-1)
1946 snear(i) = snear(i-1) + nnear(i-1)
1947 end do
1948 cfar = sfar
1949 cnear = snear
1950 do i = 1, jwork
1951 if (work(3, i) .eq. 1) then
1952 ! Far
1953 j = work(1, i)
1954 if ((j .gt. n) .or. (j .le. 0)) then
1955 write(*,*) "ALARM", j
1956 end if
1957 far(cfar(j)) = work(2, i)
1958 cfar(j) = cfar(j) + 1
1959 else if (work(3, i) .eq. 2) then
1960 ! Near
1961 j = work(1, i)
1962 if ((j .gt. n) .or. (j .le. 0)) then
1963 write(*,*) "ALARM", j
1964 end if
1965 near(cnear(j)) = work(2, i)
1966 cnear(j) = cnear(j) + 1
1967 end if
1968 end do
1969
1970 deallocate(cfar, cnear)
1971
1972end subroutine tree_get_farnear
1973
1979subroutine constants_free(constants, ddx_error)
1980 implicit none
1981 type(ddx_constants_type), intent(out) :: constants
1982 type(ddx_error_type), intent(inout) :: ddx_error
1983 integer :: istat
1984
1985 istat = 0
1986
1987 if (allocated(constants % vscales)) then
1988 deallocate(constants % vscales, stat=istat)
1989 if (istat .ne. 0) then
1990 call update_error(ddx_error, "`vscales` deallocation failed!")
1991 end if
1992 end if
1993 if (allocated(constants % v4pi2lp1)) then
1994 deallocate(constants % v4pi2lp1, stat=istat)
1995 if (istat .ne. 0) then
1996 call update_error(ddx_error, "`v4pi2lp1` deallocation failed!")
1997 end if
1998 end if
1999 if (allocated(constants % vscales_rel)) then
2000 deallocate(constants % vscales_rel, stat=istat)
2001 if (istat .ne. 0) then
2002 call update_error(ddx_error, "`vscales_rel` deallocation failed!")
2003 end if
2004 end if
2005 if (allocated(constants % vfact)) then
2006 deallocate(constants % vfact, stat=istat)
2007 if (istat .ne. 0) then
2008 call update_error(ddx_error, "`vfact` deallocation failed!")
2009 end if
2010 end if
2011 if (allocated(constants % vcnk)) then
2012 deallocate(constants % vcnk, stat=istat)
2013 if (istat .ne. 0) then
2014 call update_error(ddx_error, "`vcnk` deallocation failed!")
2015 end if
2016 end if
2017 if (allocated(constants % m2l_ztranslate_coef)) then
2018 deallocate(constants % m2l_ztranslate_coef, stat=istat)
2019 if (istat .ne. 0) then
2020 call update_error(ddx_error, &
2021 & "`m2l_ztranslate_coef` deallocation failed!")
2022 end if
2023 end if
2024 if (allocated(constants % m2l_ztranslate_adj_coef)) then
2025 deallocate(constants % m2l_ztranslate_adj_coef, stat=istat)
2026 if (istat .ne. 0) then
2027 call update_error(ddx_error, &
2028 & "`m2l_ztranslate_adj_coef` deallocation failed!")
2029 end if
2030 end if
2031 if (allocated(constants % cgrid)) then
2032 deallocate(constants % cgrid, stat=istat)
2033 if (istat .ne. 0) then
2034 call update_error(ddx_error, "`cgrid` deallocation failed!")
2035 end if
2036 end if
2037 if (allocated(constants % wgrid)) then
2038 deallocate(constants % wgrid, stat=istat)
2039 if (istat .ne. 0) then
2040 call update_error(ddx_error, "`wgrid` deallocation failed!")
2041 end if
2042 end if
2043 if (allocated(constants % vgrid)) then
2044 deallocate(constants % vgrid, stat=istat)
2045 if (istat .ne. 0) then
2046 call update_error(ddx_error, "`vgrid` deallocation failed!")
2047 end if
2048 end if
2049 if (allocated(constants % vwgrid)) then
2050 deallocate(constants % vwgrid, stat=istat)
2051 if (istat .ne. 0) then
2052 call update_error(ddx_error, "`vwgrid` deallocation failed!")
2053 end if
2054 end if
2055 if (allocated(constants % vgrid2)) then
2056 deallocate(constants % vgrid2, stat=istat)
2057 if (istat .ne. 0) then
2058 call update_error(ddx_error, "`vgrid2` deallocation failed!")
2059 end if
2060 end if
2061 if (allocated(constants % pchi)) then
2062 deallocate(constants % pchi, stat=istat)
2063 if (istat .ne. 0) then
2064 call update_error(ddx_error, "`pchi` deallocation failed!")
2065 end if
2066 end if
2067 if (allocated(constants % c_ik)) then
2068 deallocate(constants % c_ik, stat=istat)
2069 if (istat .ne. 0) then
2070 call update_error(ddx_error, "`c_ik` deallocation failed!")
2071 end if
2072 end if
2073 if (allocated(constants % si_ri)) then
2074 deallocate(constants % si_ri, stat=istat)
2075 if (istat .ne. 0) then
2076 call update_error(ddx_error, "`si_ri` deallocation failed!")
2077 end if
2078 end if
2079 if (allocated(constants % di_ri)) then
2080 deallocate(constants % di_ri, stat=istat)
2081 if (istat .ne. 0) then
2082 call update_error(ddx_error, "`di_ri` deallocation failed!")
2083 end if
2084 end if
2085 if (allocated(constants % sk_ri)) then
2086 deallocate(constants % sk_ri, stat=istat)
2087 if (istat .ne. 0) then
2088 call update_error(ddx_error, "`sk_ri` deallocation failed!")
2089 end if
2090 end if
2091 if (allocated(constants % dk_ri)) then
2092 deallocate(constants % dk_ri, stat=istat)
2093 if (istat .ne. 0) then
2094 call update_error(ddx_error, "`dk_ri` deallocation failed!")
2095 end if
2096 end if
2097 if (allocated(constants % termimat)) then
2098 deallocate(constants % termimat, stat=istat)
2099 if (istat .ne. 0) then
2100 call update_error(ddx_error, "`termimat` deallocation failed!")
2101 end if
2102 end if
2103 if (allocated(constants % b)) then
2104 deallocate(constants % b, stat=istat)
2105 if (istat .ne. 0) then
2106 call update_error(ddx_error, "`b` deallocation failed!")
2107 end if
2108 end if
2109 if (allocated(constants % l)) then
2110 deallocate(constants % l, stat=istat)
2111 if (istat .ne. 0) then
2112 call update_error(ddx_error, "`l` deallocation failed!")
2113 end if
2114 end if
2115 if (allocated(constants % inl)) then
2116 deallocate(constants % inl, stat=istat)
2117 if (istat .ne. 0) then
2118 call update_error(ddx_error, "`inl` deallocation failed!")
2119 end if
2120 end if
2121 if (allocated(constants % nl)) then
2122 deallocate(constants % nl, stat=istat)
2123 if (istat .ne. 0) then
2124 call update_error(ddx_error, "`nl` deallocation failed!")
2125 end if
2126 end if
2127 if (allocated(constants % itrnl)) then
2128 deallocate(constants % itrnl, stat=istat)
2129 if (istat .ne. 0) then
2130 call update_error(ddx_error, "`itrnl` deallocation failed!")
2131 end if
2132 end if
2133 if (allocated(constants % fi)) then
2134 deallocate(constants % fi, stat=istat)
2135 if (istat .ne. 0) then
2136 call update_error(ddx_error, "`fi` deallocation failed!")
2137 end if
2138 end if
2139 if (allocated(constants % ui)) then
2140 deallocate(constants % ui, stat=istat)
2141 if (istat .ne. 0) then
2142 call update_error(ddx_error, "`ui` deallocation failed!")
2143 end if
2144 end if
2145 if (allocated(constants % ui_cav)) then
2146 deallocate(constants % ui_cav, stat=istat)
2147 if (istat .ne. 0) then
2148 call update_error(ddx_error, "`ui_cav` deallocation failed!")
2149 end if
2150 end if
2151 if (allocated(constants % zi)) then
2152 deallocate(constants % zi, stat=istat)
2153 if (istat .ne. 0) then
2154 call update_error(ddx_error, "`zi` deallocation failed!")
2155 end if
2156 end if
2157 if (allocated(constants % zi_dr)) then
2158 deallocate(constants % zi_dr, stat=istat)
2159 if (istat .ne. 0) then
2160 call update_error(ddx_error, "`zi_dr` deallocation failed!")
2161 end if
2162 end if
2163 if (allocated(constants % ncav_sph)) then
2164 deallocate(constants % ncav_sph, stat=istat)
2165 if (istat .ne. 0) then
2166 call update_error(ddx_error, "`ncav_sph` deallocation failed!")
2167 end if
2168 end if
2169 if (allocated(constants % ccav)) then
2170 deallocate(constants % ccav, stat=istat)
2171 if (istat .ne. 0) then
2172 call update_error(ddx_error, "`ccav` deallocation failed!")
2173 end if
2174 end if
2175 if (allocated(constants % icav_ia)) then
2176 deallocate(constants % icav_ia, stat=istat)
2177 if (istat .ne. 0) then
2178 call update_error(ddx_error, "`icav_ia` deallocation failed!")
2179 end if
2180 end if
2181 if (allocated(constants % icav_ja)) then
2182 deallocate(constants % icav_ja, stat=istat)
2183 if (istat .ne. 0) then
2184 call update_error(ddx_error, "`icav_ja` deallocation failed!")
2185 end if
2186 end if
2187 if (allocated(constants % rx_prc)) then
2188 deallocate(constants % rx_prc, stat=istat)
2189 if (istat .ne. 0) then
2190 call update_error(ddx_error, "`rx_prc` deallocation failed!")
2191 end if
2192 end if
2193 if (allocated(constants % order)) then
2194 deallocate(constants % order, stat=istat)
2195 if (istat .ne. 0) then
2196 call update_error(ddx_error, "`order` deallocation failed!")
2197 end if
2198 end if
2199 if (allocated(constants % cluster)) then
2200 deallocate(constants % cluster, stat=istat)
2201 if (istat .ne. 0) then
2202 call update_error(ddx_error, "`cluster` deallocation failed!")
2203 end if
2204 end if
2205 if (allocated(constants % children)) then
2206 deallocate(constants % children, stat=istat)
2207 if (istat .ne. 0) then
2208 call update_error(ddx_error, "`children` deallocation failed!")
2209 end if
2210 end if
2211 if (allocated(constants % parent)) then
2212 deallocate(constants % parent, stat=istat)
2213 if (istat .ne. 0) then
2214 call update_error(ddx_error, "`parent` deallocation failed!")
2215 end if
2216 end if
2217 if (allocated(constants % cnode)) then
2218 deallocate(constants % cnode, stat=istat)
2219 if (istat .ne. 0) then
2220 call update_error(ddx_error, "`cnode` deallocation failed!")
2221 end if
2222 end if
2223 if (allocated(constants % rnode)) then
2224 deallocate(constants % rnode, stat=istat)
2225 if (istat .ne. 0) then
2226 call update_error(ddx_error, "`rnode` deallocation failed!")
2227 end if
2228 end if
2229 if (allocated(constants % snode)) then
2230 deallocate(constants % snode, stat=istat)
2231 if (istat .ne. 0) then
2232 call update_error(ddx_error, "`snode` deallocation failed!")
2233 end if
2234 end if
2235 if (allocated(constants % sk_rnode)) then
2236 deallocate(constants % sk_rnode, stat=istat)
2237 if (istat .ne. 0) then
2238 call update_error(ddx_error, "`sk_rnode` deallocation failed!")
2239 end if
2240 end if
2241 if (allocated(constants % si_rnode)) then
2242 deallocate(constants % si_rnode, stat=istat)
2243 if (istat .ne. 0) then
2244 call update_error(ddx_error, "`si_rnode` deallocation failed!")
2245 end if
2246 end if
2247 if (allocated(constants % nfar)) then
2248 deallocate(constants % nfar, stat=istat)
2249 if (istat .ne. 0) then
2250 call update_error(ddx_error, "`nfar` deallocation failed!")
2251 end if
2252 end if
2253 if (allocated(constants % nnear)) then
2254 deallocate(constants % nnear, stat=istat)
2255 if (istat .ne. 0) then
2256 call update_error(ddx_error, "`nnear` deallocation failed!")
2257 end if
2258 end if
2259 if (allocated(constants % far)) then
2260 deallocate(constants % far, stat=istat)
2261 if (istat .ne. 0) then
2262 call update_error(ddx_error, "`far` deallocation failed!")
2263 end if
2264 end if
2265 if (allocated(constants % near)) then
2266 deallocate(constants % near, stat=istat)
2267 if (istat .ne. 0) then
2268 call update_error(ddx_error, "`near` deallocation failed!")
2269 end if
2270 end if
2271 if (allocated(constants % sfar)) then
2272 deallocate(constants % sfar, stat=istat)
2273 if (istat .ne. 0) then
2274 call update_error(ddx_error, "`sfar` deallocation failed!")
2275 end if
2276 end if
2277 if (allocated(constants % snear)) then
2278 deallocate(constants % snear, stat=istat)
2279 if (istat .ne. 0) then
2280 call update_error(ddx_error, "`snear` deallocation failed!")
2281 end if
2282 end if
2283 if (allocated(constants % m2l_ztranslate_coef)) then
2284 deallocate(constants % m2l_ztranslate_coef, stat=istat)
2285 if (istat .ne. 0) then
2286 call update_error(ddx_error, "`m2l_ztranslate_coef` " // &
2287 & "deallocation failed!")
2288 end if
2289 end if
2290 if (allocated(constants % m2l_ztranslate_adj_coef)) then
2291 deallocate(constants % m2l_ztranslate_adj_coef, stat=istat)
2292 if (istat .ne. 0) then
2293 call update_error(ddx_error, "`m2l_ztranslate_adj_coef` " // &
2294 & "deallocation failed!")
2295 end if
2296 end if
2297
2298 call switching_free(constants%switching, ddx_error)
2299
2300end subroutine constants_free
2301
2302subroutine switching_init(params, constants, switching, ddx_error)
2303 type(ddx_params_type), intent(in) :: params
2304 type(ddx_constants_type), intent(in) :: constants
2305 type(ddx_switching_type), intent(inout) :: switching
2306 type(ddx_error_type), intent(inout) :: ddx_error
2307 integer :: info, isph, n, icav, ncav
2308 real(dp) :: d_ni, f_ni, u_ni
2309
2310 allocate(switching%u_ni(params%ngrid, params%nsph), &
2311 & switching%d_ni(params%ngrid, params%nsph), &
2312 & switching%f_ni(params%ngrid, params%nsph), &
2313 & stat=info)
2314 if (info.ne.0) then
2315 call update_error(ddx_error, &
2316 & "Allocation failed in switching_init")
2317 return
2318 end if
2319
2320 do isph = 1, params%nsph
2321 do n = 1, params%ngrid
2322 call compute_d_and_f(params, constants, isph, n, d_ni, f_ni)
2323 switching%d_ni(n, isph) = d_ni
2324 switching%f_ni(n, isph) = f_ni
2325 end do
2326 end do
2327
2328 ncav = 0
2329 do isph = 1, params%nsph
2330 do n = 1, params%ngrid
2331 u_ni = compute_u(params, constants, isph, n)
2332 if (u_ni.ne.zero) ncav = ncav + 1
2333 end do
2334 end do
2335 switching%ncav = ncav
2336
2337 allocate(switching%u_i_cav(switching%ncav), stat=info)
2338 if (info.ne.0) then
2339 call update_error(ddx_error, &
2340 & "Allocation failed in switching_init")
2341 return
2342 end if
2343
2344 icav = 0
2345 do isph = 1, params%nsph
2346 do n = 1, params%ngrid
2347 u_ni = compute_u(params, constants, isph, n)
2348 switching%u_ni(n, isph) = u_ni
2349 if (u_ni.ne.zero) then
2350 icav = icav + 1
2351 switching%u_i_cav(icav) = u_ni
2352 end if
2353 end do
2354 end do
2355
2356 if (params%force.eq.1) then
2357 allocate(switching%grad_u_ni(3, params%ngrid, params%nsph), &
2358 & switching%dr_u_ni(params%ngrid, params%nsph), &
2359 & stat=info)
2360 if (info.ne.0) then
2361 call update_error(ddx_error, &
2362 & "Allocation failed in switching_init")
2363 return
2364 end if
2365
2366 do isph = 1, params%nsph
2367 do n = 1, params%ngrid
2368 call compute_grad_i_u_i(params, constants, isph, n, &
2369 & switching%grad_u_ni(:, n, isph), &
2370 & switching%dr_u_ni(n, isph))
2371 end do
2372 end do
2373 end if
2374
2375end subroutine switching_init
2376
2377subroutine switching_free(switching, ddx_error)
2378 type(ddx_switching_type), intent(inout) :: switching
2379 type(ddx_error_type), intent(inout) :: ddx_error
2380 integer :: info
2381 if (allocated(switching%u_ni)) then
2382 deallocate(switching%u_ni, stat=info)
2383 if (info.ne.0) then
2384 call update_error(ddx_error, &
2385 & "Deallocation failed in switching_free")
2386 return
2387 end if
2388 end if
2389 if (allocated(switching%u_i_cav)) then
2390 deallocate(switching%u_i_cav, stat=info)
2391 if (info.ne.0) then
2392 call update_error(ddx_error, &
2393 & "Deallocation failed in switching_free")
2394 return
2395 end if
2396 end if
2397 if (allocated(switching%d_ni)) then
2398 deallocate(switching%d_ni, stat=info)
2399 if (info.ne.0) then
2400 call update_error(ddx_error, &
2401 & "Deallocation failed in switching_free")
2402 return
2403 end if
2404 end if
2405 if (allocated(switching%f_ni)) then
2406 deallocate(switching%f_ni, stat=info)
2407 if (info.ne.0) then
2408 call update_error(ddx_error, &
2409 & "Deallocation failed in switching_free")
2410 return
2411 end if
2412 end if
2413 if (allocated(switching%grad_u_ni)) then
2414 deallocate(switching%grad_u_ni, stat=info)
2415 if (info.ne.0) then
2416 call update_error(ddx_error, &
2417 & "Deallocation failed in switching_free")
2418 return
2419 end if
2420 end if
2421 if (allocated(switching%dr_u_ni)) then
2422 deallocate(switching%dr_u_ni, stat=info)
2423 if (info.ne.0) then
2424 call update_error(ddx_error, &
2425 & "Deallocation failed in switching_free")
2426 return
2427 end if
2428 end if
2429end subroutine switching_free
2430
2431real(dp) function compute_t(params, constants, isph, jsph, n)
2432 type(ddx_params_type), intent(in) :: params
2433 type(ddx_constants_type), intent(in) :: constants
2434 integer, intent(in) :: isph, jsph, n
2435 associate(cgrid => constants%cgrid, rsph => params%rsph, &
2436 & csph => params%csph)
2437 compute_t = sqrt((cgrid(1,n)*rsph(isph) + csph(1,isph) &
2438 & - csph(1,jsph))**2 + (cgrid(2,n)*rsph(isph) &
2439 & + csph(2,isph) - csph(2,jsph))**2 &
2440 & + (cgrid(3,n)*rsph(isph) + csph(3,isph) &
2441 & - csph(3,jsph))**2)/rsph(jsph)
2442 end associate
2443end function compute_t
2444
2445real(dp) function compute_chi(params, constants, isph, jsph, n)
2446 type(ddx_params_type), intent(in) :: params
2447 type(ddx_constants_type), intent(in) :: constants
2448 integer, intent(in) :: isph, jsph, n
2449 real(dp) :: t_ijn
2450 t_ijn = compute_t(params, constants, isph, jsph, n)
2451 compute_chi = fsw(t_ijn, params%se, params%eta)
2452end function compute_chi
2453
2454subroutine compute_d_and_f(params, constants, isph, n, d, f)
2455 type(ddx_params_type), intent(in) :: params
2456 type(ddx_constants_type), intent(in) :: constants
2457 integer, intent(in) :: isph, n
2458 real(dp), intent(out) :: d, f
2459 integer :: ij, jsph
2460 real(dp) :: chi_ijn
2461 d = one
2462 f = zero
2463 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2464 jsph = constants%nl(ij)
2465 chi_ijn = compute_chi(params, constants, isph, jsph, n)
2466 d = d*(one - chi_ijn)
2467 f = f + chi_ijn
2468 end do
2469end subroutine compute_d_and_f
2470
2471real(dp) function compute_d(params, constants, isph, n)
2472 type(ddx_params_type), intent(in) :: params
2473 type(ddx_constants_type), intent(in) :: constants
2474 integer, intent(in) :: isph, n
2475 integer :: ij, jsph
2476 real(dp) :: chi_ijn
2477 compute_d = one
2478 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2479 jsph = constants%nl(ij)
2480 chi_ijn = compute_chi(params, constants, isph, jsph, n)
2481 compute_d = compute_d*(one - chi_ijn)
2482 end do
2483end function
2484
2485real(dp) function compute_f(params, constants, isph, n)
2486 type(ddx_params_type), intent(in) :: params
2487 type(ddx_constants_type), intent(in) :: constants
2488 integer, intent(in) :: isph, n
2489 integer :: ij, jsph
2490
2491 compute_f = zero
2492 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2493 jsph = constants%nl(ij)
2494 compute_f = compute_f + compute_chi(params, constants, isph, jsph, n)
2495 end do
2496end function compute_f
2497
2498real(dp) function compute_omega(params, constants, isph, jsph, n)
2499 type(ddx_params_type), intent(in) :: params
2500 type(ddx_constants_type), intent(in) :: constants
2501 integer, intent(in) :: isph, jsph, n
2502 real(dp) :: chi_nij, d_ni, f_ni
2503
2504 chi_nij = compute_chi(params, constants, isph, jsph, n)
2505 d_ni = constants%switching%d_ni(n, isph)
2506 f_ni = constants%switching%f_ni(n, isph)
2507 compute_omega = chi_nij/(d_ni + f_ni)
2508end function compute_omega
2509
2510real(dp) function compute_u(params, constants, isph, n)
2511 type(ddx_params_type), intent(in) :: params
2512 type(ddx_constants_type), intent(in) :: constants
2513 integer, intent(in) :: isph, n
2514 integer :: ij, jsph
2515 compute_u = one
2516 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2517 jsph = constants%nl(ij)
2518 compute_u = compute_u - compute_omega(params, constants, isph, jsph, n)
2519 end do
2520 ! A little hack here: we have plenty of checks `if (u(n, isph).eq.zero)`,
2521 ! which we don't want to trigger with a noisy u. Better would be
2522 ! making the checks more robust, but this suffices for now.
2523 if (abs(compute_u) .le. 1d-15) compute_u = zero
2524end function compute_u
2525
2526subroutine compute_grad_i_u_i(params, constants, isph, n, grad, dr)
2527 type(ddx_params_type), intent(in) :: params
2528 type(ddx_constants_type), intent(in) :: constants
2529 real(dp), intent(out) :: grad(3), dr
2530 integer, intent(in) :: isph, n
2531 integer :: ij, jsph
2532 real(dp) :: vij(3), vvij, tij, sij(3), grad_tij(3), dr_tij, &
2533 & d_i, f_i, chi_ij, a, b, c, grad_p
2534
2535 grad(:) = zero
2536 dr = zero
2537 d_i = constants%switching%d_ni(n, isph)
2538 f_i = constants%switching%f_ni(n, isph)
2539 a = one/(d_i + f_i)**2
2540 b = f_i*d_i
2541 do ij = constants%inl(isph), constants%inl(isph + 1) - 1
2542 jsph = constants%nl(ij)
2543 vij = params % csph(:,isph) + &
2544 & params % rsph(isph)*constants % cgrid(:,n) - &
2545 & params % csph(:,jsph)
2546 vvij = sqrt(vij(1)*vij(1) + vij(2)*vij(2) + vij(3)*vij(3))
2547 tij = vvij/params % rsph(jsph)
2548 sij = vij/vvij
2549
2550 grad_tij(:) = sij(:)/params%rsph(jsph)
2551 dr_tij = dot_product(vij, constants%cgrid(:,n)) &
2552 & /(vvij*params%rsph(jsph))
2553
2554 grad_p = dfsw(tij, params % se, params % eta)
2555 chi_ij = fsw(tij, params % se, params % eta)
2556
2557 c = a*(d_i + b/(one - chi_ij))*grad_p
2558
2559 grad(:) = grad(:) + c*grad_tij(:)
2560 dr = dr + c*dr_tij
2561 end do
2562
2563end subroutine compute_grad_i_u_i
2564
2565end module ddx_constants
2566
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:1829
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:1503
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::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:1980
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:1594
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:1720
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:1922
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