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