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