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