ddx 0.6.0
Libary for domain-decomposition methods for polarizable continuum models
ddx_gradients.f90
1
9
11module ddx_gradients
12! Get the core-routines
13use ddx_core
14!
15contains
16
18subroutine contract_grad_l(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx)
19type(ddx_params_type), intent(in) :: params
20 type(ddx_constants_type), intent(in) :: constants
21 integer, intent(in) :: isph
22 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sigma
23 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi
24 real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm
25 real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc
26 real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin
27 real(dp), dimension(3), intent(inout) :: fx
28
29 call contract_gradi_lik(params, constants, isph, sigma, xi(:, isph), basloc, dbsloc, vplm, vcos, vsin, fx )
30 call contract_gradi_lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx )
31
32end subroutine contract_grad_l
33
35subroutine contract_gradi_lik(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx )
36 type(ddx_params_type), intent(in) :: params
37 type(ddx_constants_type), intent(in) :: constants
38 integer, intent(in) :: isph
39 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sigma
40 real(dp), dimension(params % ngrid), intent(in) :: xi
41 real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm
42 real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc
43 real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin
44 real(dp), dimension(3), intent(inout) :: fx
45 integer :: ig, ij, jsph, l, ind, m
46 real(dp) :: vvij, tij, xij, oij, t, fac, fl, f1, f2, f3, beta, tlow, thigh
47 real(dp) :: vij(3), sij(3), alp(3), va(3)
48 real(dp), external :: dnrm2
49 tlow = one - pt5*(one - params % se)*params % eta
50 thigh = one + pt5*(one + params % se)*params % eta
51
52 do ig = 1, params % ngrid
53 va = zero
54 do ij = constants % inl(isph), constants % inl(isph+1) - 1
55 jsph = constants % nl(ij)
56 vij = params % csph(:,isph) + &
57 & params % rsph(isph)*constants % cgrid(:,ig) - &
58 & params % csph(:,jsph)
59 vvij = dnrm2(3, vij, 1)
60 tij = vvij/params % rsph(jsph)
61 if (tij.ge.thigh) cycle
62 if (tij.ne.zero) then
63 sij = vij/vvij
64 else
65 sij = one
66 end if
67 call dbasis(params, constants, sij, basloc, dbsloc, vplm, vcos, vsin)
68 alp = zero
69 t = one
70 do l = 1, params % lmax
71 ind = l*l + l + 1
72 fl = dble(l)
73 fac = t/(constants % vscales(ind)**2)
74 do m = -l, l
75 f2 = fac*sigma(ind+m,jsph)
76 f1 = f2*fl*basloc(ind+m)
77 alp(:) = alp(:) + f1*sij(:) + f2*dbsloc(:,ind+m)
78 end do
79 t = t*tij
80 end do
81 beta = intmlp(params, constants, tij,sigma(:,jsph),basloc)
82 xij = fsw(tij, params % se, params % eta)
83 if (constants % fi(ig,isph).gt.one) then
84 oij = xij/constants % fi(ig,isph)
85 f2 = -oij/constants % fi(ig,isph)
86 else
87 oij = xij
88 f2 = zero
89 end if
90 f1 = oij/params % rsph(jsph)
91 va(:) = va(:) + f1*alp(:) + beta*f2*constants % zi(:,ig,isph)
92 if (tij .gt. tlow) then
93 f3 = beta*dfsw(tij,params % se,params % eta)/params % rsph(jsph)
94 if (constants % fi(ig,isph).gt.one) f3 = f3/constants % fi(ig,isph)
95 va(:) = va(:) + f3*sij(:)
96 end if
97 end do
98 fx = fx - constants % wgrid(ig)*xi(ig)*va(:)
99 end do
100end subroutine contract_gradi_lik
101
103subroutine contract_gradi_lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx )
104 type(ddx_params_type), intent(in) :: params
105 type(ddx_constants_type), intent(in) :: constants
106 integer, intent(in) :: isph
107 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sigma
108 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi
109 real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm
110 real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc
111 real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin
112 real(dp), dimension(3), intent(inout) :: fx
113
114 integer :: ig, ji, jsph, l, ind, m, jk, ksph
115 logical :: proc
116 real(dp) :: vvji, tji, xji, oji, t, fac, fl, f1, f2, beta, di, tlow, thigh
117 real(dp) :: b, g1, g2, vvjk, tjk, xjk
118 real(dp) :: vji(3), sji(3), alp(3), vb(3), vjk(3), sjk(3), vc(3)
119 real(dp) :: rho, ctheta, stheta, cphi, sphi
120 real(dp), external :: dnrm2
121
122 tlow = one - pt5*(one - params % se)*params % eta
123 thigh = one + pt5*(one + params % se)*params % eta
124
125 do ig = 1, params % ngrid
126 vb = zero
127 vc = zero
128 do ji = constants % inl(isph), constants % inl(isph+1) - 1
129 jsph = constants % nl(ji)
130 vji = params % csph(:,jsph) + &
131 & params % rsph(jsph)*constants % cgrid(:,ig) - &
132 & params % csph(:,isph)
133 vvji = dnrm2(3, vji, 1)
134 tji = vvji/params % rsph(isph)
135 if (tji.gt.thigh) cycle
136 if (tji.ne.zero) then
137 sji = vji/vvji
138 else
139 sji = one
140 end if
141 call dbasis(params, constants, sji, basloc, dbsloc, vplm, vcos, vsin)
142 alp = zero
143 t = one
144 do l = 1, params % lmax
145 ind = l*l + l + 1
146 fl = dble(l)
147 fac = t/(constants % vscales(ind)**2)
148 do m = -l, l
149 f2 = fac*sigma(ind+m,isph)
150 f1 = f2*fl*basloc(ind+m)
151 alp = alp + f1*sji + f2*dbsloc(:,ind+m)
152 end do
153 t = t*tji
154 end do
155 xji = fsw(tji, params % se, params % eta)
156 if (constants % fi(ig,jsph).gt.one) then
157 oji = xji/constants % fi(ig,jsph)
158 else
159 oji = xji
160 end if
161 f1 = oji/params % rsph(isph)
162 vb = vb + f1*alp*xi(ig,jsph)
163 if (tji .gt. tlow) then
164 beta = intmlp(params, constants, tji, sigma(:,isph), basloc)
165 if (constants % fi(ig,jsph) .gt. one) then
166 di = one/constants % fi(ig,jsph)
167 fac = di*xji
168 proc = .false.
169 b = zero
170 do jk = constants % inl(jsph), constants % inl(jsph+1) - 1
171 ksph = constants % nl(jk)
172 vjk = params % csph(:,jsph) + &
173 & params % rsph(jsph)*constants % cgrid(:,ig) - &
174 & params % csph(:,ksph)
175 !vvjk = sqrt(dot_product(vjk,vjk))
176 vvjk = dnrm2(3, vjk, 1)
177 tjk = vvjk/params % rsph(ksph)
178 if (ksph.ne.isph) then
179 if (tjk .le. thigh) then
180 proc = .true.
181 sjk = vjk/vvjk
182 !call ylmbas(sjk,basloc,vplm,vcos,vsin)
183 call ylmbas(sjk, rho, ctheta, stheta, cphi, sphi, &
184 & params % lmax, constants % vscales, basloc, vplm, &
185 & vcos, vsin)
186 g1 = intmlp(params, constants, tjk, sigma(:,ksph), basloc)
187 xjk = fsw(tjk, params % se, params % eta)
188 b = b + g1*xjk
189 end if
190 end if
191 end do
192 if (proc) then
193 g1 = di*di*dfsw(tji, params % se, params % eta)/params % rsph(isph)
194 g2 = g1*xi(ig,jsph)*b
195 vc = vc + g2*sji
196 end if
197 else
198 di = one
199 fac = zero
200 end if
201 f2 = (one-fac)*di*dfsw(tji, params % se, params % eta)/params % rsph(isph)
202 vb = vb + f2*xi(ig,jsph)*beta*sji
203 end if
204 end do
205 fx = fx + constants % wgrid(ig)*(vb - vc)
206 end do
207end subroutine contract_gradi_lji
208
210subroutine contract_grad_u(params, constants, isph, xi, phi, fx )
211 type(ddx_params_type), intent(in) :: params
212 type(ddx_constants_type), intent(in) :: constants
213 integer, intent(in) :: isph
214 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi, phi
215 real(dp), dimension(3), intent(inout) :: fx
216 integer :: ig, ji, jsph
217 real(dp) :: vvji, tji, fac, swthr
218 real(dp) :: alp(3), vji(3), sji(3)
219 real(dp), external :: dnrm2
220 do ig = 1, params % ngrid
221 alp = zero
222 if (constants % ui(ig,isph) .gt. zero .and. constants % ui(ig,isph).lt.one) then
223 alp = alp + phi(ig,isph)*xi(ig,isph)*constants % zi(:,ig,isph)
224 end if
225 do ji = constants % inl(isph), constants % inl(isph+1) - 1
226 jsph = constants % nl(ji)
227 vji = params % csph(:,jsph) + &
228 & params % rsph(jsph)*constants % cgrid(:,ig) - &
229 & params % csph(:,isph)
230 !vvji = sqrt(dot_product(vji,vji))
231 vvji = dnrm2(3, vji, 1)
232 tji = vvji/params % rsph(isph)
233 swthr = one + (params % se + 1.d0)*params % eta / 2.d0
234 if (tji.lt.swthr .and. tji.gt.swthr-params % eta .and. constants % ui(ig,jsph).gt.zero) then
235 sji = vji/vvji
236 fac = - dfsw(tji, params % se, params % eta)/params % rsph(isph)
237 alp = alp + fac*phi(ig,jsph)*xi(ig,jsph)*sji
238 end if
239 end do
240 fx = fx - constants % wgrid(ig)*alp
241 end do
242end subroutine contract_grad_u
243
252
253subroutine contract_grad_b(params, constants, isph, Xe, Xadj_e, force)
254 !! input/output
255 type(ddx_params_type), intent(in) :: params
256 type(ddx_constants_type), intent(in) :: constants
257 integer, intent(in) :: isph
258 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xe
259 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: Xadj_e
260 real(dp), dimension(3), intent(inout) :: force
261
262 call contract_gradi_bik(params, constants, isph, xe(:,:), &
263 & xadj_e(:, isph), force)
264 call contract_gradi_bji(params, constants, isph, xe(:,:), xadj_e, force)
265
266end subroutine contract_grad_b
267
284subroutine contract_grad_c(params, constants, workspace, Xr, Xe, Xadj_r_sgrid, &
285 & Xadj_e_sgrid, Xadj_r, Xadj_e, force, diff_re, ddx_error)
286 !! input/output
287 type(ddx_params_type), intent(in) :: params
288 type(ddx_constants_type), intent(in) :: constants
289 type(ddx_workspace_type), intent(inout) :: workspace
290 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xr, Xe
291 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: Xadj_r_sgrid, Xadj_e_sgrid
292 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xadj_r, Xadj_e
293 real(dp), dimension(3, params % nsph), intent(inout) :: force
294 real(dp), dimension(constants % nbasis, params % nsph), intent(out) :: diff_re
295 type(ddx_error_type), intent(inout) :: ddx_error
296
297 call contract_grad_c_worker2(params, constants, workspace, xr, xe, xadj_r_sgrid, &
298 & xadj_e_sgrid, xadj_r, xadj_e, force, diff_re, ddx_error)
299 if (ddx_error % flag .ne. 0) then
300 call update_error(ddx_error, &
301 & "contract_grad_C_worker2 returned an error, exiting")
302 return
303 end if
304 call contract_grad_c_worker1(params, constants, workspace, xadj_r_sgrid, &
305 & xadj_e_sgrid, diff_re, force, ddx_error)
306 if (ddx_error % flag .ne. 0) then
307 call update_error(ddx_error, &
308 & "contract_grad_C_worker1 returned an error, exiting")
309 return
310 end if
311
312end subroutine contract_grad_c
313
314
325subroutine contract_grad_f(params, constants, workspace, sol_adj, sol_sgrid, &
326 & gradpsi, normal_hessian_cav, force, state, ddx_error)
327 type(ddx_params_type), intent(in) :: params
328 type(ddx_constants_type), intent(in) :: constants
329 type(ddx_workspace_type), intent(inout) :: workspace
330 type(ddx_state_type), intent(inout) :: state
331 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sol_adj
332 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: sol_sgrid
333 real(dp), dimension(3, constants % ncav), intent(in) :: gradpsi
334 real(dp), dimension(3, constants % ncav), intent(in) :: normal_hessian_cav
335 real(dp), dimension(3, params % nsph), intent(inout) :: force
336 type(ddx_error_type), intent(inout) :: ddx_error
337
338 call contract_grad_f_worker1(params, constants, workspace, sol_adj, sol_sgrid, &
339 & gradpsi, force, ddx_error)
340 if (ddx_error % flag .ne. 0) then
341 call update_error(ddx_error, &
342 & "contract_grad_f_worker1 returned an error, exiting")
343 return
344 end if
345
346 call contract_grad_f_worker2(params, constants, gradpsi, &
347 & normal_hessian_cav, force, state, ddx_error)
348 if (ddx_error % flag .ne. 0) then
349 call update_error(ddx_error, &
350 & "contract_grad_f_worker2 returned an error, exiting")
351 return
352 end if
353
354end subroutine contract_grad_f
355
364subroutine contract_gradi_bik(params, constants, isph, Xe, Xadj_e, force)
365 !! input/output
366 type(ddx_params_type), intent(in) :: params
367 type(ddx_constants_type), intent(in) :: constants
368 integer, intent(in) :: isph
369 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xe
370 real(dp), dimension(params % ngrid), intent(in) :: Xadj_e
371 real(dp), dimension(3), intent(inout) :: force
372
373 ! Local Variables
374 integer :: igrid, ineigh, jsph
375 real(dp), dimension(0:params % lmax) :: SI_rijn
376 real(dp), dimension(0:params % lmax) :: DI_rijn
377 ! beta : Eq.(53) Stamm.etal.18
378 ! tlow : Lower bound for switch region
379 ! thigh : Upper bound for switch region
380 ! f1 : First factor in alpha computation
381 ! f2 : Second factor in alpha computation
382 real(dp) :: rijn, tij, beta, tlow, thigh, xij, oij, f1, f2, f3
383 ! alpha : Eq.(52) Stamm.etal.18
384 ! va : Eq.(54) Stamm.etal.18
385 real(dp) :: vij(3), sij(3), alpha(3), va(3), rj, vtij(3)
386 real(dp), external :: dnrm2
387 real(dp) :: work(params % lmax+1)
388 complex(dp) :: work_complex(params % lmax+1)
389
390 si_rijn = 0
391 di_rijn = 0
392
393 tlow = one - pt5*(one - params % se)*params % eta
394 thigh = one + pt5*(one + params % se)*params % eta
395
396 ! Loop over grid points
397 do igrid = 1, params % ngrid
398 va = zero
399 do ineigh = constants % inl(isph), constants % inl(isph+1) - 1
400 jsph = constants % nl(ineigh)
401 vij = params % csph(:,isph) + &
402 & params % rsph(isph)*constants % cgrid(:,igrid) - &
403 & params % csph(:,jsph)
404 rijn = dnrm2(3, vij, 1)
405 tij = rijn/params % rsph(jsph)
406 rj = params % rsph(jsph)
407 if (tij.ge.thigh) cycle
408 ! Computation of modified spherical Bessel function values
409 if (tij.ne.zero) then
410 sij = vij/rijn
411 else
412 sij = one
413 end if
414 vtij = vij*params % kappa
415 call fmm_l2p_bessel_grad(vtij, params % rsph(jsph)*params % kappa, &
416 & params % lmax, constants % vscales, params % kappa, &
417 & xe(:, jsph), zero, alpha)
418 call fmm_l2p_bessel_work(vtij, params % lmax, constants % vscales, &
419 & constants % SI_ri(:, jsph), one, xe(:, jsph), zero, beta, &
420 & work_complex, work)
421 xij = fsw(tij, params % se, params % eta)
422 if (constants % fi(igrid,isph).gt.one) then
423 oij = xij/constants % fi(igrid,isph)
424 f2 = -oij/constants % fi(igrid,isph)
425 else
426 oij = xij
427 f2 = zero
428 end if
429 f1 = oij
430 va(:) = va(:) + f1*alpha(:) + beta*f2*constants % zi(:,igrid,isph)
431 if (tij .gt. tlow) then
432 f3 = beta*dfsw(tij,params % se,params % eta)/params % rsph(jsph)
433 if (constants % fi(igrid,isph).gt.one) &
434 & f3 = f3/constants % fi(igrid,isph)
435 va(:) = va(:) + f3*sij(:)
436 end if
437 end do
438 force = force - constants % wgrid(igrid)*xadj_e(igrid)*va(:)
439 end do
440end subroutine contract_gradi_bik
441
442
452subroutine contract_gradi_bji(params, constants, isph, Xe, Xadj_e, force)
453 !! input/output
454 type(ddx_params_type), intent(in) :: params
455 type(ddx_constants_type), intent(in) :: constants
456 integer, intent(in) :: isph
457 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xe
458 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: Xadj_e
459 real(dp), dimension(3), intent(inout) :: force
460
461 ! Local Variables
462 ! jk : Row pointer over kth row
463 integer :: igrid, jsph, ksph, ineigh, jk
464 real(dp), dimension(0:params % lmax) :: SI_rjin, SI_rjkn
465 real(dp), dimension(0:params % lmax) :: DI_rjin, DI_rjkn
466
467 logical :: proc
468 ! fac : \delta_fj_n*\omega^\eta_ji
469 ! f1 : First factor in alpha computation
470 ! f2 : Second factor in alpha computation
471 ! beta_ji : Eq.(57) Stamm.etal.18
472 ! dj : Before Eq.(10) Stamm.etal.18
473 real(dp) :: rjin, tji, xji, oji, fac, f1, f2, beta_ji, dj, tlow, thigh
474 real(dp) :: b, beta_jk, g1, g2, rjkn, tjk, xjk
475 ! alpha : Eq.(56) Stamm.etal.18
476 ! vb : Eq.(60) Stamm.etal.18
477 ! vc : Eq.(59) Stamm.etal.18
478 real(dp) :: vji(3), sji(3), vjk(3), alpha(3), vb(3), vc(3), &
479 & vtji(3), vtjk(3)
480 ! rho : Argument for ylmbas
481 ! ctheta : Argument for ylmbas
482 ! stheta : Argument for ylmbas
483 ! cphi : Argument for ylmbas
484 ! sphi : Argument for ylmbas
485 real(dp) :: ri
486
487 real(dp), external :: dnrm2
488 real(dp) :: work(params % lmax+1)
489 complex(dp) :: work_complex(params % lmax+1)
490
491 si_rjin = 0
492 di_rjin = 0
493 si_rjkn = 0
494 di_rjkn = 0
495
496 tlow = one - pt5*(one - params % se)*params % eta
497 thigh = one + pt5*(one + params % se)*params % eta
498
499 do igrid = 1, params % ngrid
500 vb = zero
501 vc = zero
502 do ineigh = constants % inl(isph), constants % inl(isph+1) - 1
503 jsph = constants % nl(ineigh)
504 vji = params % csph(:,jsph) + &
505 & params % rsph(jsph)*constants % cgrid(:,igrid) - &
506 & params % csph(:,isph)
507 rjin = dnrm2(3, vji, 1)
508 ri = params % rsph(isph)
509 tji = rjin/ri
510 if (tji.gt.thigh) cycle
511 if (tji.ne.zero) then
512 sji = vji/rjin
513 else
514 sji = one
515 end if
516 vtji = vji*params % kappa
517 call fmm_l2p_bessel_grad(vtji, params % rsph(isph)*params % kappa, &
518 & params % lmax, constants % vscales, params % kappa, &
519 & xe(:, isph), zero, alpha)
520 xji = fsw(tji,params % se,params % eta)
521 if (constants % fi(igrid,jsph).gt.one) then
522 oji = xji/constants % fi(igrid,jsph)
523 else
524 oji = xji
525 end if
526 f1 = oji
527 vb = vb + f1*alpha*xadj_e(igrid,jsph)
528 if (tji .gt. tlow) then
529 call fmm_l2p_bessel_work(vtji, params % lmax, &
530 & constants % vscales, constants % SI_ri(:, isph), one, &
531 & xe(:, isph), zero, beta_ji, work_complex, work)
532 if (constants % fi(igrid,jsph) .gt. one) then
533 dj = one/constants % fi(igrid,jsph)
534 fac = dj*xji
535 proc = .false.
536 b = zero
537 do jk = constants % inl(jsph), constants % inl(jsph+1) - 1
538 ksph = constants % nl(jk)
539 vjk = params % csph(:,jsph) + &
540 & params % rsph(jsph)*constants % cgrid(:,igrid) - &
541 & params % csph(:,ksph)
542 rjkn = dnrm2(3, vjk, 1)
543 tjk = rjkn/params % rsph(ksph)
544 if (ksph.ne.isph) then
545 if (tjk .le. thigh) then
546 proc = .true.
547 vtjk = vjk*params % kappa
548 call fmm_l2p_bessel_work(vtjk, params % lmax, &
549 & constants % vscales, &
550 & constants % SI_ri(:, ksph), one, xe(:, ksph), &
551 & zero, beta_jk, work_complex, work)
552 xjk = fsw(tjk, params % se, params % eta)
553 b = b + beta_jk*xjk
554 end if
555 end if
556 end do
557 if (proc) then
558 g1 = dj*dj*dfsw(tji,params % se,params % eta) &
559 & /params % rsph(isph)
560 g2 = g1*xadj_e(igrid,jsph)*b
561 vc = vc + g2*sji
562 end if
563 else
564 dj = one
565 fac = zero
566 end if
567 f2 = (one-fac)*dj*dfsw(tji,params % se,params % eta) &
568 & /params % rsph(isph)
569 vb = vb + f2*xadj_e(igrid,jsph)*beta_ji*sji
570 end if
571 end do
572 force = force + constants % wgrid(igrid)*(vb - vc)
573 end do
574end subroutine contract_gradi_bji
575
576
588subroutine contract_grad_c_worker1(params, constants, workspace, Xadj_r_sgrid, &
589 & Xadj_e_sgrid, diff_re, force, ddx_error)
590 !! Inputs
591 type(ddx_params_type), intent(in) :: params
592 type(ddx_constants_type), intent(in) :: constants
593 type(ddx_workspace_type), intent(inout) :: workspace
594 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: &
595 & diff_re
596 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: &
597 & Xadj_r_sgrid, Xadj_e_sgrid
598 real(dp), dimension(3, params % nsph), intent(inout) :: force
599 type(ddx_error_type), intent(inout) :: ddx_error
600 ! Local variable
601 ! igrid0: Index for grid point n0
602 integer :: isph, jsph, igrid, l0, m0, ind0, igrid0, icav, &
603 & indl, inode, istat
604 ! term : SK_rijn/SK_rj
605 ! termi : DI_ri/SI_ri
606 ! termk : DK_ri/SK_ri
607 ! sum_int : Intermediate sum
608 ! sum_r : Intermediate sum for Laplace
609 ! sum_e : Intermediate sum for HSP
610 real(dp) :: sum_int, sum_r, sum_e
611 real(dp) :: rijn
612 real(dp) :: vij(3), sij(3), vtij(3)
613
614 ! local allocatable
615 ! phi_n_r : Phi corresponding to Laplace problem
616 ! phi_n_e : Phi corresponding to HSP problem
617 ! coefY_d : sum_{l0m0} C_ik*term*Y_l0m0^j(x_in)*Y_l0m0(s_n)
618 ! diff_re_sgrid : diff_re evaluated at grid point
619 real(dp), allocatable :: phi_n_r(:,:), phi_n_e(:,:), coefY_d(:,:,:), &
620 & diff_re_sgrid(:,:)
621
622 ! basloc : Y_lm(s_n)
623 ! vplm : Argument to call ylmbas
624 real(dp), dimension(constants % nbasis):: basloc, vplm
625 ! dbasloc : Derivative of Y_lm(s_n)
626 real(dp), dimension(3, constants % nbasis):: dbasloc
627 ! vcos : Argument to call ylmbas
628 ! vsin : Argument to call ylmbas
629 real(dp), dimension(params % lmax+1):: vcos, vsin
630 ! SK_rijn : Besssel function of first kind for rijn
631 ! DK_rijn : Derivative of Besssel function of first kind for rijn
632 real(dp), dimension(0:params % lmax) :: SK_rijn, DK_rijn
633 real(dp) :: coef(constants % nbasis0), work(constants % lmax0+1)
634 complex(dp) :: work_complex(constants % lmax0+1)
635
636 allocate(phi_n_r(params % ngrid, params % nsph), &
637 & phi_n_e(params % ngrid, params % nsph), &
638 & diff_re_sgrid(params % ngrid, params % nsph), stat=istat)
639 if (istat.ne.0) then
640 call update_error(ddx_error, "allocation ddx_error in ddx_contract_grad_C_worker1")
641 return
642 end if
643!
644 sum_int = zero
645 sum_r = zero
646 sum_e = zero
647 phi_n_r = zero
648 phi_n_e = zero
649 diff_re_sgrid = zero
650 basloc = zero
651 vplm = zero
652 dbasloc = zero
653 vcos = zero
654 vsin = zero
655 sk_rijn = zero
656 dk_rijn = zero
657
658 if (params % fmm .eq. 0) then
659 allocate(coefy_d(constants % ncav, params % ngrid, params % nsph), &
660 & stat=istat)
661 if (istat.ne.0) then
662 call update_error(ddx_error, "allocation ddx_error in fmm ddx_contract_grad_C_worker1")
663 return
664 end if
665 coefy_d = zero
666 ! Compute summation over l0, m0
667 ! Loop over the sphers j
668 do jsph = 1, params % nsph
669 ! Loop over the grid points n0
670 do igrid0 = 1, params % ngrid
671 icav = zero
672 ! Loop over spheres i
673 do isph = 1, params % nsph
674 ! Loop over grid points n
675 do igrid = 1, params % ngrid
676 ! Check for U_i^{eta}(x_in)
677 if(constants % ui(igrid, isph) .gt. zero) then
678 icav = icav + 1
679 vij = params % csph(:,isph) + &
680 & params % rsph(isph)*constants % cgrid(:,igrid) - &
681 & params % csph(:,jsph)
682 rijn = sqrt(dot_product(vij,vij))
683 sij = vij/rijn
684
685 do l0 = 0, constants % lmax0
686 do m0 = -l0,l0
687 ind0 = l0**2 + l0 + m0 + 1
688 coef(ind0) = constants % vgrid(ind0, igrid0) * &
689 & constants % C_ik(l0, jsph)
690 end do
691 end do
692 vtij = vij*params % kappa
693 call fmm_m2p_bessel_work(vtij, constants % lmax0, &
694 & constants % vscales, constants % SK_ri(:, jsph), one, &
695 & coef, zero, coefy_d(icav, igrid0, jsph), work_complex, work)
696 end if
697 end do
698 end do
699 end do
700 end do
701
702 ! Compute phi_in
703 ! Loop over spheres j
704 do jsph = 1, params % nsph
705 ! Loop over grid points n0
706 do igrid0 = 1, params % ngrid
707 icav = zero
708 sum_r = zero
709 sum_e = zero
710 ! Loop over sphers i
711 do isph = 1, params % nsph
712 ! Loop over grid points n
713 do igrid = 1, params % ngrid
714 if(constants % ui(igrid, isph) .gt. zero) then
715 icav = icav + 1
716 sum_r = sum_r + coefy_d(icav, igrid0, jsph)*xadj_r_sgrid(igrid, isph) &
717 & * constants % wgrid(igrid)*constants % ui(igrid, isph)
718 sum_e = sum_e + coefy_d(icav, igrid0, jsph)*xadj_e_sgrid(igrid, isph) &
719 & * constants % wgrid(igrid)*constants % ui(igrid, isph)
720 end if
721 end do
722 end do
723 phi_n_r(igrid0, jsph) = sum_r
724 phi_n_e(igrid0, jsph) = sum_e
725 end do
726 end do
727 else
728 ! Compute phi_n_r at first
729 ! Adjoint integration from spherical harmonics to grid points is not needed
730 ! here as ygrid already contains grid values, we just need to scale it by
731 ! weights of grid points
732 do isph = 1, params % nsph
733 workspace % tmp_grid(:, isph) = xadj_r_sgrid(:, isph) * &
734 & constants % wgrid(:) * constants % ui(:, isph)
735 end do
736 ! Adjoint FMM
737 call tree_m2p_bessel_adj(params, constants, constants % lmax0, one, &
738 & workspace % tmp_grid, zero, params % lmax, workspace % tmp_sph)
739 call tree_l2p_bessel_adj(params, constants, one, workspace % tmp_grid, &
740 & zero, workspace % tmp_node_l)
741 call tree_l2l_bessel_rotation_adj(params, constants, &
742 & workspace % tmp_node_l)
743 call tree_m2l_bessel_rotation_adj(params, constants, &
744 & workspace % tmp_node_l, workspace % tmp_node_m)
745 call tree_m2m_bessel_rotation_adj(params, constants, &
746 & workspace % tmp_node_m)
747 ! Properly load adjoint multipole harmonics into tmp_sph
748 if(constants % lmax0 .lt. params % pm) then
749 do isph = 1, params % nsph
750 inode = constants % snode(isph)
751 workspace % tmp_sph(1:constants % nbasis0, isph) = &
752 & workspace % tmp_sph(1:constants % nbasis0, isph) + &
753 & workspace % tmp_node_m(1:constants % nbasis0, inode)
754 end do
755 else
756 indl = (params % pm+1)**2
757 do isph = 1, params % nsph
758 inode = constants % snode(isph)
759 workspace % tmp_sph(1:indl, isph) = &
760 & workspace % tmp_sph(1:indl, isph) + &
761 & workspace % tmp_node_m(:, inode)
762 end do
763 end if
764 ! Scale by C_ik
765 do isph = 1, params % nsph
766 do l0 = 0, constants % lmax0
767 ind0 = l0*l0 + l0 + 1
768 workspace % tmp_sph(ind0-l0:ind0+l0, isph) = &
769 & workspace % tmp_sph(ind0-l0:ind0+l0, isph) * &
770 & constants % C_ik(l0, isph)
771 end do
772 end do
773 ! Multiply by vgrid
774 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis0, &
775 & one, constants % vgrid, constants % vgrid_nbasis, &
776 & workspace % tmp_sph, constants % nbasis, zero, phi_n_r, params % ngrid)
777 ! Compute phi_n_e now
778 ! Adjoint integration from spherical harmonics to grid points is not needed
779 ! here as ygrid already contains grid values, we just need to scale it by
780 ! weights of grid points
781 do isph = 1, params % nsph
782 workspace % tmp_grid(:, isph) = xadj_e_sgrid(:, isph) * &
783 & constants % wgrid(:) * constants % ui(:, isph)
784 end do
785 ! Adjoint FMM
786 call tree_m2p_bessel_adj(params, constants, constants % lmax0, one, &
787 & workspace % tmp_grid, zero, params % lmax, workspace % tmp_sph)
788 call tree_l2p_bessel_adj(params, constants, one, workspace % tmp_grid, &
789 & zero, workspace % tmp_node_l)
790 call tree_l2l_bessel_rotation_adj(params, constants, &
791 & workspace % tmp_node_l)
792 call tree_m2l_bessel_rotation_adj(params, constants, &
793 & workspace % tmp_node_l, workspace % tmp_node_m)
794 call tree_m2m_bessel_rotation_adj(params, constants, &
795 & workspace % tmp_node_m)
796 ! Properly load adjoint multipole harmonics into tmp_sph
797 if(constants % lmax0 .lt. params % pm) then
798 do isph = 1, params % nsph
799 inode = constants % snode(isph)
800 workspace % tmp_sph(1:constants % nbasis0, isph) = &
801 & workspace % tmp_sph(1:constants % nbasis0, isph) + &
802 & workspace % tmp_node_m(1:constants % nbasis0, inode)
803 end do
804 else
805 indl = (params % pm+1)**2
806 do isph = 1, params % nsph
807 inode = constants % snode(isph)
808 workspace % tmp_sph(1:indl, isph) = &
809 & workspace % tmp_sph(1:indl, isph) + &
810 & workspace % tmp_node_m(:, inode)
811 end do
812 end if
813 ! Scale by C_ik
814 do isph = 1, params % nsph
815 do l0 = 0, constants % lmax0
816 ind0 = l0*l0 + l0 + 1
817 workspace % tmp_sph(ind0-l0:ind0+l0, isph) = &
818 & workspace % tmp_sph(ind0-l0:ind0+l0, isph) * &
819 & constants % C_ik(l0, isph)
820 end do
821 end do
822 ! Multiply by vgrid
823 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis0, &
824 & one, constants % vgrid, constants % vgrid_nbasis, &
825 & workspace % tmp_sph, constants % nbasis, zero, phi_n_e, params % ngrid)
826 end if
827 call dgemm('T', 'N', params % ngrid, params % nsph, &
828 & constants % nbasis, one, constants % vgrid, constants % vgrid_nbasis, &
829 & diff_re , constants % nbasis, zero, diff_re_sgrid, &
830 & params % ngrid)
831 do isph = 1, params % nsph
832 call contract_grad_u(params, constants, isph, diff_re_sgrid, phi_n_r, force(:, isph))
833 call contract_grad_u(params, constants, isph, diff_re_sgrid, phi_n_e, force(:, isph))
834 end do
835
836 deallocate(phi_n_r, phi_n_e, diff_re_sgrid, stat=istat)
837 if (istat.ne.0) then
838 call update_error(ddx_error, "deallocation ddx_error in ddx_contract_grad_C_worker1")
839 return
840 end if
841 if (allocated(coefy_d)) then
842 deallocate(coefy_d, stat=istat)
843 if (istat.ne.0) then
844 call update_error(ddx_error, "deallocation ddx_error in ddx_contract_grad_C_worker1")
845 return
846 end if
847 end if
848end subroutine contract_grad_c_worker1
849
850
869subroutine contract_grad_c_worker2(params, constants, workspace, Xr, Xe, &
870 & Xadj_r_sgrid, Xadj_e_sgrid, Xadj_r, Xadj_e, force, diff_re, ddx_error)
871 !! input/output
872 type(ddx_params_type), intent(in) :: params
873 type(ddx_constants_type), intent(in) :: constants
874 type(ddx_workspace_type), intent(inout) :: workspace
875 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xr, Xe
876 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: &
877 & Xadj_r_sgrid, Xadj_e_sgrid
878 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: &
879 & Xadj_r, Xadj_e
880 real(dp), dimension(3, params % nsph), intent(inout) :: force
881 real(dp), dimension(constants % nbasis, params % nsph), intent(out) :: &
882 & diff_re
883 type(ddx_error_type), intent(inout) :: ddx_error
884 real(dp), external :: dnrm2
885 ! Local variable
886 integer :: isph, jsph, igrid, l, m, ind, l0, ind0, icav, indl, inode, &
887 & ksph, knode, jnode, knear, jsph_node, istat
888 ! val_dim3 : Intermediate value array of dimension 3
889 real(dp), dimension(3) :: vij, vtij
890 ! val : Intermediate variable to compute diff_ep
891 real(dp) :: val
892 ! large local are allocatable
893 ! phi_in : sum_{j=1}^N diff0_j * coefY_j
894 ! diff_ep_dim3 : 3 dimensional couterpart of diff_ep
895 ! sum_dim3 : Storage of sum
896 ! diff0 : dot_product([PU_j]_l0m0^l'm', l'/r_j[Xr]_jl'm' -
897 ! (i'_l'(r_j)/i_l'(r_j))[Xe]_jl'm')
898 real(dp), allocatable :: diff_ep_dim3(:,:), phi_in(:,:), sum_dim3(:,:,:), &
899 & diff0(:,:), diff1(:,:), diff1_grad(:,:,:), l2l_grad(:,:,:)
900 real(dp), dimension(0:params % lmax) :: SK_rijn, DK_rijn
901 complex(dp) :: work_complex(constants % lmax0+2)
902 real(dp) :: work(constants % lmax0+2)
903
904 allocate(diff_ep_dim3(3, constants % ncav), &
905 & phi_in(params % ngrid, params % nsph), &
906 & sum_dim3(3, constants % nbasis, params % nsph), &
907 & diff0(constants % nbasis0, params % nsph), &
908 & diff1(constants % nbasis0, params % nsph), &
909 & diff1_grad((constants % lmax0+2)**2, 3, params % nsph), &
910 & l2l_grad((params % pl+2)**2, 3, params % nsph), stat=istat)
911 if (istat.ne.0) then
912 call update_error(ddx_error, "allocation ddx_error in ddx_contract_grad_C_worker2")
913 return
914 end if
915
916 ! Setting initial values to zero
917 sk_rijn = zero
918 dk_rijn = zero
919
920 diff_re = zero
921 ! Compute l'/r_j[Xr]_jl'm' -(i'_l'(r_j)/i_l'(r_j))[Xe]_jl'm'
922 !$omp parallel do default(none) shared(params,diff_re,constants,xr,xe) &
923 !$omp private(jsph,l,m,ind) schedule(dynamic)
924 do jsph = 1, params % nsph
925 do l = 0, params % lmax
926 do m = -l,l
927 ind = l**2 + l + m + 1
928 diff_re(ind,jsph) = (params % epsp/params % eps)*(l/params % rsph(jsph)) * &
929 & xr(ind,jsph) - constants % termimat(l,jsph)*xe(ind,jsph)
930 end do
931 end do
932 end do
933
934 ! diff0 = Pchi * diff_re, linear scaling
935 diff0 = zero
936 !$omp parallel do default(none) shared(params,constants,diff_re,diff0, &
937 !$omp diff1,diff1_grad) private(jsph,l0,ind0) schedule(dynamic)
938 do jsph = 1, params % nsph
939 do l0 = 0, constants % lmax0
940 do ind0 = l0*l0+1, l0*l0+2*l0+1
941 diff0(ind0, jsph) = dot_product(diff_re(:,jsph), &
942 & constants % Pchi(:,ind0, jsph))
943 diff1(ind0, jsph) = diff0(ind0, jsph) * constants % C_ik(l0, jsph)
944 end do
945 end do
946 ! Prepare diff1_grad
947 call fmm_m2m_bessel_grad(constants % lmax0, constants % SK_ri(:, jsph), &
948 & constants % vscales, diff1(:, jsph), diff1_grad(:, :, jsph))
949 end do
950
951 if (params % fmm .eq. 0) then
952 ! phi_in = diff0 * coefY
953 ! Here, summation over j takes place
954 phi_in = zero
955 icav = 0
956 do isph = 1, params % nsph
957 do igrid = 1, params % ngrid
958 if(constants % ui(igrid, isph) .gt. zero) then
959 ! Extrenal grid point
960 icav = icav + 1
961 val = zero
962 do jsph = 1, params % nsph
963 !do ind0 = 1, constants % nbasis0
964 !!====== This place requirs coefY, that is not precomputed anymore
965 ! val = val + diff0(ind0,jsph)*constants % coefY(icav,ind0,jsph)
966 !end do
967 vij = params % csph(:, isph) + &
968 & params % rsph(isph)*constants % cgrid(:, igrid) - &
969 & params % csph(:, jsph)
970 vtij = vij * params % kappa
971 call fmm_m2p_bessel_work(vtij, constants % lmax0, &
972 & constants % vscales, constants % SK_ri(:, jsph), one, &
973 & diff1(:, jsph), one, val, work_complex, work)
974 end do
975 phi_in(igrid, isph) = val
976 end if
977 end do
978 end do
979 else
980 ! phi_in
981 ! Load input harmonics into tree data
982 workspace % tmp_sph = zero
983 workspace % tmp_sph(1:constants % nbasis0, :) = diff1(:, :)
984 if(constants % lmax0 .lt. params % pm) then
985 do isph = 1, params % nsph
986 inode = constants % snode(isph)
987 workspace % tmp_node_m(1:constants % nbasis0, inode) = &
988 & workspace % tmp_sph(1:constants % nbasis0, isph)
989 workspace % tmp_node_m(constants % nbasis0+1:, inode) = zero
990 end do
991 else
992 indl = (params % pm+1)**2
993 do isph = 1, params % nsph
994 inode = constants % snode(isph)
995 workspace % tmp_node_m(:, inode) = workspace % tmp_sph(1:indl, isph)
996 end do
997 end if
998 ! Do FMM operations
999 call tree_m2m_bessel_rotation(params, constants, workspace % tmp_node_m)
1000 call tree_m2l_bessel_rotation(params, constants, workspace % tmp_node_m, &
1001 & workspace % tmp_node_l)
1002 call tree_l2l_bessel_rotation(params, constants, workspace % tmp_node_l)
1003 call tree_l2p_bessel(params, constants, one, workspace % tmp_node_l, zero, &
1004 & phi_in)
1005 call tree_m2p_bessel(params, constants, constants % lmax0, one, &
1006 & params % lmax, workspace % tmp_sph, one, &
1007 & phi_in)
1008 ! Make phi_in zero at internal grid points
1009 !$omp parallel do default(none) shared(params,constants,phi_in) &
1010 !$omp private(isph,igrid) schedule(dynamic)
1011 do isph = 1, params % nsph
1012 do igrid = 1, params % ngrid
1013 if (constants % ui(igrid, isph) .eq. zero) then
1014 phi_in(igrid, isph) = zero
1015 end if
1016 end do
1017 end do
1018 ! Get gradients of the L2L
1019 !$omp parallel do default(none) shared(params,constants,workspace, &
1020 !$omp l2l_grad) private(isph,igrid,inode) schedule(dynamic)
1021 do isph = 1, params % nsph
1022 inode = constants % snode(isph)
1023 workspace % tmp_sph_l(:, isph) = workspace % tmp_node_l(:, inode)
1024 call fmm_l2l_bessel_grad(params % pl, &
1025 & constants % SI_ri(:, isph), constants % vscales, &
1026 & workspace % tmp_node_l(:, inode), &
1027 & l2l_grad(:, :, isph))
1028 end do
1029 workspace % tmp_sph = xadj_r + xadj_e
1030 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis, &
1031 & one, constants % vwgrid, constants % vgrid_nbasis, &
1032 & workspace % tmp_sph, constants % nbasis, zero, &
1033 & workspace % tmp_grid, params % ngrid)
1034 workspace % tmp_grid = workspace % tmp_grid * constants % ui
1035 ! Adjoint FMM with output tmp_sph2(:, :) which stores coefficients of
1036 ! harmonics of degree up to lmax+1
1037 call tree_m2p_bessel_nodiag_adj(params, constants, constants % lmax0+1, one, &
1038 & workspace % tmp_grid, zero, params % lmax+1, workspace % tmp_sph2)
1039 call tree_l2p_bessel_adj(params, constants, one, workspace % tmp_grid, zero, &
1040 & workspace % tmp_node_l)
1041 call tree_l2l_bessel_rotation_adj(params, constants, workspace % tmp_node_l)
1042 call tree_m2l_bessel_rotation_adj(params, constants, workspace % tmp_node_l, &
1043 & workspace % tmp_node_m)
1044 call tree_m2m_bessel_rotation_adj(params, constants, workspace % tmp_node_m)
1045 ! Properly load adjoint multipole harmonics into tmp_sph2 that holds
1046 ! harmonics of a degree up to lmax+1
1047 if(constants % lmax0+1 .lt. params % pm) then
1048 !$omp parallel do default(none) shared(params,constants, &
1049 !$omp workspace) private(isph,inode) schedule(dynamic)
1050 do isph = 1, params % nsph
1051 inode = constants % snode(isph)
1052 workspace % tmp_sph2(1:(constants % lmax0+2)**2, isph) = &
1053 & workspace % tmp_sph2(1:(constants % lmax0+2)**2, isph) + &
1054 & workspace % tmp_node_m(1:(constants % lmax0+2)**2, inode)
1055 end do
1056 else
1057 indl = (params % pm+1)**2
1058 !$omp parallel do default(none) shared(params,constants,indl, &
1059 !$omp workspace) private(isph,inode) schedule(dynamic)
1060 do isph = 1, params % nsph
1061 inode = constants % snode(isph)
1062 workspace % tmp_sph2(1:indl, isph) = &
1063 & workspace % tmp_sph2(1:indl, isph) + &
1064 & workspace % tmp_node_m(:, inode)
1065 end do
1066 end if
1067 end if
1068
1069 do ksph = 1, params % nsph
1070 ! Computation of derivative of U_i^e(x_in)
1071 call contract_grad_u(params, constants, ksph, xadj_r_sgrid, phi_in, force(:, ksph))
1072 call contract_grad_u(params, constants, ksph, xadj_e_sgrid, phi_in, force(:, ksph))
1073
1074 ! Aleksandr: my loop for the diff_ep_dim3
1075 diff_ep_dim3 = zero
1076 ! At first isph=ksph, jsph!=ksph
1077 icav = constants % icav_ia(ksph) - 1
1078 if (params % fmm .eq. 0) then
1079 do igrid = 1, params % ngrid
1080 if (constants % ui(igrid, ksph) .eq. zero) cycle
1081 icav = icav + 1
1082 do jsph = 1, params % nsph
1083 if (jsph .eq. ksph) cycle
1084 vij = params % csph(:,ksph) + &
1085 & params % rsph(ksph)*constants % cgrid(:,igrid) - &
1086 & params % csph(:,jsph)
1087 vtij = vij * params % kappa
1088 !call fmm_m2p_bessel_grad(vtij, &
1089 ! & params % rsph(jsph)*params % kappa, &
1090 ! & constants % lmax0, &
1091 ! & constants % vscales, params % kappa, diff1(:, jsph), one, &
1092 ! & diff_ep_dim3(:, icav))
1093 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1094 & constants % vscales, constants % SK_ri(:, jsph), &
1095 & -params % kappa, diff1_grad(:, 1, jsph), one, &
1096 & diff_ep_dim3(1, icav), work_complex, work)
1097 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1098 & constants % vscales, constants % SK_ri(:, jsph), &
1099 & -params % kappa, diff1_grad(:, 2, jsph), one, &
1100 & diff_ep_dim3(2, icav), work_complex, work)
1101 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1102 & constants % vscales, constants % SK_ri(:, jsph), &
1103 & -params % kappa, diff1_grad(:, 3, jsph), one, &
1104 & diff_ep_dim3(3, icav), work_complex, work)
1105 end do
1106 end do
1107 else
1108 knode = constants % snode(ksph)
1109 do igrid = 1, params % ngrid
1110 if (constants % ui(igrid, ksph) .eq. zero) cycle
1111 icav = icav + 1
1112 ! Far-field
1113 call dgemv('T', (params % pl+2)**2, 3, -params % kappa, &
1114 & l2l_grad(1, 1, ksph), &
1115 & (params % pl+2)**2, constants % vgrid(1, igrid), 1, &
1116 & one, diff_ep_dim3(1, icav), 1)
1117 ! Near-field
1118 do knear = constants % snear(knode), constants % snear(knode+1)-1
1119 jnode = constants % near(knear)
1120 do jsph_node = constants % cluster(1, jnode), &
1121 & constants % cluster(2, jnode)
1122 jsph = constants % order(jsph_node)
1123 if (jsph .eq. ksph) cycle
1124 vij = params % csph(:, ksph) - params % csph(:, jsph) + &
1125 & params % rsph(ksph)*constants % cgrid(:, igrid)
1126 vtij = vij * params % kappa
1127 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1128 & constants % vscales, constants % SK_ri(:, jsph), &
1129 & -params % kappa, diff1_grad(:, 1, jsph), one, &
1130 & diff_ep_dim3(1, icav), work_complex, work)
1131 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1132 & constants % vscales, constants % SK_ri(:, jsph), &
1133 & -params % kappa, diff1_grad(:, 2, jsph), one, &
1134 & diff_ep_dim3(2, icav), work_complex, work)
1135 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1136 & constants % vscales, constants % SK_ri(:, jsph), &
1137 & -params % kappa, diff1_grad(:, 3, jsph), one, &
1138 & diff_ep_dim3(3, icav), work_complex, work)
1139 end do
1140 end do
1141 end do
1142 end if
1143 sum_dim3 = zero
1144 icav = constants % icav_ia(ksph) - 1
1145 do igrid =1, params % ngrid
1146 if(constants % ui(igrid, ksph) .gt. zero) then
1147 icav = icav + 1
1148 do ind = 1, constants % nbasis
1149 sum_dim3(:,ind,ksph) = sum_dim3(:,ind,ksph) &
1150 & + diff_ep_dim3(:,icav)*constants % ui(igrid, ksph) &
1151 & *constants % vwgrid(ind, igrid)
1152 end do
1153 end if
1154 end do
1155 do ind = 1, constants % nbasis
1156 force(:, ksph) = force(:, ksph) + &
1157 & sum_dim3(:, ind, ksph)*(xadj_r(ind, ksph) + &
1158 & xadj_e(ind, ksph))
1159 end do
1160
1161 ! Now jsph=ksph and isph!=ksph
1162 if (params % fmm .eq. 0) then
1163 diff_ep_dim3 = zero
1164 sum_dim3 = zero
1165 do isph = 1, params % nsph
1166 if (isph .eq. ksph) cycle
1167 icav = constants % icav_ia(isph) - 1
1168 do igrid = 1, params % ngrid
1169 if (constants % ui(igrid, isph) .eq. zero) cycle
1170 icav = icav + 1
1171 vij = params % csph(:,isph) + &
1172 & params % rsph(isph)*constants % cgrid(:,igrid) - &
1173 & params % csph(:,ksph)
1174 vtij = vij * params % kappa
1175 !call fmm_m2p_bessel_grad(vij * params % kappa, &
1176 ! & params % rsph(ksph)*params % kappa, &
1177 ! & constants % lmax0, &
1178 ! & constants % vscales, -params % kappa, diff1(:, ksph), one, &
1179 ! & diff_ep_dim3(:, icav))
1180 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1181 & constants % vscales, constants % SK_ri(:, ksph), &
1182 & params % kappa, diff1_grad(:, 1, ksph), one, &
1183 & diff_ep_dim3(1, icav), work_complex, work)
1184 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1185 & constants % vscales, constants % SK_ri(:, ksph), &
1186 & params % kappa, diff1_grad(:, 2, ksph), one, &
1187 & diff_ep_dim3(2, icav), work_complex, work)
1188 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1189 & constants % vscales, constants % SK_ri(:, ksph), &
1190 & params % kappa, diff1_grad(:, 3, ksph), one, &
1191 & diff_ep_dim3(3, icav), work_complex, work)
1192 end do
1193 end do
1194 icav = zero
1195 do isph = 1, params % nsph
1196 do igrid =1, params % ngrid
1197 if(constants % ui(igrid, isph) .gt. zero) then
1198 icav = icav + 1
1199 do ind = 1, constants % nbasis
1200 sum_dim3(:,ind,isph) = sum_dim3(:,ind,isph) &
1201 & + diff_ep_dim3(:,icav)*constants % ui(igrid, isph) &
1202 & *constants % vwgrid(ind, igrid)
1203 end do
1204 end if
1205 end do
1206 end do
1207 ! Computation of derivative of \bf(k)_j^l0(x_in)\times Y^j_l0m0(x_in)
1208 do isph = 1, params % nsph
1209 do ind = 1, constants % nbasis
1210 force(:, ksph) = force(:, ksph) &
1211 & + sum_dim3(:, ind, isph)*(xadj_r(ind, isph) &
1212 & + xadj_e(ind, isph))
1213 end do
1214 end do
1215 else
1216 call dgemv('T', (constants % lmax0+2)**2, 3, params % kappa, &
1217 & diff1_grad(1, 1, ksph), (constants % lmax0+2)**2, &
1218 & workspace % tmp_sph2(1, ksph), 1, one, force(1, ksph), 1)
1219 end if
1220 end do
1221 deallocate(diff_ep_dim3, phi_in, sum_dim3, diff1_grad, l2l_grad, &
1222 & diff0, diff1, stat=istat)
1223 if (istat.ne.0) then
1224 call update_error(ddx_error, "deallocation ddx_error in ddx_contract_grad_C_worker2")
1225 return
1226 end if
1227end subroutine contract_grad_c_worker2
1228
1242subroutine contract_grad_f_worker1(params, constants, workspace, sol_adj, sol_sgrid, &
1243 & gradpsi, force, ddx_error)
1244 ! input/output
1245 type(ddx_params_type), intent(in) :: params
1246 type(ddx_constants_type), intent(in) :: constants
1247 type(ddx_workspace_type), intent(inout) :: workspace
1248 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: &
1249 & sol_adj
1250 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: sol_sgrid
1251 real(dp), dimension(3, constants % ncav), intent(in) :: gradpsi
1252 real(dp), dimension(3, params % nsph), intent(inout) :: force
1253 type(ddx_error_type), intent(inout) :: ddx_error
1254
1255 ! local
1256 real(dp), external :: dnrm2
1257 integer :: isph, jsph, igrid, ind, l0, ind0, icav, ksph, &
1258 & knode, jnode, knear, jsph_node, indl, inode, istat
1259 ! val_dim3 : Intermediate value array of dimension 3
1260 real(dp), dimension(3) :: vij, vtij
1261 ! val : Intermediate variable to compute diff_ep
1262 ! nderpsi : Derivative of psi on grid points
1263 real(dp) :: nderpsi, sum_int
1264
1265 ! local allocatable
1266 ! phi_in : sum_{j=1}^N diff0_j * coefY_j
1267 ! diff_ep_dim3 : 3 dimensional couterpart of diff_ep
1268 ! sum_dim3 : Storage of sum
1269 ! Debug purpose
1270 ! These variables can be taken from the subroutine update_rhs
1271 ! diff0 : dot_product([PU_j]_l0m0^l'm', l'/r_j[Xr]_jl'm' -
1272 ! (i'_l'(r_j)/i_l'(r_j))[Xe]_jl'm')
1273 ! sum_Sjin : \sum_j [S]_{jin} Eq.~(97) [QSM20.SISC]
1274 ! c0 : \sum_{n=1}^N_g w_n U_j^{x_nj}\partial_n psi_0(x_nj)Y_{l0m0}(s_n)
1275 real(dp), allocatable :: phi_in(:,:), diff_ep_dim3(:,:), &
1276 & sum_dim3(:,:,:), diff0(:,:), sum_sjin(:,:), &
1277 & c0_d(:,:), c0_d1(:,:), c0_d1_grad(:,:,:), l2l_grad(:,:,:)
1278 real(dp), dimension(0:params % lmax) :: SK_rijn, DK_rijn
1279 complex(dp) :: work_complex(constants % lmax0 + 2)
1280 real(dp) :: work(constants % lmax0 + 2)
1281
1282 allocate(phi_in(params % ngrid, params % nsph), &
1283 & diff_ep_dim3(3, constants % ncav), &
1284 & sum_dim3(3, constants % nbasis, params % nsph), &
1285 & c0_d(constants % nbasis0, params % nsph), &
1286 & c0_d1(constants % nbasis0, params % nsph), &
1287 & c0_d1_grad((constants % lmax0+2)**2, 3, params % nsph), &
1288 & sum_sjin(params % ngrid, params % nsph), &
1289 & l2l_grad((params % pl+2)**2, 3, params % nsph), &
1290 & diff0(constants % nbasis0, params % nsph), stat=istat)
1291 if (istat.ne.0) then
1292 call update_error(ddx_error, "allocation ddx_error in ddx_grad_f_worker1")
1293 return
1294 end if
1295
1296
1297 ! Setting initial values to zero
1298 sk_rijn = zero
1299 dk_rijn = zero
1300 c0_d = zero
1301 c0_d1 = zero
1302 c0_d1_grad = zero
1303
1304 icav = zero
1305 do isph = 1, params % nsph
1306 do igrid= 1, params % ngrid
1307 if ( constants % ui(igrid,isph) .gt. zero ) then
1308 icav = icav + 1
1309 nderpsi = dot_product( gradpsi(:,icav),constants % cgrid(:,igrid) )
1310 c0_d(:, isph) = c0_d(:,isph) + &
1311 & constants % wgrid(igrid)* &
1312 & constants % ui(igrid,isph)*&
1313 & nderpsi* &
1314 & constants % vgrid(1:constants % nbasis0,igrid)
1315 do l0 = 0, constants % lmax0
1316 ind0 = l0*l0 + l0 + 1
1317 c0_d1(ind0-l0:ind0+l0, isph) = c0_d(ind0-l0:ind0+l0, isph) * &
1318 & constants % C_ik(l0, isph)
1319 end do
1320 ! Prepare c0_d1_grad
1321 call fmm_m2m_bessel_grad(constants % lmax0, constants % SK_ri(:, isph), &
1322 & constants % vscales, c0_d1(:, isph), c0_d1_grad(:, :, isph))
1323 end if
1324 end do
1325 end do
1326
1327 if (params % fmm .eq. 0) then
1328 ! Compute [S]_{jin}
1329 icav = 0
1330 do isph = 1, params % nsph
1331 do igrid = 1, params % ngrid
1332 if (constants % ui(igrid,isph).gt.zero) then
1333 icav = icav + 1
1334 sum_int = zero
1335 ! Loop to compute Sijn
1336 do jsph = 1, params % nsph
1337 vij = params % csph(:, isph) + &
1338 & params % rsph(isph)*constants % cgrid(:, igrid) - &
1339 & params % csph(:, jsph)
1340 vtij = vij * params % kappa
1341 call fmm_m2p_bessel_work(vtij, constants % lmax0, &
1342 & constants % vscales, constants % SK_ri(:, jsph), one, &
1343 & c0_d1(:, jsph), one, sum_int, work_complex, work)
1344 end do ! End of loop jsph
1345 sum_sjin(igrid,isph) = -(params % epsp/params % eps)*sum_int
1346 end if
1347 end do ! End of loop igrid
1348 end do ! End of loop isph
1349 else
1350 ! Load input harmonics into tree data
1351 workspace % tmp_sph = zero
1352 workspace % tmp_sph(1:constants % nbasis0, :) = c0_d1(:, :)
1353 if(constants % lmax0 .lt. params % pm) then
1354 do isph = 1, params % nsph
1355 inode = constants % snode(isph)
1356 workspace % tmp_node_m(1:constants % nbasis0, inode) = &
1357 & workspace % tmp_sph(1:constants % nbasis0, isph)
1358 workspace % tmp_node_m(constants % nbasis0+1:, inode) = zero
1359 end do
1360 else
1361 indl = (params % pm+1)**2
1362 do isph = 1, params % nsph
1363 inode = constants % snode(isph)
1364 workspace % tmp_node_m(:, inode) = workspace % tmp_sph(1:indl, isph)
1365 end do
1366 end if
1367 ! Do FMM operations
1368 call tree_m2m_bessel_rotation(params, constants, workspace % tmp_node_m)
1369 call tree_m2l_bessel_rotation(params, constants, workspace % tmp_node_m, &
1370 & workspace % tmp_node_l)
1371 call tree_l2l_bessel_rotation(params, constants, workspace % tmp_node_l)
1372 call tree_l2p_bessel(params, constants, -params % epsp/params % eps, workspace % tmp_node_l, zero, &
1373 & sum_sjin)
1374 call tree_m2p_bessel(params, constants, constants % lmax0, -params % epsp/params % eps, &
1375 & params % lmax, workspace % tmp_sph, one, &
1376 & sum_sjin)
1377 ! Make phi_in zero at internal grid points
1378 do isph = 1, params % nsph
1379 do igrid = 1, params % ngrid
1380 if (constants % ui(igrid, isph) .eq. zero) then
1381 sum_sjin(igrid, isph) = zero
1382 end if
1383 end do
1384 end do
1385 ! Get gradients of the L2L
1386 do isph = 1, params % nsph
1387 inode = constants % snode(isph)
1388 workspace % tmp_sph_l(:, isph) = workspace % tmp_node_l(:, inode)
1389 call fmm_l2l_bessel_grad(params % pl, &
1390 & constants % SI_ri(:, isph), constants % vscales, &
1391 & workspace % tmp_node_l(:, inode), &
1392 & l2l_grad(:, :, isph))
1393 end do
1394 workspace % tmp_sph = sol_adj
1395 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis, &
1396 & one, constants % vwgrid, constants % vgrid_nbasis, &
1397 & workspace % tmp_sph, constants % nbasis, zero, &
1398 & workspace % tmp_grid, params % ngrid)
1399 workspace % tmp_grid = workspace % tmp_grid * constants % ui
1400 ! Adjoint FMM with output tmp_sph2(:, :) which stores coefficients of
1401 ! harmonics of degree up to lmax+1
1402 call tree_m2p_bessel_nodiag_adj(params, constants, constants % lmax0+1, one, &
1403 & workspace % tmp_grid, zero, params % lmax+1, workspace % tmp_sph2)
1404 call tree_l2p_bessel_adj(params, constants, one, workspace % tmp_grid, zero, &
1405 & workspace % tmp_node_l)
1406 call tree_l2l_bessel_rotation_adj(params, constants, workspace % tmp_node_l)
1407 call tree_m2l_bessel_rotation_adj(params, constants, workspace % tmp_node_l, &
1408 & workspace % tmp_node_m)
1409 call tree_m2m_bessel_rotation_adj(params, constants, workspace % tmp_node_m)
1410 ! Properly load adjoint multipole harmonics into tmp_sph2 that holds
1411 ! harmonics of a degree up to lmax+1
1412 if(constants % lmax0+1 .lt. params % pm) then
1413 do isph = 1, params % nsph
1414 inode = constants % snode(isph)
1415 workspace % tmp_sph2(1:(constants % lmax0+2)**2, isph) = &
1416 & workspace % tmp_sph2(1:(constants % lmax0+2)**2, isph) + &
1417 & workspace % tmp_node_m(1:(constants % lmax0+2)**2, inode)
1418 end do
1419 else
1420 indl = (params % pm+1)**2
1421 do isph = 1, params % nsph
1422 inode = constants % snode(isph)
1423 workspace % tmp_sph2(1:indl, isph) = &
1424 & workspace % tmp_sph2(1:indl, isph) + &
1425 & workspace % tmp_node_m(:, inode)
1426 end do
1427 end if
1428 end if
1429 do ksph = 1, params % nsph
1430 ! Computation of derivative of U_i^e(x_in)
1431 call contract_grad_u(params, constants, ksph, sol_sgrid, sum_sjin, force(:, ksph))
1432 ! Aleksandr: my loop for the diff_ep_dim3
1433 diff_ep_dim3 = zero
1434 ! At first isph=ksph, jsph!=ksph
1435 icav = constants % icav_ia(ksph) - 1
1436 if (params % fmm .eq. 0) then
1437 do igrid = 1, params % ngrid
1438 if (constants % ui(igrid, ksph) .eq. zero) cycle
1439 icav = icav + 1
1440 do jsph = 1, params % nsph
1441 if (jsph .eq. ksph) cycle
1442 vij = params % csph(:,ksph) + &
1443 & params % rsph(ksph)*constants % cgrid(:,igrid) - &
1444 & params % csph(:,jsph)
1445 vtij = vij * params % kappa
1446 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1447 & constants % vscales, constants % SK_ri(:, jsph), &
1448 & -params % kappa, c0_d1_grad(:, 1, jsph), one, &
1449 & diff_ep_dim3(1, icav), work_complex, work)
1450 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1451 & constants % vscales, constants % SK_ri(:, jsph), &
1452 & -params % kappa, c0_d1_grad(:, 2, jsph), one, &
1453 & diff_ep_dim3(2, icav), work_complex, work)
1454 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1455 & constants % vscales, constants % SK_ri(:, jsph), &
1456 & -params % kappa, c0_d1_grad(:, 3, jsph), one, &
1457 & diff_ep_dim3(3, icav), work_complex, work)
1458 end do
1459 end do
1460 else
1461 knode = constants % snode(ksph)
1462 do igrid = 1, params % ngrid
1463 if (constants % ui(igrid, ksph) .eq. zero) cycle
1464 icav = icav + 1
1465 ! Far-field
1466 call dgemv('T', (params % pl+2)**2, 3, -params % kappa, &
1467 & l2l_grad(1, 1, ksph), &
1468 & (params % pl+2)**2, constants % vgrid(1, igrid), 1, &
1469 & one, diff_ep_dim3(1, icav), 1)
1470 ! Near-field
1471 do knear = constants % snear(knode), constants % snear(knode+1)-1
1472 jnode = constants % near(knear)
1473 do jsph_node = constants % cluster(1, jnode), &
1474 & constants % cluster(2, jnode)
1475 jsph = constants % order(jsph_node)
1476 if (jsph .eq. ksph) cycle
1477 vij = params % csph(:, ksph) - params % csph(:, jsph) + &
1478 & params % rsph(ksph)*constants % cgrid(:, igrid)
1479 vtij = vij * params % kappa
1480 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1481 & constants % vscales, constants % SK_ri(:, jsph), &
1482 & -params % kappa, c0_d1_grad(:, 1, jsph), one, &
1483 & diff_ep_dim3(1, icav), work_complex, work)
1484 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1485 & constants % vscales, constants % SK_ri(:, jsph), &
1486 & -params % kappa, c0_d1_grad(:, 2, jsph), one, &
1487 & diff_ep_dim3(2, icav), work_complex, work)
1488 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1489 & constants % vscales, constants % SK_ri(:, jsph), &
1490 & -params % kappa, c0_d1_grad(:, 3, jsph), one, &
1491 & diff_ep_dim3(3, icav), work_complex, work)
1492 end do
1493 end do
1494 end do
1495 end if
1496 sum_dim3 = zero
1497 icav = constants % icav_ia(ksph) - 1
1498 do igrid =1, params % ngrid
1499 if(constants % ui(igrid, ksph) .gt. zero) then
1500 icav = icav + 1
1501 do ind = 1, constants % nbasis
1502 sum_dim3(:,ind,ksph) = sum_dim3(:,ind,ksph) &
1503 & - (params % epsp/params % eps) &
1504 & *diff_ep_dim3(:,icav)*constants % ui(igrid, ksph) &
1505 & *constants % vwgrid(ind, igrid)
1506 end do
1507 end if
1508 end do
1509 do ind = 1, constants % nbasis
1510 force(:, ksph) = force(:, ksph) + &
1511 & sum_dim3(:, ind, ksph)*sol_adj(ind, ksph)
1512 end do
1513
1514 ! Now jsph=ksph and isph!=ksph
1515 if (params % fmm .eq. 0) then
1516 diff_ep_dim3 = zero
1517 sum_dim3 = zero
1518 do isph = 1, params % nsph
1519 if (isph .eq. ksph) cycle
1520 icav = constants % icav_ia(isph) - 1
1521 do igrid = 1, params % ngrid
1522 if (constants % ui(igrid, isph) .eq. zero) cycle
1523 icav = icav + 1
1524 vij = params % csph(:,isph) + &
1525 & params % rsph(isph)*constants % cgrid(:,igrid) - &
1526 & params % csph(:,ksph)
1527 vtij = vij * params % kappa
1528 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1529 & constants % vscales, constants % SK_ri(:, ksph), &
1530 & params % kappa, c0_d1_grad(:, 1, ksph), one, &
1531 & diff_ep_dim3(1, icav), work_complex, work)
1532 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1533 & constants % vscales, constants % SK_ri(:, ksph), &
1534 & params % kappa, c0_d1_grad(:, 2, ksph), one, &
1535 & diff_ep_dim3(2, icav), work_complex, work)
1536 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1537 & constants % vscales, constants % SK_ri(:, ksph), &
1538 & params % kappa, c0_d1_grad(:, 3, ksph), one, &
1539 & diff_ep_dim3(3, icav), work_complex, work)
1540 end do
1541 end do
1542
1543 icav = zero
1544 do isph = 1, params % nsph
1545 do igrid =1, params % ngrid
1546 if(constants % ui(igrid, isph) .gt. zero) then
1547 icav = icav + 1
1548 do ind = 1, constants % nbasis
1549 sum_dim3(:,ind,isph) = sum_dim3(:,ind,isph) &
1550 & - (params % epsp/params % eps) &
1551 & *diff_ep_dim3(:,icav)*constants % ui(igrid, isph) &
1552 & *constants % vwgrid(ind, igrid)
1553 end do
1554 end if
1555 end do
1556 end do
1557
1558 ! Computation of derivative of \bf(k)_j^l0(x_in)\times Y^j_l0m0(x_in)
1559 do isph = 1, params % nsph
1560 do ind = 1, constants % nbasis
1561 force(:, ksph) = force(:, ksph) + sum_dim3(:, ind, isph)*sol_adj(ind, isph)
1562 end do
1563 end do
1564 else
1565 call dgemv('T', (constants % lmax0+2)**2, 3, -params % epsp/params % eps*params % kappa, &
1566 & c0_d1_grad(1, 1, ksph), (constants % lmax0+2)**2, &
1567 & workspace % tmp_sph2(1, ksph), 1, one, force(1, ksph), 1)
1568 end if
1569 end do
1570
1571 deallocate(phi_in, diff_ep_dim3, sum_dim3, c0_d, c0_d1, &
1572 & c0_d1_grad, sum_sjin, l2l_grad, diff0, stat=istat)
1573 if (istat.ne.0) then
1574 call update_error(ddx_error, "deallocation ddx_error in ddx_grad_f_worker1")
1575 return
1576 end if
1577
1578end subroutine contract_grad_f_worker1
1579
1580!! @param[inout] ddx_error: ddX error
1581subroutine contract_grad_f_worker2(params, constants, &
1582 & gradpsi, normal_hessian_cav, force, state, ddx_error)
1583 type(ddx_params_type), intent(in) :: params
1584 type(ddx_constants_type), intent(in) :: constants
1585 real(dp), intent(in) :: gradpsi(3, constants % ncav), &
1586 & normal_hessian_cav(3, constants % ncav)
1587 real(dp), intent(inout) :: force(3, params % nsph)
1588 type(ddx_state_type), intent(inout) :: state
1589 type(ddx_error_type), intent(inout) :: ddx_error
1590
1591 integer :: icav, isph, igrid, istat
1592 real(dp) :: nderpsi
1593 real(dp), allocatable :: gradpsi_grid(:, :)
1594
1595 allocate(gradpsi_grid(params % ngrid, params % nsph), stat=istat)
1596 if (istat.ne.0) then
1597 call update_error(ddx_error, "allocation ddx_error in ddx_grad_f_worker2")
1598 return
1599 end if
1600
1601 icav = 0
1602 do isph = 1, params % nsph
1603 do igrid = 1, params % ngrid
1604 if(constants % ui(igrid, isph) .gt. zero) then
1605 icav = icav + 1
1606 nderpsi = dot_product(gradpsi(:, icav), &
1607 & constants % cgrid(:, igrid))
1608 gradpsi_grid(igrid, isph) = nderpsi
1609 end if
1610 end do
1611 end do
1612
1613 icav = 0
1614 do isph = 1, params % nsph
1615 call contract_grad_u(params, constants, isph, gradpsi_grid, &
1616 & state % phi_n, force(:, isph))
1617 do igrid = 1, params % ngrid
1618 if(constants % ui(igrid, isph) .gt. zero) then
1619 icav = icav + 1
1620 force(:, isph) = force(:, isph) &
1621 & + constants % wgrid(igrid)*constants % ui(igrid, isph) &
1622 & *state % phi_n(igrid, isph)*normal_hessian_cav(:, icav)
1623 end if
1624 end do
1625 end do
1626
1627 deallocate(gradpsi_grid, stat=istat)
1628 if (istat.ne.0) then
1629 call update_error(ddx_error, "deallocation ddx_error in ddx_grad_f_worker2")
1630 return
1631 end if
1632
1633end subroutine contract_grad_f_worker2
1634
1636subroutine gradr_sph(params, constants, isph, vplm, vcos, vsin, basloc, &
1637 & dbsloc, g, ygrid, fx)
1638 implicit none
1639 ! Inputs
1640 type(ddx_params_type), intent(in) :: params
1641 type(ddx_constants_type), intent(in) :: constants
1642 integer, intent(in) :: isph
1643 real(dp), intent(in) :: g(constants % nbasis, params % nsph), &
1644 & ygrid(params % ngrid, params % nsph)
1645 real(dp), intent(inout) :: vplm(constants % nbasis), vcos(params % lmax+1), &
1646 & vsin(params % lmax+1), basloc(constants % nbasis), &
1647 & dbsloc(3, constants % nbasis), fx(3)
1648 ! various scratch arrays
1649 real(dp) vik(3), sik(3), vki(3), ski(3), vkj(3), skj(3), vji(3), &
1650 & sji(3), va(3), vb(3), a(3)
1651 ! jacobian matrix
1652 real(dp) sjac(3,3)
1653 ! indexes
1654 integer its, ik, ksph, l, m, ind, jsph, icomp, jcomp
1655 ! various scalar quantities
1656 real(dp) cx, cy, cz, vvki, tki, gg, fl, fac, vvkj, tkj
1657 real(dp) tt, fcl, fjj, gi, fii, vvji, tji, qji
1658 real(dp) b, vvik, tik, qik, tlow, thigh, duj
1659 real(dp) :: rho, ctheta, stheta, cphi, sphi
1660 real(dp), external :: dnrm2
1661
1662 tlow = one - pt5*(one - params % se)*params % eta
1663 thigh = one + pt5*(one + params % se)*params % eta
1664
1665 ! first set of contributions:
1666 ! diagonal block, kc and part of kb
1667
1668 fx = zero
1669 do its = 1, params % ngrid
1670 ! sum over ksph in neighbors of isph
1671 do ik = constants % inl(isph), constants % inl(isph+1) - 1
1672 ksph = constants % nl(ik)
1673 ! build geometrical quantities
1674 cx = params % csph(1,ksph) + params % rsph(ksph)*constants % cgrid(1,its)
1675 cy = params % csph(2,ksph) + params % rsph(ksph)*constants % cgrid(2,its)
1676 cz = params % csph(3,ksph) + params % rsph(ksph)*constants % cgrid(3,its)
1677 vki(1) = cx - params % csph(1,isph)
1678 vki(2) = cy - params % csph(2,isph)
1679 vki(3) = cz - params % csph(3,isph)
1680 !vvki = sqrt(vki(1)*vki(1) + vki(2)*vki(2) + &
1681 ! & vki(3)*vki(3))
1682 vvki = dnrm2(3, vki, 1)
1683 tki = vvki/params % rsph(isph)
1684
1685 ! contributions involving grad i of uk come from the switching
1686 ! region.
1687 ! note: ui avoids contributions from points that are in the
1688 ! switching between isph and ksph but are buried in a third
1689 ! sphere.
1690 if ((tki.gt.tlow).and.(tki.lt.thigh) .and. &
1691 & constants % ui(its,ksph).gt.zero) then
1692 ! other geometrical quantities
1693 ski = vki/vvki
1694
1695 ! diagonal block kk contribution, with k in n(i)
1696 gg = zero
1697 do l = 0, params % lmax
1698 ind = l*l + l + 1
1699 fl = dble(l)
1700 fac = twopi/(two*fl + one)
1701 do m = -l, l
1702 !! DEBUG comment
1703 gg = gg + fac*constants % vgrid(ind+m,its)*g(ind+m,ksph)
1704 end do
1705 end do
1706
1707 ! kc contribution
1708 do jsph = 1, params % nsph
1709 if (jsph.ne.ksph .and. jsph.ne.isph) then
1710 vkj(1) = cx - params % csph(1,jsph)
1711 vkj(2) = cy - params % csph(2,jsph)
1712 vkj(3) = cz - params % csph(3,jsph)
1713 vvkj = sqrt(vkj(1)*vkj(1) + vkj(2)*vkj(2) + &
1714 & vkj(3)*vkj(3))
1715 vvkj = dnrm2(3, vkj, 1)
1716 tkj = vvkj/params % rsph(jsph)
1717 skj = vkj/vvkj
1718 call ylmbas(skj, rho, ctheta, stheta, cphi, sphi, &
1719 & params % lmax, constants % vscales, basloc, &
1720 & vplm, vcos, vsin)
1721 tt = one/tkj
1722 do l = 0, params % lmax
1723 ind = l*l + l + 1
1724 fcl = - fourpi*dble(l)/(two*dble(l)+one)*tt
1725 do m = -l, l
1726 !! DEBUG comment
1727 gg = gg + fcl*g(ind+m,jsph)*basloc(ind+m)
1728 end do
1729 tt = tt/tkj
1730 end do
1731 !call fmm_m2p(vkj, params % rsph(jsph), &
1732 ! & params % lmax, constants % vscales_rel, -one, &
1733 ! & g(:, jsph), one, gg)
1734 end if
1735 end do
1736
1737 ! part of kb contribution
1738 call ylmbas(ski, rho, ctheta, stheta, cphi, sphi, &
1739 & params % lmax, constants % vscales, basloc, &
1740 & vplm, vcos, vsin)
1741 tt = one/tki
1742 do l = 0, params % lmax
1743 ind = l*l + l + 1
1744 fcl = - four*pi*dble(l)/(two*dble(l)+one)*tt
1745 do m = -l, l
1746 !! DEBUG comment
1747 gg = gg + fcl*g(ind+m,isph)*basloc(ind+m)
1748 end do
1749 tt = tt/tki
1750 end do
1751 !call fmm_m2p(vki, params % rsph(isph), &
1752 ! & params % lmax, constants % vscales_rel, -one, &
1753 ! & g(:, isph), one, gg)
1754
1755 ! common step, product with grad i uj
1756 duj = dfsw(tki,params % se, params % eta)/params % rsph(isph)
1757 fjj = duj*constants % wgrid(its)*gg*ygrid(its,ksph)
1758 fx(1) = fx(1) - fjj*ski(1)
1759 fx(2) = fx(2) - fjj*ski(2)
1760 fx(3) = fx(3) - fjj*ski(3)
1761 end if
1762 end do
1763
1764 ! diagonal block ii contribution
1765 if (constants % ui(its,isph).gt.zero.and.constants % ui(its,isph).lt.one) then
1766 gi = zero
1767 do l = 0, params % lmax
1768 ind = l*l + l + 1
1769 fl = dble(l)
1770 fac = twopi/(two*fl + one)
1771 do m = -l, l
1772 !! DEBUG comment
1773 gi = gi + fac*constants % vgrid(ind+m,its)*g(ind+m,isph)
1774 !gi = gi + pt5*constants % vgrid2(ind+m,its)*g(ind+m,isph)
1775 end do
1776 end do
1777 !do l = 0, (params % lmax+1)**2
1778 ! gi = gi + constants % vgrid2(l, its)*g(l, isph)
1779 !end do
1780 !gi = pt5 * gi
1781 fii = constants % wgrid(its)*gi*ygrid(its,isph)
1782 fx(1) = fx(1) + fii*constants % zi(1,its,isph)
1783 fx(2) = fx(2) + fii*constants % zi(2,its,isph)
1784 fx(3) = fx(3) + fii*constants % zi(3,its,isph)
1785 end if
1786 end do
1787
1788 ! second set of contributions:
1789 ! part of kb and ka
1790 do its = 1, params % ngrid
1791
1792 ! run over all the spheres except isph
1793 do jsph = 1, params % nsph
1794 if (constants % ui(its,jsph).gt.zero .and. jsph.ne.isph) then
1795 ! build geometrical quantities
1796 cx = params % csph(1,jsph) + params % rsph(jsph)*constants % cgrid(1,its)
1797 cy = params % csph(2,jsph) + params % rsph(jsph)*constants % cgrid(2,its)
1798 cz = params % csph(3,jsph) + params % rsph(jsph)*constants % cgrid(3,its)
1799 vji(1) = cx - params % csph(1,isph)
1800 vji(2) = cy - params % csph(2,isph)
1801 vji(3) = cz - params % csph(3,isph)
1802 !vvji = sqrt(vji(1)*vji(1) + vji(2)*vji(2) + &
1803 ! & vji(3)*vji(3))
1804 vvji = dnrm2(3, vji, 1)
1805 tji = vvji/params % rsph(isph)
1806 qji = one/vvji
1807 sji = vji/vvji
1808
1809 ! build the jacobian of sji
1810 sjac = zero
1811 sjac(1,1) = - one
1812 sjac(2,2) = - one
1813 sjac(3,3) = - one
1814 do icomp = 1, 3
1815 do jcomp = 1, 3
1816 sjac(icomp,jcomp) = qji*(sjac(icomp,jcomp) &
1817 & + sji(icomp)*sji(jcomp))
1818 end do
1819 end do
1820
1821 ! assemble the local basis and its gradient
1822 !call dbasis(sji,basloc,dbsloc,vplm,vcos,vsin)
1823 call dbasis(params, constants, sji,basloc,dbsloc,vplm,vcos,vsin)
1824
1825 ! assemble the contribution
1826 a = zero
1827 tt = one/(tji)
1828 do l = 0, params % lmax
1829 ind = l*l + l + 1
1830 fl = dble(l)
1831 fcl = - tt*fourpi*fl/(two*fl + one)
1832 do m = -l, l
1833 fac = fcl*g(ind+m,isph)
1834 b = (fl + one)*basloc(ind+m)/(params % rsph(isph)*tji)
1835
1836 ! apply the jacobian to grad y
1837 va(1) = sjac(1,1)*dbsloc(1,ind+m) + &
1838 & sjac(1,2)*dbsloc(2,ind+m) + sjac(1,3)*dbsloc(3,ind+m)
1839 va(2) = sjac(2,1)*dbsloc(1,ind+m) + &
1840 & sjac(2,2)*dbsloc(2,ind+m) + sjac(2,3)*dbsloc(3,ind+m)
1841 va(3) = sjac(3,1)*dbsloc(1,ind+m) + &
1842 & sjac(3,2)*dbsloc(2,ind+m) + sjac(3,3)*dbsloc(3,ind+m)
1843 a(1) = a(1) + fac*(sji(1)*b + va(1))
1844 a(2) = a(2) + fac*(sji(2)*b + va(2))
1845 a(3) = a(3) + fac*(sji(3)*b + va(3))
1846 end do
1847 tt = tt/tji
1848 end do
1849 fac = constants % ui(its,jsph)*constants % wgrid(its)*ygrid(its,jsph)
1850 fx(1) = fx(1) - fac*a(1)
1851 fx(2) = fx(2) - fac*a(2)
1852 fx(3) = fx(3) - fac*a(3)
1853 end if
1854 end do
1855 end do
1856
1857 ! ka contribution
1858 do its = 1, params % ngrid
1859 cx = params % csph(1,isph) + params % rsph(isph)*constants % cgrid(1,its)
1860 cy = params % csph(2,isph) + params % rsph(isph)*constants % cgrid(2,its)
1861 cz = params % csph(3,isph) + params % rsph(isph)*constants % cgrid(3,its)
1862 a = zero
1863
1864 ! iterate on all the spheres except isph
1865 do ksph = 1, params % nsph
1866 if (constants % ui(its,isph).gt.zero .and. ksph.ne.isph) then
1867 ! geometrical stuff
1868 vik(1) = cx - params % csph(1,ksph)
1869 vik(2) = cy - params % csph(2,ksph)
1870 vik(3) = cz - params % csph(3,ksph)
1871 !vvik = sqrt(vik(1)*vik(1) + vik(2)*vik(2) + &
1872 ! & vik(3)*vik(3))
1873 vvik = dnrm2(3, vik, 1)
1874 tik = vvik/params % rsph(ksph)
1875 qik = one/vvik
1876 sik = vik/vvik
1877
1878 ! build the jacobian of sik
1879 sjac = zero
1880 sjac(1,1) = one
1881 sjac(2,2) = one
1882 sjac(3,3) = one
1883 do icomp = 1, 3
1884 do jcomp = 1, 3
1885 sjac(icomp,jcomp) = qik*(sjac(icomp,jcomp) &
1886 & - sik(icomp)*sik(jcomp))
1887 end do
1888 end do
1889
1890 ! if we are in the switching region, recover grad_i u_i
1891 vb = zero
1892 if (constants % ui(its,isph).lt.one) then
1893 vb(1) = constants % zi(1,its,isph)
1894 vb(2) = constants % zi(2,its,isph)
1895 vb(3) = constants % zi(3,its,isph)
1896 end if
1897
1898 ! assemble the local basis and its gradient
1899 !call dbasis(sik,basloc,dbsloc,vplm,vcos,vsin)
1900 call dbasis(params, constants, sik,basloc,dbsloc,vplm,vcos,vsin)
1901
1902 ! assemble the contribution
1903 tt = one/(tik)
1904 do l = 0, params % lmax
1905 ind = l*l + l + 1
1906 fl = dble(l)
1907 fcl = - tt*fourpi*fl/(two*fl + one)
1908 do m = -l, l
1909 fac = fcl*g(ind+m,ksph)
1910 fac = - fac*basloc(ind+m)
1911 a(1) = a(1) + fac*vb(1)
1912 a(2) = a(2) + fac*vb(2)
1913 a(3) = a(3) + fac*vb(3)
1914
1915 fac = constants % ui(its,isph)*fcl*g(ind+m,ksph)
1916 b = - (fl + one)*basloc(ind+m)/(params % rsph(ksph)*tik)
1917
1918 ! apply the jacobian to grad y
1919 va(1) = sjac(1,1)*dbsloc(1,ind+m) + &
1920 & sjac(1,2)*dbsloc(2,ind+m) + sjac(1,3)*dbsloc(3,ind+m)
1921 va(2) = sjac(2,1)*dbsloc(1,ind+m) + &
1922 & sjac(2,2)*dbsloc(2,ind+m) + sjac(2,3)*dbsloc(3,ind+m)
1923 va(3) = sjac(3,1)*dbsloc(1,ind+m) + &
1924 & sjac(3,2)*dbsloc(2,ind+m) + sjac(3,3)*dbsloc(3,ind+m)
1925 a(1) = a(1) + fac*(sik(1)*b + va(1))
1926 a(2) = a(2) + fac*(sik(2)*b + va(2))
1927 a(3) = a(3) + fac*(sik(3)*b + va(3))
1928 end do
1929 tt = tt/tik
1930 end do
1931 end if
1932 end do
1933 fac = constants % wgrid(its)*ygrid(its,isph)
1934 fx(1) = fx(1) - fac*a(1)
1935 fx(2) = fx(2) - fac*a(2)
1936 fx(3) = fx(3) - fac*a(3)
1937 end do
1938end subroutine gradr_sph
1939
1941subroutine gradr_fmm(params, constants, workspace, g, ygrid, fx)
1942 implicit none
1943 ! Inputs
1944 type(ddx_params_type), intent(in) :: params
1945 type(ddx_constants_type), intent(in) :: constants
1946 real(dp), intent(in) :: g(constants % nbasis, params % nsph), &
1947 & ygrid(params % ngrid, params % nsph)
1948 ! Temporaries
1949 type(ddx_workspace_type), intent(inout) :: workspace
1950 ! Output
1951 real(dp), intent(out) :: fx(3, params % nsph)
1952 ! Local variables
1953 integer :: indl, indl1, l, isph, igrid, ik, ksph, &
1954 & jsph, jsph_node
1955 integer :: inear, inode, jnode
1956 real(dp) :: gg, c(3), vki(3), vvki, tki, gg3(3), tmp_gg, tmp_c(3)
1957 real(dp) :: tlow, thigh
1958 real(dp), dimension(3, 3) :: zx_coord_transform, zy_coord_transform
1959 real(dp), external :: ddot, dnrm2
1960 real(dp) :: work(params % lmax+2)
1961 !real(dp) :: l2g(params % ngrid, params % nsph)
1962 zx_coord_transform = 0
1963 zx_coord_transform(3, 2) = 1
1964 zx_coord_transform(2, 3) = 1
1965 zx_coord_transform(1, 1) = 1
1966 zy_coord_transform = 0
1967 zy_coord_transform(1, 2) = 1
1968 zy_coord_transform(2, 1) = 1
1969 zy_coord_transform(3, 3) = 1
1970 tlow = one - pt5*(one - params % se)*params % eta
1971 thigh = one + pt5*(one + params % se)*params % eta
1972 fx = zero
1973 !! Scale input harmonics at first
1974 workspace % tmp_sph(1, :) = zero
1975 indl = 2
1976 do l = 1, params % lmax
1977 indl1 = (l+1)**2
1978 workspace % tmp_sph(indl:indl1, :) = l * g(indl:indl1, :)
1979 indl = indl1 + 1
1980 end do
1981 !! Compute gradient of M2M of tmp_sph harmonics at the origin and store it
1982 !! in tmp_sph_grad. tmp_sph_grad(:, 1, :), tmp_sph_grad(:, 2, :) and
1983 !! tmp_sph_grad(:, 3, :) correspond to the OX, OY and OZ axes. Variable
1984 !! tmp_sph2 is a temporary workspace here.
1985 call tree_grad_m2m(params, constants, workspace % tmp_sph, &
1986 & workspace % tmp_sph_grad, workspace % tmp_sph2)
1987 !! Adjoint full FMM matvec to get output multipole expansions from input
1988 !! external grid points. It is used to compute R_i^B fast as a contraction
1989 !! of a gradient stored in tmp_sph_grad and a result of adjoint matvec.
1990 ! Adjoint integration from spherical harmonics to grid points is not needed
1991 ! here as ygrid already contains grid values, we just need to scale it by
1992 ! weights of grid points
1993 do isph = 1, params % nsph
1994 workspace % tmp_grid(:, isph) = ygrid(:, isph) * &
1995 & constants % wgrid(:) * constants % ui(:, isph)
1996 end do
1997 ! Adjoint FMM with output tmp_sph2(:, :) which stores coefficients of
1998 ! harmonics of degree up to lmax+1
1999 call tree_m2p_adj(params, constants, params % lmax+1, one, &
2000 & workspace % tmp_grid, zero, workspace % tmp_sph2)
2001 call tree_l2p_adj(params, constants, one, workspace % tmp_grid, zero, &
2002 & workspace % tmp_node_l, workspace % tmp_sph_l)
2003 call tree_l2l_rotation_adj(params, constants, workspace % tmp_node_l)
2004 call tree_m2l_rotation_adj(params, constants, workspace % tmp_node_l, &
2005 & workspace % tmp_node_m)
2006 call tree_m2m_rotation_adj(params, constants, workspace % tmp_node_m)
2007 ! Properly load adjoint multipole harmonics into tmp_sph2 that holds
2008 ! harmonics of a degree up to lmax+1
2009 if(params % lmax+1 .lt. params % pm) then
2010 do isph = 1, params % nsph
2011 inode = constants % snode(isph)
2012 workspace % tmp_sph2(:, isph) = workspace % tmp_sph2(:, isph) + &
2013 & workspace % tmp_node_m(1:constants % grad_nbasis, inode)
2014 end do
2015 else
2016 indl = (params % pm+1)**2
2017 do isph = 1, params % nsph
2018 inode = constants % snode(isph)
2019 workspace % tmp_sph2(1:indl, isph) = &
2020 & workspace % tmp_sph2(1:indl, isph) + &
2021 & workspace % tmp_node_m(:, inode)
2022 end do
2023 end if
2024 ! Compute second term of R_i^B as a contraction
2025 do isph = 1, params % nsph
2026 call dgemv('T', constants % grad_nbasis, 3, one, &
2027 & workspace % tmp_sph_grad(1, 1, isph), constants % grad_nbasis, &
2028 & workspace % tmp_sph2(1, isph), 1, zero, fx(1, isph), 1)
2029 end do
2030 !! Direct far-field FMM matvec to get output local expansions from input
2031 !! multipole expansions. It will be used in R_i^A.
2032 !! As of now I compute potential at all external grid points, improved
2033 !! version shall only compute it at external points in a switch region
2034 ! Load input harmonics into tree data
2035 if(params % lmax .lt. params % pm) then
2036 do isph = 1, params % nsph
2037 inode = constants % snode(isph)
2038 workspace % tmp_node_m(:constants % nbasis, inode) = &
2039 & workspace % tmp_sph(:, isph)
2040 workspace % tmp_node_m(constants % nbasis+1:, inode) = zero
2041 end do
2042 else
2043 indl = (params % pm+1)**2
2044 do isph = 1, params % nsph
2045 inode = constants % snode(isph)
2046 workspace % tmp_node_m(:, inode) = workspace % tmp_sph(1:indl, isph)
2047 end do
2048 end if
2049 ! Perform direct FMM matvec to all external grid points
2050 call tree_m2m_rotation(params, constants, workspace % tmp_node_m)
2051 call tree_m2l_rotation(params, constants, workspace % tmp_node_m, &
2052 & workspace % tmp_node_l)
2053 call tree_l2l_rotation(params, constants, workspace % tmp_node_l)
2054 call tree_l2p(params, constants, one, workspace % tmp_node_l, zero, &
2055 & workspace % tmp_grid, workspace % tmp_sph_l)
2056 call tree_m2p(params, constants, params % lmax, one, &
2057 & workspace % tmp_sph, one, workspace % tmp_grid)
2058 !! Compute gradients of L2L if pl > 0
2059 if (params % pl .gt. 0) then
2060 call tree_grad_l2l(params, constants, workspace % tmp_node_l, &
2061 & workspace % tmp_sph_l_grad, workspace % tmp_sph_l)
2062 end if
2063 !! Diagonal update of computed grid values, that is needed for R^C, a part
2064 !! of R^A and a part of R^B
2065 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis, &
2066 & pt5, constants % vgrid2, constants % vgrid_nbasis, g, &
2067 & constants % nbasis, -one, workspace % tmp_grid, params % ngrid)
2068 !! Scale temporary grid points by corresponding Lebedev weights and ygrid
2069 do igrid = 1, params % ngrid
2070 do isph = 1, params % nsph
2071 workspace % tmp_grid(igrid, isph) = &
2072 & workspace % tmp_grid(igrid, isph) * constants % wgrid(igrid) * &
2073 & ygrid(igrid, isph)
2074 end do
2075 end do
2076 !! Compute all terms of grad_i(R). The second term of R_i^B is already
2077 !! taken into account and the first term is computed together with R_i^C.
2078 do isph = 1, params % nsph
2079 do igrid = 1, params % ngrid
2080 ! Loop over all neighbouring spheres
2081 do ik = constants % inl(isph), constants % inl(isph+1) - 1
2082 ksph = constants % nl(ik)
2083 ! Only consider external grid points
2084 if(constants % ui(igrid, ksph) .eq. zero) cycle
2085 ! build geometrical quantities
2086 c = params % csph(:, ksph) + &
2087 & params % rsph(ksph)*constants % cgrid(:, igrid)
2088 vki = c - params % csph(:, isph)
2089 !vvki = sqrt(vki(1)*vki(1) + vki(2)*vki(2) + &
2090 ! & vki(3)*vki(3))
2091 vvki = dnrm2(3, vki, 1)
2092 tki = vvki / params % rsph(isph)
2093 ! Only consider such points where grad U is non-zero
2094 if((tki.le.tlow) .or. (tki.ge.thigh)) cycle
2095 ! This is entire R^C and the first R^B component (grad_i of U
2096 ! of a sum of R_kj for index inequality j!=k)
2097 ! Indexes k and j are flipped compared to the paper
2098 gg = workspace % tmp_grid(igrid, ksph)
2099 ! Compute grad_i component of forces using precomputed
2100 ! potential gg
2101 !fx(:, isph) = fx(:, isph) - &
2102 ! & dfsw(tki, params % se, params % eta)/ &
2103 ! & params % rsph(isph)*constants % wgrid(igrid)*gg* &
2104 ! & ygrid(igrid, ksph)*(vki/vvki)
2105 fx(:, isph) = fx(:, isph) - &
2106 & dfsw(tki, params % se, params % eta)/ &
2107 & params % rsph(isph)*gg*(vki/vvki)
2108 end do
2109 ! contribution from the sphere itself
2110 if((constants % ui(igrid,isph).gt.zero) .and. &
2111 & (constants % ui(igrid,isph).lt.one)) then
2112 ! R^A component (grad_i of U of a sum of R_ij for index
2113 ! inequality j!=i)
2114 ! Indexes k and j are flipped compared to the paper
2115 gg = workspace % tmp_grid(igrid, isph)
2116 ! Compute grad_i component of forces using precomputed
2117 ! potential gg
2118 !fx(:, isph) = fx(:, isph) + constants % wgrid(igrid)*gg* &
2119 ! & ygrid(igrid, isph)*constants % zi(:, igrid, isph)
2120 fx(:, isph) = fx(:, isph) + gg*constants % zi(:, igrid, isph)
2121 end if
2122 if (constants % ui(igrid, isph) .gt. zero) then
2123 ! Another R^A component (grad_i of potential of a sum of R_ij
2124 ! for index inequality j!=i)
2125 ! Indexes k and j are flipped compared to the paper
2126 ! In case pl=0 MKL does not make gg3 zero reusing old value of
2127 ! gg3, so we have to clear it manually
2128 gg3 = zero
2129 call dgemv('T', params % pl**2, 3, one, &
2130 & workspace % tmp_sph_l_grad(1, 1, isph), &
2131 & (params % pl+1)**2, constants % vgrid2(1, igrid), 1, &
2132 & zero, gg3, 1)
2133 ! Gradient of the near-field potential is a gradient of
2134 ! multipole expansion
2135 inode = constants % snode(isph)
2136 do inear = constants % snear(inode), constants % snear(inode+1)-1
2137 jnode = constants % near(inear)
2138 do jsph_node = constants % cluster(1, jnode), &
2139 & constants % cluster(2, jnode)
2140 jsph = constants % order(jsph_node)
2141 if (isph .eq. jsph) cycle
2142 c = params % csph(:, isph) + &
2143 & params % rsph(isph)*constants % cgrid(:, igrid)
2144 tmp_c = c - params % csph(:, jsph)
2145 call fmm_m2p_work(tmp_c, &
2146 & params % rsph(jsph), params % lmax+1, &
2147 & constants % vscales_rel, one, &
2148 & workspace % tmp_sph_grad(:, 1, jsph), zero, &
2149 & tmp_gg, work)
2150 gg3(1) = gg3(1) + tmp_gg
2151 call fmm_m2p_work(tmp_c, &
2152 & params % rsph(jsph), params % lmax+1, &
2153 & constants % vscales_rel, one, &
2154 & workspace % tmp_sph_grad(:, 2, jsph), zero, &
2155 & tmp_gg, work)
2156 gg3(2) = gg3(2) + tmp_gg
2157 call fmm_m2p_work(tmp_c, &
2158 & params % rsph(jsph), params % lmax+1, &
2159 & constants % vscales_rel, one, &
2160 & workspace % tmp_sph_grad(:, 3, jsph), zero, &
2161 & tmp_gg, work)
2162 gg3(3) = gg3(3) + tmp_gg
2163 end do
2164 end do
2165 ! Accumulate all computed forces
2166 fx(:, isph) = fx(:, isph) - constants % wgrid(igrid)*gg3* &
2167 & ygrid(igrid, isph)*constants % ui(igrid, isph)
2168 end if
2169 end do
2170 end do
2171end subroutine gradr_fmm
2172
2174subroutine gradr(params, constants, workspace, g, ygrid, fx)
2175 implicit none
2176 ! Inputs
2177 type(ddx_params_type), intent(in) :: params
2178 type(ddx_constants_type), intent(in) :: constants
2179 real(dp), intent(in) :: g(constants % nbasis, params % nsph), &
2180 & ygrid(params % ngrid, params % nsph)
2181 ! Temporaries
2182 type(ddx_workspace_type), intent(inout) :: workspace
2183 ! Output
2184 real(dp), intent(out) :: fx(3, params % nsph)
2185 ! Check which gradr to execute
2186 if (params % fmm .eq. 1) then
2187 call gradr_fmm(params, constants, workspace, g, ygrid, fx)
2188 else
2189 call gradr_dense(params, constants, workspace, g, ygrid, fx)
2190 end if
2191end subroutine gradr
2192
2194subroutine gradr_dense(params, constants, workspace, g, ygrid, fx)
2195 implicit none
2196 ! Inputs
2197 type(ddx_params_type), intent(in) :: params
2198 type(ddx_constants_type), intent(in) :: constants
2199 real(dp), intent(in) :: g(constants % nbasis, params % nsph), &
2200 & ygrid(params % ngrid, params % nsph)
2201 ! Temporaries
2202 type(ddx_workspace_type), intent(inout) :: workspace
2203 ! Output
2204 real(dp), intent(out) :: fx(3, params % nsph)
2205 ! Local variables
2206 integer :: isph
2207 ! Simply cycle over all spheres
2208 do isph = 1, params % nsph
2209 call gradr_sph(params, constants, isph, workspace % tmp_vplm, &
2210 & workspace % tmp_vcos, workspace % tmp_vsin, &
2211 & workspace % tmp_vylm, workspace % tmp_vdylm, &
2212 & g, ygrid, fx(:, isph))
2213 end do
2214end subroutine gradr_dense
2215
2219subroutine zeta_grad(params, constants, state, e_cav, forces)
2220 implicit none
2221 type(ddx_params_type), intent(in) :: params
2222 type(ddx_constants_type), intent(in) :: constants
2223 type(ddx_state_type), intent(inout) :: state
2224 real(dp), intent(inout) :: forces(3, params % nsph)
2225 real(dp), intent(in) :: e_cav(3, constants % ncav)
2226 ! local variables
2227 integer :: icav, isph, igrid
2228
2229 icav = 0
2230 do isph = 1, params % nsph
2231 do igrid = 1, params % ngrid
2232 if (constants % ui(igrid, isph) .eq. zero) cycle
2233 icav = icav + 1
2234 forces(:, isph) = forces(:, isph) + pt5 &
2235 & *state % zeta(icav)*e_cav(:, icav)
2236 end do
2237 end do
2238end subroutine zeta_grad
2239
2240end module ddx_gradients
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_l2p_bessel
subroutine tree_l2p_bessel(params, constants, alpha, node_l, beta, grid_v)
TODO.
Definition: ddx_core.f90:2291
ddx_core::tree_l2p_bessel_adj
subroutine tree_l2p_bessel_adj(params, constants, alpha, grid_v, beta, node_l)
TODO.
Definition: ddx_core.f90:2364
ddx_core::tree_l2l_rotation
subroutine tree_l2l_rotation(params, constants, node_l)
Transfer local coefficients over a tree.
Definition: ddx_core.f90:1866
ddx_core::tree_m2m_rotation
subroutine tree_m2m_rotation(params, constants, node_m)
Transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1656
ddx_core::tree_l2p
subroutine tree_l2p(params, constants, alpha, node_l, beta, grid_v, sph_l)
TODO.
Definition: ddx_core.f90:2252
ddx_core::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_bessel_rotation
subroutine tree_m2l_bessel_rotation(params, constants, node_m, node_l)
Transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2095
ddx_core::tree_m2m_bessel_rotation_adj
subroutine tree_m2m_bessel_rotation_adj(params, constants, node_m)
Adjoint transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1822
ddx_core::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_m2p_bessel
subroutine tree_m2p_bessel(params, constants, p, alpha, sph_p, sph_m, beta, grid_v)
TODO.
Definition: ddx_core.f90:2451
ddx_core::tree_l2l_bessel_rotation_adj
subroutine tree_l2l_bessel_rotation_adj(params, constants, node_l)
Adjoint transfer local coefficients over a tree.
Definition: ddx_core.f90:2016
ddx_core::tree_l2l_bessel_rotation
subroutine tree_l2l_bessel_rotation(params, constants, node_l)
Transfer local coefficients over a tree.
Definition: ddx_core.f90:1913
ddx_core::tree_m2p_bessel_nodiag_adj
subroutine tree_m2p_bessel_nodiag_adj(params, constants, p, alpha, grid_v, beta, sph_p, sph_m)
TODO.
Definition: ddx_core.f90:2604
ddx_core::tree_m2l_bessel_rotation_adj
subroutine tree_m2l_bessel_rotation_adj(params, constants, node_l, node_m)
Adjoint transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2146
ddx_core::tree_grad_m2m
subroutine tree_grad_m2m(params, constants, sph_m, sph_m_grad, work)
TODO.
Definition: ddx_core.f90:2650
ddx_core::tree_m2p_adj
subroutine tree_m2p_adj(params, constants, p, alpha, grid_v, beta, sph_m)
TODO.
Definition: ddx_core.f90:2503
ddx_core::tree_m2p_bessel_adj
subroutine tree_m2p_bessel_adj(params, constants, p, alpha, grid_v, beta, sph_p, sph_m)
TODO.
Definition: ddx_core.f90:2555
ddx_core::tree_m2l_rotation_adj
subroutine tree_m2l_rotation_adj(params, constants, node_l, node_m)
Adjoint transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2208
ddx_core::tree_l2l_rotation_adj
subroutine tree_l2l_rotation_adj(params, constants, node_l)
Adjoint transfer local coefficients over a tree.
Definition: ddx_core.f90:1959
ddx_core::tree_m2m_rotation_adj
subroutine tree_m2m_rotation_adj(params, constants, node_m)
Adjoint transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1778
ddx_core::tree_m2m_bessel_rotation
subroutine tree_m2m_bessel_rotation(params, constants, node_m)
Transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1717
ddx_core::dbasis
subroutine dbasis(params, constants, x, basloc, dbsloc, vplm, vcos, vsin)
Compute first derivatives of spherical harmonics.
Definition: ddx_core.f90:1161
ddx_core::intmlp
real(dp) function intmlp(params, constants, t, sigma, basloc)
TODO.
Definition: ddx_core.f90:1292
ddx_gradients
Core routines and parameters specific to gradients.
Definition: ddx_gradients.f90:11
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