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