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