ddx 1.0.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
17
19subroutine contract_grad_l(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr)
20type(ddx_params_type), intent(in) :: params
21 type(ddx_constants_type), intent(in) :: constants
22 integer, intent(in) :: isph
23 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sigma
24 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi
25 real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm
26 real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc
27 real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin
28 real(dp), dimension(3), intent(inout) :: fx
29 real(dp), optional, intent(inout) :: dr
30
31 logical :: do_dr
32 real(dp) :: dr_local
33
34 if (present(dr)) then
35 do_dr = .true.
36 else
37 do_dr = .false.
38 end if
39
40 dr_local = zero
41
42 if (params%switching.eq.0) then
43 call contract_gradi_lik(params, constants, isph, sigma, xi(:, isph), basloc, dbsloc, vplm, vcos, vsin, fx, dr_local)
44 call contract_gradi_lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr_local)
45 else
46 call contract_grad_l_new(params, constants, isph, sigma, xi, &
47 & basloc, dbsloc, vplm, vcos, vsin, fx, dr_local)
48 end if
49
50 if (do_dr) dr = dr_local
51
52end subroutine contract_grad_l
53
54subroutine contract_grad_l_new(params, constants, isph, x, xi, &
55 & basloc, dbsloc, vplm, vcos, vsin, fx, dr)
56 implicit none
57 type(ddx_params_type), intent(in) :: params
58 type(ddx_constants_type), intent(in) :: constants
59 integer, intent(in) :: isph
60 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: x
61 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi
62 real(dp), dimension(3), intent(inout) :: fx
63 real(dp), intent(inout) :: dr
64 ! passed scratch TODO: find a better solution
65 real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm
66 real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc
67 real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin
68 ! local variables
69 real(dp) :: tlow, thigh
70 integer :: ig, ij, ji, jsph, icomp, jcomp, l, m, ind, ki, ksph, kj
71 real(dp) :: vij(3), vvij, sij(3), tij, oij, qij, grad_tij(3), &
72 & dr_sij(3), dr_tij, f_i, d_i, xij
73 real(dp) :: vji(3), vvji, sji(3), tji, oji, qji, grad_tji(3), &
74 & dr_sji(3), dr_tji, f_j, d_j, xji
75 real(dp) :: fl, a, b, c, d1, d2, e, tmp_fx(3), sjac_dylm(3), &
76 & sjac(3,3), tmp_dr, dr_s_dylm, grad_p, bj
77 real(dp) :: vik(3), vvik, sik(3), tik, grad_tik(3), dr_tik, xik
78 real(dp) :: vjk(3), vvjk, sjk(3), tjk, xjk
79 real(dp) :: rho, ctheta, stheta, cphi, sphi
80
81 tlow = one - pt5*(one - params % se)*params % eta
82 thigh = one + pt5*(one + params % se)*params % eta
83
84 ! ij cases
85 do ig = 1, params % ngrid
86 do ij = constants % inl(isph), constants % inl(isph+1) - 1
87 jsph = constants % nl(ij)
88
89 ! compute the required geometrical quantities
90 oij = compute_omega(params, constants, isph, jsph, ig)
91 if (oij.eq.zero) cycle
92 vij = params % csph(:,isph) + &
93 & params % rsph(isph)*constants % cgrid(:,ig) - &
94 & params % csph(:,jsph)
95 vvij = sqrt(vij(1)*vij(1) + vij(2)*vij(2) + vij(3)*vij(3))
96 tij = vvij/params % rsph(jsph)
97 sij = vij/vvij
98
99 ! compute the differentiated geometrical quantities
100 ! wrt the nuclear position
101 grad_tij(:) = sij(:)/params%rsph(jsph)
102 sjac = zero
103 sjac(1,1) = one
104 sjac(2,2) = one
105 sjac(3,3) = one
106 qij = one/vvij
107 do icomp = 1, 3
108 do jcomp = 1, 3
109 sjac(icomp,jcomp) = qij*(sjac(icomp,jcomp) &
110 & - sij(icomp)*sij(jcomp))
111 end do
112 end do
113
114 ! compute the differentiated geometrical quantities
115 ! wrt the radius
116 dr_sij = constants%cgrid(:,ig)/vvij &
117 & - vij*dot_product(vij, constants%cgrid(:,ig))/(vvij**3)
118 dr_tij = dot_product(vij, constants%cgrid(:,ig)) &
119 & /(vvij*params%rsph(jsph))
120
121 ! compute the basis of spherical harmonics and its gradient
122 call dbasis(params, constants, sij, basloc, dbsloc, vplm, &
123 & vcos, vsin)
124
125 ! accumulate the derivative
126 e = zero
127 a = xi(ig,isph)*constants%wgrid(ig)*oij
128 tmp_fx(:) = zero
129 tmp_dr = zero
130 do l = 0, params % lmax
131 ind = l*l + l + 1
132 fl = dble(l)
133 b = a*fourpi/(two*fl + one)*tij**(l-1)
134 do m = -l, l
135
136 ! product of the jacobian of sij with the spherical
137 ! harmonics gradient
138 sjac_dylm(1) = sjac(1,1)*dbsloc(1,ind+m) + &
139 & sjac(1,2)*dbsloc(2,ind+m)+sjac(1,3) &
140 & *dbsloc(3,ind+m)
141 sjac_dylm(2) = sjac(2,1)*dbsloc(1,ind+m) + &
142 & sjac(2,2)*dbsloc(2,ind+m)+sjac(2,3)* &
143 & dbsloc(3,ind+m)
144 sjac_dylm(3) = sjac(3,1)*dbsloc(1,ind+m) + &
145 & sjac(3,2)*dbsloc(2,ind+m)+sjac(3,3) &
146 & *dbsloc(3,ind+m)
147
148 ! product of the radial derivative of sij
149 ! with the spherical harmonics gradient
150 dr_s_dylm = dr_sij(1)*dbsloc(1,ind+m) &
151 & + dr_sij(2)*dbsloc(2,ind+m) &
152 & + dr_sij(3)*dbsloc(3,ind+m)
153
154 c = b*x(ind+m,jsph)
155 d1 = c*fl*basloc(ind+m)
156 d2 = c*tij
157 tmp_fx(:) = tmp_fx(:) &
158 & - c*fl*basloc(ind+m)*grad_tij(:) &
159 & - c*tij*sjac_dylm(:)
160 dr = dr - c*fl*basloc(ind+m)*dr_tij &
161 & - c*tij*dr_s_dylm
162 e = e + b*x(ind+m,jsph)*basloc(ind+m)*tij
163 end do
164 end do
165 e = e/oij
166
167 ! add the switching derivative for the ij case
168 d_i = constants%switching%d_ni(ig,isph)
169 f_i = constants%switching%f_ni(ig,isph)
170 a = e/(d_i + f_i)
171 if (tij.ge.tlow .and. tij.le.thigh) then
172 grad_p = dfsw(tij, params % se, params % eta)
173 b = a*grad_p
174 tmp_fx(:) = tmp_fx(:) - b*grad_tij(:)
175 tmp_dr = tmp_dr - b*dr_tij
176 end if
177 xij = fsw(tij, params % se, params % eta)
178 if (xij .ne. zero) then
179 b = a*xij/(d_i + f_i)
180 do ki = constants % inl(isph), constants % inl(isph+1) - 1
181 ksph = constants % nl(ki)
182 vik = params%csph(:,isph) + &
183 & params%rsph(isph)*constants%cgrid(:,ig) - &
184 & params%csph(:,ksph)
185 vvik = sqrt(vik(1)*vik(1) + vik(2)*vik(2) &
186 & + vik(3)*vik(3))
187 tik = vvik/params % rsph(ksph)
188 if (tik.ge.tlow .and. tik.le.thigh) then
189 xik = fsw(tik, params % se, params % eta)
190 sik = vik/vvik
191 grad_tik(:) = sik(:)/params%rsph(ksph)
192 dr_tik = dot_product(vik, constants%cgrid(:,ig)) &
193 & /(vvik*params%rsph(ksph))
194 grad_p = dfsw(tik, params % se, params % eta)
195
196 if (grad_p.ne.zero) then
197 c = b*(one - d_i/(one - xik))*grad_p
198 else
199 c = zero
200 end if
201
202 tmp_fx(:) = tmp_fx(:) + c*grad_tik
203 tmp_dr = tmp_dr + c*dr_tik
204 end if
205 end do
206 end if
207
208 dr = dr + tmp_dr
209 fx(:) = fx(:) + tmp_fx(:)
210 end do
211 end do
212
213 ! ji cases
214 do ig = 1, params % ngrid
215 do ji = constants % inl(isph), constants % inl(isph+1) - 1
216 jsph = constants % nl(ji)
217
218 ! compute the geometrical quantites required
219 oji = compute_omega(params, constants, jsph, isph, ig)
220 if (oji.eq.zero) cycle
221 vji = params % csph(:,jsph) + &
222 & params % rsph(jsph)*constants % cgrid(:,ig) - &
223 & params % csph(:,isph)
224 vvji = sqrt(vji(1)*vji(1) + vji(2)*vji(2) + vji(3)*vji(3))
225 tji = vvji/params % rsph(isph)
226 sji = vji/vvji
227
228 ! computed the differentiated geometrical quantities
229 ! wrt the nuclear position
230 grad_tji(:) = -sji/params % rsph(isph)
231 sjac = zero
232 sjac(1,1) = - one
233 sjac(2,2) = - one
234 sjac(3,3) = - one
235 qji = one/vvji
236 do icomp = 1, 3
237 do jcomp = 1, 3
238 sjac(icomp,jcomp) = qji*(sjac(icomp,jcomp) &
239 & + sji(icomp)*sji(jcomp))
240 end do
241 end do
242
243 ! computed the differentiated geometrical quantities
244 ! wrt the radius
245 dr_sji = zero
246 dr_tji = -tji/params%rsph(isph)
247
248 ! compute the basis of spherical harmonics and its gradient
249 call dbasis(params, constants, sji, basloc, dbsloc, vplm, &
250 & vcos, vsin)
251
252 ! accumulate the derivative
253 e = zero
254 a = xi(ig,jsph)*constants%wgrid(ig)*oji
255 tmp_fx(:) = zero
256 tmp_dr = zero
257 do l = 0, params % lmax
258 ind = l*l + l + 1
259 fl = dble(l)
260 b = a*fourpi/(two*fl + one)*tji**(l-1)
261 do m = -l, l
262
263 ! product of the jacobian of sji with the spherical
264 ! harmonics gradient
265 sjac_dylm(1) = sjac(1,1)*dbsloc(1,ind+m) + &
266 & sjac(1,2)*dbsloc(2,ind+m)+sjac(1,3)*dbsloc(3,ind+m)
267 sjac_dylm(2) = sjac(2,1)*dbsloc(1,ind+m) + &
268 & sjac(2,2)*dbsloc(2,ind+m)+sjac(2,3)*dbsloc(3,ind+m)
269 sjac_dylm(3) = sjac(3,1)*dbsloc(1,ind+m) + &
270 & sjac(3,2)*dbsloc(2,ind+m)+sjac(3,3)*dbsloc(3,ind+m)
271
272 ! product of the radial derivative of sji
273 ! with the spherical harmonics gradient
274 dr_s_dylm = dr_sji(1)*dbsloc(1,ind+m) &
275 & + dr_sji(2)*dbsloc(2,ind+m) &
276 & + dr_sji(3)*dbsloc(3,ind+m)
277
278 c = b*x(ind+m,isph)
279 d1 = c*fl*basloc(ind+m)
280 d2 = c*tji
281 tmp_fx(:) = tmp_fx(:) &
282 & - c*fl*basloc(ind+m)*grad_tji(:) &
283 & - c*tji*sjac_dylm(:)
284 tmp_dr = tmp_dr - c*fl*basloc(ind+m)*dr_tji &
285 & - c*tji*dr_s_dylm
286 e = e + b*x(ind+m,isph)*basloc(ind+m)*tji
287 end do
288 end do
289 e = e/oji
290
291 ! add the switching derivative for the ji case
292 d_j = constants%switching%d_ni(ig,jsph)
293 f_j = constants%switching%f_ni(ig,jsph)
294 a = e/(d_j + f_j)
295 if (tji.ge.tlow .and. tji.le.thigh) then
296 grad_p = dfsw(tji, params % se, params % eta)
297 b = a*grad_p
298 tmp_fx(:) = tmp_fx(:) - b*grad_tji(:)
299 tmp_dr = tmp_dr - b*dr_tji
300 end if
301
302 fx(:) = fx(:) + tmp_fx(:)
303 dr = dr + tmp_dr
304 end do
305 end do
306
307 ! additional switching contributions, jk with j and k neighbors of i
308 tmp_fx = zero
309 tmp_dr = zero
310 do ig = 1, params % ngrid
311 do ij = constants % inl(isph), constants % inl(isph+1) - 1
312 jsph = constants % nl(ij)
313
314 ! compute ij geometrical quantities
315 vji = params % csph(:,jsph) + &
316 & params % rsph(jsph)*constants % cgrid(:,ig) - &
317 & params % csph(:,isph)
318 vvji = sqrt(vji(1)*vji(1) + vji(2)*vji(2) + vji(3)*vji(3))
319 tji = vvji/params % rsph(isph)
320 sji = vji/vvji
321
322 if (.not.(tji.ge.tlow .and. tji.le.thigh)) cycle
323
324 ! computed the differentiated geometrical quantities
325 ! wrt the nuclear position
326 grad_tji(:) = -sji/params % rsph(isph)
327
328 ! computed the differentiated geometrical quantities
329 ! wrt the radius
330 dr_tji = -tji/params%rsph(isph)
331
332 xji = fsw(tji, params % se, params % eta)
333
334 d_j = constants%switching%d_ni(ig,jsph)
335 f_j = constants%switching%f_ni(ig,jsph)
336 grad_p = dfsw(tji, params % se, params % eta)
337
338 if (grad_p.eq.zero) cycle
339
340 bj = zero
341 do kj = constants % inl(jsph), constants % inl(jsph+1) - 1
342 ksph = constants % nl(kj)
343 if (jsph.eq.ksph) cycle
344
345 ! compute the remaining geometrical quantities
346 vjk = params % csph(:,jsph) + &
347 & params % rsph(jsph)*constants % cgrid(:,ig) - &
348 & params % csph(:,ksph)
349 vvjk = sqrt(vjk(1)*vjk(1) + vjk(2)*vjk(2) &
350 & + vjk(3)*vjk(3))
351 tjk = vvjk/params % rsph(ksph)
352 sjk = vjk/vvjk
353
354 xjk = fsw(tjk, params % se, params % eta)
355 if (xjk.eq.zero) cycle
356
357 ! compute the basis of spherical harmonics
358 call ylmbas(sjk, rho, ctheta, stheta, cphi, sphi, &
359 & params % lmax, constants % vscales, basloc, &
360 & vplm, vcos, vsin)
361
362 ! accumulate the derivative
363 e = zero
364 a = xi(ig,jsph)*constants%wgrid(ig)
365 do l = 0, params % lmax
366 ind = l*l + l + 1
367 fl = dble(l)
368 b = a*fourpi/(two*fl + one)*tjk**l
369 do m = -l, l
370 e = e + b*x(ind+m,ksph)*basloc(ind+m)
371 end do
372 end do
373
374 bj = bj + e*xjk
375
376 end do
377
378 c = bj/(d_j + f_j)**2*(one - d_j/(one - xji))*grad_p
379 tmp_fx(:) = tmp_fx(:) + c*grad_tji
380 tmp_dr = tmp_dr + c*dr_tji
381 end do
382 end do
383
384 fx(:) = fx(:) + tmp_fx(:)
385 dr = dr + tmp_dr
386
387end subroutine contract_grad_l_new
388
390subroutine contract_gradi_lik(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr)
391 type(ddx_params_type), intent(in) :: params
392 type(ddx_constants_type), intent(in) :: constants
393 integer, intent(in) :: isph
394 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sigma
395 real(dp), dimension(params % ngrid), intent(in) :: xi
396 real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm
397 real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc
398 real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin
399 real(dp), dimension(3), intent(inout) :: fx
400 integer :: ig, ij, jsph, l, ind, m
401 real(dp) :: vvij, tij, xij, oij, t, fac, fl, f1, f2, f3, beta, tlow, thigh
402 real(dp) :: vij(3), sij(3), alp(3), va(3)
403 real(dp) :: dsij, dtij(3), alp1(3), alp2(3), alp_rad, dsij_rad(3), dtij_rad, alp1_rad, alp2_rad, va_rad
404 real(dp), external :: dnrm2
405 real(dp), intent(inout) :: dr
406 real(dp) :: sjac(3,3), qij
407 tlow = one - pt5*(one - params % se)*params % eta
408 thigh = one + pt5*(one + params % se)*params % eta
409
410 do ig = 1, params % ngrid
411 va = zero
412 va_rad = zero
413 do ij = constants % inl(isph), constants % inl(isph+1) - 1
414 jsph = constants % nl(ij)
415 vij = params % csph(:,isph) + &
416 & params % rsph(isph)*constants % cgrid(:,ig) - &
417 & params % csph(:,jsph)
418 vvij = dnrm2(3, vij, 1)
419 tij = vvij/params % rsph(jsph)
420 if (tij.ge.thigh) cycle
421 if (tij.ne.zero) then
422 sij = vij/vvij
423 else
424 sij = one
425 end if
426
427 dsij = 1.0_dp / (tij*params%rsph(jsph)) !scalar
428 dtij(:) = sij(:)/params%rsph(jsph) !vector
429
430 ! build the jacobian of sik
431 sjac = zero
432 sjac(1,1) = one
433 sjac(2,2) = one
434 sjac(3,3) = one
435 qij = one/vvij
436 do icomp = 1, 3
437 do jcomp = 1, 3
438 sjac(icomp,jcomp) = qij*(sjac(icomp,jcomp) &
439 & - sij(icomp)*sij(jcomp))
440 end do
441 end do
442
443 dsij_rad = constants%cgrid(:,ig)/vvij &
444 & - vij * dot_product(vij, constants%cgrid(:,ig)) / (vvij**3)
445 dtij_rad = dot_product(vij, constants%cgrid(:,ig)) &
446 & / (vvij*params%rsph(jsph))
447
448 call dbasis(params, constants, sij, basloc, dbsloc, vplm, vcos, vsin)
449 alp = zero
450 alp1 = zero
451 alp2 = zero
452
453 alp_rad = zero
454 alp1_rad = zero
455 alp2_rad = zero
456
457 t = one
458 do l = 1, params % lmax
459 ind = l*l + l + 1
460 fl = dble(l)
461 fac = t/(constants % vscales(ind)**2)
462 do m = -l, l
463 f2 = fac*sigma(ind+m,jsph)
464 f1 = f2*fl*basloc(ind+m)
465
466 alp1(:) = f1*dtij
467 alp1_rad = f1*dtij_rad
468
469 alp2(1) = sjac(1,1)*dbsloc(1,ind+m) + &
470 & sjac(1,2)*dbsloc(2,ind+m)+sjac(1,3)*dbsloc(3,ind+m)
471 alp2(2) = sjac(2,1)*dbsloc(1,ind+m) + &
472 & sjac(2,2)*dbsloc(2,ind+m)+sjac(2,3)*dbsloc(3,ind+m)
473 alp2(3) = sjac(3,1)*dbsloc(1,ind+m) + &
474 & sjac(3,2)*dbsloc(2,ind+m)+sjac(3,3)*dbsloc(3,ind+m)
475 alp2(:) = alp2(:)*tij*f2
476
477 alp2_rad = f2*tij*dot_product(dsij_rad, dbsloc(:,ind+m))
478
479 alp(:) = alp(:) + alp1(:) + alp2(:)
480 alp_rad = alp_rad + alp1_rad + alp2_rad
481
482 end do
483 t = t*tij
484 end do
485 beta = intmlp(params, constants, tij,sigma(:,jsph),basloc)
486 xij = fsw(tij, params % se, params % eta)
487 if (constants % fi(ig,isph).gt.one) then
488 oij = xij/constants % fi(ig,isph)
489 f2 = -oij/constants % fi(ig,isph)
490 else
491 oij = xij
492 f2 = zero
493 end if
494 f1 = oij
495 va(:) = va(:) + f1*alp(:) + beta*f2*constants % zi(:,ig,isph)
496 va_rad = va_rad + f1*alp_rad + beta*f2*constants % zi_dr(ig,isph)
497 if (tij .gt. tlow) then
498 f3 = beta*dfsw(tij,params % se,params % eta)
499 if (constants % fi(ig,isph).gt.one) f3 = f3/constants % fi(ig,isph)
500 va(:) = va(:) + f3*dtij
501 va_rad = va_rad + f3*dtij_rad
502 end if
503 end do
504 fx = fx - constants % wgrid(ig)*xi(ig)*va(:)
505 dr = dr - constants % wgrid(ig)*xi(ig)*va_rad
506 end do
507end subroutine contract_gradi_lik
508
510subroutine contract_gradi_lji(params, constants, isph, sigma, xi, basloc, dbsloc, vplm, vcos, vsin, fx, dr)
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) :: sigma
515 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi
516 real(dp), dimension(constants % nbasis), intent(inout) :: basloc, vplm
517 real(dp), dimension(3, constants % nbasis), intent(inout) :: dbsloc
518 real(dp), dimension(params % lmax+1), intent(inout) :: vcos, vsin
519 real(dp), dimension(3), intent(inout) :: fx
520
521 integer :: ig, ji, jsph, l, ind, m, jk, ksph
522 logical :: proc
523 real(dp) :: vvji, tji, xji, oji, t, fac, fl, f1, f2, beta, di, tlow, thigh
524 real(dp) :: b, g1, g2, vvjk, tjk, xjk
525 real(dp) :: vji(3), sji(3), alp(3), vb(3), vjk(3), sjk(3), vc(3)
526 real(dp) :: dsji, dtji(3), alp1(3), alp2(3), alp_rad, dsji_rad(3), dtji_rad, alp1_rad, alp2_rad, vb_rad, vc_rad
527 real(dp) :: rho, ctheta, stheta, cphi, sphi
528 real(dp), external :: dnrm2
529 real(dp), intent(inout) :: dr
530 real(dp) :: sjac(3,3), qji
531
532 tlow = one - pt5*(one - params % se)*params % eta
533 thigh = one + pt5*(one + params % se)*params % eta
534
535 do ig = 1, params % ngrid
536 vb = zero
537 vc = zero
538
539 vb_rad = zero
540 vc_rad = zero
541
542 do ji = constants % inl(isph), constants % inl(isph+1) - 1
543 jsph = constants % nl(ji)
544 vji = params % csph(:,jsph) + &
545 & params % rsph(jsph)*constants % cgrid(:,ig) - &
546 & params % csph(:,isph)
547 vvji = dnrm2(3, vji, 1)
548 tji = vvji/params % rsph(isph)
549 if (tji.gt.thigh) cycle
550 if (tji.ne.zero) then
551 sji = vji/vvji
552 else
553 sji = one
554 end if
555
556 ! build the jacobian of sji
557 sjac = zero
558 sjac(1,1) = - one
559 sjac(2,2) = - one
560 sjac(3,3) = - one
561 qji = one/vvji
562 do icomp = 1, 3
563 do jcomp = 1, 3
564 sjac(icomp,jcomp) = qji*(sjac(icomp,jcomp) &
565 & + sji(icomp)*sji(jcomp))
566 end do
567 end do
568
569 dsji = 1.0_dp/(tji*params % rsph(isph))
570 dtji(:) = -sji/params % rsph(isph)
571
572 dsji_rad = zero
573 dtji_rad = -tji/params%rsph(isph)
574
575 call dbasis(params, constants, sji, basloc, dbsloc, vplm, vcos, vsin)
576 alp = zero
577 alp1 = zero
578 alp2 = zero
579
580 alp_rad = zero
581 alp1_rad = zero
582 alp2_rad = zero
583
584 t = one
585 do l = 1, params % lmax
586 ind = l*l + l + 1
587 fl = dble(l)
588 fac = t/(constants % vscales(ind)**2)
589 do m = -l, l
590 f2 = fac*sigma(ind+m,isph)
591 f1 = f2*fl*basloc(ind+m)
592
593 alp1(:) = f1*dtji
594 alp1_rad = f1*dtji_rad
595
596 alp2(1) = sjac(1,1)*dbsloc(1,ind+m) + &
597 & sjac(1,2)*dbsloc(2,ind+m)+sjac(1,3)*dbsloc(3,ind+m)
598 alp2(2) = sjac(2,1)*dbsloc(1,ind+m) + &
599 & sjac(2,2)*dbsloc(2,ind+m)+sjac(2,3)*dbsloc(3,ind+m)
600 alp2(3) = sjac(3,1)*dbsloc(1,ind+m) + &
601 & sjac(3,2)*dbsloc(2,ind+m)+sjac(3,3)*dbsloc(3,ind+m)
602 alp2(:) = alp2(:)*tji*f2
603
604 alp2_rad = f2*tji*dot_product(dsji_rad,dbsloc(:,ind+m))
605
606 alp(:) = alp(:) + alp1(:) + alp2(:)
607 alp_rad = alp_rad + alp1_rad + alp2_rad
608
609 end do
610 t = t*tji
611 end do
612
613 xji = fsw(tji, params % se, params % eta)
614 if (constants % fi(ig,jsph).gt.one) then
615 oji = xji/constants % fi(ig,jsph)
616 else
617 oji = xji
618 end if
619 vb = vb + oji*alp*xi(ig,jsph)
620 vb_rad = vb_rad + oji*alp_rad*xi(ig,jsph)
621 if (tji .gt. tlow) then
622 beta = intmlp(params, constants, tji, sigma(:,isph), basloc)
623 if (constants % fi(ig,jsph) .gt. one) then
624 di = one/constants % fi(ig,jsph)
625 fac = di*xji
626 proc = .false.
627 b = zero
628 do jk = constants % inl(jsph), constants % inl(jsph+1) - 1
629 ksph = constants % nl(jk)
630 vjk = params % csph(:,jsph) + &
631 & params % rsph(jsph)*constants % cgrid(:,ig) - &
632 & params % csph(:,ksph)
633 !vvjk = sqrt(dot_product(vjk,vjk))
634 vvjk = dnrm2(3, vjk, 1)
635 tjk = vvjk/params % rsph(ksph)
636 if (ksph.ne.isph) then
637 if (tjk .le. thigh) then
638 proc = .true.
639 sjk = vjk/vvjk
640 !call ylmbas(sjk,basloc,vplm,vcos,vsin)
641 call ylmbas(sjk, rho, ctheta, stheta, cphi, sphi, &
642 & params % lmax, constants % vscales, basloc, vplm, &
643 & vcos, vsin)
644 g1 = intmlp(params, constants, tjk, sigma(:,ksph), basloc)
645 xjk = fsw(tjk, params % se, params % eta)
646 b = b + g1*xjk
647 end if
648 end if
649 end do
650 if (proc) then
651 g1 = di*di*dfsw(tji, params % se, params % eta)
652 g2 = g1*xi(ig,jsph)*b
653 vc = vc - g2*dtji
654 vc_rad = vc_rad - g2*dtji_rad
655 end if
656 else
657 di = one
658 fac = zero
659 end if
660 f2 = (one-fac)*di*dfsw(tji, params % se, params % eta)
661 vb = vb + f2*xi(ig,jsph)*beta*dtji
662 vb_rad = vb_rad + f2*xi(ig,jsph)*beta*dtji_rad
663 end if
664 end do
665 fx = fx - constants % wgrid(ig)*(vb + vc)
666 dr = dr - constants % wgrid(ig)*(vb_rad + vc_rad)
667 end do
668end subroutine contract_gradi_lji
669
670
671
673subroutine contract_grad_u(params, constants, isph, xi, phi, fx, dr)
674 implicit none
675 type(ddx_params_type), intent(in) :: params
676 type(ddx_constants_type), intent(in) :: constants
677 integer, intent(in) :: isph
678 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: xi, phi
679 real(dp), dimension(3), intent(inout) :: fx
680 real(dp), optional, intent(inout) :: dr
681 integer :: ig, ji, jsph
682 real(dp) :: vvji, tji, fac, swthr
683 real(dp) :: alp(3), vji(3), sji(3), dtji(3)
684 real(dp) :: dtji_rad, alp_rad
685 real(dp), external :: dnrm2
686 real(dp) :: dr_local
687 real(dp) :: xi_w_v_i, d_i, f_i, xi_w_v_j, d_j, f_j, grad_p, &
688 & chi_ij, chi_ji, grad_t(3), a, vij(3), vvij, tij, sij(3), dr_t
689
690 if (params%switching.eq.0) then
691 dr_local = zero
692 do ig = 1, params % ngrid
693 alp = zero
694 alp_rad = zero
695 if (constants % ui(ig,isph) .gt. zero .and. constants % ui(ig,isph).lt.one) then
696 alp = alp + phi(ig,isph)*xi(ig,isph)*constants % zi(:,ig,isph)
697 alp_rad = alp_rad + phi(ig,isph)*xi(ig,isph)*constants % zi_dr(ig,isph)
698 end if
699 do ji = constants % inl(isph), constants % inl(isph+1) - 1
700 jsph = constants % nl(ji)
701 vji = params % csph(:,jsph) + &
702 & params % rsph(jsph)*constants % cgrid(:,ig) - &
703 & params % csph(:,isph)
704 vvji = dnrm2(3, vji, 1)
705 tji = vvji/params % rsph(isph)
706 swthr = one + (params % se + 1.d0)*params % eta / 2.d0
707 if (tji.lt.swthr .and. tji.gt.swthr-params % eta .and. constants % ui(ig,jsph).gt.zero) then
708 sji = vji/vvji
709 dtji = - sji / params % rsph(isph)
710 dtji_rad = - vvji/(params%rsph(isph)**2)
711 fac = dfsw(tji, params % se, params % eta)
712 alp = alp + fac*phi(ig,jsph)*xi(ig,jsph) * dtji
713 alp_rad = alp_rad + fac*phi(ig,jsph)*xi(ig,jsph) * dtji_rad
714 end if
715 end do
716 fx = fx - constants % wgrid(ig)*alp
717 dr_local = dr_local - constants % wgrid(ig)*alp_rad
718
719 end do
720 else
721 dr_local = zero
722 do ig = 1, params % ngrid
723 xi_w_v_i = xi(ig,isph)*constants%wgrid(ig)*phi(ig,isph)
724 d_i = constants%switching%d_ni(ig,isph)
725 f_i = constants%switching%f_ni(ig,isph)
726 do ji = constants % inl(isph), constants % inl(isph+1) - 1
727 jsph = constants % nl(ji)
728 vji = params % csph(:,jsph) + &
729 & params % rsph(jsph)*constants % cgrid(:,ig) - &
730 & params % csph(:,isph)
731 vvji = sqrt(vji(1)*vji(1) + vji(2)*vji(2) &
732 & + vji(3)*vji(3))
733 sji = vji/vvji
734 tji = vvji/params % rsph(isph)
735
736 vij = params % csph(:,isph) + &
737 & params % rsph(isph)*constants % cgrid(:,ig) - &
738 & params % csph(:,jsph)
739 vvij = sqrt(vij(1)*vij(1) + vij(2)*vij(2) &
740 & + vij(3)*vij(3))
741 sij = vij/vvij
742 tij = vvij/params % rsph(jsph)
743
744 swthr = one + (params % se + one)*params % eta/two
745
746 xi_w_v_j = xi(ig,jsph)*constants%wgrid(ig)*phi(ig,jsph)
747 d_j = constants%switching%d_ni(ig,jsph)
748 f_j = constants%switching%f_ni(ig,jsph)
749
750 if (tij.lt.swthr .and. tij.gt.swthr-params % eta) then
751 ! diagonal terms (i with its neighbors j)
752 ! \partial \chi_{ij} / \partial x_i
753 grad_p = dfsw(tij, params % se, params % eta)
754 chi_ij = fsw(tij, params % se, params % eta)
755 grad_t = sij/params%rsph(jsph)
756 dr_t = (constants%cgrid(1,ig)*sij(1) &
757 & + constants%cgrid(2,ig)*sij(2) &
758 & + constants%cgrid(3,ig)*sij(3)) &
759 & /params%rsph(jsph)
760 a = xi_w_v_i/(f_i+d_i)**2 &
761 & *(-d_i - f_i*d_i/(one-chi_ij))*grad_p
762 fx = fx + a*grad_t
763 dr_local = dr_local + a*dr_t
764 end if
765
766 if (tji.lt.swthr .and. tji.gt.swthr-params % eta) then
767 ! off-diagonal terms (neighbors of i)
768 ! \partial \chi_{ij} / \partial x_j
769 grad_p = dfsw(tji, params % se, params % eta)
770 chi_ji = fsw(tji, params % se, params % eta)
771 grad_t = -sji/params%rsph(isph)
772 dr_t = -tji/params%rsph(isph)
773 a = xi_w_v_j/(f_j+d_j)**2 &
774 & *(-d_j - f_j*d_j/(one-chi_ji))*grad_p
775 fx = fx + a*grad_t
776 dr_local = dr_local + a*dr_t
777 end if
778 end do
779 end do
780 end if
781
782 if (present(dr)) dr = dr_local
783end subroutine contract_grad_u
784
793
794subroutine contract_grad_b(params, constants, isph, Xe, Xadj_e, force)
795 !! input/output
796 type(ddx_params_type), intent(in) :: params
797 type(ddx_constants_type), intent(in) :: constants
798 integer, intent(in) :: isph
799 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xe
800 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: Xadj_e
801 real(dp), dimension(3), intent(inout) :: force
802 real(dp) :: dr
803
804 if (params%switching.eq.0) then
805 call contract_gradi_bik(params, constants, isph, xe(:,:), &
806 & xadj_e(:, isph), force)
807 call contract_gradi_bji(params, constants, isph, xe(:,:), xadj_e, force)
808 else
809 call contract_grad_b_new(params, constants, isph, xe(:,:), xadj_e, force, dr)
810 end if
811
812end subroutine contract_grad_b
813
823subroutine contract_grad_b_new(params, constants, isph, Xe, Xadj_e, force, dr)
824 implicit none
825 type(ddx_params_type), intent(in) :: params
826 type(ddx_constants_type), intent(in) :: constants
827 integer, intent(in) :: isph
828 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xe
829 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: Xadj_e
830 real(dp), dimension(3), intent(inout) :: force
831 real(dp), intent(inout) :: dr
832 ! local variables
833 real(dp) :: tlow, thigh
834 integer :: ig, ij, ji, jsph, l, ind, ki, ksph, kj
835 real(dp) :: vij(3), vvij, sij(3), tij, oij, grad_tij(3), d_i, f_i, xij
836 real(dp) :: vji(3), vvji, sji(3), tji, oji, grad_tji(3), dr_tji, &
837 & d_j, f_j, xji
838 real(dp) :: vik(3), vvik, sik(3), tik, grad_tik(3), xik
839 real(dp) :: vjk(3), vvjk, tjk, xjk
840 real(dp) :: vtij(3), vtji(3), vtjk(3), alpha(3)
841 real(dp) :: beta, beta_dr, beta_jk, fac
842 real(dp) :: Xe_dr(constants % nbasis)
843 real(dp) :: a, b, c, e, bj, grad_p, tmp_fx(3), tmp_dr
844 real(dp), external :: dnrm2
845 real(dp) :: work(params % lmax+1)
846 complex(dp) :: work_complex(params % lmax+1)
847
848 tlow = one - pt5*(one - params % se)*params % eta
849 thigh = one + pt5*(one + params % se)*params % eta
850
851 do l = 0, params % lmax
852 ind = l*l + l + 1
853 fac = - params % kappa * constants % DI_ri(l, isph) &
854 & / constants % SI_ri(l, isph)
855 xe_dr(ind-l:ind+l) = fac * xe(ind-l:ind+l, isph)
856 end do
857
858 ! ij case
859 do ig = 1, params % ngrid
860 do ij = constants % inl(isph), constants % inl(isph+1) - 1
861 jsph = constants % nl(ij)
862
863 ! compute the required geometrical quantities
864 oij = compute_omega(params, constants, isph, jsph, ig)
865 if (oij.eq.zero) cycle
866 vij = params % csph(:,isph) + &
867 & params % rsph(isph)*constants % cgrid(:,ig) - &
868 & params % csph(:,jsph)
869 vvij = dnrm2(3, vij, 1)
870 tij = vvij/params % rsph(jsph)
871 sij = vij/vvij
872 vtij = vij*params % kappa
873
874 ! differentiated geometrical quantity wrt the nuclear position
875 grad_tij(:) = sij(:)/params % rsph(jsph)
876
877 ! gradient of F_ij with respect to vij (Bessel analogue of the
878 ! sum of the t^l*Y_lm(s) derivatives in the ddCOSMO case)
879 call fmm_l2p_bessel_grad(vtij, params % rsph(jsph)*params % kappa, &
880 & params % lmax, constants % vscales, params % kappa, &
881 & xe(:, jsph), zero, alpha)
882
883 ! value of F_ij
884 call fmm_l2p_bessel_work(vtij, params % lmax, constants % vscales, &
885 & constants % SI_ri(:, jsph), one, xe(:, jsph), zero, beta, &
886 & work_complex, work)
887
888 a = xadj_e(ig,isph)*constants % wgrid(ig)
889 tmp_fx(:) = - a*oij*alpha(:)
890 e = a*beta
891
892 ! add the switching derivative for the ij case
893 d_i = constants % switching % d_ni(ig,isph)
894 f_i = constants % switching % f_ni(ig,isph)
895 a = e/(d_i + f_i)
896 if (tij.ge.tlow .and. tij.le.thigh) then
897 grad_p = dfsw(tij, params % se, params % eta)
898 b = a*grad_p
899 tmp_fx(:) = tmp_fx(:) - b*grad_tij(:)
900 end if
901 xij = fsw(tij, params % se, params % eta)
902 if (xij .ne. zero) then
903 b = a*xij/(d_i + f_i)
904 do ki = constants % inl(isph), constants % inl(isph+1) - 1
905 ksph = constants % nl(ki)
906 vik = params % csph(:,isph) + &
907 & params % rsph(isph)*constants % cgrid(:,ig) - &
908 & params % csph(:,ksph)
909 vvik = dnrm2(3, vik, 1)
910 tik = vvik/params % rsph(ksph)
911 if (tik.ge.tlow .and. tik.le.thigh) then
912 xik = fsw(tik, params % se, params % eta)
913 sik = vik/vvik
914 grad_tik(:) = sik(:)/params % rsph(ksph)
915 grad_p = dfsw(tik, params % se, params % eta)
916
917 if (grad_p.ne.zero) then
918 c = b*(one - d_i/(one - xik))*grad_p
919 else
920 c = zero
921 end if
922
923 tmp_fx(:) = tmp_fx(:) + c*grad_tik
924 end if
925 end do
926 end if
927
928 force(:) = force(:) + tmp_fx(:)
929 dr = dr + dot_product(tmp_fx, constants % cgrid(:,ig))
930 end do
931 end do
932
933 ! ji cases
934 do ig = 1, params % ngrid
935 do ji = constants % inl(isph), constants % inl(isph+1) - 1
936 jsph = constants % nl(ji)
937
938 ! compute the geometrical quantites required
939 oji = compute_omega(params, constants, jsph, isph, ig)
940 if (oji.eq.zero) cycle
941 vji = params % csph(:,jsph) + &
942 & params % rsph(jsph)*constants % cgrid(:,ig) - &
943 & params % csph(:,isph)
944 vvji = dnrm2(3, vji, 1)
945 tji = vvji/params % rsph(isph)
946 sji = vji/vvji
947 vtji = vji*params % kappa
948
949 ! differentiated geometrical quantities wrt the nuclear position
950 grad_tji(:) = -sji/params % rsph(isph)
951
952 ! differentiated geometrical quantities wrt the radius
953 dr_tji = -tji/params % rsph(isph)
954
955 ! gradient of F_ji with respect to vji
956 call fmm_l2p_bessel_grad(vtji, params % rsph(isph)*params % kappa, &
957 & params % lmax, constants % vscales, params % kappa, &
958 & xe(:, isph), zero, alpha)
959
960 ! value of F_ji
961 call fmm_l2p_bessel_work(vtji, params % lmax, constants % vscales, &
962 & constants % SI_ri(:, isph), one, xe(:, isph), zero, beta, &
963 & work_complex, work)
964
965 ! derivative of F_ji with respect to r_i (normalization only)
966 call fmm_l2p_bessel_work(vtji, params % lmax, constants % vscales, &
967 & constants % SI_ri(:, isph), one, xe_dr, zero, beta_dr, &
968 & work_complex, work)
969
970 a = xadj_e(ig,jsph)*constants % wgrid(ig)
971 tmp_fx(:) = a*oji*alpha(:)
972 tmp_dr = - a*oji*beta_dr
973 e = a*beta
974
975 ! add the switching derivative for the ji case
976 d_j = constants % switching % d_ni(ig,jsph)
977 f_j = constants % switching % f_ni(ig,jsph)
978 a = e/(d_j + f_j)
979 if (tji.ge.tlow .and. tji.le.thigh) then
980 grad_p = dfsw(tji, params % se, params % eta)
981 b = a*grad_p
982 tmp_fx(:) = tmp_fx(:) - b*grad_tji(:)
983 tmp_dr = tmp_dr - b*dr_tji
984 end if
985
986 force(:) = force(:) + tmp_fx(:)
987 dr = dr + tmp_dr
988 end do
989 end do
990
991 ! additional switching contributions
992 tmp_fx = zero
993 tmp_dr = zero
994 do ig = 1, params % ngrid
995 do ij = constants % inl(isph), constants % inl(isph+1) - 1
996 jsph = constants % nl(ij)
997
998 ! compute ji geometrical quantities
999 vji = params % csph(:,jsph) + &
1000 & params % rsph(jsph)*constants % cgrid(:,ig) - &
1001 & params % csph(:,isph)
1002 vvji = dnrm2(3, vji, 1)
1003 tji = vvji/params % rsph(isph)
1004 sji = vji/vvji
1005
1006 if (.not.(tji.ge.tlow .and. tji.le.thigh)) cycle
1007
1008 ! differentiated geometrical quantities wrt the nuclear position
1009 grad_tji(:) = -sji/params % rsph(isph)
1010
1011 ! differentiated geometrical quantities wrt the radius
1012 dr_tji = -tji/params % rsph(isph)
1013
1014 xji = fsw(tji, params % se, params % eta)
1015
1016 d_j = constants % switching % d_ni(ig,jsph)
1017 f_j = constants % switching % f_ni(ig,jsph)
1018 grad_p = dfsw(tji, params % se, params % eta)
1019
1020 if (grad_p.eq.zero) cycle
1021
1022 bj = zero
1023 do kj = constants % inl(jsph), constants % inl(jsph+1) - 1
1024 ksph = constants % nl(kj)
1025 if (jsph.eq.ksph) cycle
1026
1027 ! compute the remaining geometrical quantities
1028 vjk = params % csph(:,jsph) + &
1029 & params % rsph(jsph)*constants % cgrid(:,ig) - &
1030 & params % csph(:,ksph)
1031 vvjk = dnrm2(3, vjk, 1)
1032 tjk = vvjk/params % rsph(ksph)
1033
1034 xjk = fsw(tjk, params % se, params % eta)
1035 if (xjk.eq.zero) cycle
1036
1037 ! value of F_jk (Bessel expansion of sphere k at the point
1038 ! of sphere j)
1039 vtjk = vjk*params % kappa
1040 call fmm_l2p_bessel_work(vtjk, params % lmax, &
1041 & constants % vscales, constants % SI_ri(:, ksph), one, &
1042 & xe(:, ksph), zero, beta_jk, work_complex, work)
1043
1044 e = xadj_e(ig,jsph)*constants % wgrid(ig)*beta_jk
1045 bj = bj + e*xjk
1046 end do
1047
1048 c = bj/(d_j + f_j)**2*(one - d_j/(one - xji))*grad_p
1049 tmp_fx(:) = tmp_fx(:) + c*grad_tji
1050 tmp_dr = tmp_dr + c*dr_tji
1051 end do
1052 end do
1053
1054 force(:) = force(:) + tmp_fx(:)
1055 dr = dr + tmp_dr
1056
1057end subroutine contract_grad_b_new
1058
1075subroutine contract_grad_c(params, constants, workspace, Xr, Xe, Xadj_r_sgrid, &
1076 & Xadj_e_sgrid, Xadj_r, Xadj_e, force, diff_re, ddx_error)
1077 !! input/output
1078 type(ddx_params_type), intent(in) :: params
1079 type(ddx_constants_type), intent(in) :: constants
1080 type(ddx_workspace_type), intent(inout) :: workspace
1081 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xr, Xe
1082 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: Xadj_r_sgrid, Xadj_e_sgrid
1083 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xadj_r, Xadj_e
1084 real(dp), dimension(3, params % nsph), intent(inout) :: force
1085 real(dp), dimension(constants % nbasis, params % nsph), intent(out) :: diff_re
1086 type(ddx_error_type), intent(inout) :: ddx_error
1087
1088 call contract_grad_c_worker2(params, constants, workspace, xr, xe, xadj_r_sgrid, &
1089 & xadj_e_sgrid, xadj_r, xadj_e, force, diff_re, ddx_error)
1090 if (ddx_error % flag .ne. 0) then
1091 call update_error(ddx_error, &
1092 & "contract_grad_C_worker2 returned an error, exiting")
1093 return
1094 end if
1095 call contract_grad_c_worker1(params, constants, workspace, xadj_r_sgrid, &
1096 & xadj_e_sgrid, diff_re, force, ddx_error)
1097 if (ddx_error % flag .ne. 0) then
1098 call update_error(ddx_error, &
1099 & "contract_grad_C_worker1 returned an error, exiting")
1100 return
1101 end if
1102
1103end subroutine contract_grad_c
1104
1105
1116subroutine contract_grad_f(params, constants, workspace, sol_adj, sol_sgrid, &
1117 & gradpsi, normal_hessian_cav, force, state, ddx_error)
1118 type(ddx_params_type), intent(in) :: params
1119 type(ddx_constants_type), intent(in) :: constants
1120 type(ddx_workspace_type), intent(inout) :: workspace
1121 type(ddx_state_type), intent(inout) :: state
1122 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: sol_adj
1123 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: sol_sgrid
1124 real(dp), dimension(3, constants % ncav), intent(in) :: gradpsi
1125 real(dp), dimension(3, constants % ncav), intent(in) :: normal_hessian_cav
1126 real(dp), dimension(3, params % nsph), intent(inout) :: force
1127 type(ddx_error_type), intent(inout) :: ddx_error
1128
1129 call contract_grad_f_worker1(params, constants, workspace, sol_adj, sol_sgrid, &
1130 & gradpsi, force, ddx_error)
1131 if (ddx_error % flag .ne. 0) then
1132 call update_error(ddx_error, &
1133 & "contract_grad_f_worker1 returned an error, exiting")
1134 return
1135 end if
1136
1137 call contract_grad_f_worker2(params, constants, gradpsi, &
1138 & normal_hessian_cav, force, state, ddx_error)
1139 if (ddx_error % flag .ne. 0) then
1140 call update_error(ddx_error, &
1141 & "contract_grad_f_worker2 returned an error, exiting")
1142 return
1143 end if
1144
1145end subroutine contract_grad_f
1146
1155subroutine contract_gradi_bik(params, constants, isph, Xe, Xadj_e, force)
1156 !! input/output
1157 type(ddx_params_type), intent(in) :: params
1158 type(ddx_constants_type), intent(in) :: constants
1159 integer, intent(in) :: isph
1160 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xe
1161 real(dp), dimension(params % ngrid), intent(in) :: Xadj_e
1162 real(dp), dimension(3), intent(inout) :: force
1163
1164 ! Local Variables
1165 integer :: igrid, ineigh, jsph
1166 real(dp), dimension(0:params % lmax) :: SI_rijn
1167 real(dp), dimension(0:params % lmax) :: DI_rijn
1168 ! beta : Eq.(53) Stamm.etal.18
1169 ! tlow : Lower bound for switch region
1170 ! thigh : Upper bound for switch region
1171 ! f1 : First factor in alpha computation
1172 ! f2 : Second factor in alpha computation
1173 real(dp) :: rijn, tij, beta, tlow, thigh, xij, oij, f1, f2, f3
1174 ! alpha : Eq.(52) Stamm.etal.18
1175 ! va : Eq.(54) Stamm.etal.18
1176 real(dp) :: vij(3), sij(3), alpha(3), va(3), rj, vtij(3)
1177 real(dp), external :: dnrm2
1178 real(dp) :: work(params % lmax+1)
1179 complex(dp) :: work_complex(params % lmax+1)
1180
1181 si_rijn = 0
1182 di_rijn = 0
1183
1184 tlow = one - pt5*(one - params % se)*params % eta
1185 thigh = one + pt5*(one + params % se)*params % eta
1186
1187 ! Loop over grid points
1188 do igrid = 1, params % ngrid
1189 va = zero
1190 do ineigh = constants % inl(isph), constants % inl(isph+1) - 1
1191 jsph = constants % nl(ineigh)
1192 vij = params % csph(:,isph) + &
1193 & params % rsph(isph)*constants % cgrid(:,igrid) - &
1194 & params % csph(:,jsph)
1195 rijn = dnrm2(3, vij, 1)
1196 tij = rijn/params % rsph(jsph)
1197 rj = params % rsph(jsph)
1198 if (tij.ge.thigh) cycle
1199 ! Computation of modified spherical Bessel function values
1200 if (tij.ne.zero) then
1201 sij = vij/rijn
1202 else
1203 sij = one
1204 end if
1205 vtij = vij*params % kappa
1206 call fmm_l2p_bessel_grad(vtij, params % rsph(jsph)*params % kappa, &
1207 & params % lmax, constants % vscales, params % kappa, &
1208 & xe(:, jsph), zero, alpha)
1209 call fmm_l2p_bessel_work(vtij, params % lmax, constants % vscales, &
1210 & constants % SI_ri(:, jsph), one, xe(:, jsph), zero, beta, &
1211 & work_complex, work)
1212 xij = fsw(tij, params % se, params % eta)
1213 if (constants % fi(igrid,isph).gt.one) then
1214 oij = xij/constants % fi(igrid,isph)
1215 f2 = -oij/constants % fi(igrid,isph)
1216 else
1217 oij = xij
1218 f2 = zero
1219 end if
1220 f1 = oij
1221 va(:) = va(:) + f1*alpha(:) + beta*f2*constants % zi(:,igrid,isph)
1222 if (tij .gt. tlow) then
1223 f3 = beta*dfsw(tij,params % se,params % eta)/params % rsph(jsph)
1224 if (constants % fi(igrid,isph).gt.one) &
1225 & f3 = f3/constants % fi(igrid,isph)
1226 va(:) = va(:) + f3*sij(:)
1227 end if
1228 end do
1229 force = force - constants % wgrid(igrid)*xadj_e(igrid)*va(:)
1230 end do
1231end subroutine contract_gradi_bik
1232
1233
1243subroutine contract_gradi_bji(params, constants, isph, Xe, Xadj_e, force)
1244 !! input/output
1245 type(ddx_params_type), intent(in) :: params
1246 type(ddx_constants_type), intent(in) :: constants
1247 integer, intent(in) :: isph
1248 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xe
1249 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: Xadj_e
1250 real(dp), dimension(3), intent(inout) :: force
1251
1252 ! Local Variables
1253 ! jk : Row pointer over kth row
1254 integer :: igrid, jsph, ksph, ineigh, jk
1255 real(dp), dimension(0:params % lmax) :: SI_rjin, SI_rjkn
1256 real(dp), dimension(0:params % lmax) :: DI_rjin, DI_rjkn
1257
1258 logical :: proc
1259 ! fac : \delta_fj_n*\omega^\eta_ji
1260 ! f1 : First factor in alpha computation
1261 ! f2 : Second factor in alpha computation
1262 ! beta_ji : Eq.(57) Stamm.etal.18
1263 ! dj : Before Eq.(10) Stamm.etal.18
1264 real(dp) :: rjin, tji, xji, oji, fac, f1, f2, beta_ji, dj, tlow, thigh
1265 real(dp) :: b, beta_jk, g1, g2, rjkn, tjk, xjk
1266 ! alpha : Eq.(56) Stamm.etal.18
1267 ! vb : Eq.(60) Stamm.etal.18
1268 ! vc : Eq.(59) Stamm.etal.18
1269 real(dp) :: vji(3), sji(3), vjk(3), alpha(3), vb(3), vc(3), &
1270 & vtji(3), vtjk(3)
1271 ! rho : Argument for ylmbas
1272 ! ctheta : Argument for ylmbas
1273 ! stheta : Argument for ylmbas
1274 ! cphi : Argument for ylmbas
1275 ! sphi : Argument for ylmbas
1276 real(dp) :: ri
1277
1278 real(dp), external :: dnrm2
1279 real(dp) :: work(params % lmax+1)
1280 complex(dp) :: work_complex(params % lmax+1)
1281
1282 si_rjin = 0
1283 di_rjin = 0
1284 si_rjkn = 0
1285 di_rjkn = 0
1286
1287 tlow = one - pt5*(one - params % se)*params % eta
1288 thigh = one + pt5*(one + params % se)*params % eta
1289
1290 do igrid = 1, params % ngrid
1291 vb = zero
1292 vc = zero
1293 do ineigh = constants % inl(isph), constants % inl(isph+1) - 1
1294 jsph = constants % nl(ineigh)
1295 vji = params % csph(:,jsph) + &
1296 & params % rsph(jsph)*constants % cgrid(:,igrid) - &
1297 & params % csph(:,isph)
1298 rjin = dnrm2(3, vji, 1)
1299 ri = params % rsph(isph)
1300 tji = rjin/ri
1301 if (tji.gt.thigh) cycle
1302 if (tji.ne.zero) then
1303 sji = vji/rjin
1304 else
1305 sji = one
1306 end if
1307 vtji = vji*params % kappa
1308 call fmm_l2p_bessel_grad(vtji, params % rsph(isph)*params % kappa, &
1309 & params % lmax, constants % vscales, params % kappa, &
1310 & xe(:, isph), zero, alpha)
1311 xji = fsw(tji,params % se,params % eta)
1312 if (constants % fi(igrid,jsph).gt.one) then
1313 oji = xji/constants % fi(igrid,jsph)
1314 else
1315 oji = xji
1316 end if
1317 f1 = oji
1318 vb = vb + f1*alpha*xadj_e(igrid,jsph)
1319 if (tji .gt. tlow) then
1320 call fmm_l2p_bessel_work(vtji, params % lmax, &
1321 & constants % vscales, constants % SI_ri(:, isph), one, &
1322 & xe(:, isph), zero, beta_ji, work_complex, work)
1323 if (constants % fi(igrid,jsph) .gt. one) then
1324 dj = one/constants % fi(igrid,jsph)
1325 fac = dj*xji
1326 proc = .false.
1327 b = zero
1328 do jk = constants % inl(jsph), constants % inl(jsph+1) - 1
1329 ksph = constants % nl(jk)
1330 vjk = params % csph(:,jsph) + &
1331 & params % rsph(jsph)*constants % cgrid(:,igrid) - &
1332 & params % csph(:,ksph)
1333 rjkn = dnrm2(3, vjk, 1)
1334 tjk = rjkn/params % rsph(ksph)
1335 if (ksph.ne.isph) then
1336 if (tjk .le. thigh) then
1337 proc = .true.
1338 vtjk = vjk*params % kappa
1339 call fmm_l2p_bessel_work(vtjk, params % lmax, &
1340 & constants % vscales, &
1341 & constants % SI_ri(:, ksph), one, xe(:, ksph), &
1342 & zero, beta_jk, work_complex, work)
1343 xjk = fsw(tjk, params % se, params % eta)
1344 b = b + beta_jk*xjk
1345 end if
1346 end if
1347 end do
1348 if (proc) then
1349 g1 = dj*dj*dfsw(tji,params % se,params % eta) &
1350 & /params % rsph(isph)
1351 g2 = g1*xadj_e(igrid,jsph)*b
1352 vc = vc + g2*sji
1353 end if
1354 else
1355 dj = one
1356 fac = zero
1357 end if
1358 f2 = (one-fac)*dj*dfsw(tji,params % se,params % eta) &
1359 & /params % rsph(isph)
1360 vb = vb + f2*xadj_e(igrid,jsph)*beta_ji*sji
1361 end if
1362 end do
1363 force = force + constants % wgrid(igrid)*(vb - vc)
1364 end do
1365end subroutine contract_gradi_bji
1366
1367
1379subroutine contract_grad_c_worker1(params, constants, workspace, Xadj_r_sgrid, &
1380 & Xadj_e_sgrid, diff_re, force, ddx_error)
1381 !! Inputs
1382 type(ddx_params_type), intent(in) :: params
1383 type(ddx_constants_type), intent(in) :: constants
1384 type(ddx_workspace_type), intent(inout) :: workspace
1385 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: &
1386 & diff_re
1387 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: &
1388 & Xadj_r_sgrid, Xadj_e_sgrid
1389 real(dp), dimension(3, params % nsph), intent(inout) :: force
1390 type(ddx_error_type), intent(inout) :: ddx_error
1391 ! Local variable
1392 ! igrid0: Index for grid point n0
1393 integer :: isph, jsph, igrid, l0, m0, ind0, igrid0, icav, &
1394 & indl, inode, istat
1395 ! term : SK_rijn/SK_rj
1396 ! termi : DI_ri/SI_ri
1397 ! termk : DK_ri/SK_ri
1398 ! sum_int : Intermediate sum
1399 ! sum_r : Intermediate sum for Laplace
1400 ! sum_e : Intermediate sum for HSP
1401 real(dp) :: sum_int, sum_r, sum_e
1402 real(dp) :: rijn
1403 real(dp) :: vij(3), sij(3), vtij(3)
1404
1405 ! local allocatable
1406 ! phi_n_r : Phi corresponding to Laplace problem
1407 ! phi_n_e : Phi corresponding to HSP problem
1408 ! coefY_d : sum_{l0m0} C_ik*term*Y_l0m0^j(x_in)*Y_l0m0(s_n)
1409 ! diff_re_sgrid : diff_re evaluated at grid point
1410 real(dp), allocatable :: phi_n_r(:,:), phi_n_e(:,:), coefY_d(:,:,:), &
1411 & diff_re_sgrid(:,:)
1412
1413 ! basloc : Y_lm(s_n)
1414 ! vplm : Argument to call ylmbas
1415 real(dp), dimension(constants % nbasis):: basloc, vplm
1416 ! dbasloc : Derivative of Y_lm(s_n)
1417 real(dp), dimension(3, constants % nbasis):: dbasloc
1418 ! vcos : Argument to call ylmbas
1419 ! vsin : Argument to call ylmbas
1420 real(dp), dimension(params % lmax+1):: vcos, vsin
1421 ! SK_rijn : Besssel function of first kind for rijn
1422 ! DK_rijn : Derivative of Besssel function of first kind for rijn
1423 real(dp), dimension(0:params % lmax) :: SK_rijn, DK_rijn
1424 real(dp) :: coef(constants % nbasis0), work(constants % lmax0+1)
1425 complex(dp) :: work_complex(constants % lmax0+1)
1426
1427 allocate(phi_n_r(params % ngrid, params % nsph), &
1428 & phi_n_e(params % ngrid, params % nsph), &
1429 & diff_re_sgrid(params % ngrid, params % nsph), stat=istat)
1430 if (istat.ne.0) then
1431 call update_error(ddx_error, "allocation ddx_error in ddx_contract_grad_C_worker1")
1432 return
1433 end if
1434!
1435 sum_int = zero
1436 sum_r = zero
1437 sum_e = zero
1438 phi_n_r = zero
1439 phi_n_e = zero
1440 diff_re_sgrid = zero
1441 basloc = zero
1442 vplm = zero
1443 dbasloc = zero
1444 vcos = zero
1445 vsin = zero
1446 sk_rijn = zero
1447 dk_rijn = zero
1448
1449 if (params % fmm .eq. 0) then
1450 allocate(coefy_d(constants % ncav, params % ngrid, params % nsph), &
1451 & stat=istat)
1452 if (istat.ne.0) then
1453 call update_error(ddx_error, "allocation ddx_error in fmm ddx_contract_grad_C_worker1")
1454 return
1455 end if
1456 coefy_d = zero
1457 ! Compute summation over l0, m0
1458 ! Loop over the sphers j
1459 do jsph = 1, params % nsph
1460 ! Loop over the grid points n0
1461 do igrid0 = 1, params % ngrid
1462 icav = zero
1463 ! Loop over spheres i
1464 do isph = 1, params % nsph
1465 ! Loop over grid points n
1466 do igrid = 1, params % ngrid
1467 ! Check for U_i^{eta}(x_in)
1468 if(constants % ui(igrid, isph) .gt. zero) then
1469 icav = icav + 1
1470 vij = params % csph(:,isph) + &
1471 & params % rsph(isph)*constants % cgrid(:,igrid) - &
1472 & params % csph(:,jsph)
1473 rijn = sqrt(dot_product(vij,vij))
1474 sij = vij/rijn
1475
1476 do l0 = 0, constants % lmax0
1477 do m0 = -l0,l0
1478 ind0 = l0**2 + l0 + m0 + 1
1479 coef(ind0) = constants % vgrid(ind0, igrid0) * &
1480 & constants % C_ik(l0, jsph)
1481 end do
1482 end do
1483 vtij = vij*params % kappa
1484 call fmm_m2p_bessel_work(vtij, constants % lmax0, &
1485 & constants % vscales, constants % SK_ri(:, jsph), one, &
1486 & coef, zero, coefy_d(icav, igrid0, jsph), work_complex, work)
1487 end if
1488 end do
1489 end do
1490 end do
1491 end do
1492
1493 ! Compute phi_in
1494 ! Loop over spheres j
1495 do jsph = 1, params % nsph
1496 ! Loop over grid points n0
1497 do igrid0 = 1, params % ngrid
1498 icav = zero
1499 sum_r = zero
1500 sum_e = zero
1501 ! Loop over sphers i
1502 do isph = 1, params % nsph
1503 ! Loop over grid points n
1504 do igrid = 1, params % ngrid
1505 if(constants % ui(igrid, isph) .gt. zero) then
1506 icav = icav + 1
1507 sum_r = sum_r + coefy_d(icav, igrid0, jsph)*xadj_r_sgrid(igrid, isph) &
1508 & * constants % wgrid(igrid)*constants % ui(igrid, isph)
1509 sum_e = sum_e + coefy_d(icav, igrid0, jsph)*xadj_e_sgrid(igrid, isph) &
1510 & * constants % wgrid(igrid)*constants % ui(igrid, isph)
1511 end if
1512 end do
1513 end do
1514 phi_n_r(igrid0, jsph) = sum_r
1515 phi_n_e(igrid0, jsph) = sum_e
1516 end do
1517 end do
1518 else
1519 ! Compute phi_n_r at first
1520 ! Adjoint integration from spherical harmonics to grid points is not needed
1521 ! here as ygrid already contains grid values, we just need to scale it by
1522 ! weights of grid points
1523 do isph = 1, params % nsph
1524 workspace % tmp_grid(:, isph) = xadj_r_sgrid(:, isph) * &
1525 & constants % wgrid(:) * constants % ui(:, isph)
1526 end do
1527 ! Adjoint FMM
1528 call tree_m2p_bessel_adj(params, constants, constants % lmax0, one, &
1529 & workspace % tmp_grid, zero, params % lmax, workspace % tmp_sph)
1530 call tree_l2p_bessel_adj(params, constants, one, workspace % tmp_grid, &
1531 & zero, workspace % tmp_node_l)
1532 call tree_l2l_bessel_rotation_adj(params, constants, &
1533 & workspace % tmp_node_l)
1534 call tree_m2l_bessel_rotation_adj(params, constants, &
1535 & workspace % tmp_node_l, workspace % tmp_node_m)
1536 call tree_m2m_bessel_rotation_adj(params, constants, &
1537 & workspace % tmp_node_m)
1538 ! Properly load adjoint multipole harmonics into tmp_sph
1539 if(constants % lmax0 .lt. params % pm) then
1540 do isph = 1, params % nsph
1541 inode = constants % snode(isph)
1542 workspace % tmp_sph(1:constants % nbasis0, isph) = &
1543 & workspace % tmp_sph(1:constants % nbasis0, isph) + &
1544 & workspace % tmp_node_m(1:constants % nbasis0, inode)
1545 end do
1546 else
1547 indl = (params % pm+1)**2
1548 do isph = 1, params % nsph
1549 inode = constants % snode(isph)
1550 workspace % tmp_sph(1:indl, isph) = &
1551 & workspace % tmp_sph(1:indl, isph) + &
1552 & workspace % tmp_node_m(:, inode)
1553 end do
1554 end if
1555 ! Scale by C_ik
1556 do isph = 1, params % nsph
1557 do l0 = 0, constants % lmax0
1558 ind0 = l0*l0 + l0 + 1
1559 workspace % tmp_sph(ind0-l0:ind0+l0, isph) = &
1560 & workspace % tmp_sph(ind0-l0:ind0+l0, isph) * &
1561 & constants % C_ik(l0, isph)
1562 end do
1563 end do
1564 ! Multiply by vgrid
1565 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis0, &
1566 & one, constants % vgrid, constants % vgrid_nbasis, &
1567 & workspace % tmp_sph, constants % nbasis, zero, phi_n_r, params % ngrid)
1568 ! Compute phi_n_e now
1569 ! Adjoint integration from spherical harmonics to grid points is not needed
1570 ! here as ygrid already contains grid values, we just need to scale it by
1571 ! weights of grid points
1572 do isph = 1, params % nsph
1573 workspace % tmp_grid(:, isph) = xadj_e_sgrid(:, isph) * &
1574 & constants % wgrid(:) * constants % ui(:, isph)
1575 end do
1576 ! Adjoint FMM
1577 call tree_m2p_bessel_adj(params, constants, constants % lmax0, one, &
1578 & workspace % tmp_grid, zero, params % lmax, workspace % tmp_sph)
1579 call tree_l2p_bessel_adj(params, constants, one, workspace % tmp_grid, &
1580 & zero, workspace % tmp_node_l)
1581 call tree_l2l_bessel_rotation_adj(params, constants, &
1582 & workspace % tmp_node_l)
1583 call tree_m2l_bessel_rotation_adj(params, constants, &
1584 & workspace % tmp_node_l, workspace % tmp_node_m)
1585 call tree_m2m_bessel_rotation_adj(params, constants, &
1586 & workspace % tmp_node_m)
1587 ! Properly load adjoint multipole harmonics into tmp_sph
1588 if(constants % lmax0 .lt. params % pm) then
1589 do isph = 1, params % nsph
1590 inode = constants % snode(isph)
1591 workspace % tmp_sph(1:constants % nbasis0, isph) = &
1592 & workspace % tmp_sph(1:constants % nbasis0, isph) + &
1593 & workspace % tmp_node_m(1:constants % nbasis0, inode)
1594 end do
1595 else
1596 indl = (params % pm+1)**2
1597 do isph = 1, params % nsph
1598 inode = constants % snode(isph)
1599 workspace % tmp_sph(1:indl, isph) = &
1600 & workspace % tmp_sph(1:indl, isph) + &
1601 & workspace % tmp_node_m(:, inode)
1602 end do
1603 end if
1604 ! Scale by C_ik
1605 do isph = 1, params % nsph
1606 do l0 = 0, constants % lmax0
1607 ind0 = l0*l0 + l0 + 1
1608 workspace % tmp_sph(ind0-l0:ind0+l0, isph) = &
1609 & workspace % tmp_sph(ind0-l0:ind0+l0, isph) * &
1610 & constants % C_ik(l0, isph)
1611 end do
1612 end do
1613 ! Multiply by vgrid
1614 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis0, &
1615 & one, constants % vgrid, constants % vgrid_nbasis, &
1616 & workspace % tmp_sph, constants % nbasis, zero, phi_n_e, params % ngrid)
1617 end if
1618 call dgemm('T', 'N', params % ngrid, params % nsph, &
1619 & constants % nbasis, one, constants % vgrid, constants % vgrid_nbasis, &
1620 & diff_re , constants % nbasis, zero, diff_re_sgrid, &
1621 & params % ngrid)
1622 do isph = 1, params % nsph
1623 call contract_grad_u(params, constants, isph, diff_re_sgrid, phi_n_r, force(:, isph))
1624 call contract_grad_u(params, constants, isph, diff_re_sgrid, phi_n_e, force(:, isph))
1625 end do
1626
1627 deallocate(phi_n_r, phi_n_e, diff_re_sgrid, stat=istat)
1628 if (istat.ne.0) then
1629 call update_error(ddx_error, "deallocation ddx_error in ddx_contract_grad_C_worker1")
1630 return
1631 end if
1632 if (allocated(coefy_d)) then
1633 deallocate(coefy_d, stat=istat)
1634 if (istat.ne.0) then
1635 call update_error(ddx_error, "deallocation ddx_error in ddx_contract_grad_C_worker1")
1636 return
1637 end if
1638 end if
1639end subroutine contract_grad_c_worker1
1640
1641
1660subroutine contract_grad_c_worker2(params, constants, workspace, Xr, Xe, &
1661 & Xadj_r_sgrid, Xadj_e_sgrid, Xadj_r, Xadj_e, force, diff_re, ddx_error)
1662 !! input/output
1663 type(ddx_params_type), intent(in) :: params
1664 type(ddx_constants_type), intent(in) :: constants
1665 type(ddx_workspace_type), intent(inout) :: workspace
1666 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: Xr, Xe
1667 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: &
1668 & Xadj_r_sgrid, Xadj_e_sgrid
1669 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: &
1670 & Xadj_r, Xadj_e
1671 real(dp), dimension(3, params % nsph), intent(inout) :: force
1672 real(dp), dimension(constants % nbasis, params % nsph), intent(out) :: &
1673 & diff_re
1674 type(ddx_error_type), intent(inout) :: ddx_error
1675 real(dp), external :: dnrm2
1676 ! Local variable
1677 integer :: isph, jsph, igrid, l, m, ind, l0, ind0, icav, indl, inode, &
1678 & ksph, knode, jnode, knear, jsph_node, istat
1679 ! val_dim3 : Intermediate value array of dimension 3
1680 real(dp), dimension(3) :: vij, vtij
1681 ! val : Intermediate variable to compute diff_ep
1682 real(dp) :: val
1683 ! large local are allocatable
1684 ! phi_in : sum_{j=1}^N diff0_j * coefY_j
1685 ! diff_ep_dim3 : 3 dimensional couterpart of diff_ep
1686 ! sum_dim3 : Storage of sum
1687 ! diff0 : dot_product([PU_j]_l0m0^l'm', l'/r_j[Xr]_jl'm' -
1688 ! (i'_l'(r_j)/i_l'(r_j))[Xe]_jl'm')
1689 real(dp), allocatable :: diff_ep_dim3(:,:), phi_in(:,:), sum_dim3(:,:,:), &
1690 & diff0(:,:), diff1(:,:), diff1_grad(:,:,:), l2l_grad(:,:,:)
1691 real(dp), dimension(0:params % lmax) :: SK_rijn, DK_rijn
1692 complex(dp) :: work_complex(constants % lmax0+2)
1693 real(dp) :: work(constants % lmax0+2)
1694
1695 allocate(diff_ep_dim3(3, constants % ncav), &
1696 & phi_in(params % ngrid, params % nsph), &
1697 & sum_dim3(3, constants % nbasis, params % nsph), &
1698 & diff0(constants % nbasis0, params % nsph), &
1699 & diff1(constants % nbasis0, params % nsph), &
1700 & diff1_grad((constants % lmax0+2)**2, 3, params % nsph), &
1701 & l2l_grad((params % pl+2)**2, 3, params % nsph), stat=istat)
1702 if (istat.ne.0) then
1703 call update_error(ddx_error, "allocation ddx_error in ddx_contract_grad_C_worker2")
1704 return
1705 end if
1706
1707 ! Setting initial values to zero
1708 sk_rijn = zero
1709 dk_rijn = zero
1710
1711 diff_re = zero
1712 ! Compute l'/r_j[Xr]_jl'm' -(i'_l'(r_j)/i_l'(r_j))[Xe]_jl'm'
1713 !$omp parallel do default(none) shared(params,diff_re,constants,xr,xe) &
1714 !$omp private(jsph,l,m,ind) schedule(dynamic)
1715 do jsph = 1, params % nsph
1716 do l = 0, params % lmax
1717 do m = -l,l
1718 ind = l**2 + l + m + 1
1719 diff_re(ind,jsph) = (params % epsp/params % eps)*(l/params % rsph(jsph)) * &
1720 & xr(ind,jsph) - constants % termimat(l,jsph)*xe(ind,jsph)
1721 end do
1722 end do
1723 end do
1724
1725 ! diff0 = Pchi * diff_re, linear scaling
1726 diff0 = zero
1727 !$omp parallel do default(none) shared(params,constants,diff_re,diff0, &
1728 !$omp diff1,diff1_grad) private(jsph,l0,ind0) schedule(dynamic)
1729 do jsph = 1, params % nsph
1730 do l0 = 0, constants % lmax0
1731 do ind0 = l0*l0+1, l0*l0+2*l0+1
1732 diff0(ind0, jsph) = dot_product(diff_re(:,jsph), &
1733 & constants % Pchi(:,ind0, jsph))
1734 diff1(ind0, jsph) = diff0(ind0, jsph) * constants % C_ik(l0, jsph)
1735 end do
1736 end do
1737 ! Prepare diff1_grad
1738 call fmm_m2m_bessel_grad(constants % lmax0, constants % SK_ri(:, jsph), &
1739 & constants % vscales, diff1(:, jsph), diff1_grad(:, :, jsph))
1740 end do
1741
1742 if (params % fmm .eq. 0) then
1743 ! phi_in = diff0 * coefY
1744 ! Here, summation over j takes place
1745 phi_in = zero
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 ! Extrenal grid point
1751 icav = icav + 1
1752 val = zero
1753 do jsph = 1, params % nsph
1754 !do ind0 = 1, constants % nbasis0
1755 !!====== This place requirs coefY, that is not precomputed anymore
1756 ! val = val + diff0(ind0,jsph)*constants % coefY(icav,ind0,jsph)
1757 !end do
1758 vij = params % csph(:, isph) + &
1759 & params % rsph(isph)*constants % cgrid(:, igrid) - &
1760 & params % csph(:, jsph)
1761 vtij = vij * params % kappa
1762 call fmm_m2p_bessel_work(vtij, constants % lmax0, &
1763 & constants % vscales, constants % SK_ri(:, jsph), one, &
1764 & diff1(:, jsph), one, val, work_complex, work)
1765 end do
1766 phi_in(igrid, isph) = val
1767 end if
1768 end do
1769 end do
1770 else
1771 ! phi_in
1772 ! Load input harmonics into tree data
1773 workspace % tmp_sph = zero
1774 workspace % tmp_sph(1:constants % nbasis0, :) = diff1(:, :)
1775 if(constants % lmax0 .lt. params % pm) then
1776 do isph = 1, params % nsph
1777 inode = constants % snode(isph)
1778 workspace % tmp_node_m(1:constants % nbasis0, inode) = &
1779 & workspace % tmp_sph(1:constants % nbasis0, isph)
1780 workspace % tmp_node_m(constants % nbasis0+1:, inode) = zero
1781 end do
1782 else
1783 indl = (params % pm+1)**2
1784 do isph = 1, params % nsph
1785 inode = constants % snode(isph)
1786 workspace % tmp_node_m(:, inode) = workspace % tmp_sph(1:indl, isph)
1787 end do
1788 end if
1789 ! Do FMM operations
1790 call tree_m2m_bessel_rotation(params, constants, workspace % tmp_node_m)
1791 call tree_m2l_bessel_rotation(params, constants, workspace % tmp_node_m, &
1792 & workspace % tmp_node_l)
1793 call tree_l2l_bessel_rotation(params, constants, workspace % tmp_node_l)
1794 call tree_l2p_bessel(params, constants, one, workspace % tmp_node_l, zero, &
1795 & phi_in)
1796 call tree_m2p_bessel(params, constants, constants % lmax0, one, &
1797 & params % lmax, workspace % tmp_sph, one, &
1798 & phi_in)
1799 ! Make phi_in zero at internal grid points
1800 !$omp parallel do default(none) shared(params,constants,phi_in) &
1801 !$omp private(isph,igrid) schedule(dynamic)
1802 do isph = 1, params % nsph
1803 do igrid = 1, params % ngrid
1804 if (constants % ui(igrid, isph) .eq. zero) then
1805 phi_in(igrid, isph) = zero
1806 end if
1807 end do
1808 end do
1809 ! Get gradients of the L2L
1810 !$omp parallel do default(none) shared(params,constants,workspace, &
1811 !$omp l2l_grad) private(isph,igrid,inode) schedule(dynamic)
1812 do isph = 1, params % nsph
1813 inode = constants % snode(isph)
1814 workspace % tmp_sph_l(:, isph) = workspace % tmp_node_l(:, inode)
1815 call fmm_l2l_bessel_grad(params % pl, &
1816 & constants % SI_ri(:, isph), constants % vscales, &
1817 & workspace % tmp_node_l(:, inode), &
1818 & l2l_grad(:, :, isph))
1819 end do
1820 workspace % tmp_sph = xadj_r + xadj_e
1821 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis, &
1822 & one, constants % vwgrid, constants % vgrid_nbasis, &
1823 & workspace % tmp_sph, constants % nbasis, zero, &
1824 & workspace % tmp_grid, params % ngrid)
1825 workspace % tmp_grid = workspace % tmp_grid * constants % ui
1826 ! Adjoint FMM with output tmp_sph2(:, :) which stores coefficients of
1827 ! harmonics of degree up to lmax+1
1828 call tree_m2p_bessel_nodiag_adj(params, constants, constants % lmax0+1, one, &
1829 & workspace % tmp_grid, zero, params % lmax+1, workspace % tmp_sph2)
1830 call tree_l2p_bessel_adj(params, constants, one, workspace % tmp_grid, zero, &
1831 & workspace % tmp_node_l)
1832 call tree_l2l_bessel_rotation_adj(params, constants, workspace % tmp_node_l)
1833 call tree_m2l_bessel_rotation_adj(params, constants, workspace % tmp_node_l, &
1834 & workspace % tmp_node_m)
1835 call tree_m2m_bessel_rotation_adj(params, constants, workspace % tmp_node_m)
1836 ! Properly load adjoint multipole harmonics into tmp_sph2 that holds
1837 ! harmonics of a degree up to lmax+1
1838 if(constants % lmax0+1 .lt. params % pm) then
1839 !$omp parallel do default(none) shared(params,constants, &
1840 !$omp workspace) private(isph,inode) schedule(dynamic)
1841 do isph = 1, params % nsph
1842 inode = constants % snode(isph)
1843 workspace % tmp_sph2(1:(constants % lmax0+2)**2, isph) = &
1844 & workspace % tmp_sph2(1:(constants % lmax0+2)**2, isph) + &
1845 & workspace % tmp_node_m(1:(constants % lmax0+2)**2, inode)
1846 end do
1847 else
1848 indl = (params % pm+1)**2
1849 !$omp parallel do default(none) shared(params,constants,indl, &
1850 !$omp workspace) private(isph,inode) schedule(dynamic)
1851 do isph = 1, params % nsph
1852 inode = constants % snode(isph)
1853 workspace % tmp_sph2(1:indl, isph) = &
1854 & workspace % tmp_sph2(1:indl, isph) + &
1855 & workspace % tmp_node_m(:, inode)
1856 end do
1857 end if
1858 end if
1859
1860 do ksph = 1, params % nsph
1861 ! Computation of derivative of U_i^e(x_in)
1862 call contract_grad_u(params, constants, ksph, xadj_r_sgrid, phi_in, force(:, ksph))
1863 call contract_grad_u(params, constants, ksph, xadj_e_sgrid, phi_in, force(:, ksph))
1864
1865 ! Aleksandr: my loop for the diff_ep_dim3
1866 diff_ep_dim3 = zero
1867 ! At first isph=ksph, jsph!=ksph
1868 icav = constants % icav_ia(ksph) - 1
1869 if (params % fmm .eq. 0) then
1870 do igrid = 1, params % ngrid
1871 if (constants % ui(igrid, ksph) .eq. zero) cycle
1872 icav = icav + 1
1873 do jsph = 1, params % nsph
1874 if (jsph .eq. ksph) cycle
1875 vij = params % csph(:,ksph) + &
1876 & params % rsph(ksph)*constants % cgrid(:,igrid) - &
1877 & params % csph(:,jsph)
1878 vtij = vij * params % kappa
1879 !call fmm_m2p_bessel_grad(vtij, &
1880 ! & params % rsph(jsph)*params % kappa, &
1881 ! & constants % lmax0, &
1882 ! & constants % vscales, params % kappa, diff1(:, jsph), one, &
1883 ! & diff_ep_dim3(:, icav))
1884 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1885 & constants % vscales, constants % SK_ri(:, jsph), &
1886 & -params % kappa, diff1_grad(:, 1, jsph), one, &
1887 & diff_ep_dim3(1, icav), work_complex, work)
1888 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1889 & constants % vscales, constants % SK_ri(:, jsph), &
1890 & -params % kappa, diff1_grad(:, 2, jsph), one, &
1891 & diff_ep_dim3(2, icav), work_complex, work)
1892 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1893 & constants % vscales, constants % SK_ri(:, jsph), &
1894 & -params % kappa, diff1_grad(:, 3, jsph), one, &
1895 & diff_ep_dim3(3, icav), work_complex, work)
1896 end do
1897 end do
1898 else
1899 knode = constants % snode(ksph)
1900 do igrid = 1, params % ngrid
1901 if (constants % ui(igrid, ksph) .eq. zero) cycle
1902 icav = icav + 1
1903 ! Far-field
1904 call dgemv('T', (params % pl+2)**2, 3, -params % kappa, &
1905 & l2l_grad(1, 1, ksph), &
1906 & (params % pl+2)**2, constants % vgrid(1, igrid), 1, &
1907 & one, diff_ep_dim3(1, icav), 1)
1908 ! Near-field
1909 do knear = constants % snear(knode), constants % snear(knode+1)-1
1910 jnode = constants % near(knear)
1911 do jsph_node = constants % cluster(1, jnode), &
1912 & constants % cluster(2, jnode)
1913 jsph = constants % order(jsph_node)
1914 if (jsph .eq. ksph) cycle
1915 vij = params % csph(:, ksph) - params % csph(:, jsph) + &
1916 & params % rsph(ksph)*constants % cgrid(:, igrid)
1917 vtij = vij * params % kappa
1918 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1919 & constants % vscales, constants % SK_ri(:, jsph), &
1920 & -params % kappa, diff1_grad(:, 1, jsph), one, &
1921 & diff_ep_dim3(1, icav), work_complex, work)
1922 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1923 & constants % vscales, constants % SK_ri(:, jsph), &
1924 & -params % kappa, diff1_grad(:, 2, jsph), one, &
1925 & diff_ep_dim3(2, icav), work_complex, work)
1926 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1927 & constants % vscales, constants % SK_ri(:, jsph), &
1928 & -params % kappa, diff1_grad(:, 3, jsph), one, &
1929 & diff_ep_dim3(3, icav), work_complex, work)
1930 end do
1931 end do
1932 end do
1933 end if
1934 sum_dim3 = zero
1935 icav = constants % icav_ia(ksph) - 1
1936 do igrid =1, params % ngrid
1937 if(constants % ui(igrid, ksph) .gt. zero) then
1938 icav = icav + 1
1939 do ind = 1, constants % nbasis
1940 sum_dim3(:,ind,ksph) = sum_dim3(:,ind,ksph) &
1941 & + diff_ep_dim3(:,icav)*constants % ui(igrid, ksph) &
1942 & *constants % vwgrid(ind, igrid)
1943 end do
1944 end if
1945 end do
1946 do ind = 1, constants % nbasis
1947 force(:, ksph) = force(:, ksph) + &
1948 & sum_dim3(:, ind, ksph)*(xadj_r(ind, ksph) + &
1949 & xadj_e(ind, ksph))
1950 end do
1951
1952 ! Now jsph=ksph and isph!=ksph
1953 if (params % fmm .eq. 0) then
1954 diff_ep_dim3 = zero
1955 sum_dim3 = zero
1956 do isph = 1, params % nsph
1957 if (isph .eq. ksph) cycle
1958 icav = constants % icav_ia(isph) - 1
1959 do igrid = 1, params % ngrid
1960 if (constants % ui(igrid, isph) .eq. zero) cycle
1961 icav = icav + 1
1962 vij = params % csph(:,isph) + &
1963 & params % rsph(isph)*constants % cgrid(:,igrid) - &
1964 & params % csph(:,ksph)
1965 vtij = vij * params % kappa
1966 !call fmm_m2p_bessel_grad(vij * params % kappa, &
1967 ! & params % rsph(ksph)*params % kappa, &
1968 ! & constants % lmax0, &
1969 ! & constants % vscales, -params % kappa, diff1(:, ksph), one, &
1970 ! & diff_ep_dim3(:, icav))
1971 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1972 & constants % vscales, constants % SK_ri(:, ksph), &
1973 & params % kappa, diff1_grad(:, 1, ksph), one, &
1974 & diff_ep_dim3(1, icav), work_complex, work)
1975 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1976 & constants % vscales, constants % SK_ri(:, ksph), &
1977 & params % kappa, diff1_grad(:, 2, ksph), one, &
1978 & diff_ep_dim3(2, icav), work_complex, work)
1979 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
1980 & constants % vscales, constants % SK_ri(:, ksph), &
1981 & params % kappa, diff1_grad(:, 3, ksph), one, &
1982 & diff_ep_dim3(3, icav), work_complex, work)
1983 end do
1984 end do
1985 icav = zero
1986 do isph = 1, params % nsph
1987 do igrid =1, params % ngrid
1988 if(constants % ui(igrid, isph) .gt. zero) then
1989 icav = icav + 1
1990 do ind = 1, constants % nbasis
1991 sum_dim3(:,ind,isph) = sum_dim3(:,ind,isph) &
1992 & + diff_ep_dim3(:,icav)*constants % ui(igrid, isph) &
1993 & *constants % vwgrid(ind, igrid)
1994 end do
1995 end if
1996 end do
1997 end do
1998 ! Computation of derivative of \bf(k)_j^l0(x_in)\times Y^j_l0m0(x_in)
1999 do isph = 1, params % nsph
2000 do ind = 1, constants % nbasis
2001 force(:, ksph) = force(:, ksph) &
2002 & + sum_dim3(:, ind, isph)*(xadj_r(ind, isph) &
2003 & + xadj_e(ind, isph))
2004 end do
2005 end do
2006 else
2007 call dgemv('T', (constants % lmax0+2)**2, 3, params % kappa, &
2008 & diff1_grad(1, 1, ksph), (constants % lmax0+2)**2, &
2009 & workspace % tmp_sph2(1, ksph), 1, one, force(1, ksph), 1)
2010 end if
2011 end do
2012 deallocate(diff_ep_dim3, phi_in, sum_dim3, diff1_grad, l2l_grad, &
2013 & diff0, diff1, stat=istat)
2014 if (istat.ne.0) then
2015 call update_error(ddx_error, "deallocation ddx_error in ddx_contract_grad_C_worker2")
2016 return
2017 end if
2018end subroutine contract_grad_c_worker2
2019
2033subroutine contract_grad_f_worker1(params, constants, workspace, sol_adj, sol_sgrid, &
2034 & gradpsi, force, ddx_error)
2035 ! input/output
2036 type(ddx_params_type), intent(in) :: params
2037 type(ddx_constants_type), intent(in) :: constants
2038 type(ddx_workspace_type), intent(inout) :: workspace
2039 real(dp), dimension(constants % nbasis, params % nsph), intent(in) :: &
2040 & sol_adj
2041 real(dp), dimension(params % ngrid, params % nsph), intent(in) :: sol_sgrid
2042 real(dp), dimension(3, constants % ncav), intent(in) :: gradpsi
2043 real(dp), dimension(3, params % nsph), intent(inout) :: force
2044 type(ddx_error_type), intent(inout) :: ddx_error
2045
2046 ! local
2047 real(dp), external :: dnrm2
2048 integer :: isph, jsph, igrid, ind, l0, ind0, icav, ksph, &
2049 & knode, jnode, knear, jsph_node, indl, inode, istat
2050 ! val_dim3 : Intermediate value array of dimension 3
2051 real(dp), dimension(3) :: vij, vtij
2052 ! val : Intermediate variable to compute diff_ep
2053 ! nderpsi : Derivative of psi on grid points
2054 real(dp) :: nderpsi, sum_int
2055
2056 ! local allocatable
2057 ! phi_in : sum_{j=1}^N diff0_j * coefY_j
2058 ! diff_ep_dim3 : 3 dimensional couterpart of diff_ep
2059 ! sum_dim3 : Storage of sum
2060 ! Debug purpose
2061 ! These variables can be taken from the subroutine update_rhs
2062 ! diff0 : dot_product([PU_j]_l0m0^l'm', l'/r_j[Xr]_jl'm' -
2063 ! (i'_l'(r_j)/i_l'(r_j))[Xe]_jl'm')
2064 ! sum_Sjin : \sum_j [S]_{jin} Eq.~(97) [QSM20.SISC]
2065 ! c0 : \sum_{n=1}^N_g w_n U_j^{x_nj}\partial_n psi_0(x_nj)Y_{l0m0}(s_n)
2066 real(dp), allocatable :: phi_in(:,:), diff_ep_dim3(:,:), &
2067 & sum_dim3(:,:,:), diff0(:,:), sum_sjin(:,:), &
2068 & c0_d(:,:), c0_d1(:,:), c0_d1_grad(:,:,:), l2l_grad(:,:,:)
2069 real(dp), dimension(0:params % lmax) :: SK_rijn, DK_rijn
2070 complex(dp) :: work_complex(constants % lmax0 + 2)
2071 real(dp) :: work(constants % lmax0 + 2)
2072
2073 allocate(phi_in(params % ngrid, params % nsph), &
2074 & diff_ep_dim3(3, constants % ncav), &
2075 & sum_dim3(3, constants % nbasis, params % nsph), &
2076 & c0_d(constants % nbasis0, params % nsph), &
2077 & c0_d1(constants % nbasis0, params % nsph), &
2078 & c0_d1_grad((constants % lmax0+2)**2, 3, params % nsph), &
2079 & sum_sjin(params % ngrid, params % nsph), &
2080 & l2l_grad((params % pl+2)**2, 3, params % nsph), &
2081 & diff0(constants % nbasis0, params % nsph), stat=istat)
2082 if (istat.ne.0) then
2083 call update_error(ddx_error, "allocation ddx_error in ddx_grad_f_worker1")
2084 return
2085 end if
2086
2087
2088 ! Setting initial values to zero
2089 sk_rijn = zero
2090 dk_rijn = zero
2091 c0_d = zero
2092 c0_d1 = zero
2093 c0_d1_grad = zero
2094
2095 icav = zero
2096 do isph = 1, params % nsph
2097 do igrid= 1, params % ngrid
2098 if ( constants % ui(igrid,isph) .gt. zero ) then
2099 icav = icav + 1
2100 nderpsi = dot_product( gradpsi(:,icav),constants % cgrid(:,igrid) )
2101 c0_d(:, isph) = c0_d(:,isph) + &
2102 & constants % wgrid(igrid)* &
2103 & constants % ui(igrid,isph)*&
2104 & nderpsi* &
2105 & constants % vgrid(1:constants % nbasis0,igrid)
2106 do l0 = 0, constants % lmax0
2107 ind0 = l0*l0 + l0 + 1
2108 c0_d1(ind0-l0:ind0+l0, isph) = c0_d(ind0-l0:ind0+l0, isph) * &
2109 & constants % C_ik(l0, isph)
2110 end do
2111 ! Prepare c0_d1_grad
2112 call fmm_m2m_bessel_grad(constants % lmax0, constants % SK_ri(:, isph), &
2113 & constants % vscales, c0_d1(:, isph), c0_d1_grad(:, :, isph))
2114 end if
2115 end do
2116 end do
2117
2118 if (params % fmm .eq. 0) then
2119 ! Compute [S]_{jin}
2120 icav = 0
2121 do isph = 1, params % nsph
2122 do igrid = 1, params % ngrid
2123 if (constants % ui(igrid,isph).gt.zero) then
2124 icav = icav + 1
2125 sum_int = zero
2126 ! Loop to compute Sijn
2127 do jsph = 1, params % nsph
2128 vij = params % csph(:, isph) + &
2129 & params % rsph(isph)*constants % cgrid(:, igrid) - &
2130 & params % csph(:, jsph)
2131 vtij = vij * params % kappa
2132 call fmm_m2p_bessel_work(vtij, constants % lmax0, &
2133 & constants % vscales, constants % SK_ri(:, jsph), one, &
2134 & c0_d1(:, jsph), one, sum_int, work_complex, work)
2135 end do ! End of loop jsph
2136 sum_sjin(igrid,isph) = -(params % epsp/params % eps)*sum_int
2137 end if
2138 end do ! End of loop igrid
2139 end do ! End of loop isph
2140 else
2141 ! Load input harmonics into tree data
2142 workspace % tmp_sph = zero
2143 workspace % tmp_sph(1:constants % nbasis0, :) = c0_d1(:, :)
2144 if(constants % lmax0 .lt. params % pm) then
2145 do isph = 1, params % nsph
2146 inode = constants % snode(isph)
2147 workspace % tmp_node_m(1:constants % nbasis0, inode) = &
2148 & workspace % tmp_sph(1:constants % nbasis0, isph)
2149 workspace % tmp_node_m(constants % nbasis0+1:, inode) = zero
2150 end do
2151 else
2152 indl = (params % pm+1)**2
2153 do isph = 1, params % nsph
2154 inode = constants % snode(isph)
2155 workspace % tmp_node_m(:, inode) = workspace % tmp_sph(1:indl, isph)
2156 end do
2157 end if
2158 ! Do FMM operations
2159 call tree_m2m_bessel_rotation(params, constants, workspace % tmp_node_m)
2160 call tree_m2l_bessel_rotation(params, constants, workspace % tmp_node_m, &
2161 & workspace % tmp_node_l)
2162 call tree_l2l_bessel_rotation(params, constants, workspace % tmp_node_l)
2163 call tree_l2p_bessel(params, constants, -params % epsp/params % eps, workspace % tmp_node_l, zero, &
2164 & sum_sjin)
2165 call tree_m2p_bessel(params, constants, constants % lmax0, -params % epsp/params % eps, &
2166 & params % lmax, workspace % tmp_sph, one, &
2167 & sum_sjin)
2168 ! Make phi_in zero at internal grid points
2169 do isph = 1, params % nsph
2170 do igrid = 1, params % ngrid
2171 if (constants % ui(igrid, isph) .eq. zero) then
2172 sum_sjin(igrid, isph) = zero
2173 end if
2174 end do
2175 end do
2176 ! Get gradients of the L2L
2177 do isph = 1, params % nsph
2178 inode = constants % snode(isph)
2179 workspace % tmp_sph_l(:, isph) = workspace % tmp_node_l(:, inode)
2180 call fmm_l2l_bessel_grad(params % pl, &
2181 & constants % SI_ri(:, isph), constants % vscales, &
2182 & workspace % tmp_node_l(:, inode), &
2183 & l2l_grad(:, :, isph))
2184 end do
2185 workspace % tmp_sph = sol_adj
2186 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis, &
2187 & one, constants % vwgrid, constants % vgrid_nbasis, &
2188 & workspace % tmp_sph, constants % nbasis, zero, &
2189 & workspace % tmp_grid, params % ngrid)
2190 workspace % tmp_grid = workspace % tmp_grid * constants % ui
2191 ! Adjoint FMM with output tmp_sph2(:, :) which stores coefficients of
2192 ! harmonics of degree up to lmax+1
2193 call tree_m2p_bessel_nodiag_adj(params, constants, constants % lmax0+1, one, &
2194 & workspace % tmp_grid, zero, params % lmax+1, workspace % tmp_sph2)
2195 call tree_l2p_bessel_adj(params, constants, one, workspace % tmp_grid, zero, &
2196 & workspace % tmp_node_l)
2197 call tree_l2l_bessel_rotation_adj(params, constants, workspace % tmp_node_l)
2198 call tree_m2l_bessel_rotation_adj(params, constants, workspace % tmp_node_l, &
2199 & workspace % tmp_node_m)
2200 call tree_m2m_bessel_rotation_adj(params, constants, workspace % tmp_node_m)
2201 ! Properly load adjoint multipole harmonics into tmp_sph2 that holds
2202 ! harmonics of a degree up to lmax+1
2203 if(constants % lmax0+1 .lt. params % pm) then
2204 do isph = 1, params % nsph
2205 inode = constants % snode(isph)
2206 workspace % tmp_sph2(1:(constants % lmax0+2)**2, isph) = &
2207 & workspace % tmp_sph2(1:(constants % lmax0+2)**2, isph) + &
2208 & workspace % tmp_node_m(1:(constants % lmax0+2)**2, inode)
2209 end do
2210 else
2211 indl = (params % pm+1)**2
2212 do isph = 1, params % nsph
2213 inode = constants % snode(isph)
2214 workspace % tmp_sph2(1:indl, isph) = &
2215 & workspace % tmp_sph2(1:indl, isph) + &
2216 & workspace % tmp_node_m(:, inode)
2217 end do
2218 end if
2219 end if
2220 do ksph = 1, params % nsph
2221 ! Computation of derivative of U_i^e(x_in)
2222 call contract_grad_u(params, constants, ksph, sol_sgrid, sum_sjin, force(:, ksph))
2223 ! Aleksandr: my loop for the diff_ep_dim3
2224 diff_ep_dim3 = zero
2225 ! At first isph=ksph, jsph!=ksph
2226 icav = constants % icav_ia(ksph) - 1
2227 if (params % fmm .eq. 0) then
2228 do igrid = 1, params % ngrid
2229 if (constants % ui(igrid, ksph) .eq. zero) cycle
2230 icav = icav + 1
2231 do jsph = 1, params % nsph
2232 if (jsph .eq. ksph) cycle
2233 vij = params % csph(:,ksph) + &
2234 & params % rsph(ksph)*constants % cgrid(:,igrid) - &
2235 & params % csph(:,jsph)
2236 vtij = vij * params % kappa
2237 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
2238 & constants % vscales, constants % SK_ri(:, jsph), &
2239 & -params % kappa, c0_d1_grad(:, 1, jsph), one, &
2240 & diff_ep_dim3(1, icav), work_complex, work)
2241 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
2242 & constants % vscales, constants % SK_ri(:, jsph), &
2243 & -params % kappa, c0_d1_grad(:, 2, jsph), one, &
2244 & diff_ep_dim3(2, icav), work_complex, work)
2245 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
2246 & constants % vscales, constants % SK_ri(:, jsph), &
2247 & -params % kappa, c0_d1_grad(:, 3, jsph), one, &
2248 & diff_ep_dim3(3, icav), work_complex, work)
2249 end do
2250 end do
2251 else
2252 knode = constants % snode(ksph)
2253 do igrid = 1, params % ngrid
2254 if (constants % ui(igrid, ksph) .eq. zero) cycle
2255 icav = icav + 1
2256 ! Far-field
2257 call dgemv('T', (params % pl+2)**2, 3, -params % kappa, &
2258 & l2l_grad(1, 1, ksph), &
2259 & (params % pl+2)**2, constants % vgrid(1, igrid), 1, &
2260 & one, diff_ep_dim3(1, icav), 1)
2261 ! Near-field
2262 do knear = constants % snear(knode), constants % snear(knode+1)-1
2263 jnode = constants % near(knear)
2264 do jsph_node = constants % cluster(1, jnode), &
2265 & constants % cluster(2, jnode)
2266 jsph = constants % order(jsph_node)
2267 if (jsph .eq. ksph) cycle
2268 vij = params % csph(:, ksph) - params % csph(:, jsph) + &
2269 & params % rsph(ksph)*constants % cgrid(:, igrid)
2270 vtij = vij * params % kappa
2271 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
2272 & constants % vscales, constants % SK_ri(:, jsph), &
2273 & -params % kappa, c0_d1_grad(:, 1, jsph), one, &
2274 & diff_ep_dim3(1, icav), work_complex, work)
2275 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
2276 & constants % vscales, constants % SK_ri(:, jsph), &
2277 & -params % kappa, c0_d1_grad(:, 2, jsph), one, &
2278 & diff_ep_dim3(2, icav), work_complex, work)
2279 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
2280 & constants % vscales, constants % SK_ri(:, jsph), &
2281 & -params % kappa, c0_d1_grad(:, 3, jsph), one, &
2282 & diff_ep_dim3(3, icav), work_complex, work)
2283 end do
2284 end do
2285 end do
2286 end if
2287 sum_dim3 = zero
2288 icav = constants % icav_ia(ksph) - 1
2289 do igrid =1, params % ngrid
2290 if(constants % ui(igrid, ksph) .gt. zero) then
2291 icav = icav + 1
2292 do ind = 1, constants % nbasis
2293 sum_dim3(:,ind,ksph) = sum_dim3(:,ind,ksph) &
2294 & - (params % epsp/params % eps) &
2295 & *diff_ep_dim3(:,icav)*constants % ui(igrid, ksph) &
2296 & *constants % vwgrid(ind, igrid)
2297 end do
2298 end if
2299 end do
2300 do ind = 1, constants % nbasis
2301 force(:, ksph) = force(:, ksph) + &
2302 & sum_dim3(:, ind, ksph)*sol_adj(ind, ksph)
2303 end do
2304
2305 ! Now jsph=ksph and isph!=ksph
2306 if (params % fmm .eq. 0) then
2307 diff_ep_dim3 = zero
2308 sum_dim3 = zero
2309 do isph = 1, params % nsph
2310 if (isph .eq. ksph) cycle
2311 icav = constants % icav_ia(isph) - 1
2312 do igrid = 1, params % ngrid
2313 if (constants % ui(igrid, isph) .eq. zero) cycle
2314 icav = icav + 1
2315 vij = params % csph(:,isph) + &
2316 & params % rsph(isph)*constants % cgrid(:,igrid) - &
2317 & params % csph(:,ksph)
2318 vtij = vij * params % kappa
2319 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
2320 & constants % vscales, constants % SK_ri(:, ksph), &
2321 & params % kappa, c0_d1_grad(:, 1, ksph), one, &
2322 & diff_ep_dim3(1, icav), work_complex, work)
2323 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
2324 & constants % vscales, constants % SK_ri(:, ksph), &
2325 & params % kappa, c0_d1_grad(:, 2, ksph), one, &
2326 & diff_ep_dim3(2, icav), work_complex, work)
2327 call fmm_m2p_bessel_work(vtij, constants % lmax0+1, &
2328 & constants % vscales, constants % SK_ri(:, ksph), &
2329 & params % kappa, c0_d1_grad(:, 3, ksph), one, &
2330 & diff_ep_dim3(3, icav), work_complex, work)
2331 end do
2332 end do
2333
2334 icav = zero
2335 do isph = 1, params % nsph
2336 do igrid =1, params % ngrid
2337 if(constants % ui(igrid, isph) .gt. zero) then
2338 icav = icav + 1
2339 do ind = 1, constants % nbasis
2340 sum_dim3(:,ind,isph) = sum_dim3(:,ind,isph) &
2341 & - (params % epsp/params % eps) &
2342 & *diff_ep_dim3(:,icav)*constants % ui(igrid, isph) &
2343 & *constants % vwgrid(ind, igrid)
2344 end do
2345 end if
2346 end do
2347 end do
2348
2349 ! Computation of derivative of \bf(k)_j^l0(x_in)\times Y^j_l0m0(x_in)
2350 do isph = 1, params % nsph
2351 do ind = 1, constants % nbasis
2352 force(:, ksph) = force(:, ksph) + sum_dim3(:, ind, isph)*sol_adj(ind, isph)
2353 end do
2354 end do
2355 else
2356 call dgemv('T', (constants % lmax0+2)**2, 3, -params % epsp/params % eps*params % kappa, &
2357 & c0_d1_grad(1, 1, ksph), (constants % lmax0+2)**2, &
2358 & workspace % tmp_sph2(1, ksph), 1, one, force(1, ksph), 1)
2359 end if
2360 end do
2361
2362 deallocate(phi_in, diff_ep_dim3, sum_dim3, c0_d, c0_d1, &
2363 & c0_d1_grad, sum_sjin, l2l_grad, diff0, stat=istat)
2364 if (istat.ne.0) then
2365 call update_error(ddx_error, "deallocation ddx_error in ddx_grad_f_worker1")
2366 return
2367 end if
2368
2369end subroutine contract_grad_f_worker1
2370
2371!! @param[inout] ddx_error: ddX error
2372subroutine contract_grad_f_worker2(params, constants, &
2373 & gradpsi, normal_hessian_cav, force, state, ddx_error)
2374 type(ddx_params_type), intent(in) :: params
2375 type(ddx_constants_type), intent(in) :: constants
2376 real(dp), intent(in) :: gradpsi(3, constants % ncav), &
2377 & normal_hessian_cav(3, constants % ncav)
2378 real(dp), intent(inout) :: force(3, params % nsph)
2379 type(ddx_state_type), intent(inout) :: state
2380 type(ddx_error_type), intent(inout) :: ddx_error
2381
2382 integer :: icav, isph, igrid, istat
2383 real(dp) :: nderpsi
2384 real(dp), allocatable :: gradpsi_grid(:, :)
2385
2386 allocate(gradpsi_grid(params % ngrid, params % nsph), stat=istat)
2387 if (istat.ne.0) then
2388 call update_error(ddx_error, "allocation ddx_error in ddx_grad_f_worker2")
2389 return
2390 end if
2391
2392 icav = 0
2393 do isph = 1, params % nsph
2394 do igrid = 1, params % ngrid
2395 if(constants % ui(igrid, isph) .gt. zero) then
2396 icav = icav + 1
2397 nderpsi = dot_product(gradpsi(:, icav), &
2398 & constants % cgrid(:, igrid))
2399 gradpsi_grid(igrid, isph) = nderpsi
2400 end if
2401 end do
2402 end do
2403
2404 icav = 0
2405 do isph = 1, params % nsph
2406 call contract_grad_u(params, constants, isph, gradpsi_grid, &
2407 & state % phi_n, force(:, isph))
2408 do igrid = 1, params % ngrid
2409 if(constants % ui(igrid, isph) .gt. zero) then
2410 icav = icav + 1
2411 force(:, isph) = force(:, isph) &
2412 & + constants % wgrid(igrid)*constants % ui(igrid, isph) &
2413 & *state % phi_n(igrid, isph)*normal_hessian_cav(:, icav)
2414 end if
2415 end do
2416 end do
2417
2418 deallocate(gradpsi_grid, stat=istat)
2419 if (istat.ne.0) then
2420 call update_error(ddx_error, "deallocation ddx_error in ddx_grad_f_worker2")
2421 return
2422 end if
2423
2424end subroutine contract_grad_f_worker2
2425
2427subroutine gradr_sph(params, constants, isph, vplm, vcos, vsin, basloc, &
2428 & dbsloc, g, ygrid, fx)
2429 implicit none
2430 ! Inputs
2431 type(ddx_params_type), intent(in) :: params
2432 type(ddx_constants_type), intent(in) :: constants
2433 integer, intent(in) :: isph
2434 real(dp), intent(in) :: g(constants % nbasis, params % nsph), &
2435 & ygrid(params % ngrid, params % nsph)
2436 real(dp), intent(inout) :: vplm(constants % nbasis), vcos(params % lmax+1), &
2437 & vsin(params % lmax+1), basloc(constants % nbasis), &
2438 & dbsloc(3, constants % nbasis), fx(3)
2439 ! various scratch arrays
2440 real(dp) vik(3), sik(3), vki(3), ski(3), vkj(3), skj(3), vji(3), &
2441 & sji(3), va(3), vb(3), a(3)
2442 ! jacobian matrix
2443 real(dp) sjac(3,3)
2444 ! indexes
2445 integer its, ik, ksph, l, m, ind, jsph, icomp, jcomp
2446 ! various scalar quantities
2447 real(dp) cx, cy, cz, vvki, tki, gg, fl, fac, vvkj, tkj
2448 real(dp) tt, fcl, fjj, gi, fii, vvji, tji, qji
2449 real(dp) b, vvik, tik, qik, tlow, thigh, duj
2450 real(dp) :: rho, ctheta, stheta, cphi, sphi
2451 real(dp), external :: dnrm2
2452 real(dp) :: d_k, f_k, chi_ki, grad_p
2453
2454 tlow = one - pt5*(one - params % se)*params % eta
2455 thigh = one + pt5*(one + params % se)*params % eta
2456
2457 ! first set of contributions:
2458 ! diagonal block, kc and part of kb
2459
2460 fx = zero
2461 do its = 1, params % ngrid
2462 ! sum over ksph in neighbors of isph
2463 do ik = constants % inl(isph), constants % inl(isph+1) - 1
2464 ksph = constants % nl(ik)
2465 ! build geometrical quantities
2466 cx = params % csph(1,ksph) + params % rsph(ksph)*constants % cgrid(1,its)
2467 cy = params % csph(2,ksph) + params % rsph(ksph)*constants % cgrid(2,its)
2468 cz = params % csph(3,ksph) + params % rsph(ksph)*constants % cgrid(3,its)
2469 vki(1) = cx - params % csph(1,isph)
2470 vki(2) = cy - params % csph(2,isph)
2471 vki(3) = cz - params % csph(3,isph)
2472 !vvki = sqrt(vki(1)*vki(1) + vki(2)*vki(2) + &
2473 ! & vki(3)*vki(3))
2474 vvki = dnrm2(3, vki, 1)
2475 tki = vvki/params % rsph(isph)
2476
2477 ! contributions involving grad i of uk come from the switching
2478 ! region.
2479 ! note: ui avoids contributions from points that are in the
2480 ! switching between isph and ksph but are buried in a third
2481 ! sphere.
2482 if ((tki.gt.tlow).and.(tki.lt.thigh) .and. &
2483 & constants % ui(its,ksph).gt.zero) then
2484 ! other geometrical quantities
2485 ski = vki/vvki
2486
2487 ! diagonal block kk contribution, with k in n(i)
2488 gg = zero
2489 do l = 0, params % lmax
2490 ind = l*l + l + 1
2491 fl = dble(l)
2492 fac = twopi/(two*fl + one)
2493 do m = -l, l
2494 !! DEBUG comment
2495 gg = gg + fac*constants % vgrid(ind+m,its)*g(ind+m,ksph)
2496 end do
2497 end do
2498
2499 ! kc contribution
2500 do jsph = 1, params % nsph
2501 if (jsph.ne.ksph .and. jsph.ne.isph) then
2502 vkj(1) = cx - params % csph(1,jsph)
2503 vkj(2) = cy - params % csph(2,jsph)
2504 vkj(3) = cz - params % csph(3,jsph)
2505 vvkj = sqrt(vkj(1)*vkj(1) + vkj(2)*vkj(2) + &
2506 & vkj(3)*vkj(3))
2507 vvkj = dnrm2(3, vkj, 1)
2508 tkj = vvkj/params % rsph(jsph)
2509 skj = vkj/vvkj
2510 call ylmbas(skj, rho, ctheta, stheta, cphi, sphi, &
2511 & params % lmax, constants % vscales, basloc, &
2512 & vplm, vcos, vsin)
2513 tt = one/tkj
2514 do l = 0, params % lmax
2515 ind = l*l + l + 1
2516 fcl = - fourpi*dble(l)/(two*dble(l)+one)*tt
2517 do m = -l, l
2518 !! DEBUG comment
2519 gg = gg + fcl*g(ind+m,jsph)*basloc(ind+m)
2520 end do
2521 tt = tt/tkj
2522 end do
2523 !call fmm_m2p(vkj, params % rsph(jsph), &
2524 ! & params % lmax, constants % vscales_rel, -one, &
2525 ! & g(:, jsph), one, gg)
2526 end if
2527 end do
2528
2529 ! part of kb contribution
2530 call ylmbas(ski, rho, ctheta, stheta, cphi, sphi, &
2531 & params % lmax, constants % vscales, basloc, &
2532 & vplm, vcos, vsin)
2533 tt = one/tki
2534 do l = 0, params % lmax
2535 ind = l*l + l + 1
2536 fcl = - four*pi*dble(l)/(two*dble(l)+one)*tt
2537 do m = -l, l
2538 !! DEBUG comment
2539 gg = gg + fcl*g(ind+m,isph)*basloc(ind+m)
2540 end do
2541 tt = tt/tki
2542 end do
2543 !call fmm_m2p(vki, params % rsph(isph), &
2544 ! & params % lmax, constants % vscales_rel, -one, &
2545 ! & g(:, isph), one, gg)
2546
2547 ! common step, product with grad i uj
2548 if (params%switching.eq.0) then
2549 duj = dfsw(tki,params % se, params % eta)/params % rsph(isph)
2550 else
2551 d_k = constants%switching%d_ni(its, ksph)
2552 f_k = constants%switching%f_ni(its, ksph)
2553 chi_ki = fsw(tki,params%se,params%eta)
2554 grad_p = dfsw(tki,params%se,params%eta)
2555 duj = one/(d_k + f_k)**2 &
2556 & *(d_k + d_k*f_k/(one - chi_ki)) &
2557 & *grad_p/params%rsph(isph)
2558 end if
2559 fjj = duj*constants % wgrid(its)*gg*ygrid(its,ksph)
2560 fx(1) = fx(1) - fjj*ski(1)
2561 fx(2) = fx(2) - fjj*ski(2)
2562 fx(3) = fx(3) - fjj*ski(3)
2563 end if
2564 end do
2565
2566 ! diagonal block ii contribution
2567 if (constants % ui(its,isph).gt.zero.and.constants % ui(its,isph).lt.one) then
2568 gi = zero
2569 do l = 0, params % lmax
2570 ind = l*l + l + 1
2571 fl = dble(l)
2572 fac = twopi/(two*fl + one)
2573 do m = -l, l
2574 !! DEBUG comment
2575 gi = gi + fac*constants % vgrid(ind+m,its)*g(ind+m,isph)
2576 !gi = gi + pt5*constants % vgrid2(ind+m,its)*g(ind+m,isph)
2577 end do
2578 end do
2579 !do l = 0, (params % lmax+1)**2
2580 ! gi = gi + constants % vgrid2(l, its)*g(l, isph)
2581 !end do
2582 !gi = pt5 * gi
2583 fii = constants % wgrid(its)*gi*ygrid(its,isph)
2584 fx(1) = fx(1) + fii*constants % zi(1,its,isph)
2585 fx(2) = fx(2) + fii*constants % zi(2,its,isph)
2586 fx(3) = fx(3) + fii*constants % zi(3,its,isph)
2587 end if
2588 end do
2589
2590 ! second set of contributions:
2591 ! part of kb and ka
2592 do its = 1, params % ngrid
2593
2594 ! run over all the spheres except isph
2595 do jsph = 1, params % nsph
2596 if (constants % ui(its,jsph).gt.zero .and. jsph.ne.isph) then
2597 ! build geometrical quantities
2598 cx = params % csph(1,jsph) + params % rsph(jsph)*constants % cgrid(1,its)
2599 cy = params % csph(2,jsph) + params % rsph(jsph)*constants % cgrid(2,its)
2600 cz = params % csph(3,jsph) + params % rsph(jsph)*constants % cgrid(3,its)
2601 vji(1) = cx - params % csph(1,isph)
2602 vji(2) = cy - params % csph(2,isph)
2603 vji(3) = cz - params % csph(3,isph)
2604 !vvji = sqrt(vji(1)*vji(1) + vji(2)*vji(2) + &
2605 ! & vji(3)*vji(3))
2606 vvji = dnrm2(3, vji, 1)
2607 tji = vvji/params % rsph(isph)
2608 qji = one/vvji
2609 sji = vji/vvji
2610
2611 ! build the jacobian of sji
2612 sjac = zero
2613 sjac(1,1) = - one
2614 sjac(2,2) = - one
2615 sjac(3,3) = - one
2616 do icomp = 1, 3
2617 do jcomp = 1, 3
2618 sjac(icomp,jcomp) = qji*(sjac(icomp,jcomp) &
2619 & + sji(icomp)*sji(jcomp))
2620 end do
2621 end do
2622
2623 ! assemble the local basis and its gradient
2624 !call dbasis(sji,basloc,dbsloc,vplm,vcos,vsin)
2625 call dbasis(params, constants, sji,basloc,dbsloc,vplm,vcos,vsin)
2626
2627 ! assemble the contribution
2628 a = zero
2629 tt = one/(tji)
2630 do l = 0, params % lmax
2631 ind = l*l + l + 1
2632 fl = dble(l)
2633 fcl = - tt*fourpi*fl/(two*fl + one)
2634 do m = -l, l
2635 fac = fcl*g(ind+m,isph)
2636 b = (fl + one)*basloc(ind+m)/(params % rsph(isph)*tji)
2637
2638 ! apply the jacobian to grad y
2639 va(1) = sjac(1,1)*dbsloc(1,ind+m) + &
2640 & sjac(1,2)*dbsloc(2,ind+m) + sjac(1,3)*dbsloc(3,ind+m)
2641 va(2) = sjac(2,1)*dbsloc(1,ind+m) + &
2642 & sjac(2,2)*dbsloc(2,ind+m) + sjac(2,3)*dbsloc(3,ind+m)
2643 va(3) = sjac(3,1)*dbsloc(1,ind+m) + &
2644 & sjac(3,2)*dbsloc(2,ind+m) + sjac(3,3)*dbsloc(3,ind+m)
2645 a(1) = a(1) + fac*(sji(1)*b + va(1))
2646 a(2) = a(2) + fac*(sji(2)*b + va(2))
2647 a(3) = a(3) + fac*(sji(3)*b + va(3))
2648 end do
2649 tt = tt/tji
2650 end do
2651 fac = constants % ui(its,jsph)*constants % wgrid(its)*ygrid(its,jsph)
2652 fx(1) = fx(1) - fac*a(1)
2653 fx(2) = fx(2) - fac*a(2)
2654 fx(3) = fx(3) - fac*a(3)
2655 end if
2656 end do
2657 end do
2658
2659 ! ka contribution
2660 do its = 1, params % ngrid
2661 cx = params % csph(1,isph) + params % rsph(isph)*constants % cgrid(1,its)
2662 cy = params % csph(2,isph) + params % rsph(isph)*constants % cgrid(2,its)
2663 cz = params % csph(3,isph) + params % rsph(isph)*constants % cgrid(3,its)
2664 a = zero
2665
2666 ! iterate on all the spheres except isph
2667 do ksph = 1, params % nsph
2668 if (constants % ui(its,isph).gt.zero .and. ksph.ne.isph) then
2669 ! geometrical stuff
2670 vik(1) = cx - params % csph(1,ksph)
2671 vik(2) = cy - params % csph(2,ksph)
2672 vik(3) = cz - params % csph(3,ksph)
2673 !vvik = sqrt(vik(1)*vik(1) + vik(2)*vik(2) + &
2674 ! & vik(3)*vik(3))
2675 vvik = dnrm2(3, vik, 1)
2676 tik = vvik/params % rsph(ksph)
2677 qik = one/vvik
2678 sik = vik/vvik
2679
2680 ! build the jacobian of sik
2681 sjac = zero
2682 sjac(1,1) = one
2683 sjac(2,2) = one
2684 sjac(3,3) = one
2685 do icomp = 1, 3
2686 do jcomp = 1, 3
2687 sjac(icomp,jcomp) = qik*(sjac(icomp,jcomp) &
2688 & - sik(icomp)*sik(jcomp))
2689 end do
2690 end do
2691
2692 ! if we are in the switching region, recover grad_i u_i
2693 vb = zero
2694 if (constants % ui(its,isph).lt.one) then
2695 vb(1) = constants % zi(1,its,isph)
2696 vb(2) = constants % zi(2,its,isph)
2697 vb(3) = constants % zi(3,its,isph)
2698 end if
2699
2700 ! assemble the local basis and its gradient
2701 !call dbasis(sik,basloc,dbsloc,vplm,vcos,vsin)
2702 call dbasis(params, constants, sik,basloc,dbsloc,vplm,vcos,vsin)
2703
2704 ! assemble the contribution
2705 tt = one/(tik)
2706 do l = 0, params % lmax
2707 ind = l*l + l + 1
2708 fl = dble(l)
2709 fcl = - tt*fourpi*fl/(two*fl + one)
2710 do m = -l, l
2711 fac = fcl*g(ind+m,ksph)
2712 fac = - fac*basloc(ind+m)
2713 a(1) = a(1) + fac*vb(1)
2714 a(2) = a(2) + fac*vb(2)
2715 a(3) = a(3) + fac*vb(3)
2716
2717 fac = constants % ui(its,isph)*fcl*g(ind+m,ksph)
2718 b = - (fl + one)*basloc(ind+m)/(params % rsph(ksph)*tik)
2719
2720 ! apply the jacobian to grad y
2721 va(1) = sjac(1,1)*dbsloc(1,ind+m) + &
2722 & sjac(1,2)*dbsloc(2,ind+m) + sjac(1,3)*dbsloc(3,ind+m)
2723 va(2) = sjac(2,1)*dbsloc(1,ind+m) + &
2724 & sjac(2,2)*dbsloc(2,ind+m) + sjac(2,3)*dbsloc(3,ind+m)
2725 va(3) = sjac(3,1)*dbsloc(1,ind+m) + &
2726 & sjac(3,2)*dbsloc(2,ind+m) + sjac(3,3)*dbsloc(3,ind+m)
2727 a(1) = a(1) + fac*(sik(1)*b + va(1))
2728 a(2) = a(2) + fac*(sik(2)*b + va(2))
2729 a(3) = a(3) + fac*(sik(3)*b + va(3))
2730 end do
2731 tt = tt/tik
2732 end do
2733 end if
2734 end do
2735 fac = constants % wgrid(its)*ygrid(its,isph)
2736 fx(1) = fx(1) - fac*a(1)
2737 fx(2) = fx(2) - fac*a(2)
2738 fx(3) = fx(3) - fac*a(3)
2739 end do
2740end subroutine gradr_sph
2741
2743subroutine gradr_fmm(params, constants, workspace, g, ygrid, fx)
2744 implicit none
2745 ! Inputs
2746 type(ddx_params_type), intent(in) :: params
2747 type(ddx_constants_type), intent(in) :: constants
2748 real(dp), intent(in) :: g(constants % nbasis, params % nsph), &
2749 & ygrid(params % ngrid, params % nsph)
2750 ! Temporaries
2751 type(ddx_workspace_type), intent(inout) :: workspace
2752 ! Output
2753 real(dp), intent(out) :: fx(3, params % nsph)
2754 ! Local variables
2755 integer :: indl, indl1, l, isph, igrid, ik, ksph, &
2756 & jsph, jsph_node
2757 integer :: inear, inode, jnode
2758 real(dp) :: gg, c(3), vki(3), vvki, tki, gg3(3), tmp_gg, tmp_c(3)
2759 real(dp) :: d_k, f_k, chi_ki, grad_p, duj
2760 real(dp) :: tlow, thigh
2761 real(dp), dimension(3, 3) :: zx_coord_transform, zy_coord_transform
2762 real(dp), external :: ddot, dnrm2
2763 real(dp) :: work(params % lmax+2)
2764 !real(dp) :: l2g(params % ngrid, params % nsph)
2765 zx_coord_transform = 0
2766 zx_coord_transform(3, 2) = 1
2767 zx_coord_transform(2, 3) = 1
2768 zx_coord_transform(1, 1) = 1
2769 zy_coord_transform = 0
2770 zy_coord_transform(1, 2) = 1
2771 zy_coord_transform(2, 1) = 1
2772 zy_coord_transform(3, 3) = 1
2773 tlow = one - pt5*(one - params % se)*params % eta
2774 thigh = one + pt5*(one + params % se)*params % eta
2775 fx = zero
2776 !! Scale input harmonics at first
2777 workspace % tmp_sph(1, :) = zero
2778 indl = 2
2779 do l = 1, params % lmax
2780 indl1 = (l+1)**2
2781 workspace % tmp_sph(indl:indl1, :) = l * g(indl:indl1, :)
2782 indl = indl1 + 1
2783 end do
2784 !! Compute gradient of M2M of tmp_sph harmonics at the origin and store it
2785 !! in tmp_sph_grad. tmp_sph_grad(:, 1, :), tmp_sph_grad(:, 2, :) and
2786 !! tmp_sph_grad(:, 3, :) correspond to the OX, OY and OZ axes. Variable
2787 !! tmp_sph2 is a temporary workspace here.
2788 call tree_grad_m2m(params, constants, workspace % tmp_sph, &
2789 & workspace % tmp_sph_grad, workspace % tmp_sph2)
2790 !! Adjoint full FMM matvec to get output multipole expansions from input
2791 !! external grid points. It is used to compute R_i^B fast as a contraction
2792 !! of a gradient stored in tmp_sph_grad and a result of adjoint matvec.
2793 ! Adjoint integration from spherical harmonics to grid points is not needed
2794 ! here as ygrid already contains grid values, we just need to scale it by
2795 ! weights of grid points
2796 do isph = 1, params % nsph
2797 workspace % tmp_grid(:, isph) = ygrid(:, isph) * &
2798 & constants % wgrid(:) * constants % ui(:, isph)
2799 end do
2800 ! Adjoint FMM with output tmp_sph2(:, :) which stores coefficients of
2801 ! harmonics of degree up to lmax+1
2802 call tree_m2p_adj(params, constants, params % lmax+1, one, &
2803 & workspace % tmp_grid, zero, workspace % tmp_sph2)
2804 call tree_l2p_adj(params, constants, one, workspace % tmp_grid, zero, &
2805 & workspace % tmp_node_l, workspace % tmp_sph_l)
2806 call tree_l2l_rotation_adj(params, constants, workspace % tmp_node_l)
2807 call tree_m2l_rotation_adj(params, constants, workspace % tmp_node_l, &
2808 & workspace % tmp_node_m)
2809 call tree_m2m_rotation_adj(params, constants, workspace % tmp_node_m)
2810 ! Properly load adjoint multipole harmonics into tmp_sph2 that holds
2811 ! harmonics of a degree up to lmax+1
2812 if(params % lmax+1 .lt. params % pm) then
2813 do isph = 1, params % nsph
2814 inode = constants % snode(isph)
2815 workspace % tmp_sph2(:, isph) = workspace % tmp_sph2(:, isph) + &
2816 & workspace % tmp_node_m(1:constants % grad_nbasis, inode)
2817 end do
2818 else
2819 indl = (params % pm+1)**2
2820 do isph = 1, params % nsph
2821 inode = constants % snode(isph)
2822 workspace % tmp_sph2(1:indl, isph) = &
2823 & workspace % tmp_sph2(1:indl, isph) + &
2824 & workspace % tmp_node_m(:, inode)
2825 end do
2826 end if
2827 ! Compute second term of R_i^B as a contraction
2828 do isph = 1, params % nsph
2829 call dgemv('T', constants % grad_nbasis, 3, one, &
2830 & workspace % tmp_sph_grad(1, 1, isph), constants % grad_nbasis, &
2831 & workspace % tmp_sph2(1, isph), 1, zero, fx(1, isph), 1)
2832 end do
2833 !! Direct far-field FMM matvec to get output local expansions from input
2834 !! multipole expansions. It will be used in R_i^A.
2835 !! As of now I compute potential at all external grid points, improved
2836 !! version shall only compute it at external points in a switch region
2837 ! Load input harmonics into tree data
2838 if(params % lmax .lt. params % pm) then
2839 do isph = 1, params % nsph
2840 inode = constants % snode(isph)
2841 workspace % tmp_node_m(:constants % nbasis, inode) = &
2842 & workspace % tmp_sph(:, isph)
2843 workspace % tmp_node_m(constants % nbasis+1:, inode) = zero
2844 end do
2845 else
2846 indl = (params % pm+1)**2
2847 do isph = 1, params % nsph
2848 inode = constants % snode(isph)
2849 workspace % tmp_node_m(:, inode) = workspace % tmp_sph(1:indl, isph)
2850 end do
2851 end if
2852 ! Perform direct FMM matvec to all external grid points
2853 call tree_m2m_rotation(params, constants, workspace % tmp_node_m)
2854 call tree_m2l_rotation(params, constants, workspace % tmp_node_m, &
2855 & workspace % tmp_node_l)
2856 call tree_l2l_rotation(params, constants, workspace % tmp_node_l)
2857 call tree_l2p(params, constants, one, workspace % tmp_node_l, zero, &
2858 & workspace % tmp_grid, workspace % tmp_sph_l)
2859 call tree_m2p(params, constants, params % lmax, one, &
2860 & workspace % tmp_sph, one, workspace % tmp_grid)
2861 !! Compute gradients of L2L if pl > 0
2862 if (params % pl .gt. 0) then
2863 call tree_grad_l2l(params, constants, workspace % tmp_node_l, &
2864 & workspace % tmp_sph_l_grad, workspace % tmp_sph_l)
2865 end if
2866 !! Diagonal update of computed grid values, that is needed for R^C, a part
2867 !! of R^A and a part of R^B
2868 call dgemm('T', 'N', params % ngrid, params % nsph, constants % nbasis, &
2869 & pt5, constants % vgrid2, constants % vgrid_nbasis, g, &
2870 & constants % nbasis, -one, workspace % tmp_grid, params % ngrid)
2871 !! Scale temporary grid points by corresponding Lebedev weights and ygrid
2872 do igrid = 1, params % ngrid
2873 do isph = 1, params % nsph
2874 workspace % tmp_grid(igrid, isph) = &
2875 & workspace % tmp_grid(igrid, isph) * constants % wgrid(igrid) * &
2876 & ygrid(igrid, isph)
2877 end do
2878 end do
2879 !! Compute all terms of grad_i(R). The second term of R_i^B is already
2880 !! taken into account and the first term is computed together with R_i^C.
2881 do isph = 1, params % nsph
2882 do igrid = 1, params % ngrid
2883 ! Loop over all neighbouring spheres
2884 do ik = constants % inl(isph), constants % inl(isph+1) - 1
2885 ksph = constants % nl(ik)
2886 ! Only consider external grid points
2887 if(constants % ui(igrid, ksph) .eq. zero) cycle
2888 ! build geometrical quantities
2889 c = params % csph(:, ksph) + &
2890 & params % rsph(ksph)*constants % cgrid(:, igrid)
2891 vki = c - params % csph(:, isph)
2892 !vvki = sqrt(vki(1)*vki(1) + vki(2)*vki(2) + &
2893 ! & vki(3)*vki(3))
2894 vvki = dnrm2(3, vki, 1)
2895 tki = vvki / params % rsph(isph)
2896 ! Only consider such points where grad U is non-zero
2897 if((tki.le.tlow) .or. (tki.ge.thigh)) cycle
2898 ! This is entire R^C and the first R^B component (grad_i of U
2899 ! of a sum of R_kj for index inequality j!=k)
2900 ! Indexes k and j are flipped compared to the paper
2901 gg = workspace % tmp_grid(igrid, ksph)
2902 ! Compute grad_i component of forces using precomputed
2903 ! potential gg
2904 !fx(:, isph) = fx(:, isph) - &
2905 ! & dfsw(tki, params % se, params % eta)/ &
2906 ! & params % rsph(isph)*constants % wgrid(igrid)*gg* &
2907 ! & ygrid(igrid, ksph)*(vki/vvki)
2908 if (params%switching.eq.0) then
2909 fx(:, isph) = fx(:, isph) - &
2910 & dfsw(tki, params % se, params % eta)/ &
2911 & params % rsph(isph)*gg*(vki/vvki)
2912 else
2913 d_k = constants%switching%d_ni(igrid, ksph)
2914 f_k = constants%switching%f_ni(igrid, ksph)
2915 chi_ki = fsw(tki,params%se,params%eta)
2916 grad_p = dfsw(tki,params%se,params%eta)
2917 duj = one/(d_k + f_k)**2 &
2918 & *(d_k + d_k*f_k/(one - chi_ki)) &
2919 & *grad_p/params%rsph(isph)
2920 fx(:, isph) = fx(:, isph) - gg*duj*(vki/vvki)
2921 end if
2922 end do
2923 ! contribution from the sphere itself
2924 if((constants % ui(igrid,isph).gt.zero) .and. &
2925 & (constants % ui(igrid,isph).lt.one)) then
2926 ! R^A component (grad_i of U of a sum of R_ij for index
2927 ! inequality j!=i)
2928 ! Indexes k and j are flipped compared to the paper
2929 gg = workspace % tmp_grid(igrid, isph)
2930 ! Compute grad_i component of forces using precomputed
2931 ! potential gg
2932 !fx(:, isph) = fx(:, isph) + constants % wgrid(igrid)*gg* &
2933 ! & ygrid(igrid, isph)*constants % zi(:, igrid, isph)
2934 fx(:, isph) = fx(:, isph) + gg*constants % zi(:, igrid, isph)
2935 end if
2936 if (constants % ui(igrid, isph) .gt. zero) then
2937 ! Another R^A component (grad_i of potential of a sum of R_ij
2938 ! for index inequality j!=i)
2939 ! Indexes k and j are flipped compared to the paper
2940 ! In case pl=0 MKL does not make gg3 zero reusing old value of
2941 ! gg3, so we have to clear it manually
2942 gg3 = zero
2943 call dgemv('T', params % pl**2, 3, one, &
2944 & workspace % tmp_sph_l_grad(1, 1, isph), &
2945 & (params % pl+1)**2, constants % vgrid2(1, igrid), 1, &
2946 & zero, gg3, 1)
2947 ! Gradient of the near-field potential is a gradient of
2948 ! multipole expansion
2949 inode = constants % snode(isph)
2950 do inear = constants % snear(inode), constants % snear(inode+1)-1
2951 jnode = constants % near(inear)
2952 do jsph_node = constants % cluster(1, jnode), &
2953 & constants % cluster(2, jnode)
2954 jsph = constants % order(jsph_node)
2955 if (isph .eq. jsph) cycle
2956 c = params % csph(:, isph) + &
2957 & params % rsph(isph)*constants % cgrid(:, igrid)
2958 tmp_c = c - params % csph(:, jsph)
2959 call fmm_m2p_work(tmp_c, &
2960 & params % rsph(jsph), params % lmax+1, &
2961 & constants % vscales_rel, one, &
2962 & workspace % tmp_sph_grad(:, 1, jsph), zero, &
2963 & tmp_gg, work)
2964 gg3(1) = gg3(1) + tmp_gg
2965 call fmm_m2p_work(tmp_c, &
2966 & params % rsph(jsph), params % lmax+1, &
2967 & constants % vscales_rel, one, &
2968 & workspace % tmp_sph_grad(:, 2, jsph), zero, &
2969 & tmp_gg, work)
2970 gg3(2) = gg3(2) + tmp_gg
2971 call fmm_m2p_work(tmp_c, &
2972 & params % rsph(jsph), params % lmax+1, &
2973 & constants % vscales_rel, one, &
2974 & workspace % tmp_sph_grad(:, 3, jsph), zero, &
2975 & tmp_gg, work)
2976 gg3(3) = gg3(3) + tmp_gg
2977 end do
2978 end do
2979 ! Accumulate all computed forces
2980 fx(:, isph) = fx(:, isph) - constants % wgrid(igrid)*gg3* &
2981 & ygrid(igrid, isph)*constants % ui(igrid, isph)
2982 end if
2983 end do
2984 end do
2985end subroutine gradr_fmm
2986
2988subroutine gradr(params, constants, workspace, g, ygrid, fx)
2989 implicit none
2990 ! Inputs
2991 type(ddx_params_type), intent(in) :: params
2992 type(ddx_constants_type), intent(in) :: constants
2993 real(dp), intent(in) :: g(constants % nbasis, params % nsph), &
2994 & ygrid(params % ngrid, params % nsph)
2995 ! Temporaries
2996 type(ddx_workspace_type), intent(inout) :: workspace
2997 ! Output
2998 real(dp), intent(out) :: fx(3, params % nsph)
2999 ! Check which gradr to execute
3000 if (params % fmm .eq. 1) then
3001 call gradr_fmm(params, constants, workspace, g, ygrid, fx)
3002 else
3003 call gradr_dense(params, constants, workspace, g, ygrid, fx)
3004 end if
3005end subroutine gradr
3006
3008subroutine gradr_dense(params, constants, workspace, g, ygrid, fx)
3009 implicit none
3010 ! Inputs
3011 type(ddx_params_type), intent(in) :: params
3012 type(ddx_constants_type), intent(in) :: constants
3013 real(dp), intent(in) :: g(constants % nbasis, params % nsph), &
3014 & ygrid(params % ngrid, params % nsph)
3015 ! Temporaries
3016 type(ddx_workspace_type), intent(inout) :: workspace
3017 ! Output
3018 real(dp), intent(out) :: fx(3, params % nsph)
3019 ! Local variables
3020 integer :: isph
3021 ! Simply cycle over all spheres
3022 do isph = 1, params % nsph
3023 call gradr_sph(params, constants, isph, workspace % tmp_vplm, &
3024 & workspace % tmp_vcos, workspace % tmp_vsin, &
3025 & workspace % tmp_vylm, workspace % tmp_vdylm, &
3026 & g, ygrid, fx(:, isph))
3027 end do
3028end subroutine gradr_dense
3029
3033subroutine zeta_grad(params, constants, state, e_cav, forces)
3034 implicit none
3035 type(ddx_params_type), intent(in) :: params
3036 type(ddx_constants_type), intent(in) :: constants
3037 type(ddx_state_type), intent(inout) :: state
3038 real(dp), intent(inout) :: forces(3, params % nsph)
3039 real(dp), intent(in) :: e_cav(3, constants % ncav)
3040 ! local variables
3041 integer :: icav, isph, igrid
3042
3043 icav = 0
3044 do isph = 1, params % nsph
3045 do igrid = 1, params % ngrid
3046 if (constants % ui(igrid, isph) .eq. zero) cycle
3047 icav = icav + 1
3048 forces(:, isph) = forces(:, isph) + pt5 &
3049 & *state % zeta(icav)*e_cav(:, icav)
3050 end do
3051 end do
3052end subroutine zeta_grad
3053
3057subroutine zeta_grad_dr(params, constants, state, e_cav, dr)
3058 implicit none
3059 type(ddx_params_type), intent(in) :: params
3060 type(ddx_constants_type), intent(in) :: constants
3061 type(ddx_state_type), intent(inout) :: state
3062 real(dp), intent(inout) :: dr(params % nsph)
3063 real(dp), intent(in) :: e_cav(3, constants % ncav)
3064 ! local variables
3065 integer :: icav, isph, igrid
3066
3067 icav = 0
3068 do isph = 1, params % nsph
3069 do igrid = 1, params % ngrid
3070 if (constants % ui(igrid, isph) .eq. zero) cycle
3071 icav = icav + 1
3072 dr(isph) = dr(isph) + pt5 &
3073 & * dot_product(constants % cgrid(:,igrid), e_cav(:, icav))*state % zeta(icav)
3074 end do
3075 end do
3076end subroutine zeta_grad_dr
3077
3078end 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:2785
ddx_core::tree_l2p_bessel
subroutine tree_l2p_bessel(params, constants, alpha, node_l, beta, grid_v)
TODO.
Definition: ddx_core.f90:2362
ddx_core::tree_l2p_bessel_adj
subroutine tree_l2p_bessel_adj(params, constants, alpha, grid_v, beta, node_l)
TODO.
Definition: ddx_core.f90:2435
ddx_core::tree_l2l_rotation
subroutine tree_l2l_rotation(params, constants, node_l)
Transfer local coefficients over a tree.
Definition: ddx_core.f90:1937
ddx_core::tree_m2m_rotation
subroutine tree_m2m_rotation(params, constants, node_m)
Transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1727
ddx_core::tree_l2p
subroutine tree_l2p(params, constants, alpha, node_l, beta, grid_v, sph_l)
TODO.
Definition: ddx_core.f90:2323
ddx_core::tree_l2p_adj
subroutine tree_l2p_adj(params, constants, alpha, grid_v, beta, node_l, sph_l)
TODO.
Definition: ddx_core.f90:2397
ddx_core::tree_m2l_bessel_rotation
subroutine tree_m2l_bessel_rotation(params, constants, node_m, node_l)
Transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2166
ddx_core::tree_m2m_bessel_rotation_adj
subroutine tree_m2m_bessel_rotation_adj(params, constants, node_m)
Adjoint transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1893
ddx_core::tree_m2l_rotation
subroutine tree_m2l_rotation(params, constants, node_m, node_l)
Transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2117
ddx_core::tree_m2p
subroutine tree_m2p(params, constants, p, alpha, sph_m, beta, grid_v)
TODO.
Definition: ddx_core.f90:2472
ddx_core::tree_m2p_bessel
subroutine tree_m2p_bessel(params, constants, p, alpha, sph_p, sph_m, beta, grid_v)
TODO.
Definition: ddx_core.f90:2522
ddx_core::tree_l2l_bessel_rotation_adj
subroutine tree_l2l_bessel_rotation_adj(params, constants, node_l)
Adjoint transfer local coefficients over a tree.
Definition: ddx_core.f90:2087
ddx_core::tree_l2l_bessel_rotation
subroutine tree_l2l_bessel_rotation(params, constants, node_l)
Transfer local coefficients over a tree.
Definition: ddx_core.f90:1984
ddx_core::tree_m2p_bessel_nodiag_adj
subroutine tree_m2p_bessel_nodiag_adj(params, constants, p, alpha, grid_v, beta, sph_p, sph_m)
TODO.
Definition: ddx_core.f90:2675
ddx_core::tree_m2l_bessel_rotation_adj
subroutine tree_m2l_bessel_rotation_adj(params, constants, node_l, node_m)
Adjoint transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2217
ddx_core::tree_grad_m2m
subroutine tree_grad_m2m(params, constants, sph_m, sph_m_grad, work)
TODO.
Definition: ddx_core.f90:2721
ddx_core::tree_m2p_adj
subroutine tree_m2p_adj(params, constants, p, alpha, grid_v, beta, sph_m)
TODO.
Definition: ddx_core.f90:2574
ddx_core::tree_m2p_bessel_adj
subroutine tree_m2p_bessel_adj(params, constants, p, alpha, grid_v, beta, sph_p, sph_m)
TODO.
Definition: ddx_core.f90:2626
ddx_core::tree_m2l_rotation_adj
subroutine tree_m2l_rotation_adj(params, constants, node_l, node_m)
Adjoint transfer multipole local coefficients into local over a tree.
Definition: ddx_core.f90:2279
ddx_core::tree_l2l_rotation_adj
subroutine tree_l2l_rotation_adj(params, constants, node_l)
Adjoint transfer local coefficients over a tree.
Definition: ddx_core.f90:2030
ddx_core::tree_m2m_rotation_adj
subroutine tree_m2m_rotation_adj(params, constants, node_m)
Adjoint transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1849
ddx_core::tree_m2m_bessel_rotation
subroutine tree_m2m_bessel_rotation(params, constants, node_m)
Transfer multipole coefficients over a tree.
Definition: ddx_core.f90:1788
ddx_core::dbasis
subroutine dbasis(params, constants, x, basloc, dbsloc, vplm, vcos, vsin)
Compute first derivatives of spherical harmonics.
Definition: ddx_core.f90:1179
ddx_core::intmlp
real(dp) function intmlp(params, constants, t, sigma, basloc)
TODO.
Definition: ddx_core.f90:1310
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