ddx 0.6.8
Libary for domain-decomposition methods for polarizable continuum models
ddx_core.f90
Go to the documentation of this file.
1
11
13module ddx_core
14! Get ddx_params_type and all compile-time definitions
15use ddx_parameters
16! Get ddx_constants_type and all run-time constants
17use ddx_constants
18! Get ddx_workspace_type for temporary buffers
19use ddx_workspace
20! Get harmonics-related functions
21use ddx_harmonics
22! Enable OpenMP
23use omp_lib
24implicit none
25
27
33type ddx_state_type
34 !! High-level entities to be stored and accessed by users
35
36 !!
37 !! Quantities common to all models.
38 !!
42 real(dp), allocatable :: psi(:, :)
45 real(dp), allocatable :: phi_cav(:)
48 real(dp), allocatable :: gradphi_cav(:, :)
50 real(dp), allocatable :: phi_grid(:, :)
52 real(dp), allocatable :: zeta(:)
54 real(dp) :: force_time
55
56 !!
57 !! ddCOSMO quantities (used also by ddPCM).
58 !!
61 real(dp), allocatable :: phi(:, :)
64 real(dp), allocatable :: xs(:, :)
67 integer :: xs_niter
71 real(dp), allocatable :: xs_rel_diff(:)
74 real(dp) :: xs_time
77 real(dp), allocatable :: sgrid(:, :)
80 real(dp), allocatable :: s(:, :)
83 integer :: s_niter
87 real(dp), allocatable :: s_rel_diff(:)
90 real(dp) :: s_time
91
92 !!
93 !! ddPCM specific quantities.
94 !!
97 real(dp), allocatable :: phiinf(:, :)
100 real(dp), allocatable :: phieps(:, :)
103 integer :: phieps_niter
107 real(dp), allocatable :: phieps_rel_diff(:)
110 real(dp) :: phieps_time
114 real(dp), allocatable :: g(:, :)
117 real(dp), allocatable :: y(:, :)
120 integer :: y_niter
123 real(dp), allocatable :: y_rel_diff(:)
125 real(dp) :: y_time
128 real(dp), allocatable :: ygrid(:, :)
132 real(dp), allocatable :: q(:, :)
135 real(dp), allocatable :: qgrid(:, :)
136
137 !!
138 !! ddLPB quantities
139 !!
142 real(dp), allocatable :: rhs_lpb(:,:,:)
145 real(dp), allocatable :: rhs_adj_lpb(:,:,:)
148 real(dp), allocatable :: x_lpb(:,:,:)
152 real(dp), allocatable :: x_adj_lpb(:,:,:)
155 real(dp), allocatable :: g_lpb(:,:)
158 real(dp), allocatable :: f_lpb(:,:)
161 integer :: x_lpb_niter
164 integer :: x_adj_lpb_niter
167 real(dp) :: x_lpb_time
170 real(dp) :: x_adj_lpb_time
174 real(dp), allocatable :: x_lpb_rel_diff(:)
178 real(dp), allocatable :: x_adj_lpb_rel_diff(:)
181 real(dp), allocatable :: zeta_dip(:, :)
182 real(dp), allocatable :: x_adj_re_grid(:, :)
183 real(dp), allocatable :: x_adj_r_grid(:, :)
184 real(dp), allocatable :: x_adj_e_grid(:, :)
185 real(dp), allocatable :: phi_n(:, :)
186 real(dp) :: hsp_time
187 real(dp) :: hsp_adj_time
188
189 !! Some flags to see if quantities are properly initialized
190 logical :: rhs_done
191 logical :: adjoint_rhs_done
192 logical :: guess_done
193 logical :: solved
194 logical :: adjoint_guess_done
195 logical :: adjoint_solved
196
197end type ddx_state_type
198
203type ddx_electrostatics_type
205 real(dp), allocatable :: phi_cav(:)
208 real(dp), allocatable :: e_cav(:, :)
211 real(dp), allocatable :: g_cav(:, :, :)
213 logical :: do_phi = .false.
215 logical :: do_e = .false.
217 logical :: do_g = .false.
218end type ddx_electrostatics_type
219
222type ddx_type
223 !! New types inside the old one for an easier shift to the new design
224 type(ddx_params_type) :: params
225 type(ddx_constants_type) :: constants
226 type(ddx_workspace_type) :: workspace
227end type ddx_type
228
229contains
230
231!------------------------------------------------------------------------------
265!------------------------------------------------------------------------------
266subroutine allocate_model(nsph, x, y, z, rvdw, model, lmax, ngrid, force, fmm, pm, &
267 & pl, se, eta, eps, kappa, matvecmem, maxiter, jacobi_ndiis, nproc, &
268 & output_filename, ddx_data, ddx_error)
269 ! Inputs
270 implicit none
271 integer, intent(in) :: nsph, model, lmax, force, fmm, pm, pl, &
272 & matvecmem, maxiter, jacobi_ndiis, &
273 & ngrid, nproc
274 real(dp), intent(in):: x(nsph), y(nsph), z(nsph), &
275 & rvdw(nsph), se, eta, eps, kappa
276 character(len=255), intent(in) :: output_filename
277 ! Output
278 type(ddx_type), target, intent(out) :: ddx_data
279 type(ddx_error_type), intent(inout) :: ddx_error
280 ! Local variables
281 real(dp), allocatable :: csph(:, :)
282 integer :: istatus
283 if (ddx_error % flag .ne. 0) then
284 call update_error(ddx_error, "allocate_model received input in error state, " // &
285 & "exiting")
286 return
287 end if
288 ! Init underlying objects
289 allocate(csph(3, nsph), stat=istatus)
290 if (istatus .ne. 0) then
291 call update_error(ddx_error, "Allocation failed in allocate_model")
292 return
293 end if
294 csph(1, :) = x
295 csph(2, :) = y
296 csph(3, :) = z
297 call params_init(model, force, eps, kappa, eta, se, lmax, ngrid, &
298 & matvecmem, maxiter, jacobi_ndiis, fmm, pm, pl, nproc, nsph, &
299 & csph, rvdw, output_filename, ddx_data % params, ddx_error)
300 if (ddx_error % flag .ne. 0) then
301 call update_error(ddx_error, "params_init returned an error, exiting")
302 return
303 end if
304 call constants_init(ddx_data % params, ddx_data % constants, ddx_error)
305 if (ddx_error % flag .ne. 0) then
306 call update_error(ddx_error, "constants_init returned an error, exiting")
307 return
308 end if
309 call workspace_init(ddx_data % params, ddx_data % constants, &
310 & ddx_data % workspace, ddx_error)
311 if (ddx_error % flag .ne. 0) then
312 call update_error(ddx_error, "workspace_init returned an error, exiting")
313 return
314 end if
315 deallocate(csph, stat=istatus)
316 if (istatus .ne. 0) then
317 call update_error(ddx_error, "Deallocation failed in allocate_model")
318 return
319 end if
320end subroutine allocate_model
321
322!------------------------------------------------------------------------------
327!------------------------------------------------------------------------------
328subroutine deallocate_model(ddx_data, ddx_error)
329 implicit none
330 ! Input/output
331 type(ddx_type), intent(inout) :: ddx_data
332 type(ddx_error_type), intent(inout) :: ddx_error
333 ! Local variables
334 call workspace_free(ddx_data % workspace, ddx_error)
335 call constants_free(ddx_data % constants, ddx_error)
336 call params_free(ddx_data % params, ddx_error)
337end subroutine deallocate_model
338
348subroutine allocate_electrostatics(params, constants, electrostatics, ddx_error)
349 implicit none
350 type(ddx_params_type), intent(in) :: params
351 type(ddx_constants_type), intent(in) :: constants
352 type(ddx_electrostatics_type), intent(out) :: electrostatics
353 type(ddx_error_type), intent(inout) :: ddx_error
354
355 ! local variables
356 integer :: info
357
358 ! Figure out the required properties
359 if (params % model .eq. 3) then
360 if (params % force .eq. 1) then
361 electrostatics % do_phi = .true.
362 electrostatics % do_e = .true.
363 electrostatics % do_g = .true.
364 else
365 electrostatics % do_phi = .true.
366 electrostatics % do_e = .true.
367 electrostatics % do_g = .false.
368 end if
369 else
370 if (params % force .eq. 1) then
371 electrostatics % do_phi = .true.
372 electrostatics % do_e = .true.
373 electrostatics % do_g = .false.
374 else
375 electrostatics % do_phi = .true.
376 electrostatics % do_e = .false.
377 electrostatics % do_g = .false.
378 end if
379 end if
380
381 ! Allocate the arrays for the required electrostatic properties
382 if (electrostatics % do_phi) then
383 allocate(electrostatics % phi_cav(constants % ncav), stat=info)
384 if (info .ne. 0) then
385 call update_error(ddx_error, &
386 & "Allocation failed in allocate_electrostatics")
387 return
388 end if
389 end if
390 if (electrostatics % do_e) then
391 allocate(electrostatics % e_cav(3, constants % ncav), stat=info)
392 if (info .ne. 0) then
393 call update_error(ddx_error, &
394 & "Allocation failed in allocate_electrostatics")
395 return
396 end if
397 end if
398 if (electrostatics % do_g) then
399 allocate(electrostatics % g_cav(3, 3, constants % ncav), stat=info)
400 if (info .ne. 0) then
401 call update_error(ddx_error, &
402 & "Allocation failed in allocate_electrostatics")
403 return
404 end if
405 end if
406end subroutine allocate_electrostatics
407
414subroutine deallocate_electrostatics(electrostatics, ddx_error)
415 implicit none
416 type(ddx_electrostatics_type), intent(inout) :: electrostatics
417 type(ddx_error_type), intent(inout) :: ddx_error
418
419 ! local variables
420 integer :: info
421
422 ! Allocate the arrays for the required electrostatic properties
423 if (allocated(electrostatics % phi_cav)) then
424 deallocate(electrostatics % phi_cav, stat=info)
425 if (info .ne. 0) then
426 call update_error(ddx_error, &
427 & "Deallocation failed in deallocate_electrostatics")
428 end if
429 end if
430 if (allocated(electrostatics % e_cav)) then
431 deallocate(electrostatics % e_cav, stat=info)
432 if (info .ne. 0) then
433 call update_error(ddx_error, &
434 & "Deallocation failed in deallocate_electrostatics")
435 end if
436 end if
437 if (allocated(electrostatics % g_cav)) then
438 deallocate(electrostatics % g_cav, stat=info)
439 if (info .ne. 0) then
440 call update_error(ddx_error, &
441 & "Deallocation failed in deallocate_electrostatics")
442 end if
443 end if
444end subroutine deallocate_electrostatics
445
454subroutine allocate_state(params, constants, state, ddx_error)
455 implicit none
456 type(ddx_params_type), intent(in) :: params
457 type(ddx_constants_type), intent(in) :: constants
458 type(ddx_state_type), intent(out) :: state
459 type(ddx_error_type), intent(inout) :: ddx_error
460 integer :: istatus
461
462 state % rhs_done = .false.
463 state % adjoint_rhs_done = .false.
464 state % guess_done = .false.
465 state % solved = .false.
466 state % adjoint_guess_done = .false.
467 state % adjoint_solved = .false.
468
469 allocate(state % psi(constants % nbasis, params % nsph), stat=istatus, source=0.0_dp)
470 if (istatus .ne. 0) then
471 call update_error(ddx_error, "allocate_state: `psi` allocation failed")
472 return
473 end if
474 allocate(state % phi_cav(constants % ncav), stat=istatus, source=0.0_dp)
475 if (istatus .ne. 0) then
476 call update_error(ddx_error, "allocate_state: `phi_cav` allocation failed")
477 return
478 end if
479 allocate(state % gradphi_cav(3, constants % ncav), stat=istatus, source=0.0_dp)
480 if (istatus .ne. 0) then
481 call update_error(ddx_error, &
482 & "allocate_state: `gradphi_cav` allocation failed")
483 return
484 end if
485 allocate(state % q(constants % nbasis, &
486 & params % nsph), stat=istatus, source=0.0_dp)
487 if (istatus .ne. 0) then
488 call update_error(ddx_error, "allocate_state: `q` " // &
489 & "allocation failed")
490 return
491 end if
492
493 ! COSMO model
494 if (params % model .eq. 1) then
495 allocate(state % phi_grid(params % ngrid, &
496 & params % nsph), stat=istatus, source=0.0_dp)
497 if (istatus .ne. 0) then
498 call update_error(ddx_error, "allocate_state: `phi_grid` " // &
499 & "allocation failed")
500 return
501 end if
502 allocate(state % phi(constants % nbasis, &
503 & params % nsph), stat=istatus, source=0.0_dp)
504 if (istatus .ne. 0) then
505 call update_error(ddx_error, "allocate_state: `phi` " // &
506 & "allocation failed")
507 return
508 end if
509 allocate(state % xs(constants % nbasis, &
510 & params % nsph), stat=istatus, source=0.0_dp)
511 if (istatus .ne. 0) then
512 call update_error(ddx_error, "allocate_state: `xs` " // &
513 & "allocation failed")
514 return
515 end if
516 allocate(state % xs_rel_diff(params % maxiter), &
517 & stat=istatus, source=0.0_dp)
518 if (istatus .ne. 0) then
519 call update_error(ddx_error, "allocate_state: `xs_rel_diff` " // &
520 & "allocation failed")
521 return
522 end if
523 allocate(state % s(constants % nbasis, &
524 & params % nsph), stat=istatus, source=0.0_dp)
525 if (istatus .ne. 0) then
526 call update_error(ddx_error, "allocate_state: `s` " // &
527 & "allocation failed")
528 return
529 end if
530 allocate(state % s_rel_diff(params % maxiter), &
531 & stat=istatus, source=0.0_dp)
532 if (istatus .ne. 0) then
533 call update_error(ddx_error, "allocate_state: `s_rel_diff` " // &
534 & "allocation failed")
535 return
536 end if
537 allocate(state % sgrid(params % ngrid, &
538 & params % nsph), stat=istatus, source=0.0_dp)
539 if (istatus .ne. 0) then
540 call update_error(ddx_error, "allocate_state: `sgrid` " // &
541 & "allocation failed")
542 return
543 end if
544 allocate(state % zeta(constants % ncav), stat=istatus, source=0.0_dp)
545 if (istatus .ne. 0) then
546 call update_error(ddx_error, "allocate_state: `zeta` " // &
547 & "allocation failed")
548 return
549 end if
550 ! PCM model
551 else if (params % model .eq. 2) then
552 allocate(state % phi_grid(params % ngrid, &
553 & params % nsph), stat=istatus, source=0.0_dp)
554 if (istatus .ne. 0) then
555 call update_error(ddx_error, "allocate_state: `phi_grid` " // &
556 & "allocation failed")
557 return
558 end if
559 allocate(state % phi(constants % nbasis, &
560 & params % nsph), stat=istatus, source=0.0_dp)
561 if (istatus .ne. 0) then
562 call update_error(ddx_error, "allocate_state: `phi` " // &
563 & "allocation failed")
564 return
565 end if
566 allocate(state % phiinf(constants % nbasis, &
567 & params % nsph), stat=istatus, source=0.0_dp)
568 if (istatus .ne. 0) then
569 call update_error(ddx_error, "allocate_state: `phiinf` " // &
570 & "allocation failed")
571 return
572 end if
573 allocate(state % phieps(constants % nbasis, &
574 & params % nsph), stat=istatus, source=0.0_dp)
575 if (istatus .ne. 0) then
576 call update_error(ddx_error, "allocate_state: `phieps` " // &
577 & "allocation failed")
578 return
579 end if
580 allocate(state % phieps_rel_diff(params % maxiter), &
581 & stat=istatus, source=0.0_dp)
582 if (istatus .ne. 0) then
583 call update_error(ddx_error, "allocate_state: `xs_rel_diff` " // &
584 & "allocation failed")
585 return
586 end if
587 allocate(state % xs(constants % nbasis, &
588 & params % nsph), stat=istatus, source=0.0_dp)
589 if (istatus .ne. 0) then
590 call update_error(ddx_error, "allocate_state: `xs` " // &
591 & "allocation failed")
592 return
593 end if
594 allocate(state % xs_rel_diff(params % maxiter), &
595 & stat=istatus, source=0.0_dp)
596 if (istatus .ne. 0) then
597 call update_error(ddx_error, "allocate_state: `xs_rel_diff` " // &
598 & "allocation failed")
599 return
600 end if
601 allocate(state % s(constants % nbasis, &
602 & params % nsph), stat=istatus, source=0.0_dp)
603 if (istatus .ne. 0) then
604 call update_error(ddx_error, "allocate_state: `s` " // &
605 & "allocation failed")
606 return
607 end if
608 allocate(state % s_rel_diff(params % maxiter), &
609 & stat=istatus, source=0.0_dp)
610 if (istatus .ne. 0) then
611 call update_error(ddx_error, "allocate_state: `xs_rel_diff` " // &
612 & "allocation failed")
613 return
614 end if
615 allocate(state % sgrid(params % ngrid, &
616 & params % nsph), stat=istatus, source=0.0_dp)
617 if (istatus .ne. 0) then
618 call update_error(ddx_error, "allocate_state: `sgrid` " // &
619 & "allocation failed")
620 return
621 end if
622 allocate(state % y(constants % nbasis, &
623 & params % nsph), stat=istatus, source=0.0_dp)
624 if (istatus .ne. 0) then
625 call update_error(ddx_error, "allocate_state: `y` " // &
626 & "allocation failed")
627 return
628 end if
629 allocate(state % y_rel_diff(params % maxiter), &
630 & stat=istatus, source=0.0_dp)
631 if (istatus .ne. 0) then
632 call update_error(ddx_error, "allocate_state: `y_rel_diff` " // &
633 & "allocation failed")
634 return
635 end if
636 allocate(state % ygrid(params % ngrid, &
637 & params % nsph), stat=istatus, source=0.0_dp)
638 if (istatus .ne. 0) then
639 call update_error(ddx_error, "allocate_state: `ygrid` " // &
640 & "allocation failed")
641 return
642 end if
643 allocate(state % g(constants % nbasis, &
644 & params % nsph), stat=istatus, source=0.0_dp)
645 if (istatus .ne. 0) then
646 call update_error(ddx_error, "allocate_state: `g` " // &
647 & "allocation failed")
648 return
649 end if
650 allocate(state % qgrid(params % ngrid, &
651 & params % nsph), stat=istatus, source=0.0_dp)
652 if (istatus .ne. 0) then
653 call update_error(ddx_error, "allocate_state: `qgrid` " // &
654 & "allocation failed")
655 return
656 end if
657 allocate(state % zeta(constants % ncav), stat=istatus, source=0.0_dp)
658 if (istatus .ne. 0) then
659 call update_error(ddx_error, "allocate_state: `zeta` " // &
660 & "allocation failed")
661 return
662 end if
663 ! LPB model
664 else if (params % model .eq. 3) then
665 allocate(state % phi_grid(params % ngrid, &
666 & params % nsph), stat=istatus, source=0.0_dp)
667 if (istatus .ne. 0) then
668 call update_error(ddx_error, "allocate_state: `phi_grid` " // &
669 & "allocation failed")
670 return
671 end if
672 allocate(state % phi(constants % nbasis, &
673 & params % nsph), stat=istatus, source=0.0_dp)
674 if (istatus .ne. 0) then
675 call update_error(ddx_error, "allocate_state: `phi` " // &
676 & "allocation failed")
677 return
678 end if
679 allocate(state % zeta(constants % ncav), stat=istatus, source=0.0_dp)
680 if (istatus .ne. 0) then
681 call update_error(ddx_error, "allocate_state: `zeta` " // &
682 & "allocation failed")
683 !write(*, *) "ddx_Error in allocation of M2P matrices"
684 return
685 end if
686 allocate(state % x_lpb_rel_diff(params % maxiter), &
687 & stat=istatus, source=0.0_dp)
688 if (istatus .ne. 0) then
689 call update_error(ddx_error, "allocate_state: `x_lpb_rel_diff` " // &
690 & "allocation failed")
691 return
692 end if
693 allocate(state % rhs_lpb(constants % nbasis, &
694 & params % nsph, 2), &
695 & stat=istatus, source=0.0_dp)
696 if (istatus .ne. 0) then
697 call update_error(ddx_error, "allocate_state: `rhs_lpb` " // &
698 & "allocation failed")
699 return
700 end if
701 allocate(state % rhs_adj_lpb(constants % nbasis, &
702 & params % nsph, 2), &
703 & stat=istatus, source=0.0_dp)
704 if (istatus .ne. 0) then
705 call update_error(ddx_error, "allocate_state: `rhs_adj_lpb` " // &
706 & "allocation failed")
707 return
708 end if
709 allocate(state % x_lpb(constants % nbasis, &
710 & params % nsph, 2), &
711 & stat=istatus, source=0.0_dp)
712 if (istatus .ne. 0) then
713 call update_error(ddx_error, "allocate_state: `x_lpb` " // &
714 & "allocation failed")
715 return
716 end if
717 allocate(state % x_adj_lpb(constants % nbasis, &
718 & params % nsph, 2), stat=istatus, source=0.0_dp)
719 if (istatus .ne. 0) then
720 call update_error(ddx_error, "allocate_state: `x_adj_lpb` " // &
721 & "allocation failed")
722 return
723 end if
724 allocate(state % x_adj_lpb_rel_diff(params % maxiter), &
725 & stat=istatus, source=0.0_dp)
726 if (istatus .ne. 0) then
727 call update_error(ddx_error, &
728 & "allocate_state: `x_adj_lpb_rel_diff` allocation failed")
729 return
730 end if
731 allocate(state % g_lpb(params % ngrid, &
732 & params % nsph), stat=istatus, source=0.0_dp)
733 if (istatus .ne. 0) then
734 call update_error(ddx_error, "allocate_state: `g_lpb` " // &
735 & "allocation failed")
736 return
737 end if
738 allocate(state % f_lpb(params % ngrid, &
739 & params % nsph), stat=istatus, source=0.0_dp)
740 if (istatus .ne. 0) then
741 call update_error(ddx_error, "allocate_state: `f_lpb` " // &
742 & "allocation failed")
743 return
744 end if
745 allocate(state % zeta_dip(3, constants % ncav), stat=istatus, source=0.0_dp)
746 if (istatus .ne. 0) then
747 call update_error(ddx_error, "allocate_state: `zeta_dip` " // &
748 & "allocation failed")
749 return
750 end if
751 allocate(state % x_adj_re_grid(params % ngrid, params % nsph), &
752 & stat=istatus, source=0.0_dp)
753 if (istatus .ne. 0) then
754 call update_error(ddx_error, "allocate_state: `x_adj_re_grid` " // &
755 & "allocation failed")
756 return
757 end if
758 allocate(state % x_adj_r_grid(params % ngrid, params % nsph), &
759 & stat=istatus, source=0.0_dp)
760 if (istatus .ne. 0) then
761 call update_error(ddx_error, "allocate_state: `x_adj_r_grid` " // &
762 & "allocation failed")
763 return
764 end if
765 allocate(state % x_adj_e_grid(params % ngrid, params % nsph), &
766 & stat=istatus, source=0.0_dp)
767 if (istatus .ne. 0) then
768 call update_error(ddx_error, "allocate_state: `x_adj_e_grid` " // &
769 & "allocation failed")
770 return
771 end if
772 allocate(state % phi_n(params % ngrid, params % nsph), &
773 & stat=istatus, source=0.0_dp)
774 if (istatus .ne. 0) then
775 call update_error(ddx_error, "allocate_state: `phi_n` " // &
776 & "allocation failed")
777 return
778 end if
779 end if
780end subroutine allocate_state
781
782
789subroutine deallocate_state(state, ddx_error)
790 implicit none
791 type(ddx_state_type), intent(inout) :: state
792 type(ddx_error_type), intent(inout) :: ddx_error
793 integer :: istatus
794
795 if (allocated(state % phi_cav)) then
796 deallocate(state % phi_cav, stat=istatus)
797 if (istatus .ne. 0) then
798 call update_error(ddx_error, "`phi_cav` deallocation failed!")
799 return
800 endif
801 end if
802 if (allocated(state % gradphi_cav)) then
803 deallocate(state % gradphi_cav, stat=istatus)
804 if (istatus .ne. 0) then
805 call update_error(ddx_error, "`gradphi_cav` deallocation failed!")
806 return
807 endif
808 end if
809 if (allocated(state % psi)) then
810 deallocate(state % psi, stat=istatus)
811 if (istatus .ne. 0) then
812 call update_error(ddx_error, "`psi` deallocation failed!")
813 return
814 endif
815 end if
816 if (allocated(state % phi_grid)) then
817 deallocate(state % phi_grid, stat=istatus)
818 if (istatus .ne. 0) then
819 call update_error(ddx_error, "`phi_grid` deallocation failed!")
820 return
821 endif
822 end if
823 if (allocated(state % phi)) then
824 deallocate(state % phi, stat=istatus)
825 if (istatus .ne. 0) then
826 call update_error(ddx_error, "`phi` deallocation failed!")
827 endif
828 end if
829 if (allocated(state % phiinf)) then
830 deallocate(state % phiinf, stat=istatus)
831 if (istatus .ne. 0) then
832 call update_error(ddx_error, "`phiinf` deallocation failed!")
833 endif
834 end if
835 if (allocated(state % phieps)) then
836 deallocate(state % phieps, stat=istatus)
837 if (istatus .ne. 0) then
838 call update_error(ddx_error, "`phieps` deallocation failed!")
839 endif
840 end if
841 if (allocated(state % phieps_rel_diff)) then
842 deallocate(state % phieps_rel_diff, stat=istatus)
843 if (istatus .ne. 0) then
844 call update_error(ddx_error, "`phieps_rel_diff` deallocation failed!")
845 endif
846 end if
847 if (allocated(state % xs)) then
848 deallocate(state % xs, stat=istatus)
849 if (istatus .ne. 0) then
850 call update_error(ddx_error, "`xs` deallocation failed!")
851 endif
852 end if
853 if (allocated(state % xs_rel_diff)) then
854 deallocate(state % xs_rel_diff, stat=istatus)
855 if (istatus .ne. 0) then
856 call update_error(ddx_error, "`xs_rel_diff` deallocation failed!")
857 endif
858 end if
859 if (allocated(state % s)) then
860 deallocate(state % s, stat=istatus)
861 if (istatus .ne. 0) then
862 call update_error(ddx_error, "`s` deallocation failed!")
863 endif
864 end if
865 if (allocated(state % s_rel_diff)) then
866 deallocate(state % s_rel_diff, stat=istatus)
867 if (istatus .ne. 0) then
868 call update_error(ddx_error, "`s_rel_diff` deallocation failed!")
869 endif
870 end if
871 if (allocated(state % sgrid)) then
872 deallocate(state % sgrid, stat=istatus)
873 if (istatus .ne. 0) then
874 call update_error(ddx_error, "`sgrid` deallocation failed!")
875 endif
876 end if
877 if (allocated(state % y)) then
878 deallocate(state % y, stat=istatus)
879 if (istatus .ne. 0) then
880 call update_error(ddx_error, "`y` deallocation failed!")
881 endif
882 end if
883 if (allocated(state % y_rel_diff)) then
884 deallocate(state % y_rel_diff, stat=istatus)
885 if (istatus .ne. 0) then
886 call update_error(ddx_error, "`y_rel_diff` deallocation failed!")
887 endif
888 end if
889 if (allocated(state % ygrid)) then
890 deallocate(state % ygrid, stat=istatus)
891 if (istatus .ne. 0) then
892 call update_error(ddx_error, "`ygrid` deallocation failed!")
893 endif
894 end if
895 if (allocated(state % g)) then
896 deallocate(state % g, stat=istatus)
897 if (istatus .ne. 0) then
898 call update_error(ddx_error, "`g` deallocation failed!")
899 endif
900 end if
901 if (allocated(state % q)) then
902 deallocate(state % q, stat=istatus)
903 if (istatus .ne. 0) then
904 call update_error(ddx_error, "`q` deallocation failed!")
905 endif
906 end if
907 if (allocated(state % qgrid)) then
908 deallocate(state % qgrid, stat=istatus)
909 if (istatus .ne. 0) then
910 call update_error(ddx_error, "`qgrid` deallocation failed!")
911 endif
912 end if
913 if (allocated(state % zeta)) then
914 deallocate(state % zeta, stat=istatus)
915 if (istatus .ne. 0) then
916 call update_error(ddx_error, "`zeta` deallocation failed!")
917 endif
918 end if
919 if (allocated(state % x_lpb)) then
920 deallocate(state % x_lpb, stat=istatus)
921 if (istatus .ne. 0) then
922 call update_error(ddx_error, "`x_lpb` deallocation failed!")
923 endif
924 end if
925 if (allocated(state % x_lpb_rel_diff)) then
926 deallocate(state % x_lpb_rel_diff, stat=istatus)
927 if (istatus .ne. 0) then
928 call update_error(ddx_error, "`x_lpb_rel_diff` deallocation failed!")
929 end if
930 end if
931 if (allocated(state % x_adj_lpb)) then
932 deallocate(state % x_adj_lpb, stat=istatus)
933 if (istatus .ne. 0) then
934 call update_error(ddx_error, "x_adj_lpb deallocation failed!")
935 endif
936 end if
937 if (allocated(state % x_adj_lpb_rel_diff)) then
938 deallocate(state % x_adj_lpb_rel_diff, stat=istatus)
939 if (istatus .ne. 0) then
940 call update_error(ddx_error, "`x_adj_lpb_rel_diff` deallocation failed!")
941 end if
942 end if
943 if (allocated(state % g_lpb)) then
944 deallocate(state % g_lpb, stat=istatus)
945 if (istatus .ne. 0) then
946 call update_error(ddx_error, "`g_lpb` deallocation failed!")
947 endif
948 end if
949 if (allocated(state % f_lpb)) then
950 deallocate(state % f_lpb, stat=istatus)
951 if (istatus .ne. 0) then
952 call update_error(ddx_error, "`f_lpb` deallocation failed!")
953 endif
954 end if
955 if (allocated(state % rhs_lpb)) then
956 deallocate(state % rhs_lpb, stat=istatus)
957 if (istatus .ne. 0) then
958 call update_error(ddx_error, "`rhs_lpb` deallocation failed!")
959 endif
960 end if
961 if (allocated(state % rhs_adj_lpb)) then
962 deallocate(state % rhs_adj_lpb, stat=istatus)
963 if (istatus .ne. 0) then
964 call update_error(ddx_error, "`rhs_adj_lpb` deallocation failed!")
965 endif
966 end if
967 if (allocated(state % zeta_dip)) then
968 deallocate(state % zeta_dip, stat=istatus)
969 if (istatus .ne. 0) then
970 call update_error(ddx_error, "`zeta_dip` deallocation failed!")
971 endif
972 end if
973 if (allocated(state % x_adj_re_grid)) then
974 deallocate(state % x_adj_re_grid, stat=istatus)
975 if (istatus .ne. 0) then
976 call update_error(ddx_error, "`x_adj_re_grid` deallocation failed!")
977 endif
978 end if
979 if (allocated(state % x_adj_r_grid)) then
980 deallocate(state % x_adj_r_grid, stat=istatus)
981 if (istatus .ne. 0) then
982 call update_error(ddx_error, "`x_adj_r_grid` deallocation failed!")
983 endif
984 end if
985 if (allocated(state % x_adj_e_grid)) then
986 deallocate(state % x_adj_e_grid, stat=istatus)
987 if (istatus .ne. 0) then
988 call update_error(ddx_error, "`x_adj_e_grid` deallocation failed!")
989 endif
990 end if
991 if (allocated(state % phi_n)) then
992 deallocate(state % phi_n, stat=istatus)
993 if (istatus .ne. 0) then
994 call update_error(ddx_error, "`phi_n` deallocation failed!")
995 endif
996 end if
997end subroutine deallocate_state
998
999!------------------------------------------------------------------------------
1011subroutine print_spherical(iunit, label, nbasis, lmax, ncol, icol, x)
1012 implicit none
1013 ! Inputs
1014 character (len=*), intent(in) :: label
1015 integer, intent(in) :: nbasis, lmax, ncol, icol, iunit
1016 real(dp), intent(in) :: x(nbasis, ncol)
1017 ! Local variables
1018 integer :: l, m, ind, noff, nprt, ic, j
1019 ! Print header:
1020 if (ncol .eq. 1) then
1021 write (iunit,'(3x,a,1x,a,i4,a)') label, "(column ", icol, ")"
1022 else
1023 write (iunit,'(3x,a)') label
1024 endif
1025 ! Print entries:
1026 if (ncol .eq. 1) then
1027 do l = 0, lmax
1028 ind = l*l + l + 1
1029 do m = -l, l
1030 write(iunit,1000) l, m, x(ind+m, 1)
1031 end do
1032 end do
1033 else
1034 noff = mod(ncol, 5)
1035 nprt = max(ncol-noff, 0)
1036 do ic = 1, nprt, 5
1037 write(iunit,1010) (j, j = ic, ic+4)
1038 do l = 0, lmax
1039 ind = l*l + l + 1
1040 do m = -l, l
1041 write(iunit,1020) l, m, x(ind+m, ic:ic+4)
1042 end do
1043 end do
1044 end do
1045 write (iunit,1010) (j, j = nprt+1, nprt+noff)
1046 do l = 0, lmax
1047 ind = l*l + l + 1
1048 do m = -l, l
1049 write(iunit,1020) l, m, x(ind+m, nprt+1:nprt+noff)
1050 end do
1051 end do
1052 end if
1053 1000 format(1x,i3,i4,f14.8)
1054 1010 format(8x,5i14)
1055 1020 format(1x,i3,i4,5f14.8)
1056end subroutine print_spherical
1057
1058!------------------------------------------------------------------------------
1068subroutine print_nodes(iunit, label, ngrid, ncol, icol, x)
1069 implicit none
1070 ! Inputs
1071 character (len=*), intent(in) :: label
1072 integer, intent(in) :: ngrid, ncol, icol, iunit
1073 real(dp), intent(in) :: x(ngrid, ncol)
1074 ! Local variables
1075 integer :: ig, noff, nprt, ic, j
1076 ! Print header :
1077 if (ncol .eq. 1) then
1078 write (iunit,'(3x,a,1x,a,i4,a)') label, "(column ", icol, ")"
1079 else
1080 write (iunit,'(3x,a)') label
1081 endif
1082 ! Print entries :
1083 if (ncol .eq. 1) then
1084 do ig = 1, ngrid
1085 write(iunit,1000) ig, x(ig, 1)
1086 enddo
1087 else
1088 noff = mod(ncol, 5)
1089 nprt = max(ncol-noff, 0)
1090 do ic = 1, nprt, 5
1091 write(iunit,1010) (j, j = ic, ic+4)
1092 do ig = 1, ngrid
1093 write(iunit,1020) ig, x(ig, ic:ic+4)
1094 end do
1095 end do
1096 write (iunit,1010) (j, j = nprt+1, nprt+noff)
1097 do ig = 1, ngrid
1098 write(iunit,1020) ig, x(ig, nprt+1:nprt+noff)
1099 end do
1100 end if
1101 !
1102 1000 format(1x,i5,f14.8)
1103 1010 format(6x,5i14)
1104 1020 format(1x,i5,5f14.8)
1105 !
1106end subroutine print_nodes
1107
1108!------------------------------------------------------------------------------
1113!------------------------------------------------------------------------------
1114subroutine print_ddvector(ddx_data, label, vector)
1115 implicit none
1116 ! Inputs
1117 type(ddx_type), intent(in) :: ddx_data
1118 character(len=*) :: label
1119 real(dp) :: vector(ddx_data % constants % nbasis, ddx_data % params % nsph)
1120 ! Local variables
1121 integer :: isph, lm
1122 write(6,*) label
1123 do isph = 1, ddx_data % params % nsph
1124 do lm = 1, ddx_data % constants % nbasis
1125 write(6,'(F15.8)') vector(lm,isph)
1126 end do
1127 end do
1128 return
1129end subroutine print_ddvector
1130
1131
1132!------------------------------------------------------------------------------
1147subroutine ddintegrate(nsph, nbasis, ngrid, vwgrid, ldvwgrid, x_grid, x_lm)
1148 implicit none
1149 !! Inputs
1150 integer, intent(in) :: nsph, nbasis, ngrid, ldvwgrid
1151 real(dp), intent(in) :: vwgrid(ldvwgrid, ngrid)
1152 real(dp), intent(in) :: x_grid(ngrid, nsph)
1153 !! Output
1154 real(dp), intent(out) :: x_lm(nbasis, nsph)
1155 !! Just call a single dgemm to do the job
1156 call dgemm('N', 'N', nbasis, nsph, ngrid, one, vwgrid, ldvwgrid, x_grid, &
1157 & ngrid, zero, x_lm, nbasis)
1158end subroutine ddintegrate
1159
1160!------------------------------------------------------------------------------
1174!------------------------------------------------------------------------------
1175subroutine dbasis(params, constants, x, basloc, dbsloc, vplm, vcos, vsin)
1176 implicit none
1177 type(ddx_params_type), intent(in) :: params
1178 type(ddx_constants_type), intent(in) :: constants
1179 real(dp), dimension(3), intent(in) :: x
1180 real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm
1181 real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc
1182 real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin
1183 integer :: l, m, ind
1184 real(dp) :: cthe, sthe, cphi, sphi, plm, fln, pp1, pm1, pp, VC, VS
1185 real(dp) :: et(3), ep(3)
1186
1187 ! get cos(\theta), sin(\theta), cos(\phi) and sin(\phi) from the cartesian
1188 ! coordinates of x.
1189 cthe = x(3)
1190 sthe = sqrt(one - cthe*cthe)
1191 ! not ( NORTH or SOUTH pole )
1192 if ( sthe.ne.zero ) then
1193 cphi = x(1)/sthe
1194 sphi = x(2)/sthe
1195 ! NORTH or SOUTH pole
1196 else
1197 cphi = one
1198 sphi = zero
1199 end if
1200 ! evaluate the derivatives of theta and phi:
1201 et(1) = cthe*cphi
1202 et(2) = cthe*sphi
1203 et(3) = -sthe
1204 ! not ( NORTH or SOUTH pole )
1205 if( sthe.ne.zero ) then
1206 ep(1) = -sphi/sthe
1207 ep(2) = cphi/sthe
1208 ep(3) = zero
1209 ! NORTH or SOUTH pole
1210 else
1211 ep(1) = zero
1212 ep(2) = one
1213 ep(3) = zero
1214 end if
1215 vc=zero
1216 vs=cthe
1217 ! evaluate the generalized legendre polynomials. Temporary workspace
1218 ! is of size (p+1) here, so we use vcos for that purpose
1219 call polleg_work( cthe, sthe, params % lmax, vplm, vcos )
1220 !
1221 ! evaluate cos(m*phi) and sin(m*phi) arrays. notice that this is
1222 ! pointless if z = 1, as the only non vanishing terms will be the
1223 ! ones with m=0.
1224 !
1225 ! not ( NORTH or SOUTH pole )
1226 if ( sthe.ne.zero ) then
1227 call trgev( cphi, sphi, params % lmax, vcos, vsin )
1228 ! NORTH or SOUTH pole
1229 else
1230 vcos = one
1231 vsin = zero
1232 end if
1233 ! now build the spherical harmonics. we will distinguish m=0,
1234 ! m>0 and m<0:
1235 !
1236 basloc = zero
1237 dbsloc = zero
1238 do l = 0, params % lmax
1239 ind = l*l + l + 1
1240 ! m = 0
1241 fln = constants % vscales(ind)
1242 basloc(ind) = fln*vplm(ind)
1243 if (l.gt.0) then
1244 dbsloc(:,ind) = fln*vplm(ind+1)*et(:)
1245 else
1246 dbsloc(:,ind) = zero
1247 end if
1248 do m = 1, l
1249 fln = constants % vscales(ind+m)
1250 plm = fln*vplm(ind+m)
1251 pp1 = zero
1252 if (m.lt.l) pp1 = -pt5*vplm(ind+m+1)
1253 pm1 = pt5*(dble(l+m)*dble(l-m+1)*vplm(ind+m-1))
1254 pp = pp1 + pm1
1255 !
1256 ! m > 0
1257 ! -----
1258 !
1259 basloc(ind+m) = plm*vcos(m+1)
1260 !
1261 ! not ( NORTH or SOUTH pole )
1262 if ( sthe.ne.zero ) then
1263 !
1264 dbsloc(:,ind+m) = -fln*pp*vcos(m+1)*et(:) - dble(m)*plm*vsin(m+1)*ep(:)
1265 !
1266 !
1267 ! NORTH or SOUTH pole
1268 else
1269 !
1270 dbsloc(:,ind+m) = -fln*pp*vcos(m+1)*et(:) - fln*pp*ep(:)*vc
1271 !
1272 !
1273 endif
1274 !
1275 ! m < 0
1276 ! -----
1277 !
1278 basloc(ind-m) = plm*vsin(m+1)
1279 !
1280 ! not ( NORTH or SOUTH pole )
1281 if ( sthe.ne.zero ) then
1282 !
1283 dbsloc(:,ind-m) = -fln*pp*vsin(m+1)*et(:) + dble(m)*plm*vcos(m+1)*ep(:)
1284 !
1285 ! NORTH or SOUTH pole
1286 else
1287 !
1288 dbsloc(:,ind-m) = -fln*pp*vsin(m+1)*et(:) - fln*pp*ep(:)*vs
1289 !
1290 endif
1291 !
1292 enddo
1293 enddo
1294end subroutine dbasis
1295
1296! Purpose : compute
1297!
1298! l
1299! sum 4pi/(2l+1) t * Y_l^m( s ) * sigma_l^m
1300! l,m
1301!
1302! which is need to compute action of COSMO matrix L.
1303!------------------------------------------------------------------------------------------------
1305!------------------------------------------------------------------------------
1306real(dp) function intmlp(params, constants, t, sigma, basloc )
1307 implicit none
1308!
1309 type(ddx_params_type), intent(in) :: params
1310 type(ddx_constants_type), intent(in) :: constants
1311 real(dp), intent(in) :: t
1312 real(dp), dimension(constants % nbasis), intent(in) :: sigma, basloc
1313!
1314 integer :: l, ind
1315 real(dp) :: tt, ss, fac
1316!
1317!------------------------------------------------------------------------------------------------
1318!
1319! initialize t^l
1320 tt = one
1321!
1322! initialize
1323 ss = zero
1324!
1325! loop over l
1326 do l = 0, params % lmax
1327!
1328 ind = l*l + l + 1
1329!
1330! update factor 4pi / (2l+1) * t^l
1331 fac = tt / constants % vscales(ind)**2
1332!
1333! contract over l,m and accumulate
1334 ss = ss + fac * dot_product( basloc(ind-l:ind+l), &
1335 sigma( ind-l:ind+l) )
1336!
1337! update t^l
1338 tt = tt*t
1339!
1340 enddo
1341!
1342! redirect
1343 intmlp = ss
1344!
1345!
1346end function intmlp
1347
1348!------------------------------------------------------------------------------
1351!------------------------------------------------------------------------------
1352subroutine wghpot(ncav, phi_cav, nsph, ngrid, ui, phi_grid, g)
1353 implicit none
1354 !! Inputs
1355 integer, intent(in) :: ncav, nsph, ngrid
1356 real(dp), intent(in) :: phi_cav(ncav), ui(ngrid, nsph)
1357 !! Outputs
1358 real(dp), intent(out) :: g(ngrid, nsph), phi_grid(ngrid, nsph)
1359 !! Local variables
1360 integer isph, igrid, icav
1361 !! Code
1362 ! Initialize
1363 icav = 0
1364 g = zero
1365 phi_grid = zero
1366 ! Loop over spheres
1367 do isph = 1, nsph
1368 ! Loop over points
1369 do igrid = 1, ngrid
1370 ! Non-zero contribution from point
1371 if (ui(igrid, isph) .ne. zero) then
1372 ! Advance cavity point counter
1373 icav = icav + 1
1374 phi_grid(igrid, isph) = phi_cav(icav)
1375 ! Weigh by (negative) characteristic function
1376 g(igrid, isph) = -ui(igrid, isph) * phi_cav(icav)
1377 endif
1378 enddo
1379 enddo
1380end subroutine wghpot
1381
1383subroutine hsnorm(lmax, nbasis, u, unorm)
1384 implicit none
1385 integer, intent(in) :: lmax, nbasis
1386 real(dp), dimension(nbasis), intent(in) :: u
1387 real(dp), intent(inout) :: unorm
1388 integer :: l, m, ind
1389 real(dp) :: fac
1390 ! initialize
1391 unorm = zero
1392 do l = 0, lmax
1393 ind = l*l + l + 1
1394 fac = one/(one + dble(l))
1395 do m = -l, l
1396 unorm = unorm + fac*u(ind+m)*u(ind+m)
1397 end do
1398 end do
1399 ! the much neglected square root
1400 unorm = sqrt(unorm)
1401end subroutine hsnorm
1402
1404real(dp) function hnorm(lmax, nbasis, nsph, x)
1405 implicit none
1406 integer, intent(in) :: lmax, nbasis, nsph
1407 real(dp), dimension(nbasis, nsph), intent(in) :: x
1408 integer :: isph
1409 real(dp) :: vrms, fac
1410
1411 vrms = 0.0_dp
1412 !$omp parallel do default(none) shared(nsph,lmax,nbasis,x) &
1413 !$omp private(isph,fac) schedule(dynamic) reduction(+:vrms)
1414 do isph = 1, nsph
1415 call hsnorm(lmax, nbasis, x(:,isph), fac)
1416 vrms = vrms + fac*fac
1417 enddo
1418 hnorm = sqrt(vrms/dble(nsph))
1419end function hnorm
1420
1421real(dp) function rmsnorm(lmax, nbasis, nsph, x)
1422 implicit none
1423 ! TODO: nbasis and lmax are redundant, remove one
1424 integer, intent(in) :: lmax, nbasis, nsph
1425 real(dp), dimension(nbasis, nsph), intent(in) :: x
1426 integer :: n
1427 real(dp) :: vrms, vmax
1428 ! lmax is not actually used, it is here to comply to the interface
1429 if (lmax.eq.0) continue
1430 n = nbasis*nsph
1431 call rmsvec(n,x,vrms,vmax)
1432 rmsnorm = vrms
1433end function rmsnorm
1434
1436subroutine rmsvec(n, v, vrms, vmax)
1437 implicit none
1438 integer, intent(in) :: n
1439 real(dp), dimension(n), intent(in) :: v
1440 real(dp), intent(inout) :: vrms, vmax
1441 integer :: i
1442
1443 vrms = zero
1444 vmax = zero
1445 do i = 1,n
1446 vmax = max(vmax,abs(v(i)))
1447 vrms = vrms + v(i)*v(i)
1448 end do
1449 ! the much neglected square root
1450 vrms = sqrt(vrms/dble(n))
1451endsubroutine rmsvec
1452
1453subroutine adjrhs(params, constants, isph, xi, vlm, work)
1454 implicit none
1455 type(ddx_params_type), intent(in) :: params
1456 type(ddx_constants_type), intent(in) :: constants
1457 integer, intent(in) :: isph
1458 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi
1459 real(dp), dimension(constants % nbasis), intent(inout) :: vlm
1460 real(dp), dimension(params % lmax+1), intent(inout) :: work
1461
1462 integer :: ij, jsph, ig
1463 real(dp) :: vji(3), vvji, tji, xji, oji, fac
1464
1465 do ij = constants % inl(isph), constants % inl(isph+1)-1
1466 jsph = constants % nl(ij)
1467 do ig = 1, params % ngrid
1468 vji = params % csph(:,jsph) + params % rsph(jsph)* &
1469 & constants % cgrid(:,ig) - params % csph(:,isph)
1470 vvji = sqrt(dot_product(vji,vji))
1471 tji = vvji/params % rsph(isph)
1472 if ( tji.lt.( one + (params % se+one)/two*params % eta ) ) then
1473 xji = fsw( tji, params % se, params % eta )
1474 if ( constants % fi(ig,jsph).gt.one ) then
1475 oji = xji/ constants % fi(ig,jsph)
1476 else
1477 oji = xji
1478 endif
1479 fac = constants % wgrid(ig) * xi(ig,jsph) * oji
1480 call fmm_l2p_adj_work(vji, fac, params % rsph(isph), &
1481 & params % lmax, constants % vscales_rel, one, vlm, work)
1482 endif
1483 enddo
1484 enddo
1485end subroutine adjrhs
1486
1487subroutine calcv(params, constants, isph, pot, sigma, work)
1488 implicit none
1489 type(ddx_params_type), intent(in) :: params
1490 type(ddx_constants_type), intent(in) :: constants
1491 integer, intent(in) :: isph
1492 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sigma
1493 real(dp), dimension(params % ngrid), intent(inout) :: pot
1494 real(dp), dimension(params % lmax+1), intent(inout) :: work
1495
1496 integer :: its, ij, jsph
1497 real(dp) :: vij(3)
1498 real(dp) :: vvij, tij, xij, oij, thigh
1499
1500 thigh = one + pt5*(params % se + one)*params % eta
1501 pot(:) = zero
1502 ! loop over grid points
1503 do its = 1, params % ngrid
1504 ! contribution from integration point present
1505 if (constants % ui(its,isph).lt.one) then
1506 ! loop over neighbors of i-sphere
1507 do ij = constants % inl(isph), constants % inl(isph+1)-1
1508 jsph = constants % nl(ij)
1509 ! compute t_n^ij = | r_i + \rho_i s_n - r_j | / \rho_j
1510 vij = params % csph(:,isph) + params % rsph(isph)* &
1511 & constants % cgrid(:,its) - params % csph(:,jsph)
1512 vvij = sqrt( dot_product( vij, vij ) )
1513 tij = vvij / params % rsph(jsph)
1514 ! point is INSIDE j-sphere
1515 if (tij.lt.thigh) then
1516 xij = fsw(tij, params % se, params % eta)
1517 if (constants % fi(its,isph).gt.one) then
1518 oij = xij / constants % fi(its,isph)
1519 else
1520 oij = xij
1521 end if
1522 call fmm_l2p_work(vij, params % rsph(jsph), params % lmax, &
1523 & constants % vscales_rel, oij, sigma(:, jsph), one, &
1524 & pot(its), work)
1525 end if
1526 end do
1527 end if
1528 end do
1529end subroutine calcv
1530
1531!------------------------------------------------------------------------------
1533!------------------------------------------------------------------------------
1534subroutine ddeval_grid(params, constants, alpha, x_sph, beta, x_grid)
1535 implicit none
1536 !! Inputs
1537 type(ddx_params_type), intent(in) :: params
1538 type(ddx_constants_type), intent(in) :: constants
1539 real(dp), intent(in) :: alpha, x_sph(constants % nbasis, params % nsph), &
1540 & beta
1541 !! Output
1542 real(dp), intent(inout) :: x_grid(params % ngrid, params % nsph)
1543 !! The code
1544 call ddeval_grid_work(constants % nbasis, params % ngrid, params % nsph, &
1545 & constants % vgrid, constants % vgrid_nbasis, alpha, x_sph, beta, &
1546 & x_grid)
1547end subroutine ddeval_grid
1548
1549!------------------------------------------------------------------------------
1551!------------------------------------------------------------------------------
1552subroutine ddeval_grid_work(nbasis, ngrid, nsph, vgrid, ldvgrid, alpha, &
1553 & x_sph, beta, x_grid)
1554 implicit none
1555 !! Inputs
1556 integer, intent(in) :: nbasis, ngrid, nsph, ldvgrid
1557 real(dp), intent(in) :: vgrid(ldvgrid, ngrid), alpha, &
1558 & x_sph(nbasis, nsph), beta
1559 !! Output
1560 real(dp), intent(inout) :: x_grid(ngrid, nsph)
1561 !! Local variables
1562 external :: dgemm
1563 !! The code
1564 call dgemm('T', 'N', ngrid, nsph, nbasis, alpha, vgrid, ldvgrid, x_sph, &
1565 & nbasis, beta, x_grid, ngrid)
1566end subroutine ddeval_grid_work
1567
1569subroutine ddintegrate_sph(params, constants, alpha, x_grid, beta, x_sph)
1570 implicit none
1571 !! Inputs
1572 type(ddx_params_type), intent(in) :: params
1573 type(ddx_constants_type), intent(in) :: constants
1574 real(dp), intent(in) :: alpha, x_grid(params % ngrid, params % nsph), &
1575 & beta
1576 !! Output
1577 real(dp), intent(inout) :: x_sph(constants % nbasis, params % nsph)
1578 !! The code
1579 call ddintegrate_sph_work(constants % nbasis, params % ngrid, &
1580 & params % nsph, constants % vwgrid, constants % vgrid_nbasis, alpha, &
1581 & x_grid, beta, x_sph)
1582end subroutine ddintegrate_sph
1583
1585subroutine ddintegrate_sph_work(nbasis, ngrid, nsph, vwgrid, ldvwgrid, alpha, &
1586 & x_grid, beta, x_sph)
1587 implicit none
1588 !! Inputs
1589 integer, intent(in) :: nbasis, ngrid, nsph, ldvwgrid
1590 real(dp), intent(in) :: vwgrid(ldvwgrid, ngrid), alpha, &
1591 & x_grid(ngrid, nsph), beta
1592 !! Outputs
1593 real(dp), intent(inout) :: x_sph(nbasis, nsph)
1594 !! Local variables
1595 !! The code
1596 call dgemm('N', 'N', nbasis, nsph, ngrid, alpha, vwgrid, ldvwgrid, &
1597 & x_grid, ngrid, beta, x_sph, nbasis)
1598end subroutine ddintegrate_sph_work
1599
1600!------------------------------------------------------------------------------
1602!------------------------------------------------------------------------------
1603subroutine ddcav_to_grid(params, constants, x_cav, x_grid)
1604 implicit none
1605 !! Inputs
1606 type(ddx_params_type), intent(in) :: params
1607 type(ddx_constants_type), intent(in) :: constants
1608 real(dp), intent(in) :: x_cav(constants % ncav)
1609 !! Output
1610 real(dp), intent(inout) :: x_grid(params % ngrid, params % nsph)
1611 !! The code
1612 call ddcav_to_grid_work(params % ngrid, params % nsph, constants % ncav, &
1613 & constants % icav_ia, constants % icav_ja, x_cav, x_grid)
1614end subroutine ddcav_to_grid
1615
1616!------------------------------------------------------------------------------
1618!------------------------------------------------------------------------------
1619subroutine ddcav_to_grid_work(ngrid, nsph, ncav, icav_ia, icav_ja, x_cav, &
1620 & x_grid)
1621 implicit none
1622 !! Inputs
1623 integer, intent(in) :: ngrid, nsph, ncav
1624 integer, intent(in) :: icav_ia(nsph+1), icav_ja(ncav)
1625 real(dp), intent(in) :: x_cav(ncav)
1626 !! Output
1627 real(dp), intent(out) :: x_grid(ngrid, nsph)
1628 !! Local variables
1629 integer :: isph, icav, igrid, igrid_old
1630 !! The code
1631 do isph = 1, nsph
1632 igrid_old = 0
1633 igrid = 0
1634 do icav = icav_ia(isph), icav_ia(isph+1)-1
1635 igrid = icav_ja(icav)
1636 x_grid(igrid_old+1:igrid-1, isph) = zero
1637 x_grid(igrid, isph) = x_cav(icav)
1638 igrid_old = igrid
1639 end do
1640 x_grid(igrid+1:ngrid, isph) = zero
1641 end do
1642end subroutine ddcav_to_grid_work
1643
1648subroutine ddproject_cav(params, constants, s, xi)
1649 implicit none
1650 type(ddx_params_type), intent(in) :: params
1651 type(ddx_constants_type), intent(in) :: constants
1652 real(dp), intent(in) :: s(constants%nbasis, params%nsph)
1653 real(dp), intent(out) :: xi(constants%ncav)
1654 integer :: its, isph, ii
1655 ii = 0
1656 do isph = 1, params%nsph
1657 do its = 1, params%ngrid
1658 if (constants%ui(its, isph) .gt. zero) then
1659 ii = ii + 1
1660 xi(ii) = constants%ui(its, isph)*dot_product( &
1661 & constants%vwgrid(:constants % nbasis, its), s(:, isph))
1662 end if
1663 end do
1664 end do
1665end subroutine ddproject_cav
1666
1667!------------------------------------------------------------------------------
1669!------------------------------------------------------------------------------
1670subroutine tree_m2m_rotation(params, constants, node_m)
1671 implicit none
1672 ! Inputs
1673 type(ddx_params_type), intent(in) :: params
1674 type(ddx_constants_type), intent(in) :: constants
1675 ! Output
1676 real(dp), intent(inout) :: node_m((params % pm+1)**2, constants % nclusters)
1677 ! Temporary workspace
1678 real(dp) :: work(6*params % pm**2 + 19*params % pm + 8)
1679 ! Call corresponding work routine
1680 call tree_m2m_rotation_work(params, constants, node_m, work)
1681end subroutine tree_m2m_rotation
1682
1683!------------------------------------------------------------------------------
1685!------------------------------------------------------------------------------
1686subroutine tree_m2m_rotation_work(params, constants, node_m, work)
1687 implicit none
1688 ! Inputs
1689 type(ddx_params_type), intent(in) :: params
1690 type(ddx_constants_type), intent(in) :: constants
1691 ! Output
1692 real(dp), intent(inout) :: node_m((params % pm+1)**2, constants % nclusters)
1693 ! Temporary workspace
1694 real(dp), intent(out) :: work(6*params % pm**2 + 19*params % pm + 8)
1695 ! Local variables
1696 integer :: i, j
1697 real(dp) :: c1(3), c(3), r1, r
1698 ! Bottom-to-top pass
1699 !!$omp parallel do default(none) shared(constants,params,node_m) &
1700 !!$omp private(i,j,c1,c,r1,r,work)
1701 do i = constants % nclusters, 1, -1
1702 ! Leaf node does not need any update
1703 if (constants % children(1, i) == 0) cycle
1704 c = constants % cnode(:, i)
1705 r = constants % rnode(i)
1706 ! First child initializes output
1707 j = constants % children(1, i)
1708 c1 = constants % cnode(:, j)
1709 r1 = constants % rnode(j)
1710 c1 = c1 - c
1711 call fmm_m2m_rotation_work(c1, r1, r, &
1712 & params % pm, &
1713 & constants % vscales, &
1714 & constants % vcnk, one, &
1715 & node_m(:, j), zero, node_m(:, i), work)
1716 ! All other children update the same output
1717 do j = constants % children(1, i)+1, constants % children(2, i)
1718 c1 = constants % cnode(:, j)
1719 r1 = constants % rnode(j)
1720 c1 = c1 - c
1721 call fmm_m2m_rotation_work(c1, r1, r, params % pm, &
1722 & constants % vscales, constants % vcnk, one, &
1723 & node_m(:, j), one, node_m(:, i), work)
1724 end do
1725 end do
1726end subroutine tree_m2m_rotation_work
1727
1728!------------------------------------------------------------------------------
1730!------------------------------------------------------------------------------
1731subroutine tree_m2m_bessel_rotation(params, constants, node_m)
1732 implicit none
1733 ! Inputs
1734 type(ddx_params_type), intent(in) :: params
1735 type(ddx_constants_type), intent(in) :: constants
1736 ! Output
1737 real(dp), intent(inout) :: node_m((params % pm+1)**2, constants % nclusters)
1738 ! Temporary workspace
1739 !real(dp) :: work(6*params % pm**2 + 19*params % pm + 8)
1740 ! Call corresponding work routine
1741 call tree_m2m_bessel_rotation_work(params, constants, node_m)
1742end subroutine tree_m2m_bessel_rotation
1743
1744!------------------------------------------------------------------------------
1746!------------------------------------------------------------------------------
1747subroutine tree_m2m_bessel_rotation_work(params, constants, node_m)
1748 use complex_bessel
1749 implicit none
1750 ! Inputs
1751 type(ddx_params_type), intent(in) :: params
1752 type(ddx_constants_type), intent(in) :: constants
1753 ! Output
1754 real(dp), intent(inout) :: node_m((params % pm+1)**2, constants % nclusters)
1755 ! Temporary workspace
1756 real(dp) :: work(6*params % pm**2 + 19*params % pm + 8)
1757 complex(dp) :: work_complex(2*params % pm+1)
1758 ! Local variables
1759 integer :: i, j
1760 real(dp) :: c1(3), c(3), r1, r
1761 ! Bottom-to-top pass
1762 do i = constants % nclusters, 1, -1
1763 ! Leaf node does not need any update
1764 if (constants % children(1, i) == 0) cycle
1765 c = constants % cnode(:, i)
1766 r = constants % rnode(i)
1767 ! First child initializes output
1768 j = constants % children(1, i)
1769 c1 = constants % cnode(:, j)
1770 r1 = constants % rnode(j)
1771 c1 = params % kappa*(c1 - c)
1772 call fmm_m2m_bessel_rotation_work(c1, &
1773 & constants % SK_rnode(:, j), constants % SK_rnode(:, i), &
1774 & params % pm, constants % vscales, one, node_m(:, j), zero, &
1775 & node_m(:, i), work, work_complex)
1776 ! All other children update the same output
1777 do j = constants % children(1, i)+1, constants % children(2, i)
1778 c1 = constants % cnode(:, j)
1779 r1 = constants % rnode(j)
1780 c1 = params % kappa*(c1 - c)
1781 call fmm_m2m_bessel_rotation_work(c1, &
1782 & constants % SK_rnode(:, j), constants % SK_rnode(:, i), &
1783 & params % pm, constants % vscales, one, node_m(:, j), one, &
1784 & node_m(:, i), work, work_complex)
1785 end do
1786 end do
1787end subroutine tree_m2m_bessel_rotation_work
1788
1789!------------------------------------------------------------------------------
1791!------------------------------------------------------------------------------
1792subroutine tree_m2m_rotation_adj(params, constants, node_m)
1793 implicit none
1794 ! Inputs
1795 type(ddx_params_type), intent(in) :: params
1796 type(ddx_constants_type), intent(in) :: constants
1797 ! Output
1798 real(dp), intent(inout) :: node_m((params % pm+1)**2, constants % nclusters)
1799 ! Temporary workspace
1800 real(dp) :: work(6*params % pm**2 + 19*params % pm + 8)
1801 ! Call corresponding work routine
1802 call tree_m2m_rotation_adj_work(params, constants, node_m, work)
1803end subroutine tree_m2m_rotation_adj
1804
1805!------------------------------------------------------------------------------
1807!------------------------------------------------------------------------------
1808subroutine tree_m2m_rotation_adj_work(params, constants, node_m, work)
1809 implicit none
1810 ! Inputs
1811 type(ddx_params_type), intent(in) :: params
1812 type(ddx_constants_type), intent(in) :: constants
1813 ! Output
1814 real(dp), intent(inout) :: node_m((params % pm+1)**2, constants % nclusters)
1815 ! Temporary workspace
1816 real(dp), intent(out) :: work(6*params % pm**2 + 19*params % pm + 8)
1817 ! Local variables
1818 integer :: i, j
1819 real(dp) :: c1(3), c(3), r1, r
1820 ! Top-to-bottom pass
1821 do i = 2, constants % nclusters
1822 j = constants % parent(i)
1823 c = constants % cnode(:, j)
1824 r = constants % rnode(j)
1825 c1 = constants % cnode(:, i)
1826 r1 = constants % rnode(i)
1827 c1 = c - c1
1828 call fmm_m2m_rotation_adj_work(c1, r, r1, params % pm, &
1829 & constants % vscales, constants % vcnk, one, node_m(:, j), one, &
1830 & node_m(:, i), work)
1831 end do
1832end subroutine tree_m2m_rotation_adj_work
1833!------------------------------------------------------------------------------
1835!------------------------------------------------------------------------------
1836subroutine tree_m2m_bessel_rotation_adj(params, constants, node_m)
1837 implicit none
1838 ! Inputs
1839 type(ddx_params_type), intent(in) :: params
1840 type(ddx_constants_type), intent(in) :: constants
1841 ! Output
1842 real(dp), intent(inout) :: node_m((params % pm+1)**2, constants % nclusters)
1843 ! Temporary workspace
1844 !real(dp), intent(out) :: work(6*params % pm**2 + 19*params % pm + 8)
1845 call tree_m2m_bessel_rotation_adj_work(params, constants, node_m)
1846end subroutine tree_m2m_bessel_rotation_adj
1847
1848!------------------------------------------------------------------------------
1850!------------------------------------------------------------------------------
1851subroutine tree_m2m_bessel_rotation_adj_work(params, constants, node_m)
1852 implicit none
1853 ! Inputs
1854 type(ddx_params_type), intent(in) :: params
1855 type(ddx_constants_type), intent(in) :: constants
1856 ! Output
1857 real(dp), intent(inout) :: node_m((params % pm+1)**2, constants % nclusters)
1858 ! Temporary workspace
1859 real(dp) :: work(6*params % pm**2 + 19*params % pm + 8)
1860 complex(dp) :: work_complex(2*params % pm+1)
1861 ! Local variables
1862 integer :: i, j
1863 real(dp) :: c1(3), c(3)
1864 ! Top-to-bottom pass
1865 do i = 2, constants % nclusters
1866 j = constants % parent(i)
1867 c = constants % cnode(:, j)
1868 c1 = constants % cnode(:, i)
1869 c1 = params % kappa*(c1 - c)
1870 ! Little bit confusion about the i and j indices
1871 call fmm_m2m_bessel_rotation_adj_work(c1, constants % SK_rnode(:, i), &
1872 & constants % SK_rnode(:, j), params % pm, constants % vscales, &
1873 & one, node_m(:, j), one, node_m(:, i), work, work_complex)
1874 end do
1875end subroutine tree_m2m_bessel_rotation_adj_work
1876
1877!------------------------------------------------------------------------------
1879!------------------------------------------------------------------------------
1880subroutine tree_l2l_rotation(params, constants, node_l)
1881 implicit none
1882 ! Inputs
1883 type(ddx_params_type), intent(in) :: params
1884 type(ddx_constants_type), intent(in) :: constants
1885 ! Output
1886 real(dp), intent(inout) :: node_l((params % pl+1)**2, constants % nclusters)
1887 ! Temporary workspace
1888 real(dp) :: work(6*params % pl**2 + 19*params % pl + 8)
1889 ! Call corresponding work routine
1890 call tree_l2l_rotation_work(params, constants, node_l, work)
1891end subroutine tree_l2l_rotation
1892
1893!------------------------------------------------------------------------------
1895!------------------------------------------------------------------------------
1896subroutine tree_l2l_rotation_work(params, constants, node_l, work)
1897 implicit none
1898 ! Inputs
1899 type(ddx_params_type), intent(in) :: params
1900 type(ddx_constants_type), intent(in) :: constants
1901 ! Output
1902 real(dp), intent(inout) :: node_l((params % pl+1)**2, constants % nclusters)
1903 ! Temporary workspace
1904 real(dp), intent(out) :: work(6*params % pl**2 + 19*params % pl + 8)
1905 ! Local variables
1906 integer :: i, j
1907 real(dp) :: c1(3), c(3), r1, r
1908 ! Top-to-bottom pass
1909 !!$omp parallel do default(none) shared(constants,params,node_l) &
1910 !!$omp private(i,j,c1,c,r1,r,work)
1911 do i = 2, constants % nclusters
1912 j = constants % parent(i)
1913 c = constants % cnode(:, j)
1914 r = constants % rnode(j)
1915 c1 = constants % cnode(:, i)
1916 r1 = constants % rnode(i)
1917 c1 = c - c1
1918 call fmm_l2l_rotation_work(c1, r, r1, params % pl, &
1919 & constants % vscales, constants % vfact, one, &
1920 & node_l(:, j), one, node_l(:, i), work)
1921 end do
1922end subroutine tree_l2l_rotation_work
1923
1924!------------------------------------------------------------------------------
1926!------------------------------------------------------------------------------
1927subroutine tree_l2l_bessel_rotation(params, constants, node_l)
1928 implicit none
1929 ! Inputs
1930 type(ddx_params_type), intent(in) :: params
1931 type(ddx_constants_type), intent(in) :: constants
1932 ! Output
1933 real(dp), intent(inout) :: node_l((params % pl+1)**2, constants % nclusters)
1934 ! Temporary workspace
1935 !real(dp) :: work(6*params % pl**2 + 19*params % pl + 8)
1936 ! Call corresponding work routine
1937 call tree_l2l_bessel_rotation_work(params, constants, node_l)
1938end subroutine tree_l2l_bessel_rotation
1939
1940!------------------------------------------------------------------------------
1942!------------------------------------------------------------------------------
1943subroutine tree_l2l_bessel_rotation_work(params, constants, node_l)
1944 use complex_bessel
1945 implicit none
1946 ! Inputs
1947 type(ddx_params_type), intent(in) :: params
1948 type(ddx_constants_type), intent(in) :: constants
1949 ! Output
1950 real(dp), intent(inout) :: node_l((params % pl+1)**2, constants % nclusters)
1951 ! Temporary workspace
1952 real(dp) :: work(6*params % pl**2 + 19*params % pl + 8)
1953 complex(dp) :: work_complex(2*params % pl+1)
1954 ! Local variables
1955 integer :: i, j
1956 real(dp) :: c_child(3), c_parent(3), c_diff(3)
1957 ! Top-to-bottom pass
1958 do i = 2, constants % nclusters
1959 j = constants % parent(i)
1960 c_child = constants % cnode(:, j)
1961 c_parent = constants % cnode(:, i)
1962 c_diff = params % kappa*(c_child - c_parent)
1963 call fmm_l2l_bessel_rotation_work(c_diff, &
1964 & constants % si_rnode(:, j), constants % si_rnode(:, i), &
1965 & params % pl, constants % vscales, one, &
1966 & node_l(:, j), one, node_l(:, i), work, work_complex)
1967 end do
1968end subroutine tree_l2l_bessel_rotation_work
1969
1970!------------------------------------------------------------------------------
1972!------------------------------------------------------------------------------
1973subroutine tree_l2l_rotation_adj(params, constants, node_l)
1974 implicit none
1975 ! Inputs
1976 type(ddx_params_type), intent(in) :: params
1977 type(ddx_constants_type), intent(in) :: constants
1978 ! Output
1979 real(dp), intent(inout) :: node_l((params % pl+1)**2, constants % nclusters)
1980 ! Temporary workspace
1981 real(dp) :: work(6*params % pl**2 + 19*params % pl + 8)
1982 ! Call corresponding work routine
1983 call tree_l2l_rotation_adj_work(params, constants, node_l, work)
1984end subroutine tree_l2l_rotation_adj
1985
1986!------------------------------------------------------------------------------
1988!------------------------------------------------------------------------------
1989subroutine tree_l2l_rotation_adj_work(params, constants, node_l, work)
1990 implicit none
1991 ! Inputs
1992 type(ddx_params_type), intent(in) :: params
1993 type(ddx_constants_type), intent(in) :: constants
1994 ! Output
1995 real(dp), intent(inout) :: node_l((params % pl+1)**2, constants % nclusters)
1996 ! Temporary workspace
1997 real(dp), intent(out) :: work(6*params % pl**2 + 19*params % pl + 8)
1998 ! Local variables
1999 integer :: i, j
2000 real(dp) :: c1(3), c(3), r1, r
2001 ! Bottom-to-top pass
2002 do i = constants % nclusters, 1, -1
2003 ! Leaf node does not need any update
2004 if (constants % children(1, i) == 0) cycle
2005 c = constants % cnode(:, i)
2006 r = constants % rnode(i)
2007 ! First child initializes output
2008 j = constants % children(1, i)
2009 c1 = constants % cnode(:, j)
2010 r1 = constants % rnode(j)
2011 c1 = c1 - c
2012 call fmm_l2l_rotation_adj_work(c1, r1, r, params % pl, &
2013 & constants % vscales, constants % vfact, one, &
2014 & node_l(:, j), zero, node_l(:, i), work)
2015 ! All other children update the same output
2016 do j = constants % children(1, i)+1, constants % children(2, i)
2017 c1 = constants % cnode(:, j)
2018 r1 = constants % rnode(j)
2019 c1 = c1 - c
2020 call fmm_l2l_rotation_adj_work(c1, r1, r, params % pl, &
2021 & constants % vscales, constants % vfact, one, &
2022 & node_l(:, j), one, node_l(:, i), work)
2023 end do
2024 end do
2025end subroutine tree_l2l_rotation_adj_work
2026
2027!------------------------------------------------------------------------------
2029!------------------------------------------------------------------------------
2030subroutine tree_l2l_bessel_rotation_adj(params, constants, node_l)
2031 implicit none
2032 ! Inputs
2033 type(ddx_params_type), intent(in) :: params
2034 type(ddx_constants_type), intent(in) :: constants
2035 ! Output
2036 real(dp), intent(inout) :: node_l((params % pl+1)**2, constants % nclusters)
2037 ! Temporary workspace
2038 real(dp) :: work(6*params % pl**2 + 19*params % pl + 8)
2039 complex(dp) :: work_complex(2*params % pl+1)
2040 ! Local variables
2041 integer :: i, j
2042 real(dp) :: c_parent(3), c_child(3), c_diff(3)
2043 ! Bottom-to-top pass
2044 do i = constants % nclusters, 1, -1
2045 c_child = constants % cnode(:, i)
2046 j = constants % parent(i)
2047 if (j == 0) cycle
2048 c_parent = constants % cnode(:, j)
2049 c_diff = params % kappa*(c_parent - c_child)
2050 call fmm_l2l_bessel_rotation_adj_work(c_diff, &
2051 & constants % si_rnode(:, i), constants % si_rnode(:, j), &
2052 & params % pl, constants % vscales, one, &
2053 & node_l(:, i), one, node_l(:, j), work, work_complex)
2054 end do
2055end subroutine tree_l2l_bessel_rotation_adj
2056
2057!------------------------------------------------------------------------------
2059!------------------------------------------------------------------------------
2060subroutine tree_m2l_rotation(params, constants, node_m, node_l)
2061 implicit none
2062 ! Inputs
2063 type(ddx_params_type), intent(in) :: params
2064 type(ddx_constants_type), intent(in) :: constants
2065 real(dp), intent(in) :: node_m((params % pm+1)**2, constants % nclusters)
2066 ! Output
2067 real(dp), intent(out) :: node_l((params % pl+1)**2, constants % nclusters)
2068 ! Temporary workspace
2069 real(dp) :: work(6*max(params % pm, params % pl)**2 + &
2070 & 19*max(params % pm, params % pl) + 8)
2071 ! Local variables
2072 integer :: i, j, k
2073 real(dp) :: c1(3), c(3), r1, r
2074 ! Any order of this cycle is OK
2075 !$omp parallel do default(none) shared(constants,params,node_m,node_l) &
2076 !$omp private(i,c,r,k,c1,r1,work) schedule(dynamic)
2077 do i = 1, constants % nclusters
2078 ! If no far admissible pairs just set output to zero
2079 if (constants % nfar(i) .eq. 0) then
2080 node_l(:, i) = zero
2081 cycle
2082 end if
2083 c = constants % cnode(:, i)
2084 r = constants % rnode(i)
2085 ! Use the first far admissible pair to initialize output
2086 k = constants % far(constants % sfar(i))
2087 c1 = constants % cnode(:, k)
2088 r1 = constants % rnode(k)
2089 c1 = c1 - c
2090 call fmm_m2l_rotation_work(c1, r1, r, params % pm, params % pl, &
2091 & constants % vscales, constants % m2l_ztranslate_coef, one, &
2092 & node_m(:, k), zero, node_l(:, i), work)
2093 do j = constants % sfar(i)+1, constants % sfar(i+1)-1
2094 k = constants % far(j)
2095 c1 = constants % cnode(:, k)
2096 r1 = constants % rnode(k)
2097 c1 = c1 - c
2098 call fmm_m2l_rotation_work(c1, r1, r, params % pm, &
2099 & params % pl, constants % vscales, &
2100 & constants % m2l_ztranslate_coef, one, node_m(:, k), one, &
2101 & node_l(:, i), work)
2102 end do
2103 end do
2104end subroutine tree_m2l_rotation
2105
2106!------------------------------------------------------------------------------
2108!------------------------------------------------------------------------------
2109subroutine tree_m2l_bessel_rotation(params, constants, node_m, node_l)
2110 use complex_bessel
2111 implicit none
2112 ! Inputs
2113 type(ddx_params_type), intent(in) :: params
2114 type(ddx_constants_type), intent(in) :: constants
2115 real(dp), intent(in) :: node_m((params % pm+1)**2, constants % nclusters)
2116 ! Output
2117 real(dp), intent(out) :: node_l((params % pl+1)**2, constants % nclusters)
2118 ! Temporary workspace
2119 real(dp) :: work(6*params % pm**2 + 19*params % pm + 8)
2120 complex(dp) :: work_complex(2*params % pm+1)
2121 ! Local variables
2122 integer :: i, j, k
2123 real(dp) :: c1(3), c(3), r1, r
2124 ! Any order of this cycle is OK
2125 !$omp parallel do default(none) shared(constants,params,node_m,node_l) &
2126 !$omp private(i,c,r,k,c1,r1,work,work_complex) schedule(dynamic)
2127 do i = 1, constants % nclusters
2128 ! If no far admissible pairs just set output to zero
2129 if (constants % nfar(i) .eq. 0) then
2130 node_l(:, i) = zero
2131 cycle
2132 end if
2133 c = constants % cnode(:, i)
2134 r = constants % rnode(i)
2135 ! Use the first far admissible pair to initialize output
2136 k = constants % far(constants % sfar(i))
2137 c1 = constants % cnode(:, k)
2138 r1 = constants % rnode(k)
2139 c1 = params % kappa*(c1 - c)
2140 call fmm_m2l_bessel_rotation_work(c1, &
2141 & constants % SK_rnode(:, k), constants % SI_rnode(:, i), &
2142 & params % pm, &
2143 & constants % vscales, one, &
2144 & node_m(:, k), zero, node_l(:, i), work, work_complex)
2145 do j = constants % sfar(i)+1, constants % sfar(i+1)-1
2146 k = constants % far(j)
2147 c1 = constants % cnode(:, k)
2148 r1 = constants % rnode(k)
2149 c1 = params % kappa*(c1 - c)
2150 call fmm_m2l_bessel_rotation_work(c1, constants % SK_rnode(:, k), &
2151 & constants % SI_rnode(:, i), params % pm, &
2152 & constants % vscales, one, &
2153 & node_m(:, k), one, node_l(:, i), work, work_complex)
2154 end do
2155 end do
2156end subroutine tree_m2l_bessel_rotation
2157!------------------------------------------------------------------------------
2159!------------------------------------------------------------------------------
2160subroutine tree_m2l_bessel_rotation_adj(params, constants, node_l, node_m)
2161 implicit none
2162 ! Inputs
2163 type(ddx_params_type), intent(in) :: params
2164 type(ddx_constants_type), intent(in) :: constants
2165 real(dp), intent(in) :: node_l((params % pl+1)**2, constants % nclusters)
2166 ! Output
2167 real(dp), intent(inout) :: node_m((params % pm+1)**2, constants % nclusters)
2168 !Call corresponding work routine
2169 call tree_m2l_bessel_rotation_adj_work(params, constants, node_l, node_m)
2170end subroutine tree_m2l_bessel_rotation_adj
2171
2172
2173!------------------------------------------------------------------------------
2175!------------------------------------------------------------------------------
2176subroutine tree_m2l_bessel_rotation_adj_work(params, constants, node_l, node_m)
2177 implicit none
2178 ! Inputs
2179 type(ddx_params_type), intent(in) :: params
2180 type(ddx_constants_type), intent(in) :: constants
2181 real(dp), intent(in) :: node_l((params % pl+1)**2, constants % nclusters)
2182 ! Output
2183 real(dp), intent(out) :: node_m((params % pm+1)**2, constants % nclusters)
2184 ! Temporary workspace
2185 real(dp) :: work(6*params % pm**2 + 19*params % pm + 8)
2186 complex(dp) :: work_complex(2*params % pm+1)
2187 ! Local variables
2188 integer :: i, j, k
2189 real(dp) :: c1(3), c(3), r
2190 node_m = zero
2191 !$omp parallel do shared(constants,params,node_l,node_m) &
2192 !$omp private(i,c,r,k,c1,work,work_complex,j) schedule(dynamic)
2193 do i = 1, constants % nclusters
2194 ! If no far admissible pairs just set output to zero
2195 if (constants % nfar(i) .eq. 0) then
2196 cycle
2197 end if
2198 c = constants % cnode(:, i)
2199 r = constants % rnode(i)
2200 ! Use the first far admissible pair to initialize output
2201 k = constants % far(constants % sfar(i))
2202 c1 = constants % cnode(:, k)
2203 c1 = params % kappa*(c - c1)
2204 call fmm_m2l_bessel_rotation_adj_work(c1, constants % SI_rnode(:, k), &
2205 & constants % SK_rnode(:, i), params % pm, constants % vscales, &
2206 & one, node_l(:, k), one, node_m(:, i), work, work_complex)
2207 do j = constants % sfar(i)+1, constants % sfar(i+1)-1
2208 k = constants % far(j)
2209 c1 = constants % cnode(:, k)
2210 c1 = params % kappa*(c - c1)
2211 call fmm_m2l_bessel_rotation_adj_work(c1, &
2212 & constants % SI_rnode(:, k), constants % SK_rnode(:, i), &
2213 & params % pm, constants % vscales, one, node_l(:, k), &
2214 & one, node_m(:, i), work, work_complex)
2215 end do
2216 end do
2217end subroutine tree_m2l_bessel_rotation_adj_work
2218
2219!------------------------------------------------------------------------------
2221!------------------------------------------------------------------------------
2222subroutine tree_m2l_rotation_adj(params, constants, node_l, node_m)
2223 implicit none
2224 ! Inputs
2225 type(ddx_params_type), intent(in) :: params
2226 type(ddx_constants_type), intent(in) :: constants
2227 real(dp), intent(in) :: node_l((params % pl+1)**2, constants % nclusters)
2228 ! Output
2229 real(dp), intent(out) :: node_m((params % pm+1)**2, constants % nclusters)
2230 ! Local variables
2231 integer :: i, j, k
2232 real(dp) :: c1(3), c(3), r1, r
2233 node_m = zero
2234 !$omp parallel do default(none) shared(constants,params,node_m,node_l) &
2235 !$omp private(i,c,r,k,c1,r1,j) schedule(dynamic)
2236 do i = 1, constants % nclusters
2237 ! If no far admissible pairs just set output to zero
2238 if (constants % nfar(i) .eq. 0) then
2239 cycle
2240 end if
2241 c = constants % cnode(:, i)
2242 r = constants % rnode(i)
2243 ! Use the first far admissible pair to initialize output
2244 k = constants % far(constants % sfar(i))
2245 c1 = constants % cnode(:, k)
2246 r1 = constants % rnode(k)
2247 c1 = c1 - c
2248 call fmm_m2l_rotation_adj(c1, r1, r, params % pl, params % pm, &
2249 & constants % vscales, constants % m2l_ztranslate_adj_coef, one, &
2250 & node_l(:, k), one, node_m(:, i))
2251 do j = constants % sfar(i)+1, constants % sfar(i+1)-1
2252 k = constants % far(j)
2253 c1 = constants % cnode(:, k)
2254 r1 = constants % rnode(k)
2255 c1 = c1 - c
2256 call fmm_m2l_rotation_adj(c1, r1, r, params % pl, params % pm, &
2257 & constants % vscales, constants % m2l_ztranslate_adj_coef, one, &
2258 & node_l(:, k), one, node_m(:, i))
2259 end do
2260 end do
2261end subroutine tree_m2l_rotation_adj
2262
2263!------------------------------------------------------------------------------
2265!------------------------------------------------------------------------------
2266subroutine tree_l2p(params, constants, alpha, node_l, beta, grid_v, sph_l)
2267 implicit none
2268 ! Inputs
2269 type(ddx_params_type), intent(in) :: params
2270 type(ddx_constants_type), intent(in) :: constants
2271 real(dp), intent(in) :: node_l((params % pl+1)**2, constants % nclusters), &
2272 & alpha, beta
2273 ! Output
2274 real(dp), intent(inout) :: grid_v(params % ngrid, params % nsph)
2275 ! Scratch
2276 real(dp), intent(out) :: sph_l((params % pl+1)**2, params % nsph)
2277 ! Local variables
2278 integer :: isph
2279 external :: dgemm
2280
2281
2282 ! Init output
2283 if (beta .eq. zero) then
2284 grid_v = zero
2285 else
2286 grid_v = beta * grid_v
2287 end if
2288 ! Get data from all clusters to spheres
2289 !$omp parallel do default(none) shared(params,constants,node_l,sph_l) &
2290 !$omp private(isph) schedule(dynamic)
2291 do isph = 1, params % nsph
2292 sph_l(:, isph) = node_l(:, constants % snode(isph))
2293 end do
2294 ! Get values at grid points
2295 call dgemm('T', 'N', params % ngrid, params % nsph, &
2296 & (params % pl+1)**2, alpha, constants % vgrid2, &
2297 & constants % vgrid_nbasis, sph_l, (params % pl+1)**2, beta, grid_v, &
2298 & params % ngrid)
2299
2300end subroutine tree_l2p
2301
2302!------------------------------------------------------------------------------
2304!------------------------------------------------------------------------------
2305subroutine tree_l2p_bessel(params, constants, alpha, node_l, beta, grid_v)
2306 implicit none
2307 ! Inputs
2308 type(ddx_params_type), intent(in) :: params
2309 type(ddx_constants_type), intent(in) :: constants
2310 real(dp), intent(in) :: node_l((params % pl+1)**2, constants % nclusters), &
2311 & alpha, beta
2312 ! Output
2313 real(dp), intent(inout) :: grid_v(params % ngrid, params % nsph)
2314 ! Local variables
2315 real(dp) :: sph_l((params % pl+1)**2, params % nsph)
2316 integer :: isph
2317 external :: dgemm
2318 ! Init output
2319 if (beta .eq. zero) then
2320 grid_v = zero
2321 else
2322 grid_v = beta * grid_v
2323 end if
2324 ! Get data from all clusters to spheres
2325 !$omp parallel do default(none) shared(params,constants,node_l,sph_l) &
2326 !$omp private(isph) schedule(dynamic)
2327 do isph = 1, params % nsph
2328 sph_l(:, isph) = node_l(:, constants % snode(isph))
2329 end do
2330 ! Get values at grid points
2331 call dgemm('T', 'N', params % ngrid, params % nsph, &
2332 & (params % pl+1)**2, alpha, constants % vgrid, &
2333 & constants % vgrid_nbasis, sph_l, (params % pl+1)**2, beta, grid_v, &
2334 & params % ngrid)
2335end subroutine tree_l2p_bessel
2336
2337!------------------------------------------------------------------------------
2339!------------------------------------------------------------------------------
2340subroutine tree_l2p_adj(params, constants, alpha, grid_v, beta, node_l, sph_l)
2341 implicit none
2342 ! Inputs
2343 type(ddx_params_type), intent(in) :: params
2344 type(ddx_constants_type), intent(in) :: constants
2345 real(dp), intent(in) :: grid_v(params % ngrid, params % nsph), alpha, &
2346 & beta
2347 ! Output
2348 real(dp), intent(inout) :: node_l((params % pl+1)**2, &
2349 & constants % nclusters)
2350 ! Scractch
2351 real(dp), intent(out) :: sph_l((params % pl+1)**2, params % nsph)
2352 ! Local variables
2353 integer :: isph, inode
2354 external :: dgemm
2355 ! Init output
2356 if (beta .eq. zero) then
2357 node_l = zero
2358 else
2359 node_l = beta * node_l
2360 end if
2361 ! Get weights of spherical harmonics at each sphere
2362 call dgemm('N', 'N', (params % pl+1)**2, params % nsph, &
2363 & params % ngrid, one, constants % vgrid2, constants % vgrid_nbasis, &
2364 & grid_v, params % ngrid, zero, sph_l, &
2365 & (params % pl+1)**2)
2366 ! Get data from all clusters to spheres
2367 !$omp parallel do default(none) shared(params,constants,node_l,sph_l, &
2368 !$omp alpha) private(isph,inode) schedule(dynamic)
2369 do isph = 1, params % nsph
2370 inode = constants % snode(isph)
2371 node_l(:, inode) = node_l(:, inode) + alpha*sph_l(:, isph)
2372 end do
2373end subroutine tree_l2p_adj
2374
2375!------------------------------------------------------------------------------
2377!------------------------------------------------------------------------------
2378subroutine tree_l2p_bessel_adj(params, constants, alpha, grid_v, beta, node_l)
2379 implicit none
2380 ! Inputs
2381 type(ddx_params_type), intent(in) :: params
2382 type(ddx_constants_type), intent(in) :: constants
2383 real(dp), intent(in) :: grid_v(params % ngrid, params % nsph), alpha, &
2384 & beta
2385 ! Output
2386 real(dp), intent(inout) :: node_l((params % pl+1)**2, &
2387 & constants % nclusters)
2388 ! Local variables
2389 real(dp) :: sph_l((params % pl+1)**2, params % nsph)
2390 integer :: isph, inode
2391 external :: dgemm
2392 ! Init output
2393 if (beta .eq. zero) then
2394 node_l = zero
2395 else
2396 node_l = beta * node_l
2397 end if
2398 ! Get weights of spherical harmonics at each sphere
2399 call dgemm('N', 'N', (params % pl+1)**2, params % nsph, &
2400 & params % ngrid, one, constants % vgrid, constants % vgrid_nbasis, &
2401 & grid_v, params % ngrid, zero, sph_l, &
2402 & (params % pl+1)**2)
2403 ! Get data from all clusters to spheres
2404 !$omp parallel do default(none) shared(params,constants,node_l,sph_l, &
2405 !$omp alpha) private(isph,inode) schedule(dynamic)
2406 do isph = 1, params % nsph
2407 inode = constants % snode(isph)
2408 node_l(:, inode) = node_l(:, inode) + alpha*sph_l(:, isph)
2409 end do
2410end subroutine tree_l2p_bessel_adj
2411
2412!------------------------------------------------------------------------------
2414!------------------------------------------------------------------------------
2415subroutine tree_m2p(params, constants, p, alpha, sph_m, beta, grid_v)
2416 implicit none
2417 ! Inputs
2418 type(ddx_params_type), intent(in) :: params
2419 type(ddx_constants_type), intent(in) :: constants
2420 integer, intent(in) :: p
2421 real(dp), intent(in) :: sph_m((p+1)**2, params % nsph), alpha, beta
2422 ! Output
2423 real(dp), intent(inout) :: grid_v(params % ngrid, params % nsph)
2424 ! Local variables
2425 integer :: isph, inode, jnear, jnode, jsph, igrid
2426 real(dp) :: c(3)
2427 ! Temporary workspace
2428 real(dp) :: work(p+1)
2429 ! Init output
2430 if (beta .eq. zero) then
2431 grid_v = zero
2432 else
2433 grid_v = beta * grid_v
2434 end if
2435 ! Cycle over all spheres
2436 !$omp parallel do default(none) shared(params,constants,grid_v,p, &
2437 !$omp alpha,sph_m), private(isph,inode,jnear,jnode,jsph,igrid,c,work) &
2438 !$omp schedule(dynamic)
2439 do isph = 1, params % nsph
2440 ! Cycle over all near-field admissible pairs of spheres
2441 inode = constants % snode(isph)
2442 do jnear = constants % snear(inode), constants % snear(inode+1)-1
2443 ! Near-field interactions are possible only between leaf nodes,
2444 ! which must contain only a single input sphere
2445 jnode = constants % near(jnear)
2446 jsph = constants % order(constants % cluster(1, jnode))
2447 ! Ignore self-interaction
2448 if(isph .eq. jsph) cycle
2449 ! Accumulate interaction for external grid points only
2450 do igrid = 1, params % ngrid
2451 if(constants % ui(igrid, isph) .eq. zero) cycle
2452 c = constants % cgrid(:, igrid)*params % rsph(isph) - &
2453 & params % csph(:, jsph) + params % csph(:, isph)
2454 call fmm_m2p_work(c, params % rsph(jsph), p, &
2455 & constants % vscales_rel, alpha, sph_m(:, jsph), one, &
2456 & grid_v(igrid, isph), work)
2457 end do
2458 end do
2459 end do
2460end subroutine tree_m2p
2461
2462!------------------------------------------------------------------------------
2464!------------------------------------------------------------------------------
2465subroutine tree_m2p_bessel(params, constants, p, alpha, sph_p, sph_m, beta, grid_v)
2466 implicit none
2467 ! Inputs
2468 type(ddx_params_type), intent(in) :: params
2469 type(ddx_constants_type), intent(in) :: constants
2470 integer, intent(in) :: p, sph_p
2471 real(dp), intent(in) :: sph_m((sph_p+1)**2, params % nsph), alpha, beta
2472 ! Output
2473 real(dp), intent(inout) :: grid_v(params % ngrid, params % nsph)
2474 ! Local variables
2475 integer :: isph, inode, jnear, jnode, jsph, igrid
2476 real(dp) :: c(3)
2477 ! Temporary workspace
2478 real(dp) :: work(p+1)
2479 complex(dp) :: work_complex(p+1)
2480 ! Init output
2481 if (beta .eq. zero) then
2482 grid_v = zero
2483 else
2484 grid_v = beta * grid_v
2485 end if
2486 ! Cycle over all spheres
2487 !$omp parallel do default(none) shared(params,constants,grid_v,p, &
2488 !$omp alpha,sph_m), private(isph,inode,jnear,jnode,jsph,igrid,c,work, &
2489 !$omp work_complex) schedule(dynamic)
2490 do isph = 1, params % nsph
2491 ! Cycle over all near-field admissible pairs of spheres
2492 inode = constants % snode(isph)
2493 do jnear = constants % snear(inode), constants % snear(inode+1)-1
2494 ! Near-field interactions are possible only between leaf nodes,
2495 ! which must contain only a single input sphere
2496 jnode = constants % near(jnear)
2497 jsph = constants % order(constants % cluster(1, jnode))
2498 ! Ignore self-interaction
2499 !if(isph .eq. jsph) cycle
2500 ! Accumulate interaction for external grid points only
2501 do igrid = 1, params % ngrid
2502 if(constants % ui(igrid, isph) .eq. zero) cycle
2503 c = constants % cgrid(:, igrid)*params % rsph(isph) - &
2504 & params % csph(:, jsph) + params % csph(:, isph)
2505 c = c * params % kappa
2506 call fmm_m2p_bessel_work(c, p, constants % vscales, &
2507 & constants % SK_ri(:, jsph), alpha, sph_m(:, jsph), one, &
2508 & grid_v(igrid, isph), work_complex, work)
2509 end do
2510 end do
2511 end do
2512end subroutine tree_m2p_bessel
2513
2514!------------------------------------------------------------------------------
2516!------------------------------------------------------------------------------
2517subroutine tree_m2p_adj(params, constants, p, alpha, grid_v, beta, sph_m)
2518 implicit none
2519 ! Inputs
2520 type(ddx_params_type), intent(in) :: params
2521 type(ddx_constants_type), intent(in) :: constants
2522 integer, intent(in) :: p
2523 real(dp), intent(in) :: grid_v(params % ngrid, params % nsph), alpha, &
2524 & beta
2525 ! Output
2526 real(dp), intent(inout) :: sph_m((p+1)**2, params % nsph)
2527 ! Temporary workspace
2528 real(dp) :: work(p+1)
2529 ! Local variables
2530 integer :: isph, inode, jnear, jnode, jsph, igrid
2531 real(dp) :: c(3)
2532 ! Init output
2533 if (beta .eq. zero) then
2534 sph_m = zero
2535 else
2536 sph_m = beta * sph_m
2537 end if
2538 ! Cycle over all spheres
2539 !$omp parallel do default(none) shared(params,constants,grid_v,p, &
2540 !$omp alpha,sph_m), private(isph,inode,jnear,jnode,jsph,igrid,c,work) &
2541 !$omp schedule(dynamic)
2542 do isph = 1, params % nsph
2543 ! Cycle over all near-field admissible pairs of spheres
2544 inode = constants % snode(isph)
2545 do jnear = constants % snear(inode), constants % snear(inode+1)-1
2546 ! Near-field interactions are possible only between leaf nodes,
2547 ! which must contain only a single input sphere
2548 jnode = constants % near(jnear)
2549 jsph = constants % order(constants % cluster(1, jnode))
2550 ! Ignore self-interaction
2551 if(isph .eq. jsph) cycle
2552 ! Accumulate interaction for external grid points only
2553 do igrid = 1, params % ngrid
2554 if(constants % ui(igrid, jsph) .eq. zero) cycle
2555 c = constants % cgrid(:, igrid)*params % rsph(jsph) - &
2556 & params % csph(:, isph) + params % csph(:, jsph)
2557 call fmm_m2p_adj_work(c, alpha*grid_v(igrid, jsph), &
2558 & params % rsph(isph), p, constants % vscales_rel, one, &
2559 & sph_m(:, isph), work)
2560 end do
2561 end do
2562 end do
2563end subroutine tree_m2p_adj
2564
2565!------------------------------------------------------------------------------
2567!------------------------------------------------------------------------------
2568subroutine tree_m2p_bessel_adj(params, constants, p, alpha, grid_v, beta, sph_p, &
2569 & sph_m)
2570 implicit none
2571 ! Inputs
2572 type(ddx_params_type), intent(in) :: params
2573 type(ddx_constants_type), intent(in) :: constants
2574 integer, intent(in) :: p, sph_p
2575 real(dp), intent(in) :: grid_v(params % ngrid, params % nsph), alpha, &
2576 & beta
2577 ! Output
2578 real(dp), intent(inout) :: sph_m((sph_p+1)**2, params % nsph)
2579 ! Local variables
2580 !real(dp), allocatable :: tmp(:, :, :)
2581 integer :: isph, inode, jnear, jnode, jsph, igrid
2582 real(dp) :: c(3)
2583 ! Init output
2584 if (beta .eq. zero) then
2585 sph_m = zero
2586 else
2587 sph_m = beta * sph_m
2588 end if
2589 ! Cycle over all spheres
2590 !$omp parallel do default(none) shared(params,constants,p,sph_m, &
2591 !$omp alpha,grid_v) private(isph,inode,jnear,jnode,jsph,igrid,c) &
2592 !$omp schedule(dynamic)
2593 do isph = 1, params % nsph
2594 ! Cycle over all near-field admissible pairs of spheres
2595 inode = constants % snode(isph)
2596 do jnear = constants % snear(inode), constants % snear(inode+1)-1
2597 ! Near-field interactions are possible only between leaf nodes,
2598 ! which must contain only a single input sphere
2599 jnode = constants % near(jnear)
2600 jsph = constants % order(constants % cluster(1, jnode))
2601 ! Accumulate interaction for external grid points only
2602 do igrid = 1, params % ngrid
2603 if(constants % ui(igrid, jsph) .eq. zero) cycle
2604 c = constants % cgrid(:, igrid)*params % rsph(jsph) - &
2605 & params % csph(:, isph) + params % csph(:, jsph)
2606 call fmm_m2p_bessel_adj(c, alpha*grid_v(igrid, jsph), &
2607 & params % rsph(isph), params % kappa, p, &
2608 & constants % vscales, one, sph_m(:, isph))
2609 end do
2610 end do
2611 end do
2612end subroutine tree_m2p_bessel_adj
2613
2614!------------------------------------------------------------------------------
2616!------------------------------------------------------------------------------
2617subroutine tree_m2p_bessel_nodiag_adj(params, constants, p, alpha, grid_v, beta, sph_p, &
2618 & sph_m)
2619 implicit none
2620 ! Inputs
2621 type(ddx_params_type), intent(in) :: params
2622 type(ddx_constants_type), intent(in) :: constants
2623 integer, intent(in) :: p, sph_p
2624 real(dp), intent(in) :: grid_v(params % ngrid, params % nsph), alpha, &
2625 & beta
2626 ! Output
2627 real(dp), intent(inout) :: sph_m((sph_p+1)**2, params % nsph)
2628 ! Local variables
2629 integer :: isph, inode, jnear, jnode, jsph, igrid
2630 real(dp) :: c(3)
2631 ! Init output
2632 if (beta .eq. zero) then
2633 sph_m = zero
2634 else
2635 sph_m = beta * sph_m
2636 end if
2637 ! Cycle over all spheres
2638 do isph = 1, params % nsph
2639 ! Cycle over all near-field admissible pairs of spheres
2640 inode = constants % snode(isph)
2641 do jnear = constants % snear(inode), constants % snear(inode+1)-1
2642 ! Near-field interactions are possible only between leaf nodes,
2643 ! which must contain only a single input sphere
2644 jnode = constants % near(jnear)
2645 jsph = constants % order(constants % cluster(1, jnode))
2646 ! Ignore self-interaction
2647 if(isph .eq. jsph) cycle
2648 ! Accumulate interaction for external grid points only
2649 do igrid = 1, params % ngrid
2650 if(constants % ui(igrid, isph) .eq. zero) cycle
2651 c = constants % cgrid(:, igrid)*params % rsph(isph) - &
2652 & params % csph(:, jsph) + params % csph(:, isph)
2653 call fmm_m2p_bessel_adj(c, alpha*grid_v(igrid, isph), &
2654 & params % rsph(jsph), params % kappa, p, constants % vscales, one, &
2655 & sph_m(:, jsph))
2656 end do
2657 end do
2658 end do
2659end subroutine tree_m2p_bessel_nodiag_adj
2660
2661!------------------------------------------------------------------------------
2663!------------------------------------------------------------------------------
2664subroutine tree_grad_m2m(params, constants, sph_m, sph_m_grad, work)
2665 implicit none
2666 ! Inputs
2667 type(ddx_params_type), intent(in) :: params
2668 type(ddx_constants_type), intent(in) :: constants
2669 real(dp), intent(in) :: sph_m(constants % nbasis, params % nsph)
2670 ! Output
2671 real(dp), intent(inout) :: sph_m_grad((params % lmax+2)**2, 3, &
2672 & params % nsph)
2673 ! Temporary workspace
2674 real(dp), intent(inout) :: work((params % lmax+2)**2, params % nsph)
2675 ! Local variables
2676 integer :: isph, l, indi, indj, m
2677 real(dp) :: tmp1, tmp2
2678 real(dp), dimension(3, 3) :: zx_coord_transform, zy_coord_transform
2679 ! Set coordinate transformations
2680 zx_coord_transform = zero
2681 zx_coord_transform(3, 2) = one
2682 zx_coord_transform(2, 3) = one
2683 zx_coord_transform(1, 1) = one
2684 zy_coord_transform = zero
2685 zy_coord_transform(1, 2) = one
2686 zy_coord_transform(2, 1) = one
2687 zy_coord_transform(3, 3) = one
2688 ! At first reflect harmonics of a degree up to lmax
2689 sph_m_grad(1:constants % nbasis, 3, :) = sph_m
2690 do isph = 1, params % nsph
2691 call fmm_sph_transform(params % lmax, zx_coord_transform, one, &
2692 & sph_m(:, isph), zero, sph_m_grad(1:constants % nbasis, 1, isph))
2693 call fmm_sph_transform(params % lmax, zy_coord_transform, one, &
2694 & sph_m(:, isph), zero, sph_m_grad(1:constants % nbasis, 2, isph))
2695 end do
2696 ! Derivative of M2M translation over OZ axis at the origin consists of 2
2697 ! steps:
2698 ! 1) increase degree l and scale by sqrt((2*l+1)*(l*l-m*m)) / sqrt(2*l-1)
2699 ! 2) scale by 1/rsph(isph)
2700 do l = params % lmax+1, 1, -1
2701 indi = l*l + l + 1
2702 indj = indi - 2*l
2703 tmp1 = sqrt(dble(2*l+1)) / sqrt(dble(2*l-1))
2704 do m = 1-l, l-1
2705 tmp2 = sqrt(dble(l*l-m*m)) * tmp1
2706 sph_m_grad(indi+m, :, :) = tmp2 * sph_m_grad(indj+m, :, :)
2707 end do
2708 sph_m_grad(indi+l, :, :) = zero
2709 sph_m_grad(indi-l, :, :) = zero
2710 end do
2711 sph_m_grad(1, :, :) = zero
2712 ! Scale by 1/rsph(isph) and rotate harmonics of degree up to lmax+1 back to
2713 ! the initial axis. Coefficient of 0-th degree is zero so we ignore it.
2714 do isph = 1, params % nsph
2715 sph_m_grad(:, 3, isph) = sph_m_grad(:, 3, isph) / params % rsph(isph)
2716 work(:, isph) = sph_m_grad(:, 1, isph) / params % rsph(isph)
2717 call fmm_sph_transform(params % lmax+1, zx_coord_transform, one, &
2718 & work(:, isph), zero, sph_m_grad(:, 1, isph))
2719 work(:, isph) = sph_m_grad(:, 2, isph) / params % rsph(isph)
2720 call fmm_sph_transform(params % lmax+1, zy_coord_transform, one, &
2721 & work(:, isph), zero, sph_m_grad(:, 2, isph))
2722 end do
2723end subroutine tree_grad_m2m
2724
2725!------------------------------------------------------------------------------
2727!------------------------------------------------------------------------------
2728subroutine tree_grad_l2l(params, constants, node_l, sph_l_grad, work)
2729 implicit none
2730 ! Inputs
2731 type(ddx_params_type), intent(in) :: params
2732 type(ddx_constants_type), intent(in) :: constants
2733 real(dp), intent(in) :: node_l((params % pl+1)**2, constants % nclusters)
2734 ! Output
2735 real(dp), intent(out) :: sph_l_grad((params % pl+1)**2, 3, params % nsph)
2736 ! Temporary workspace
2737 real(dp), intent(out) :: work((params % pl+1)**2, params % nsph)
2738 ! Local variables
2739 integer :: isph, inode, l, indi, indj, m
2740 real(dp) :: tmp1, tmp2
2741 real(dp), dimension(3, 3) :: zx_coord_transform, zy_coord_transform
2742 ! Gradient of L2L reduces degree by 1, so exit if degree of harmonics is 0
2743 ! or -1 (which means no FMM at all)
2744 if (params % pl .le. 0) return
2745 ! Set coordinate transformations
2746 zx_coord_transform = zero
2747 zx_coord_transform(3, 2) = one
2748 zx_coord_transform(2, 3) = one
2749 zx_coord_transform(1, 1) = one
2750 zy_coord_transform = zero
2751 zy_coord_transform(1, 2) = one
2752 zy_coord_transform(2, 1) = one
2753 zy_coord_transform(3, 3) = one
2754 ! At first reflect harmonics of a degree up to pl
2755 do isph = 1, params % nsph
2756 inode = constants % snode(isph)
2757 sph_l_grad(:, 3, isph) = node_l(:, inode)
2758 call fmm_sph_transform(params % pl, zx_coord_transform, one, &
2759 & node_l(:, inode), zero, sph_l_grad(:, 1, isph))
2760 call fmm_sph_transform(params % pl, zy_coord_transform, one, &
2761 & node_l(:, inode), zero, sph_l_grad(:, 2, isph))
2762 end do
2763 ! Derivative of L2L translation over OZ axis at the origin consists of 2
2764 ! steps:
2765 ! 1) decrease degree l and scale by sqrt((2*l-1)*(l*l-m*m)) / sqrt(2*l+1)
2766 ! 2) scale by 1/rsph(isph)
2767 do l = 1, params % pl
2768 indi = l*l + l + 1
2769 indj = indi - 2*l
2770 tmp1 = -sqrt(dble(2*l-1)) / sqrt(dble(2*l+1))
2771 do m = 1-l, l-1
2772 tmp2 = sqrt(dble(l*l-m*m)) * tmp1
2773 sph_l_grad(indj+m, :, :) = tmp2 * sph_l_grad(indi+m, :, :)
2774 end do
2775 end do
2776 ! Scale by 1/rsph(isph) and rotate harmonics of degree up to pl-1 back to
2777 ! the initial axis. Coefficient of pl-th degree is zero so we ignore it.
2778 do isph = 1, params % nsph
2779 sph_l_grad(1:params % pl**2, 3, isph) = &
2780 & sph_l_grad(1:params % pl**2, 3, isph) / params % rsph(isph)
2781 work(1:params % pl**2, isph) = &
2782 & sph_l_grad(1:params % pl**2, 1, isph) / params % rsph(isph)
2783 call fmm_sph_transform(params % pl-1, zx_coord_transform, one, &
2784 & work(1:params % pl**2, isph), zero, &
2785 & sph_l_grad(1:params % pl**2, 1, isph))
2786 work(1:params % pl**2, isph) = &
2787 & sph_l_grad(1:params % pl**2, 2, isph) / params % rsph(isph)
2788 call fmm_sph_transform(params % pl-1, zy_coord_transform, one, &
2789 & work(1:params % pl**2, isph), zero, &
2790 & sph_l_grad(1:params % pl**2, 2, isph))
2791 end do
2792 ! Set degree pl to zero to avoid problems if user actually uses it
2793 l = params % pl
2794 indi = l*l + l + 1
2795 sph_l_grad(indi-l:indi+l, :, :) = zero
2796end subroutine tree_grad_l2l
2797
2803subroutine get_banner(string)
2804 implicit none
2805 character (len=2047), intent(out) :: string
2806 character (len=10) :: vstr
2807 write(vstr, *) "0.6.8"
2808 write(string, *) &
2809 & " +----------------------------------------------------------------+", &
2810 & new_line('a'), &
2811 & " | |", &
2812 & new_line('a'), &
2813 & " | 888 888 Y8b d8Y |", &
2814 & new_line('a'), &
2815 & " | 888 888 Y8b d8Y |", &
2816 & new_line('a'), &
2817 & " | 888 888 Y8888Y |", &
2818 & new_line('a'), &
2819 & " | .d88888 .d88888 Y88Y |", &
2820 & new_line('a'), &
2821 & " | d88 888 d88 888 d88b |", &
2822 & new_line('a'), &
2823 & " | 888 888 888 888 d8888b |", &
2824 & new_line('a'), &
2825 & " | Y88b 888 Y88b 888 d8Y Y8b |", &
2826 & new_line('a'), &
2827 & " | Y88888 Y88888 d8Y Y8b |", &
2828 & new_line('a'), &
2829 & " | |", &
2830 & new_line('a'), &
2831 & " | https://ddsolvation.github.io/ddX/ |", &
2832 & new_line('a'), &
2833 & " | Version:", vstr, " |", &
2834 & new_line('a'), &
2835 & " | |", &
2836 & new_line('a'), &
2837 & " +----------------------------------------------------------------+"
2838end subroutine get_banner
2839
2851subroutine cav_to_spherical(params, constants, workspace, property_cav, &
2852 & property_sph)
2853 implicit none
2854 type(ddx_params_type), intent(in) :: params
2855 type(ddx_constants_type), intent(in) :: constants
2856 type(ddx_workspace_type), intent(inout) :: workspace
2857 real(dp), intent(in) :: property_cav(constants % ncav)
2858 real(dp), intent(out) :: property_sph(constants % nbasis, params % nsph)
2859
2860 ! multiply by the characteristic function U
2861 workspace % tmp_cav = property_cav * constants % ui_cav
2862
2863 ! extend the function to the sphere intersection with zeros
2864 call ddcav_to_grid_work(params % ngrid, params % nsph, constants % ncav, &
2865 & constants % icav_ia, constants % icav_ja, workspace % tmp_cav, &
2866 & workspace % tmp_grid)
2867
2868 ! integrate against spherical harmonics
2869 call ddintegrate(params % nsph, constants % nbasis, &
2870 & params % ngrid, constants % vwgrid, constants % vgrid_nbasis, &
2871 & workspace % tmp_grid, property_sph)
2872end subroutine cav_to_spherical
2873
2874end module ddx_core
ddx_core::allocate_electrostatics
subroutine allocate_electrostatics(params, constants, electrostatics, ddx_error)
Given the chosen model, find the required electrostatic properties, and allocate the arrays for them ...
Definition: ddx_core.f90:349
ddx_core::deallocate_electrostatics
subroutine deallocate_electrostatics(electrostatics, ddx_error)
Deallocate the electrostatic properties.
Definition: ddx_core.f90:415
ddx_core::allocate_state
subroutine allocate_state(params, constants, state, ddx_error)
Initialize the ddx_state object.
Definition: ddx_core.f90:455
ddx_core::deallocate_state
subroutine deallocate_state(state, ddx_error)
Deallocate the ddx_state object.
Definition: ddx_core.f90:790
ddx_core::get_banner
subroutine get_banner(string)
Print the ddX logo.
Definition: ddx_core.f90:2804
ddx_constants
Run-time constants.
Definition: ddx_constants.f90:17
ddx_constants::constants_init
subroutine constants_init(params, constants, ddx_error)
Compute all necessary constants.
Definition: ddx_constants.f90:237
ddx_constants::constants_free
subroutine constants_free(constants, ddx_error)
Deallocate the constants.
Definition: ddx_constants.f90:1836
ddx_constants::fsw
real(dp) function fsw(t, se, eta)
Switching function.
Definition: ddx_constants.f90:1277
ddx_core
Core routines and parameters of the ddX software.
Definition: ddx_core.f90:13
ddx_core::tree_l2l_rotation_adj_work
subroutine tree_l2l_rotation_adj_work(params, constants, node_l, work)
Adjoint transfer local coefficients over a tree.
Definition: ddx_core.f90:1990
ddx_core::deallocate_model
subroutine deallocate_model(ddx_data, ddx_error)
Deallocate object with corresponding data.
Definition: ddx_core.f90:329
ddx_core::tree_grad_l2l
subroutine tree_grad_l2l(params, constants, node_l, sph_l_grad, work)
TODO.
Definition: ddx_core.f90:2729
ddx_core::tree_l2p_bessel
subroutine tree_l2p_bessel(params, constants, alpha, node_l, beta, grid_v)
TODO.
Definition: ddx_core.f90:2306
ddx_core::tree_l2p_bessel_adj
subroutine tree_l2p_bessel_adj(params, constants, alpha, grid_v, beta, node_l)
TODO.
Definition: ddx_core.f90:2379
ddx_core::tree_l2l_rotation
subroutine tree_l2l_rotation(params, constants, node_l)
Transfer local coefficients over a tree.
Definition: ddx_core.f90:1881
ddx_core::tree_m2m_rotation
subroutine tree_m2m_rotation(params, constants, node_m)
Transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1671
ddx_core::tree_l2p
subroutine tree_l2p(params, constants, alpha, node_l, beta, grid_v, sph_l)
TODO.
Definition: ddx_core.f90:2267
ddx_core::wghpot
subroutine wghpot(ncav, phi_cav, nsph, ngrid, ui, phi_grid, g)
Weigh potential at cavity points by characteristic function TODO use cavity points in CSR format.
Definition: ddx_core.f90:1353
ddx_core::hnorm
real(dp) function hnorm(lmax, nbasis, nsph, x)
Compute the global Sobolev H^(-1/2)-norm of x.
Definition: ddx_core.f90:1405
ddx_core::tree_l2p_adj
subroutine tree_l2p_adj(params, constants, alpha, grid_v, beta, node_l, sph_l)
TODO.
Definition: ddx_core.f90:2341
ddx_core::cav_to_spherical
subroutine cav_to_spherical(params, constants, workspace, property_cav, property_sph)
Transform a function defined at the exposed cavity points (cav) to a spherical harmonics expansion....
Definition: ddx_core.f90:2853
ddx_core::tree_l2l_rotation_work
subroutine tree_l2l_rotation_work(params, constants, node_l, work)
Transfer local coefficients over a tree.
Definition: ddx_core.f90:1897
ddx_core::tree_m2l_bessel_rotation
subroutine tree_m2l_bessel_rotation(params, constants, node_m, node_l)
Transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2110
ddx_core::tree_m2m_bessel_rotation_adj
subroutine tree_m2m_bessel_rotation_adj(params, constants, node_m)
Adjoint transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1837
ddx_core::rmsvec
subroutine rmsvec(n, v, vrms, vmax)
compute root-mean-square and max norm
Definition: ddx_core.f90:1437
ddx_core::print_nodes
subroutine print_nodes(iunit, label, ngrid, ncol, icol, x)
Print array of quadrature points.
Definition: ddx_core.f90:1069
ddx_core::tree_m2m_bessel_rotation_work
subroutine tree_m2m_bessel_rotation_work(params, constants, node_m)
Transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1748
ddx_core::ddeval_grid_work
subroutine ddeval_grid_work(nbasis, ngrid, nsph, vgrid, ldvgrid, alpha, x_sph, beta, x_grid)
Evaluate values of spherical harmonics at Lebedev grid points.
Definition: ddx_core.f90:1554
ddx_core::tree_m2l_rotation
subroutine tree_m2l_rotation(params, constants, node_m, node_l)
Transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2061
ddx_core::ddintegrate_sph
subroutine ddintegrate_sph(params, constants, alpha, x_grid, beta, x_sph)
Integrate values at grid points into spherical harmonics.
Definition: ddx_core.f90:1570
ddx_core::tree_m2p
subroutine tree_m2p(params, constants, p, alpha, sph_m, beta, grid_v)
TODO.
Definition: ddx_core.f90:2416
ddx_core::ddintegrate_sph_work
subroutine ddintegrate_sph_work(nbasis, ngrid, nsph, vwgrid, ldvwgrid, alpha, x_grid, beta, x_sph)
Integrate values at grid points into spherical harmonics.
Definition: ddx_core.f90:1587
ddx_core::tree_m2p_bessel
subroutine tree_m2p_bessel(params, constants, p, alpha, sph_p, sph_m, beta, grid_v)
TODO.
Definition: ddx_core.f90:2466
ddx_core::tree_l2l_bessel_rotation_adj
subroutine tree_l2l_bessel_rotation_adj(params, constants, node_l)
Adjoint transfer local coefficients over a tree.
Definition: ddx_core.f90:2031
ddx_core::ddintegrate
subroutine ddintegrate(nsph, nbasis, ngrid, vwgrid, ldvwgrid, x_grid, x_lm)
Integrate against spherical harmonics.
Definition: ddx_core.f90:1148
ddx_core::tree_l2l_bessel_rotation
subroutine tree_l2l_bessel_rotation(params, constants, node_l)
Transfer local coefficients over a tree.
Definition: ddx_core.f90:1928
ddx_core::print_ddvector
subroutine print_ddvector(ddx_data, label, vector)
Print dd Solution vector.
Definition: ddx_core.f90:1115
ddx_core::tree_l2l_bessel_rotation_work
subroutine tree_l2l_bessel_rotation_work(params, constants, node_l)
Transfer local coefficients over a tree.
Definition: ddx_core.f90:1944
ddx_core::tree_m2p_bessel_nodiag_adj
subroutine tree_m2p_bessel_nodiag_adj(params, constants, p, alpha, grid_v, beta, sph_p, sph_m)
TODO.
Definition: ddx_core.f90:2619
ddx_core::tree_m2l_bessel_rotation_adj
subroutine tree_m2l_bessel_rotation_adj(params, constants, node_l, node_m)
Adjoint transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2161
ddx_core::tree_m2m_rotation_work
subroutine tree_m2m_rotation_work(params, constants, node_m, work)
Transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1687
ddx_core::tree_grad_m2m
subroutine tree_grad_m2m(params, constants, sph_m, sph_m_grad, work)
TODO.
Definition: ddx_core.f90:2665
ddx_core::tree_m2p_adj
subroutine tree_m2p_adj(params, constants, p, alpha, grid_v, beta, sph_m)
TODO.
Definition: ddx_core.f90:2518
ddx_core::tree_m2p_bessel_adj
subroutine tree_m2p_bessel_adj(params, constants, p, alpha, grid_v, beta, sph_p, sph_m)
TODO.
Definition: ddx_core.f90:2570
ddx_core::hsnorm
subroutine hsnorm(lmax, nbasis, u, unorm)
Compute the local Sobolev H^(-1/2)-norm on one sphere of u.
Definition: ddx_core.f90:1384
ddx_core::tree_m2l_bessel_rotation_adj_work
subroutine tree_m2l_bessel_rotation_adj_work(params, constants, node_l, node_m)
Adjoint transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2177
ddx_core::tree_m2l_rotation_adj
subroutine tree_m2l_rotation_adj(params, constants, node_l, node_m)
Adjoint transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2223
ddx_core::tree_l2l_rotation_adj
subroutine tree_l2l_rotation_adj(params, constants, node_l)
Adjoint transfer local coefficients over a tree.
Definition: ddx_core.f90:1974
ddx_core::tree_m2m_rotation_adj
subroutine tree_m2m_rotation_adj(params, constants, node_m)
Adjoint transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1793
ddx_core::ddproject_cav
subroutine ddproject_cav(params, constants, s, xi)
Integrate by a characteristic function at Lebedev grid points \xi(n,i) = sum w_n U_n^i Y_l^m(s_n) [S_...
Definition: ddx_core.f90:1649
ddx_core::ddeval_grid
subroutine ddeval_grid(params, constants, alpha, x_sph, beta, x_grid)
Evaluate values of spherical harmonics at Lebedev grid points.
Definition: ddx_core.f90:1535
ddx_core::ddcav_to_grid_work
subroutine ddcav_to_grid_work(ngrid, nsph, ncav, icav_ia, icav_ja, x_cav, x_grid)
Unwrap values at cavity points into values at all grid points.
Definition: ddx_core.f90:1621
ddx_core::tree_m2m_bessel_rotation
subroutine tree_m2m_bessel_rotation(params, constants, node_m)
Transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1732
ddx_core::tree_m2m_bessel_rotation_adj_work
subroutine tree_m2m_bessel_rotation_adj_work(params, constants, node_m)
Adjoint transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1852
ddx_core::ddcav_to_grid
subroutine ddcav_to_grid(params, constants, x_cav, x_grid)
Unwrap values at cavity points into values at all grid points.
Definition: ddx_core.f90:1604
ddx_core::allocate_model
subroutine allocate_model(nsph, x, y, z, rvdw, model, lmax, ngrid, force, fmm, pm, pl, se, eta, eps, kappa, matvecmem, maxiter, jacobi_ndiis, nproc, output_filename, ddx_data, ddx_error)
Initialize ddX input with a full set of parameters.
Definition: ddx_core.f90:269
ddx_core::tree_m2m_rotation_adj_work
subroutine tree_m2m_rotation_adj_work(params, constants, node_m, work)
Adjoint transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1809
ddx_core::dbasis
subroutine dbasis(params, constants, x, basloc, dbsloc, vplm, vcos, vsin)
Compute first derivatives of spherical harmonics.
Definition: ddx_core.f90:1176
ddx_core::print_spherical
subroutine print_spherical(iunit, label, nbasis, lmax, ncol, icol, x)
Print array of spherical harmonics.
Definition: ddx_core.f90:1012
ddx_core::intmlp
real(dp) function intmlp(params, constants, t, sigma, basloc)
TODO.
Definition: ddx_core.f90:1307
ddx_harmonics
Harmonics-related core routines.
Definition: ddx_harmonics.f90:16
ddx_harmonics::fmm_l2p_adj_work
subroutine fmm_l2p_adj_work(c, src_q, dst_r, p, vscales_rel, beta, dst_l, work)
Adjoint L2P.
Definition: ddx_harmonics.f90:2673
ddx_harmonics::fmm_m2m_bessel_rotation_adj_work
subroutine fmm_m2m_bessel_rotation_adj_work(c, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:6014
ddx_harmonics::fmm_m2p_bessel_adj
subroutine fmm_m2p_bessel_adj(c, src_q, dst_r, kappa, p, vscales, beta, dst_m)
Adjoint M2P operation.
Definition: ddx_harmonics.f90:1704
ddx_harmonics::polleg_work
subroutine polleg_work(ctheta, stheta, p, vplm, work)
Compute associated Legendre polynomials.
Definition: ddx_harmonics.f90:483
ddx_harmonics::trgev
subroutine trgev(cphi, sphi, p, vcos, vsin)
Compute arrays of and .
Definition: ddx_harmonics.f90:349
ddx_harmonics::fmm_m2p_adj_work
subroutine fmm_m2p_adj_work(c, src_q, dst_r, p, vscales_rel, beta, dst_m, work)
Adjoint M2P operation.
Definition: ddx_harmonics.f90:1405
ddx_harmonics::fmm_l2l_bessel_rotation_adj_work
subroutine fmm_l2l_bessel_rotation_adj_work(c, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct L2L translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:7424
ddx_harmonics::fmm_m2m_bessel_rotation_work
subroutine fmm_m2m_bessel_rotation_work(c, src_sk, dst_sk, p, vscales, alpha, src_m, beta, dst_m, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:5878
ddx_harmonics::fmm_l2l_rotation_work
subroutine fmm_l2l_rotation_work(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Direct L2L translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:7156
ddx_harmonics::fmm_m2p_work
subroutine fmm_m2p_work(c, src_r, p, vscales_rel, alpha, src_m, beta, dst_v, work)
Accumulate potential, induced by multipole spherical harmonics.
Definition: ddx_harmonics.f90:893
ddx_harmonics::fmm_m2l_rotation_adj
subroutine fmm_m2l_rotation_adj(c, src_r, dst_r, pl, pm, vscales, m2l_ztranslate_adj_coef, alpha, src_l, beta, dst_m)
Adjoint M2L translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:8930
ddx_harmonics::fmm_m2m_rotation_work
subroutine fmm_m2m_rotation_work(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Direct M2M translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:5743
ddx_harmonics::fmm_sph_transform
subroutine fmm_sph_transform(p, r1, alpha, src, beta, dst)
Transform coefficients of spherical harmonics to a new cartesion system.
Definition: ddx_harmonics.f90:2999
ddx_harmonics::fmm_m2l_rotation_work
subroutine fmm_m2l_rotation_work(c, src_r, dst_r, pm, pl, vscales, m2l_ztranslate_coef, alpha, src_m, beta, dst_l, work)
Direct M2L translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:8100
ddx_harmonics::fmm_l2l_rotation_adj_work
subroutine fmm_l2l_rotation_adj_work(c, src_r, dst_r, p, vscales, vfact, alpha, src_l, beta, dst_l, work)
Adjoint L2L translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:7582
ddx_harmonics::fmm_m2l_bessel_rotation_work
subroutine fmm_m2l_bessel_rotation_work(c, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:8237
ddx_harmonics::fmm_l2l_bessel_rotation_work
subroutine fmm_l2l_bessel_rotation_work(c, src_si, dst_si, p, vscales, alpha, src_l, beta, dst_l, work, work_complex)
Direct L2L translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:7290
ddx_harmonics::fmm_m2l_bessel_rotation_adj_work
subroutine fmm_m2l_bessel_rotation_adj_work(c, src_sk, dst_si, p, vscales, alpha, src_m, beta, dst_l, work, work_complex)
Direct M2M translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:8373
ddx_harmonics::fmm_m2m_rotation_adj_work
subroutine fmm_m2m_rotation_adj_work(c, src_r, dst_r, p, vscales, vcnk, alpha, src_m, beta, dst_m, work)
Adjoint M2M translation by 4 rotations and 1 translation.
Definition: ddx_harmonics.f90:6129
ddx_harmonics::fmm_m2p_bessel_work
subroutine fmm_m2p_bessel_work(c, p, vscales, SK_ri, alpha, src_m, beta, dst_v, work_complex, work)
Accumulate Bessel potential, induced by multipole spherical harmonics.
Definition: ddx_harmonics.f90:1134
ddx_harmonics::fmm_l2p_work
subroutine fmm_l2p_work(c, src_r, p, vscales_rel, alpha, src_l, beta, dst_v, work)
Accumulate potential, induced by local spherical harmonics.
Definition: ddx_harmonics.f90:2167
ddx_parameters
Module to treat properly user input parameters.
Definition: ddx_parameters.f90:15
ddx_parameters::params_free
subroutine params_free(params, ddx_error)
Deallocate the parameter object.
Definition: ddx_parameters.f90:345
ddx_parameters::params_init
subroutine params_init(model, force, eps, kappa, eta, se, lmax, ngrid, matvecmem, maxiter, jacobi_ndiis, fmm, pm, pl, nproc, nsph, csph, rsph, output_filename, params, ddx_error)
Initialize and check input parameters.
Definition: ddx_parameters.f90:123
ddx_workspace
Workspace for temporary buffers.
Definition: ddx_workspace.f90:15
ddx_workspace::workspace_init
subroutine workspace_init(params, constants, workspace, ddx_error)
Initialize and allocate the temporary workspaces.
Definition: ddx_workspace.f90:115
ddx_workspace::workspace_free
subroutine workspace_free(workspace, ddx_error)
Deallocate the temporary workspaces.
Definition: ddx_workspace.f90:290
ddx_constants::ddx_constants_type
Container for precomputed constants.
Definition: ddx_constants.f90:28
ddx_core::ddx_electrostatics_type
Container for the electrostatic properties. Since different methods require different electrostatic p...
Definition: ddx_core.f90:203
ddx_core::ddx_state_type
This defined type contains the primal and adjoint RHSs, the solution of the primal and adjoint linear...
Definition: ddx_core.f90:33
ddx_core::ddx_type
Main ddX type that stores all the required information. Container for the params, contants and worksp...
Definition: ddx_core.f90:222
ddx_parameters::ddx_params_type
Type to check and store user input parameters.
Definition: ddx_parameters.f90:25
ddx_workspace::ddx_workspace_type
Container for temporary arrays.
Definition: ddx_workspace.f90:25