ddx 0.6.0
Libary for domain-decomposition methods for polarizable continuum models
ddx_multipolar_solutes.f90
1
2module ddx_multipolar_solutes
3
4use ddx_definitions
5use ddx_core
6use ddx_gradients
7
8implicit none
9
10contains
11
14
28subroutine multipole_electrostatics(params, constants, workspace, multipoles, &
29 & mmax, electrostatics, ddx_error)
30 implicit none
31 type(ddx_params_type), intent(in) :: params
32 type(ddx_constants_type), intent(in) :: constants
33 type(ddx_workspace_type), intent(inout) :: workspace
34 integer, intent(in) :: mmax
35 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
36 type(ddx_electrostatics_type), intent(out) :: electrostatics
37 type(ddx_error_type), intent(inout) :: ddx_error
38
39 call allocate_electrostatics(params, constants, electrostatics, ddx_error)
40 if (ddx_error % flag .ne. 0) then
41 call update_error(ddx_error, &
42 & "allocate_electrostatics returned an error, exiting")
43 return
44 end if
45
46 ! Compute the required electrostatic properties
47 if (electrostatics % do_phi .and. electrostatics % do_e &
48 & .and. electrostatics % do_g) then
49 call multipole_electrostatics_2(params, constants, workspace, &
50 & multipoles, 0, electrostatics % phi_cav, &
51 & electrostatics % e_cav, electrostatics % g_cav, ddx_error)
52 else if (electrostatics % do_phi .and. electrostatics % do_e) then
53 call multipole_electrostatics_1(params, constants, workspace, &
54 & multipoles, 0, electrostatics % phi_cav, &
55 & electrostatics % e_cav, ddx_error)
56 else
57 call multipole_electrostatics_0(params, constants, workspace, &
58 & multipoles, 0, electrostatics % phi_cav, ddx_error)
59 end if
60
61end subroutine multipole_electrostatics
62
80subroutine multipole_electrostatics_2(params, constants, workspace, multipoles, &
81 & mmax, phi_cav, e_cav, g_cav, ddx_error)
82 implicit none
83 type(ddx_params_type), intent(in) :: params
84 type(ddx_workspace_type), intent(inout) :: workspace
85 type(ddx_constants_type), intent(in) :: constants
86 type(ddx_error_type), intent(inout) :: ddx_error
87 integer, intent(in) :: mmax
88 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
89 real(dp), intent(out) :: phi_cav(constants % ncav)
90 real(dp), intent(out) :: e_cav(3, constants % ncav)
91 real(dp), intent(out) :: g_cav(3, 3, constants % ncav)
92
93 if (params % fmm .eq. 0) then
94 call build_g_dense(multipoles, params % csph, mmax, params % nsph, &
95 & phi_cav, constants % ccav, constants % ncav, e_cav, g_cav, &
96 & ddx_error)
97 else if (params % fmm .eq. 1) then
98 call build_g_fmm(params, constants, workspace, multipoles, &
99 & mmax, phi_cav, e_cav, g_cav, ddx_error)
100 end if
101end subroutine multipole_electrostatics_2
102
121subroutine build_g_dense(multipoles, cm, mmax, nm, phi_cav, ccav, ncav, &
122 & e_cav, g_cav, ddx_error)
123 implicit none
124 integer, intent(in) :: mmax, nm, ncav
125 real(dp), intent(in) :: multipoles((mmax + 1)**2, nm)
126 real(dp), intent(in) :: cm(3, nm), ccav(3, ncav)
127 real(dp), intent(out) :: phi_cav(ncav), e_cav(3, ncav), g_cav(3, 3, ncav)
128 type(ddx_error_type), intent(inout) :: ddx_error
129 real(dp), allocatable :: tmp_m_grad(:, :, :), vscales(:), &
130 & vscales_rel(:), v4pi2lp1(:), tmp_m_hess(:, :, :, :), tmp(:, :)
131 integer icav, im, info
132 real(dp) :: v, ex, ey, ez, c(3), gxx, gxy, gxz, gyy, gyz, gzz
133
134 ! allocate some space for the M2M gradients and precompute
135 ! the quantities for the m2p
136 allocate(tmp_m_grad((mmax + 2)**2, 3, nm), vscales((mmax + 3)**2), &
137 & vscales_rel((mmax + 3)**2), v4pi2lp1(mmax + 3), &
138 & tmp_m_hess((mmax + 3)**2, 3, nm, 3), tmp((mmax + 2)**2, 3), &
139 & stat=info)
140 if (info .ne. 0) then
141 call update_error(ddx_error, "Allocation failed in build_g_dense!")
142 return
143 end if
144 call ylmscale(mmax + 2, vscales, v4pi2lp1, vscales_rel)
145
146 ! call the helper routine for the M2M gradient and hessian
147 call grad_m2m(multipoles, mmax, nm, tmp_m_grad, ddx_error)
148 tmp = tmp_m_grad(:, 1, :)
149 call grad_m2m(tmp, mmax + 1, nm, tmp_m_hess(:, :, :, 1), ddx_error)
150 tmp = tmp_m_grad(:, 2, :)
151 call grad_m2m(tmp, mmax + 1, nm, tmp_m_hess(:, :, :, 2), ddx_error)
152 tmp = tmp_m_grad(:, 3, :)
153 call grad_m2m(tmp, mmax + 1, nm, tmp_m_hess(:, :, :, 3), ddx_error)
154 if (ddx_error % flag .ne. 0) then
155 call update_error(ddx_error, "grad_m2m returned an error, exiting")
156 return
157 end if
158
159 ! loop over the targets and the sources and assemble the electric
160 ! potential and field
161 do icav = 1, ncav
162 v = zero
163 ex = zero
164 ey = zero
165 ez = zero
166 gxx = zero
167 gxy = zero
168 gxz = zero
169 gyy = zero
170 gyz = zero
171 gzz = zero
172 do im = 1, nm
173 c(:) = ccav(:, icav) - cm(:, im)
174 call fmm_m2p(c, one, mmax, vscales_rel, one, multipoles(:, im), &
175 & one, v)
176
177 call fmm_m2p(c, one, mmax + 1, vscales_rel, one, &
178 & tmp_m_grad(:, 1, im), one, ex)
179 call fmm_m2p(c, one, mmax + 1, vscales_rel, one, &
180 & tmp_m_grad(:, 2, im), one, ey)
181 call fmm_m2p(c, one, mmax + 1, vscales_rel, one, &
182 & tmp_m_grad(:, 3, im), one, ez)
183
184 call fmm_m2p(c, one, mmax + 2, vscales_rel, one, &
185 & tmp_m_hess(:, 1, im, 1), one, gxx)
186 call fmm_m2p(c, one, mmax + 2, vscales_rel, one, &
187 & tmp_m_hess(:, 2, im, 1), one, gxy)
188 call fmm_m2p(c, one, mmax + 2, vscales_rel, one, &
189 & tmp_m_hess(:, 3, im, 1), one, gxz)
190
191 call fmm_m2p(c, one, mmax + 2, vscales_rel, one, &
192 & tmp_m_hess(:, 2, im, 2), one, gyy)
193 call fmm_m2p(c, one, mmax + 2, vscales_rel, one, &
194 & tmp_m_hess(:, 3, im, 2), one, gyz)
195
196 call fmm_m2p(c, one, mmax + 2, vscales_rel, one, &
197 & tmp_m_hess(:, 3, im, 3), one, gzz)
198 end do
199 phi_cav(icav) = v
200 e_cav(1, icav) = ex
201 e_cav(2, icav) = ey
202 e_cav(3, icav) = ez
203 g_cav(1, 1, icav) = gxx
204 g_cav(1, 2, icav) = gxy
205 g_cav(1, 3, icav) = gxz
206 g_cav(2, 1, icav) = gxy
207 g_cav(2, 2, icav) = gyy
208 g_cav(2, 3, icav) = gyz
209 g_cav(3, 1, icav) = gxz
210 g_cav(3, 2, icav) = gyz
211 g_cav(3, 3, icav) = gzz
212 end do
213
214 deallocate(tmp_m_grad, vscales, vscales_rel, v4pi2lp1, tmp_m_hess, &
215 & tmp, stat=info)
216 if (info .ne. 0) then
217 call update_error(ddx_error, "Deallocation failed in build_e_dense!")
218 return
219 end if
220end subroutine build_g_dense
221
237subroutine multipole_electrostatics_1(params, constants, workspace, multipoles, &
238 & mmax, phi_cav, e_cav, ddx_error)
239 implicit none
240 type(ddx_params_type), intent(in) :: params
241 type(ddx_workspace_type), intent(inout) :: workspace
242 type(ddx_constants_type), intent(in) :: constants
243 type(ddx_error_type), intent(inout) :: ddx_error
244 integer, intent(in) :: mmax
245 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
246 real(dp), intent(out) :: phi_cav(constants % ncav)
247 real(dp), intent(out) :: e_cav(3, constants % ncav)
248 if (params % fmm .eq. 0) then
249 call build_e_dense(multipoles, params % csph, mmax, params % nsph, &
250 & phi_cav, constants % ccav, constants % ncav, e_cav, ddx_error)
251 else if (params % fmm .eq. 1) then
252 call build_e_fmm(params, constants, workspace, multipoles, &
253 & mmax, phi_cav, e_cav, ddx_error)
254 end if
255end subroutine multipole_electrostatics_1
256
272subroutine build_e_dense(multipoles, cm, mmax, nm, phi_cav, ccav, ncav, &
273 & e_cav, ddx_error)
274 implicit none
275 integer, intent(in) :: mmax, nm, ncav
276 real(dp), intent(in) :: cm(3, nm), ccav(3, ncav), &
277 & multipoles((mmax + 1)**2, nm)
278 real(dp), intent(out) :: phi_cav(ncav), e_cav(3, ncav)
279 type(ddx_error_type), intent(inout) :: ddx_error
280 ! local variables
281 real(dp), allocatable :: tmp_m_grad(:, :, :), vscales(:), &
282 & vscales_rel(:), v4pi2lp1(:)
283 integer icav, im, info
284 real(dp) :: v, ex, ey, ez, c(3)
285
286 ! allocate some space for the M2M gradients and precompute
287 ! the quantities for the m2p
288 allocate(tmp_m_grad((mmax + 2)**2, 3, nm), vscales((mmax + 2)**2), &
289 & vscales_rel((mmax + 2)**2), v4pi2lp1(mmax + 2), stat=info)
290 if (info .ne. 0) then
291 call update_error(ddx_error, "Allocation failed in build_e_dense!")
292 return
293 end if
294 call ylmscale(mmax + 1, vscales, v4pi2lp1, vscales_rel)
295
296 ! call the helper routine for the M2M gradients
297 call grad_m2m(multipoles, mmax, nm, tmp_m_grad, ddx_error)
298 if (ddx_error % flag .ne. 0) then
299 call update_error(ddx_error, "grad_m2m returned an error, exiting")
300 return
301 end if
302
303
304 ! loop over the targets and the sources and assemble the electric
305 ! potential and field
306 do icav = 1, ncav
307 v = zero
308 ex = zero
309 ey = zero
310 ez = zero
311 do im = 1, nm
312 c(:) = ccav(:, icav) - cm(:, im)
313 call fmm_m2p(c, one, mmax, vscales_rel, one, multipoles(:, im), &
314 & one, v)
315
316 call fmm_m2p(c, one, mmax + 1, vscales_rel, one, &
317 & tmp_m_grad(:, 1, im), one, ex)
318 call fmm_m2p(c, one, mmax + 1, vscales_rel, one, &
319 & tmp_m_grad(:, 2, im), one, ey)
320 call fmm_m2p(c, one, mmax + 1, vscales_rel, one, &
321 & tmp_m_grad(:, 3, im), one, ez)
322 end do
323 phi_cav(icav) = v
324 e_cav(1, icav) = ex
325 e_cav(2, icav) = ey
326 e_cav(3, icav) = ez
327 end do
328
329 deallocate(tmp_m_grad, vscales, vscales_rel, v4pi2lp1, stat=info)
330 if (info .ne. 0) then
331 call update_error(ddx_error, "Deallocation failed in build_e_dense!")
332 return
333 end if
334end subroutine build_e_dense
335
348subroutine multipole_electrostatics_0(params, constants, workspace, multipoles, &
349 & mmax, phi_cav, ddx_error)
350 implicit none
351 type(ddx_params_type), intent(in) :: params
352 type(ddx_workspace_type), intent(inout) :: workspace
353 type(ddx_constants_type), intent(in) :: constants
354 type(ddx_error_type), intent(inout) :: ddx_error
355 integer, intent(in) :: mmax
356 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
357 real(dp), intent(out) :: phi_cav(constants % ncav)
358 if (params % fmm .eq. 0) then
359 call build_phi_dense(multipoles, params % csph, mmax, params % nsph, &
360 & phi_cav, constants % ccav, constants % ncav, ddx_error)
361 else if (params % fmm .eq. 1) then
362 call build_phi_fmm(params, constants, workspace, multipoles, mmax, &
363 & phi_cav)
364 end if
365end subroutine multipole_electrostatics_0
366
381subroutine build_phi_dense(multipoles, cm, mmax, nm, phi_cav, ccav, ncav, &
382 & ddx_error)
383 implicit none
384 integer, intent(in) :: mmax, nm, ncav
385 real(dp), intent(in) :: cm(3, nm), ccav(3, ncav), &
386 & multipoles((mmax + 1)**2, nm)
387 real(dp), intent(out) :: phi_cav(ncav)
388 type(ddx_error_type), intent(inout) :: ddx_error
389 ! local variables
390 integer icav, im, info
391 real(dp) :: v, c(3)
392 real(dp), allocatable :: vscales(:), vscales_rel(:), v4pi2lp1(:)
393
394 ! precompute the quantities for the m2p
395 allocate(vscales((mmax + 1)**2), vscales_rel((mmax + 1)**2), &
396 & v4pi2lp1(mmax + 1), stat=info)
397 if (info .ne. 0) then
398 call update_error(ddx_error, 'Allocation failed in build_phi_dense!')
399 return
400 end if
401 call ylmscale(mmax, vscales, v4pi2lp1, vscales_rel)
402
403 do icav = 1, ncav
404 v = zero
405 do im = 1, nm
406 c(:) = ccav(:, icav) - cm(:, im)
407 call fmm_m2p(c, one, mmax, vscales_rel, one, &
408 & multipoles(:, im), one, v)
409 end do
410 phi_cav(icav) = v
411 end do
412
413 deallocate(vscales, vscales_rel, v4pi2lp1, stat=info)
414 if (info .ne. 0) then
415 call update_error(ddx_error, 'Deallocation failed in build_phi_dense!')
416 return
417 end if
418
419end subroutine build_phi_dense
420
434subroutine build_e_fmm(params, constants, workspace, multipoles, &
435 & mmax, phi_cav, e_cav, ddx_error)
436 implicit none
437 type(ddx_params_type), intent(in) :: params
438 type(ddx_workspace_type), intent(inout) :: workspace
439 type(ddx_constants_type), intent(in) :: constants
440 type(ddx_error_type), intent(inout) :: ddx_error
441 integer, intent(in) :: mmax
442 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
443 real(dp), intent(out) :: phi_cav(constants % ncav), &
444 & e_cav(3, constants % ncav)
445 ! local variables
446 integer :: info, isph, igrid, inode, jnode, jsph, jnear, icav
447 real(dp) :: ex, ey, ez, c(3)
448 real(dp), allocatable :: tmp_m_grad(:, :, :), grid_grad(:, :, :)
449
450 allocate(tmp_m_grad((mmax + 2)**2, 3, params % nsph), &
451 & grid_grad(params % ngrid, 3, params % nsph), stat=info)
452 if (info .ne. 0) then
453 call update_error(ddx_error, "Allocation failed in build_e_fmm!")
454 return
455 end if
456
457 ! compute the gradient of the m2m trasformation
458 call grad_m2m(multipoles, mmax, params % nsph, tmp_m_grad, ddx_error)
459 if (ddx_error % flag .ne. 0) then
460 call update_error(ddx_error, "grad_m2m returned an error, exiting")
461 return
462 end if
463
464 ! copy the multipoles in the right places
465 call load_m(params, constants, workspace, multipoles, mmax)
466
467 ! perform the m2m, m2l and l2l steps
468 call do_fmm(params, constants, workspace)
469
470 ! near field potential (m2p)
471 call tree_m2p(params, constants, params % lmax, one, &
472 & workspace % tmp_sph, one, workspace % tmp_grid)
473
474 ! far field potential, each sphere at its own points (l2p)
475 call dgemm('T', 'N', params % ngrid, params % nsph, &
476 & constants % nbasis, one, constants % vgrid2, &
477 & constants % vgrid_nbasis, workspace % tmp_sph, &
478 & constants % nbasis, one, workspace % tmp_grid, &
479 & params % ngrid)
480
481 ! near field gradients
482 do isph = 1, params % nsph
483 do igrid = 1, params % ngrid
484 if(constants % ui(igrid, isph) .eq. zero) cycle
485 inode = constants % snode(isph)
486 ex = zero
487 ey = zero
488 ez = zero
489 do jnear = constants % snear(inode), constants % snear(inode+1) - 1
490 jnode = constants % near(jnear)
491 jsph = constants % order(constants % cluster(1, jnode))
492 c = params % csph(:, isph) + constants % cgrid(:, igrid) &
493 & *params % rsph(isph) - params % csph(:, jsph)
494 call fmm_m2p(c, one, mmax + 1, constants % vscales_rel, &
495 & one, tmp_m_grad(:, 1, jsph), one, ex)
496 call fmm_m2p(c, one, mmax + 1, constants % vscales_rel, &
497 & one, tmp_m_grad(:, 2, jsph), one, ey)
498 call fmm_m2p(c, one, mmax + 1, constants % vscales_rel, &
499 & one, tmp_m_grad(:, 3, jsph), one, ez)
500 end do
501 grid_grad(igrid, 1, isph) = ex
502 grid_grad(igrid, 2, isph) = ey
503 grid_grad(igrid, 3, isph) = ez
504 end do
505 end do
506
507 ! far-field FMM gradients (only if pl > 0)
508 if (params % pl .gt. 0) then
509 call tree_grad_l2l(params, constants, workspace % tmp_node_l, &
510 & workspace % tmp_sph_l_grad, workspace % tmp_sph_l)
511 call dgemm('T', 'N', params % ngrid, 3*params % nsph, &
512 & (params % pl)**2, one, constants % vgrid2, &
513 & constants % vgrid_nbasis, workspace % tmp_sph_l_grad, &
514 & (params % pl+1)**2, one, grid_grad, params % ngrid)
515 end if
516
517 ! discard the internal points
518 icav = 0
519 do isph = 1, params % nsph
520 do igrid = 1, params % ngrid
521 if(constants % ui(igrid, isph) .eq. zero) cycle
522 icav = icav + 1
523 phi_cav(icav) = workspace % tmp_grid(igrid, isph)
524 e_cav(:, icav) = grid_grad(igrid, :, isph)
525 end do
526 end do
527
528 deallocate(tmp_m_grad, grid_grad, stat=info)
529 if (info .ne. 0) then
530 call update_error(ddx_error, "Deallocation failed in build_e_fmm!")
531 return
532 end if
533
534end subroutine build_e_fmm
535
551subroutine build_g_fmm(params, constants, workspace, multipoles, &
552 & mmax, phi_cav, e_cav, g_cav, ddx_error)
553 implicit none
554 type(ddx_params_type), intent(in) :: params
555 type(ddx_workspace_type), intent(inout) :: workspace
556 type(ddx_constants_type), intent(in) :: constants
557 type(ddx_error_type), intent(inout) :: ddx_error
558 integer, intent(in) :: mmax
559 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
560 real(dp), intent(out) :: phi_cav(constants % ncav), &
561 & e_cav(3, constants % ncav), g_cav(3, 3, constants % ncav)
562 ! local variables
563 integer :: info, isph, igrid, inode, jnode, jsph, jnear, icav
564 real(dp) :: ex, ey, ez, c(3), gxx, gxy, gxz, gyy, gyz, gzz
565 real(dp), allocatable :: m_grad(:, :, :), grid_grad(:, :, :), &
566 & m_hess_x(:, :, :), m_hess_y(:, :, :), m_hess_z(:, :, :), &
567 & m_grad_component(:, :), grid_hess_x(:, :, :), &
568 & grid_hess_y(:, :, :), grid_hess_z(:, :, :), grid_hessian2(:, :, :)
569
570 allocate(m_hess_x((mmax+3)**2, 3, params % nsph), &
571 & m_hess_y((mmax+3)**2, 3, params % nsph), &
572 & m_hess_z((mmax+3)**2, 3, params % nsph), &
573 & m_grad_component((mmax+2)**2, params % nsph), &
574 & m_grad((mmax + 2)**2, 3, params % nsph), &
575 & grid_grad(params % ngrid, 3, params % nsph), &
576 & grid_hess_x(params % ngrid, 3, params % nsph), &
577 & grid_hess_y(params % ngrid, 3, params % nsph), &
578 & grid_hess_z(params % ngrid, 3, params % nsph), &
579 & grid_hessian2(params % ngrid, 3, params % nsph), &
580 & stat=info)
581 if (info .ne. 0) then
582 call update_error(ddx_error, "Allocation failed in build_g_fmm!")
583 return
584 end if
585
586 ! compute the gradient of the m2m trasformation
587 call grad_m2m(multipoles, mmax, params % nsph, m_grad, ddx_error)
588
589 ! starting from the previously computed gradient, compute the
590 ! hessian of the target multipoles
591 m_grad_component = m_grad(:, 1, :)
592 call grad_m2m(m_grad_component, mmax + 1, params % nsph, m_hess_x, ddx_error)
593 m_grad_component = m_grad(:, 2, :)
594 call grad_m2m(m_grad_component, mmax + 1, params % nsph, m_hess_y, ddx_error)
595 m_grad_component = m_grad(:, 3, :)
596 call grad_m2m(m_grad_component, mmax + 1, params % nsph, m_hess_z, ddx_error)
597 if (ddx_error % flag .ne. 0) then
598 call update_error(ddx_error, "grad_m2m returned an error, exiting")
599 return
600 end if
601
602
603 ! copy the multipoles to the right places
604 call load_m(params, constants, workspace, multipoles, mmax)
605
606 ! perform the m2m, m2l and l2l steps
607 call do_fmm(params, constants, workspace)
608
609 ! near field potential (m2p)
610 call tree_m2p(params, constants, params % lmax, one, &
611 & workspace % tmp_sph, one, workspace % tmp_grid)
612
613 ! far field potential, each sphere at its own points (l2p)
614 call dgemm('T', 'N', params % ngrid, params % nsph, &
615 & constants % nbasis, one, constants % vgrid2, &
616 & constants % vgrid_nbasis, workspace % tmp_sph, &
617 & constants % nbasis, one, workspace % tmp_grid, &
618 & params % ngrid)
619
620 ! near field gradients
621 do isph = 1, params % nsph
622 do igrid = 1, params % ngrid
623 if(constants % ui(igrid, isph) .eq. zero) cycle
624 inode = constants % snode(isph)
625 ex = zero
626 ey = zero
627 ez = zero
628 do jnear = constants % snear(inode), constants % snear(inode+1) - 1
629 jnode = constants % near(jnear)
630 jsph = constants % order(constants % cluster(1, jnode))
631 c = params % csph(:, isph) + constants % cgrid(:, igrid) &
632 & *params % rsph(isph) - params % csph(:, jsph)
633 call fmm_m2p(c, one, mmax + 1, constants % vscales_rel, &
634 & one, m_grad(:, 1, jsph), one, ex)
635 call fmm_m2p(c, one, mmax + 1, constants % vscales_rel, &
636 & one, m_grad(:, 2, jsph), one, ey)
637 call fmm_m2p(c, one, mmax + 1, constants % vscales_rel, &
638 & one, m_grad(:, 3, jsph), one, ez)
639 end do
640 grid_grad(igrid, 1, isph) = ex
641 grid_grad(igrid, 2, isph) = ey
642 grid_grad(igrid, 3, isph) = ez
643 end do
644 end do
645
646 ! near field hessians
647 do isph = 1, params % nsph
648 do igrid = 1, params % ngrid
649 if(constants % ui(igrid, isph) .eq. zero) cycle
650 inode = constants % snode(isph)
651 gxx = zero
652 gxy = zero
653 gxz = zero
654 gyy = zero
655 gyz = zero
656 gzz = zero
657 do jnear = constants % snear(inode), constants % snear(inode+1) - 1
658 jnode = constants % near(jnear)
659 jsph = constants % order(constants % cluster(1, jnode))
660 c = params % csph(:, isph) + constants % cgrid(:, igrid) &
661 & *params % rsph(isph) - params % csph(:, jsph)
662 call fmm_m2p(c, one, mmax + 2, constants % vscales_rel, &
663 & one, m_hess_x(:, 1, jsph), one, gxx)
664 call fmm_m2p(c, one, mmax + 2, constants % vscales_rel, &
665 & one, m_hess_x(:, 2, jsph), one, gxy)
666 call fmm_m2p(c, one, mmax + 2, constants % vscales_rel, &
667 & one, m_hess_x(:, 3, jsph), one, gxz)
668 call fmm_m2p(c, one, mmax + 2, constants % vscales_rel, &
669 & one, m_hess_y(:, 2, jsph), one, gyy)
670 call fmm_m2p(c, one, mmax + 2, constants % vscales_rel, &
671 & one, m_hess_y(:, 3, jsph), one, gyz)
672 call fmm_m2p(c, one, mmax + 2, constants % vscales_rel, &
673 & one, m_hess_z(:, 3, jsph), one, gzz)
674 end do
675 grid_hess_x(igrid, 1, isph) = gxx
676 grid_hess_x(igrid, 2, isph) = gxy
677 grid_hess_x(igrid, 3, isph) = gxz
678 grid_hess_y(igrid, 1, isph) = gxy
679 grid_hess_y(igrid, 2, isph) = gyy
680 grid_hess_y(igrid, 3, isph) = gyz
681 grid_hess_z(igrid, 1, isph) = gxz
682 grid_hess_z(igrid, 2, isph) = gyz
683 grid_hess_z(igrid, 3, isph) = gzz
684 end do
685 end do
686
687 ! far-field FMM gradients (only if pl > 0)
688 if (params % pl .gt. 0) then
689 call tree_grad_l2l(params, constants, workspace % tmp_node_l, &
690 & workspace % tmp_sph_l_grad, workspace % tmp_sph_l)
691 call dgemm('T', 'N', params % ngrid, 3*params % nsph, &
692 & (params % pl)**2, one, constants % vgrid2, &
693 & constants % vgrid_nbasis, workspace % tmp_sph_l_grad, &
694 & (params % pl+1)**2, one, grid_grad, params % ngrid)
695 call tree_grad_l2l(params, constants, workspace % tmp_node_l, &
696 & workspace % tmp_sph_l_grad, workspace % tmp_sph_l)
697 end if
698
699 ! Far field FMM hessians
700 if (params % pl .gt. 1) then
701 ! Load previously computed gradient into leaves, since
702 ! tree_grad_l2l currently takes local expansions of entire
703 ! tree. In future it might be changed.
704 do isph = 1, params % nsph
705 inode = constants % snode(isph)
706 workspace % tmp_node_l(:, inode) = &
707 & workspace % tmp_sph_l_grad(:, 1, isph)
708 end do
709 ! Get gradient of a gradient of L2L. Currently this uses input
710 ! pl maximal degree of local harmonics but in reality we need
711 ! only pl-1 maximal degree since coefficients of harmonics of a
712 ! degree pl are zeros.
713 call tree_grad_l2l(params, constants, workspace % tmp_node_l, &
714 & workspace % tmp_sph_l_grad2, workspace % tmp_sph_l)
715 ! Apply L2P for every axis
716 call dgemm('T', 'N', params % ngrid, 3*params % nsph, &
717 & (params % pl-1)**2, one, constants % vgrid2, &
718 & constants % vgrid_nbasis, workspace % tmp_sph_l_grad2, &
719 & (params % pl+1)**2, zero, grid_hessian2, params % ngrid)
720 ! Properly copy hessian
721 grid_hess_x = grid_hess_x + grid_hessian2
722
723 ! same for y
724 do isph = 1, params % nsph
725 inode = constants % snode(isph)
726 workspace % tmp_node_l(:, inode) = &
727 & workspace % tmp_sph_l_grad(:, 2, isph)
728 end do
729 call tree_grad_l2l(params, constants, workspace % tmp_node_l, &
730 & workspace % tmp_sph_l_grad2, workspace % tmp_sph_l)
731 call dgemm('T', 'N', params % ngrid, 3*params % nsph, &
732 & (params % pl-1)**2, one, constants % vgrid2, &
733 & constants % vgrid_nbasis, workspace % tmp_sph_l_grad2, &
734 & (params % pl+1)**2, zero, grid_hessian2, params % ngrid)
735 grid_hess_y = grid_hess_y + grid_hessian2
736
737 ! same for z
738 do isph = 1, params % nsph
739 inode = constants % snode(isph)
740 workspace % tmp_node_l(:, inode) = &
741 & workspace % tmp_sph_l_grad(:, 3, isph)
742 end do
743 call tree_grad_l2l(params, constants, workspace % tmp_node_l, &
744 & workspace % tmp_sph_l_grad2, workspace % tmp_sph_l)
745 call dgemm('T', 'N', params % ngrid, 3*params % nsph, &
746 & (params % pl-1)**2, one, constants % vgrid2, &
747 & constants % vgrid_nbasis, workspace % tmp_sph_l_grad2, &
748 & (params % pl+1)**2, zero, grid_hessian2, params % ngrid)
749 grid_hess_z = grid_hess_z + grid_hessian2
750
751 end if
752
753 ! discard the internal points
754 icav = 0
755 do isph = 1, params % nsph
756 do igrid = 1, params % ngrid
757 if(constants % ui(igrid, isph) .eq. zero) cycle
758 icav = icav + 1
759 phi_cav(icav) = workspace % tmp_grid(igrid, isph)
760 e_cav(:, icav) = grid_grad(igrid, :, isph)
761 g_cav(1, :, icav) = grid_hess_x(igrid, :, isph)
762 g_cav(2, :, icav) = grid_hess_y(igrid, :, isph)
763 g_cav(3, :, icav) = grid_hess_z(igrid, :, isph)
764 end do
765 end do
766
767 deallocate(m_hess_x, m_hess_y, m_hess_z, m_grad_component, m_grad, &
768 & grid_grad, grid_hess_x, grid_hess_y, grid_hess_z, grid_hessian2, &
769 & stat=info)
770 if (info .ne. 0) then
771 call update_error(ddx_error, "Deallocation failed in build_g_fmm!")
772 return
773 end if
774
775end subroutine build_g_fmm
776
787subroutine build_phi_fmm(params, constants, workspace, multipoles, mmax, &
788 & phi_cav)
789 implicit none
790 type(ddx_params_type), intent(in) :: params
791 type(ddx_workspace_type), intent(inout) :: workspace
792 type(ddx_constants_type), intent(in) :: constants
793 integer, intent(in) :: mmax
794 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
795 real(dp), intent(out) :: phi_cav(constants % ncav)
796 ! local variables
797 integer isph, igrid, icav
798
799 ! copy the multipoles in the right places
800 call load_m(params, constants, workspace, multipoles, mmax)
801
802 ! perform the m2m, m2l and l2l steps
803 call do_fmm(params, constants, workspace)
804
805 ! near field
806 call tree_m2p(params, constants, params % lmax, one, &
807 & workspace % tmp_sph, one, workspace % tmp_grid)
808
809 ! Potential from each sphere to its own grid points (l2p)
810 call dgemm('T', 'N', params % ngrid, params % nsph, &
811 & constants % nbasis, one, constants % vgrid2, &
812 & constants % vgrid_nbasis, workspace % tmp_sph, &
813 & constants % nbasis, one, workspace % tmp_grid, &
814 & params % ngrid)
815
816 ! discard the internal points
817 icav = 0
818 do isph = 1, params % nsph
819 do igrid = 1, params % ngrid
820 if(constants % ui(igrid, isph) .eq. zero) cycle
821 icav = icav + 1
822 phi_cav(icav) = workspace % tmp_grid(igrid, isph)
823 end do
824 end do
825
826end subroutine build_phi_fmm
827
838subroutine multipole_psi(params, multipoles, mmax, psi)
839 implicit none
840 type(ddx_params_type), intent(in) :: params
841 integer, intent(in) :: mmax
842 real(dp), intent(in) :: multipoles((mmax+1)**2, params % nsph)
843 real(dp), intent(out) :: psi((params % lmax+1)**2, params % nsph)
844 integer :: isph, l, m, i
845 real(dp) :: v
846
847 psi = zero
848 do isph = 1, params % nsph
849 do l = 0, mmax
850 v = fourpi/((two*dble(l) + one)*(params % rsph(isph)**(l)))
851 i = l*l + l + 1
852 do m = -l, l
853 psi(i + m, isph) = v*multipoles(i + m, isph)
854 end do
855 end do
856 end do
857end subroutine multipole_psi
858
868subroutine load_m(params, constants, workspace, multipoles, mmax)
869 implicit none
870 type(ddx_params_type), intent(in) :: params
871 type(ddx_workspace_type), intent(inout) :: workspace
872 type(ddx_constants_type), intent(in) :: constants
873 integer, intent(in) :: mmax
874 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
875 integer :: isph, l, ind, m, inode, tmp_max
876
877 tmp_max = min(params % pm, mmax)
878
879 ! P2M is not needed as we have already multipolar distributions
880 ! we just have to copy the multipoles in the proper places
881 do isph = 1, params % nsph
882 inode = constants % snode(isph)
883 workspace % tmp_sph(:, isph) = zero
884 workspace % tmp_node_m(:, inode) = zero
885 do l = 0, tmp_max
886 ind = l*l + l + 1
887 do m = -l, l
888 workspace % tmp_sph(ind+m, isph) = &
889 & multipoles(ind+m, isph)/(params%rsph(isph)**(l+1))
890 workspace % tmp_node_m(ind+m, inode) = &
891 & multipoles(ind+m, isph)/(params%rsph(isph)**(l+1))
892 end do
893 end do
894 end do
895
896end subroutine load_m
897
905subroutine do_fmm(params, constants, workspace)
906 implicit none
907 type(ddx_params_type), intent(in) :: params
908 type(ddx_workspace_type), intent(inout) :: workspace
909 type(ddx_constants_type), intent(in) :: constants
910 call tree_m2m_rotation(params, constants, workspace % tmp_node_m)
911 call tree_m2l_rotation(params, constants, workspace % tmp_node_m, &
912 & workspace % tmp_node_l)
913 call tree_l2l_rotation(params, constants, workspace % tmp_node_l)
914 call tree_l2p(params, constants, one, workspace % tmp_node_l, zero, &
915 & workspace % tmp_grid, workspace % tmp_sph_l)
916end subroutine do_fmm
917
928subroutine grad_m2m(multipoles, mmax, nm, tmp_m_grad, ddx_error)
929 implicit none
930 integer, intent(in) :: mmax, nm
931 real(dp), intent(in) :: multipoles((mmax + 1)**2, nm)
932 real(dp), intent(out) :: tmp_m_grad((mmax + 2)**2, 3, nm)
933 type(ddx_error_type), intent(inout) :: ddx_error
934 ! local variables
935 real(dp), dimension(3, 3) :: zx_coord_transform, zy_coord_transform
936 real(dp), allocatable :: tmp(:, :)
937 integer :: info, im, l, indi, indj, m
938 real(dp) :: tmp1, tmp2
939
940 allocate(tmp((mmax + 2)**2, nm), stat=info)
941 if (info .ne. 0) then
942 call update_error(ddx_error, "Allocation failed in grad_m2m!")
943 return
944 end if
945
946 zx_coord_transform = zero
947 zx_coord_transform(3, 2) = one
948 zx_coord_transform(2, 3) = one
949 zx_coord_transform(1, 1) = one
950 zy_coord_transform = zero
951 zy_coord_transform(1, 2) = one
952 zy_coord_transform(2, 1) = one
953 zy_coord_transform(3, 3) = one
954
955 tmp_m_grad(1:(mmax + 1)**2, 3, :) = multipoles
956 do im = 1, nm
957 call fmm_sph_transform(mmax, zx_coord_transform, one, &
958 & multipoles(:, im), zero, tmp_m_grad(1:(mmax + 1)**2, 1, im))
959 call fmm_sph_transform(mmax, zy_coord_transform, one, &
960 & multipoles(:, im), zero, tmp_m_grad(1:(mmax + 1)**2, 2, im))
961 end do
962
963 do l = mmax+1, 1, -1
964 indi = l*l + l + 1
965 indj = indi - 2*l
966 tmp1 = sqrt(dble(2*l+1)) / sqrt(dble(2*l-1))
967 do m = 1-l, l-1
968 tmp2 = sqrt(dble(l*l-m*m)) * tmp1
969 tmp_m_grad(indi+m, :, :) = tmp2 * tmp_m_grad(indj+m, :, :)
970 end do
971 tmp_m_grad(indi+l, :, :) = zero
972 tmp_m_grad(indi-l, :, :) = zero
973 end do
974 tmp_m_grad(1, :, :) = zero
975
976 do im = 1, nm
977 tmp_m_grad(:, 3, im) = tmp_m_grad(:, 3, im)
978 tmp(:, im) = tmp_m_grad(:, 1, im)
979 call fmm_sph_transform(mmax+1, zx_coord_transform, one, &
980 & tmp(:, im), zero, tmp_m_grad(:, 1, im))
981 tmp(:, im) = tmp_m_grad(:, 2, im)
982 call fmm_sph_transform(mmax+1, zy_coord_transform, one, &
983 & tmp(:, im), zero, tmp_m_grad(:, 2, im))
984 end do
985
986 deallocate(tmp, stat=info)
987 if (info .ne. 0) then
988 call update_error(ddx_error, "Deallocation failed in grad_m2m!")
989 return
990 end if
991
992end subroutine grad_m2m
993
1003subroutine grad_phi_for_charges(params, constants, workspace, state, &
1004 & charges, forces, ddx_error)
1005 implicit none
1006 type(ddx_params_type), intent(in) :: params
1007 type(ddx_workspace_type), intent(inout) :: workspace
1008 type(ddx_constants_type), intent(in) :: constants
1009 type(ddx_state_type), intent(inout) :: state
1010 type(ddx_error_type), intent(inout) :: ddx_error
1011 real(dp), intent(in) :: charges(params % nsph)
1012 real(dp), intent(inout) :: forces(3, params % nsph)
1013 ! local variables
1014 integer :: info
1015 real(dp), allocatable :: multipoles(:, :)
1016
1017 allocate(multipoles(1, params % nsph), stat=info)
1018 if (info .ne. 0) then
1019 call update_error(ddx_error, "Allocation failed in grad_phi_for_charges")
1020 return
1021 end if
1022 ! convert the charges to multipoles
1023 multipoles(1, :) = charges/sqrt4pi
1024 call grad_phi(params, constants, workspace, state, 0, multipoles, forces, &
1025 & ddx_error)
1026 deallocate(multipoles, stat=info)
1027 if (info .ne. 0) then
1028 call update_error(ddx_error, "Deallocation failed in grad_phi_for_charges")
1029 return
1030 end if
1031
1032end subroutine grad_phi_for_charges
1033
1047subroutine grad_phi(params, constants, workspace, state, mmax, &
1048 & multipoles, forces, ddx_error)
1049 implicit none
1050 type(ddx_params_type), intent(in) :: params
1051 type(ddx_workspace_type), intent(inout) :: workspace
1052 type(ddx_constants_type), intent(in) :: constants
1053 type(ddx_state_type), intent(inout) :: state
1054 type(ddx_error_type), intent(inout) :: ddx_error
1055 integer, intent(in) :: mmax
1056 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
1057 real(dp), intent(inout) :: forces(3, params % nsph)
1058 ! local variables
1059 integer :: info, im, lm, ind, l, m
1060 real(dp), allocatable :: adj_phi(:, :), m_grad(:, :, :)
1061 real(dp) :: fac
1062
1063 ! get some space for the adjoint potential, note that we need it
1064 ! up to mmax + 1 as we are doing derivatives
1065 allocate(adj_phi((mmax + 2)**2, params % nsph), &
1066 & m_grad((mmax + 2)**2, 3, params % nsph), stat=info)
1067 if (info .ne. 0) then
1068 call update_error(ddx_error, "Allocation failed in grad_phi!")
1069 return
1070 end if
1071
1072 ! build the adjoint potential
1073 call build_adj_phi(params, constants, workspace, state % zeta, &
1074 & mmax + 1, adj_phi, ddx_error)
1075 if (ddx_error % flag .ne. 0) then
1076 call update_error(ddx_error, "build_adj_phi returned an error, exiting")
1077 return
1078 end if
1079
1080 ! build the gradient of the M2M transformation
1081 call grad_m2m(multipoles, mmax, params % nsph, m_grad, ddx_error)
1082 if (ddx_error % flag .ne. 0) then
1083 call update_error(ddx_error, "grad_m2m returned an error, exiting")
1084 return
1085 end if
1086
1087 ! contract the two ingredients to build the second contribution
1088 do im = 1, params % nsph
1089 do l = 1, mmax + 1
1090 ind = l*l + l + 1
1091 fac = (-one)**(l+1)
1092 do m = -l, l
1093 lm = ind + m
1094 forces(1, im) = forces(1, im) &
1095 & + fac*pt5*m_grad(lm, 1, im)*adj_phi(lm, im)
1096 forces(2, im) = forces(2, im) &
1097 & + fac*pt5*m_grad(lm, 2, im)*adj_phi(lm, im)
1098 forces(3, im) = forces(3, im) &
1099 & + fac*pt5*m_grad(lm, 3, im)*adj_phi(lm, im)
1100 end do
1101 end do
1102 end do
1103
1104 deallocate(adj_phi, m_grad, stat=info)
1105 if (info .ne. 0) then
1106 call update_error(ddx_error, "Deallocation failed in grad_phi!")
1107 return
1108 end if
1109
1110end subroutine grad_phi
1111
1125subroutine build_adj_phi(params, constants, workspace, charges, mmax, &
1126 & adj_phi, ddx_error)
1127 implicit none
1128 type(ddx_params_type), intent(in) :: params
1129 type(ddx_workspace_type), intent(inout) :: workspace
1130 type(ddx_constants_type), intent(in) :: constants
1131 type(ddx_error_type), intent(inout) :: ddx_error
1132 integer, intent(in) :: mmax
1133 real(dp), intent(in) :: charges(constants % ncav)
1134 real(dp), intent(out) :: adj_phi((mmax + 1)**2, params % nsph)
1135 if (params % fmm .eq. 0) then
1136 call build_adj_phi_dense(charges, constants % ccav, constants % ncav, &
1137 & params % csph, mmax, params % nsph, adj_phi, ddx_error)
1138 else if (params % fmm .eq. 1) then
1139 call build_adj_phi_fmm(params, constants, workspace, charges, mmax, &
1140 & adj_phi, ddx_error)
1141 end if
1142end subroutine build_adj_phi
1143
1157subroutine build_adj_phi_dense(qcav, ccav, ncav, cm, mmax, nm, adj_phi, &
1158 & ddx_error)
1159 implicit none
1160 integer, intent (in) :: ncav, mmax, nm
1161 real(dp), intent(in) :: qcav(ncav), ccav(3, ncav), cm(3, nm)
1162 real(dp), intent(out) :: adj_phi((mmax + 1)**2, nm)
1163 type(ddx_error_type), intent(inout) :: ddx_error
1164 ! local variables
1165 integer :: icav, im, info
1166 real(dp) :: c(3)
1167 real(dp), allocatable :: vscales(:), vscales_rel(:), v4pi2lp1(:), &
1168 & work(:)
1169
1170 allocate(vscales((mmax + 1)**2), vscales_rel((mmax + 1)**2), &
1171 & v4pi2lp1(mmax + 1), work((mmax + 1)**2 + 3*mmax), stat=info)
1172 if (info .ne. 0) then
1173 call update_error(ddx_error, 'Allocation failed in build_adj_phi_dense!')
1174 return
1175 end if
1176 call ylmscale(mmax, vscales, v4pi2lp1, vscales_rel)
1177
1178 do im = 1, nm
1179 adj_phi(:, im) = zero
1180 do icav = 1, ncav
1181 c(:) = cm(:, im) - ccav(:, icav)
1182 call fmm_m2p_adj_work(c, qcav(icav), one, mmax, vscales_rel, &
1183 & one, adj_phi(:, im), work)
1184 end do
1185 end do
1186
1187 deallocate(vscales, vscales_rel, v4pi2lp1, work, stat=info)
1188 if (info .ne. 0) then
1189 call update_error(ddx_error, 'Deallocation failed in build_adj_phi_dense!')
1190 return
1191 end if
1192
1193end subroutine build_adj_phi_dense
1194
1208subroutine build_adj_phi_fmm(params, constants, workspace, charges, mmax, &
1209 & adj_phi, ddx_error)
1210 implicit none
1211 type(ddx_params_type), intent(in) :: params
1212 type(ddx_workspace_type), intent(inout) :: workspace
1213 type(ddx_constants_type), intent(in) :: constants
1214 type(ddx_error_type), intent(inout) :: ddx_error
1215 integer, intent(in) :: mmax
1216 real(dp), intent(in) :: charges(constants % ncav)
1217 real(dp), intent(out) :: adj_phi((mmax + 1)**2, params % nsph)
1218 ! local variables
1219 real(dp), allocatable :: tmp_grid(:, :), sph_m(:, :), work(:)
1220 integer :: info, icav, isph, igrid, inode, indl, l, jnear, &
1221 & jnode, jsph, m, ind
1222 real(dp) :: c(3), fac
1223
1224 allocate(tmp_grid(params % ngrid, params % nsph), &
1225 & sph_m((mmax + 1)**2, params % nsph), &
1226 & work((mmax + 1)**2 + 3*mmax), stat=info)
1227 if (info .ne. 0) then
1228 call update_error(ddx_error, "Allocation failed in build_adj_phi_fmm!")
1229 return
1230 end if
1231
1232 tmp_grid = zero
1233
1234 ! expand the input from cavity points to all the Lebedev points
1235 icav = 0
1236 do isph = 1, params % nsph
1237 do igrid = 1, params % ngrid
1238 if (constants % ui(igrid, isph) .eq. zero) cycle
1239 icav = icav + 1
1240 tmp_grid(igrid, isph) = charges(icav)
1241 end do
1242 end do
1243 adj_phi = zero
1244
1245 ! near field. We do not call tree_m2p_adj as this routine skips the
1246 ! isph = jsph case which instead is needed.
1247 do isph = 1, params % nsph
1248 inode = constants % snode(isph)
1249 do jnear = constants % snear(inode), constants % snear(inode+1)-1
1250 jnode = constants % near(jnear)
1251 jsph = constants % order(constants % cluster(1, jnode))
1252 do igrid = 1, params % ngrid
1253 c = constants % cgrid(:, igrid)*params % rsph(isph) - &
1254 & params % csph(:, jsph) + params % csph(:, isph)
1255 call fmm_m2p_adj_work(c, tmp_grid(igrid, isph), &
1256 & params % rsph(jsph), mmax, constants % vscales_rel, one, &
1257 & adj_phi(:, jsph), work)
1258 end do
1259 end do
1260 end do
1261
1262 ! compute the far field and store it in tmp_node_m
1263 call tree_l2p_adj(params, constants, one, tmp_grid, zero, &
1264 & workspace % tmp_node_l, workspace % tmp_sph_l)
1265 call tree_l2l_rotation_adj(params, constants, workspace % tmp_node_l)
1266 call tree_m2l_rotation_adj(params, constants, workspace % tmp_node_l, &
1267 & workspace % tmp_node_m)
1268 call tree_m2m_rotation_adj(params, constants, workspace % tmp_node_m)
1269
1270 ! move the far field from tmp_node_m to adj_phi
1271 if(mmax .lt. params % pm) then
1272 do isph = 1, params % nsph
1273 inode = constants % snode(isph)
1274 adj_phi(:, isph) = adj_phi(:, isph) &
1275 & + workspace % tmp_node_m(1:(mmax + 1)**2, inode)
1276 end do
1277 else
1278 indl = (params % pm+1)**2
1279 do isph = 1, params % nsph
1280 inode = constants % snode(isph)
1281 adj_phi(1:indl, isph) = adj_phi(1:indl, isph) &
1282 & + workspace % tmp_node_m(:, inode)
1283 end do
1284 end if
1285
1286 ! scale by a factor to make it consistent with the non FMMized adjoint
1287 ! potential
1288 do isph = 1, params % nsph
1289 do l = 0, mmax
1290 ind = l*l + l + 1
1291 fac = - (-one/params % rsph(isph))**(l + 1)
1292 do m = -l, l
1293 adj_phi(ind + m, isph) = fac*adj_phi(ind + m, isph)
1294 end do
1295 end do
1296 end do
1297
1298 deallocate(tmp_grid, work, stat=info)
1299 if (info .ne. 0) then
1300 call update_error(ddx_error, "Deallocation failed in build_adj_phi_fmm")
1301 return
1302 end if
1303
1304end subroutine build_adj_phi_fmm
1305
1317subroutine grad_e_for_charges(params, constants, workspace, state, &
1318 & charges, force, ddx_error)
1319 implicit none
1320 type(ddx_params_type), intent(in) :: params
1321 type(ddx_workspace_type), intent(inout) :: workspace
1322 type(ddx_constants_type), intent(in) :: constants
1323 type(ddx_state_type), intent(inout) :: state
1324 type(ddx_error_type), intent(inout) :: ddx_error
1325 real(dp), intent(in) :: charges(params % nsph)
1326 real(dp), intent(inout) :: force(3, params % nsph)
1327 real(dp), allocatable :: multipoles(:, :)
1328 integer :: info
1329
1330 allocate(multipoles(1, params % nsph), stat=info)
1331 if (info .ne. 0) then
1332 call update_error(ddx_error, "Allocation failed in grad_e_for_charges")
1333 return
1334 end if
1335 ! convert the charges to multipoles
1336 multipoles(1, :) = charges/sqrt4pi
1337 call grad_e(params, constants, workspace, state, 0, multipoles, force, &
1338 & ddx_error)
1339 deallocate(multipoles, stat=info)
1340 if (info .ne. 0) then
1341 call update_error(ddx_error, "Deallocation failed in grad_e_for_charges")
1342 return
1343 end if
1344
1345end subroutine grad_e_for_charges
1346
1360subroutine grad_e(params, constants, workspace, state, mmax, &
1361 & multipoles, force, ddx_error)
1362 implicit none
1363 type(ddx_params_type), intent(in) :: params
1364 type(ddx_workspace_type), intent(inout) :: workspace
1365 type(ddx_constants_type), intent(in) :: constants
1366 type(ddx_state_type), intent(inout) :: state
1367 type(ddx_error_type), intent(inout) :: ddx_error
1368 integer, intent(in) :: mmax
1369 real(dp), intent(in) :: multipoles((mmax+1)**2, params % nsph)
1370 real(dp), intent(inout) :: force(3, params % nsph)
1371
1372 real(dp), allocatable :: zeta_comp(:), m_grad(:, :, :), &
1373 & m_hess_x(:, :, :), m_hess_y(:, :, :), m_hess_z(:, :, :), &
1374 & m_grad_component(:, :), adj_phi_x(:, :), adj_phi_y(:, :), &
1375 & adj_phi_z(:, :)
1376 integer :: info, isph, ind, l, m, lm
1377 real(dp) :: fac
1378
1379 allocate(zeta_comp(constants % ncav), &
1380 & adj_phi_x((mmax+3)**2, params % nsph), &
1381 & adj_phi_y((mmax+3)**2, params % nsph), &
1382 & adj_phi_z((mmax+3)**2, params % nsph), &
1383 & m_hess_x((mmax+3)**2, 3, params % nsph), &
1384 & m_hess_y((mmax+3)**2, 3, params % nsph), &
1385 & m_hess_z((mmax+3)**2, 3, params % nsph), &
1386 & m_grad((mmax+2)**2, 3, params % nsph), &
1387 & m_grad_component((mmax+2)**2, params % nsph), &
1388 & stat=info)
1389 if (info.ne.0) then
1390 call update_error(ddx_error, "Allocation failed in grad_e")
1391 return
1392 end if
1393
1394 ! compute the gradient of the target charges
1395 call grad_m2m(multipoles, mmax, params % nsph, m_grad, ddx_error)
1396
1397 ! starting from the previously computed gradient, compute the
1398 ! hessian of the target multipoles
1399 m_grad_component = m_grad(:, 1, :)
1400 call grad_m2m(m_grad_component, mmax + 1, params % nsph, m_hess_x, ddx_error)
1401 m_grad_component = m_grad(:, 2, :)
1402 call grad_m2m(m_grad_component, mmax + 1, params % nsph, m_hess_y, ddx_error)
1403 m_grad_component = m_grad(:, 3, :)
1404 call grad_m2m(m_grad_component, mmax + 1, params % nsph, m_hess_z, ddx_error)
1405 if (ddx_error % flag .ne. 0) then
1406 call update_error(ddx_error, "grad_m2m returned an error, exiting")
1407 return
1408 end if
1409
1410
1411 ! now we compute the adjoint potential of the x, y and z components
1412 ! of the input dipoles separately
1413 zeta_comp = state % zeta_dip(1, :)
1414 call build_adj_phi(params, constants, workspace, zeta_comp, &
1415 & mmax + 2, adj_phi_x, ddx_error)
1416 zeta_comp = state % zeta_dip(2, :)
1417 call build_adj_phi(params, constants, workspace, zeta_comp, &
1418 & mmax + 2, adj_phi_y, ddx_error)
1419 zeta_comp = state % zeta_dip(3, :)
1420 call build_adj_phi(params, constants, workspace, zeta_comp, &
1421 & mmax + 2, adj_phi_z, ddx_error)
1422 if (ddx_error % flag .ne. 0) then
1423 call update_error(ddx_error, "build_adj_phi returned an error, exiting")
1424 return
1425 end if
1426
1427
1428 ! contract the hessian of the multipoles with the three adjoint
1429 ! potentials to get the forces contributions
1430 do isph = 1, params % nsph
1431 do l = 2, mmax + 2
1432 ind = l*l + l + 1
1433 fac = -(-one)**(l+1)/2.0d0
1434 do m = -l, l
1435 lm = ind + m
1436 force(1, isph) = force(1, isph) &
1437 & + fac*m_hess_x(lm, 1, isph)*adj_phi_x(lm, isph) &
1438 & + fac*m_hess_y(lm, 1, isph)*adj_phi_y(lm, isph) &
1439 & + fac*m_hess_z(lm, 1, isph)*adj_phi_z(lm, isph)
1440 force(2, isph) = force(2, isph) &
1441 & + fac*m_hess_x(lm, 2, isph)*adj_phi_x(lm, isph) &
1442 & + fac*m_hess_y(lm, 2, isph)*adj_phi_y(lm, isph) &
1443 & + fac*m_hess_z(lm, 2, isph)*adj_phi_z(lm, isph)
1444 force(3, isph) = force(3, isph) &
1445 & + fac*m_hess_x(lm, 3, isph)*adj_phi_x(lm, isph) &
1446 & + fac*m_hess_y(lm, 3, isph)*adj_phi_y(lm, isph) &
1447 & + fac*m_hess_z(lm, 3, isph)*adj_phi_z(lm, isph)
1448 end do
1449 end do
1450 end do
1451
1452 deallocate(zeta_comp, adj_phi_x, adj_phi_y, adj_phi_z, m_hess_x, &
1453 & m_hess_y, m_hess_z, m_grad, m_grad_component, stat=info)
1454 if (info.ne.0) then
1455 call update_error(ddx_error, "Allocation failed in grad_e")
1456 return
1457 end if
1458
1459end subroutine grad_e
1460
1475subroutine multipole_force_terms(params, constants, workspace, state, mmax, &
1476 & multipoles, forces, ddx_error)
1477 implicit none
1478 type(ddx_params_type), intent(in) :: params
1479 type(ddx_workspace_type), intent(inout) :: workspace
1480 type(ddx_constants_type), intent(in) :: constants
1481 type(ddx_state_type), intent(inout) :: state
1482 type(ddx_error_type), intent(inout) :: ddx_error
1483 integer, intent(in) :: mmax
1484 real(dp), intent(in) :: multipoles((mmax + 1)**2, params % nsph)
1485 real(dp), intent(inout) :: forces(3, params % nsph)
1486
1487 call grad_phi(params, constants, workspace, state, mmax, multipoles, &
1488 & forces, ddx_error)
1489 if (ddx_error % flag .ne. 0) then
1490 call update_error(ddx_error, &
1491 & "multipole_force_terms: grad_phi returned an error, exiting")
1492 return
1493 end if
1494
1495 if (params % model .eq. 3) then
1496 call grad_e(params, constants, workspace, state, mmax, multipoles, &
1497 & forces, ddx_error)
1498 if (ddx_error % flag .ne. 0) then
1499 call update_error(ddx_error, &
1500 & "multipole_force_terms: grad_e returned an error, exiting")
1501 return
1502 end if
1503 end if
1504end subroutine multipole_force_terms
1505
1506end module ddx_multipolar_solutes
ddx_core::allocate_electrostatics
subroutine allocate_electrostatics(params, constants, electrostatics, ddx_error)
Given the chosen model, find the required electrostatic properties, and allocate the arrays for them ...
Definition: ddx_core.f90:341
ddx_multipolar_solutes::multipole_electrostatics
subroutine multipole_electrostatics(params, constants, workspace, multipoles, mmax, electrostatics, ddx_error)
Given a multipolar distribution, compute the required electrostatic properties for the chosen model.
Definition: ddx_multipolar_solutes.f90:30
ddx_multipolar_solutes::multipole_force_terms
subroutine multipole_force_terms(params, constants, workspace, state, mmax, multipoles, forces, ddx_error)
Given a multipolar distribution in real spherical harmonics and centered on the spheres,...
Definition: ddx_multipolar_solutes.f90:1477
ddx_multipolar_solutes::multipole_psi
subroutine multipole_psi(params, multipoles, mmax, psi)
Given a multipolar distribution, assemble the RHS psi. The multipoles must be centered on the ddx sph...
Definition: ddx_multipolar_solutes.f90:839
ddx_core
Core routines and parameters of the ddX software.
Definition: ddx_core.f90:13
ddx_core::tree_grad_l2l
subroutine tree_grad_l2l(params, constants, node_l, sph_l_grad, work)
TODO.
Definition: ddx_core.f90:2714
ddx_core::tree_l2l_rotation
subroutine tree_l2l_rotation(params, constants, node_l)
Transfer local coefficients over a tree.
Definition: ddx_core.f90:1866
ddx_core::tree_m2m_rotation
subroutine tree_m2m_rotation(params, constants, node_m)
Transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1656
ddx_core::tree_l2p
subroutine tree_l2p(params, constants, alpha, node_l, beta, grid_v, sph_l)
TODO.
Definition: ddx_core.f90:2252
ddx_core::tree_l2p_adj
subroutine tree_l2p_adj(params, constants, alpha, grid_v, beta, node_l, sph_l)
TODO.
Definition: ddx_core.f90:2326
ddx_core::tree_m2l_rotation
subroutine tree_m2l_rotation(params, constants, node_m, node_l)
Transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2046
ddx_core::tree_m2p
subroutine tree_m2p(params, constants, p, alpha, sph_m, beta, grid_v)
TODO.
Definition: ddx_core.f90:2401
ddx_core::tree_m2l_rotation_adj
subroutine tree_m2l_rotation_adj(params, constants, node_l, node_m)
Adjoint transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2208
ddx_core::tree_l2l_rotation_adj
subroutine tree_l2l_rotation_adj(params, constants, node_l)
Adjoint transfer local coefficients over a tree.
Definition: ddx_core.f90:1959
ddx_core::tree_m2m_rotation_adj
subroutine tree_m2m_rotation_adj(params, constants, node_m)
Adjoint transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1778
ddx_definitions
Compile-time constants and definitions.
Definition: ddx_definitions.f90:13
ddx_gradients
Core routines and parameters specific to gradients.
Definition: ddx_gradients.f90:11
ddx_multipolar_solutes
Routines to build rhs (phi and psi)
Definition: ddx_multipolar_solutes.f90:2
ddx_core::ddx_electrostatics_type
Container for the electrostatic properties. Since different methods require different electrostatic p...
Definition: ddx_core.f90:195
ddx_core::ddx_state_type
This defined type contains the primal and adjoint RHSs, the solution of the primal and adjoint linear...
Definition: ddx_core.f90:33