Octopus
math.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
22module math_oct_m
23 use debug_oct_m
24 use global_oct_m
25 use, intrinsic :: iso_fortran_env
30
31 implicit none
32
33 private
34 public :: &
35 is_close, &
37 ylmr_cmplx, &
38 ylmr_real, &
39 weights, &
40 factorial, &
41 hermite, &
46 ddelta, &
47 member, &
51 even, &
52 odd, &
56 pad, &
57 pad_pow2, &
58 log2, &
60 phi1, &
61 phi2, &
63 is_prime, &
72
73
74 interface is_close
75 module procedure dis_close_scalar, zis_close_scalar
76 end interface
77
78 interface diagonal_matrix
80 end interface diagonal_matrix
81
85
89
90
91 ! ------------------------------------------------------------------------------
94 interface interpolate
97 end interface interpolate
98 ! ------------------------------------------------------------------------------
99
100 interface log2
101 module procedure dlog2, ilog2, llog2
102 end interface log2
103
104 interface pad
105 module procedure pad4, pad8, pad48, pad88
106 end interface pad
107
108contains
109
128 elemental logical function dis_close_scalar(x, y, rtol, atol)
129 real(real64), intent(in) :: x, y
130 real(real64), optional, intent(in) :: rtol
131 real(real64), optional, intent(in) :: atol
132 real(real64) :: atol_, rtol_
133
134 atol_ = optional_default(atol, 1.e-8_real64)
135 rtol_ = optional_default(rtol, 1.e-5_real64)
136
137 dis_close_scalar = abs(x - y) <= (atol_ + rtol_ * abs(y))
138 end function dis_close_scalar
139
141 elemental logical function zis_close_scalar(x, y, rtol, atol)
142 complex(real64), intent(in) :: x, y
143 real(real64), optional, intent(in) :: rtol
144 real(real64), optional, intent(in) :: atol
145 real(real64) :: atol_, rtol_
146
147 atol_ = optional_default(atol, 1.e-8_real64)
148 rtol_ = optional_default(rtol, 1.e-5_real64)
149
150 zis_close_scalar = abs(x - y) <= (atol_ + rtol_ * abs(y))
151 end function zis_close_scalar
152
153
154 ! ---------------------------------------------------------
157 pure function idiagonal_matrix(dim, diag) result(matrix)
158 integer, intent(in) :: dim
159 integer, intent(in) :: diag
160 integer :: matrix(dim, dim)
161
162 integer :: ii
163
164 matrix = 0
165 do ii = 1, dim
166 matrix(ii, ii) = diag
167 end do
168
169 end function idiagonal_matrix
170
171 ! ---------------------------------------------------------
172 recursive function hermite(n, x) result (h)
173 integer, intent(in) :: n
174 real(real64), intent(in) :: x
176 real(real64) :: h
177
178 ! no push_sub, called too frequently
180 if (n <= 0) then
181 h = m_one
182 elseif (n == 1) then
183 h = m_two*x
184 else
185 h = m_two*x*hermite(n-1,x) - m_two*(n-1)*hermite(n-2,x)
186 end if
188 end function hermite
189
190
191 ! ---------------------------------------------------------
192 recursive function factorial (n) RESULT (fac)
193 integer, intent(in) :: n
194
195 integer :: fac
196
197 ! no push_sub for recursive function
198
199 if (n <= 1) then
200 fac = 1
201 else
202 fac = n * factorial(n-1)
203 end if
204 end function factorial
205
206 ! ---------------------------------------------------------
208 subroutine ylmr_cmplx(xx, li, mi, ylm)
209 real(real64), intent(in) :: xx(3)
210 integer, intent(in) :: li, mi
211 complex(real64), intent(out) :: ylm
212
213 integer :: i
214 real(real64) :: dx, dy, dz, r, plm, cosm, sinm, cosmm1, sinmm1, cosphi, sinphi
215 real(real64), parameter :: tiny = 1.e-15_real64
216
217 ! no push_sub: this routine is called too many times
218
219 ! if l=0, no calculations are required
220 if (li == 0) then
221 ylm = cmplx(sqrt(m_one/(m_four*m_pi)), m_zero, real64)
222 return
223 end if
224
225 r = sqrt(xx(1)**2 + xx(2)**2 + xx(3)**2)
226
227 ! if r=0, direction is undefined => make ylm=0 except for l=0
228 if (r <= tiny) then
229 ylm = m_z0
230 return
231 end if
232 dx = xx(1)/r
233 dy = xx(2)/r
234 dz = xx(3)/r
235
236 ! get the associated Legendre polynomial (including the normalization factor)
237 ! We need to avoid getting outside [-1:1] due to rounding errors
238 if(abs(dz) > m_one) dz = sign(m_one, dz)
239 plm = loct_legendre_sphplm(li, abs(mi), dz)
240
241 ! compute sin(|m|*phi) and cos(|m|*phi)
242 r = sqrt(dx*dx + dy*dy)
243 if(r > tiny) then
244 cosphi = dx/r
245 sinphi = dy/r
246 else ! In this case the values are ill defined so we choose \phi=0
247 cosphi = m_zero
248 sinphi = m_zero
249 end if
251 cosm = m_one
252 sinm = m_zero
253 do i = 1, abs(mi)
254 cosmm1 = cosm
255 sinmm1 = sinm
256 cosm = cosmm1*cosphi - sinmm1*sinphi
257 sinm = cosmm1*sinphi + sinmm1*cosphi
258 end do
259
260 !And now ylm
261 ylm = plm*cmplx(cosm, sinm, real64)
262
263 if (mi < 0) then
264 ylm = conjg(ylm)
265 do i = 1, abs(mi)
266 ylm = -ylm
267 end do
268 end if
269
270 end subroutine ylmr_cmplx
271
272
273 ! ---------------------------------------------------------
279 subroutine ylmr_real(xx, li, mi, ylm)
280 real(real64), intent(in) :: xx(3)
281 integer, intent(in) :: li, mi
282 real(real64), intent(out) :: ylm
283
284 integer, parameter :: lmaxd = 20
285 real(real64), parameter :: tiny = 1.e-15_real64
286 integer :: i, ilm0, l, m, mabs
287 integer, save :: lmax = -1
288
289 real(real64) :: cmi, cosm, cosmm1, cosphi, dphase, dplg, fac, &
290 fourpi, plgndr, phase, pll, pmm, pmmp1, sinm, &
291 sinmm1, sinphi, rsize, rx, ry, rz, xysize
292 real(real64), save :: c(0:(lmaxd+1)*(lmaxd+1))
293
294 ! no push_sub: called too frequently
295
296 ! evaluate normalization constants once and for all
297 if (li > lmax) then
298 fourpi = m_four*m_pi
299 do l = 0, li
300 ilm0 = l*l + l
301 do m = 0, l
302 fac = (2*l+1)/fourpi
303 do i = l - m + 1, l + m
304 fac = fac/i
305 end do
306 c(ilm0 + m) = sqrt(fac)
307 ! next line because ylm are real combinations of m and -m
308 if (m /= 0) c(ilm0 + m) = c(ilm0 + m)*sqrt(m_two)
309 c(ilm0 - m) = c(ilm0 + m)
310 end do
311 end do
312 lmax = li
313 end if
314
315 ! if l=0, no calculations are required
316 if (li == 0) then
317 ylm = c(0)
318 return
319 end if
320
321 ! if r=0, direction is undefined => make ylm=0 except for l=0
322 rsize = sqrt(xx(1)**2 + xx(2)**2 + xx(3)**2)
323 if (rsize < tiny) then
324 ylm = m_zero
325 return
326 end if
327
328 rx = xx(1)/rsize
329 ry = xx(2)/rsize
330 rz = xx(3)/rsize
331
332 ! explicit formulas for l=1 and l=2
333 if (li == 1) then
334 select case (mi)
335 case (-1)
336 ylm = (-c(1))*ry
337 case (0)
338 ylm = c(2)*rz
339 case (1)
340 ylm = (-c(3))*rx
341 end select
342 return
343 end if
344
345 if (li == 2) then
346 select case (mi)
347 case (-2)
348 ylm = c(4)*6.0_real64*rx*ry
349 case (-1)
350 ylm = (-c(5))*m_three*ry*rz
351 case (0)
352 ylm = c(6)*m_half*(m_three*rz*rz - m_one)
353 case (1)
354 ylm = (-c(7))*m_three*rx*rz
355 case (2)
356 ylm = c(8)*m_three*(rx*rx - ry*ry)
357 end select
358 return
359 end if
360
361 ! general algorithm based on routine plgndr of numerical recipes
362 mabs = abs(mi)
363 xysize = sqrt(max(rx*rx + ry*ry, tiny))
364 cosphi = rx/xysize
365 sinphi = ry/xysize
366 cosm = m_one
367 sinm = m_zero
368 do m = 1, mabs
369 cosmm1 = cosm
370 sinmm1 = sinm
371 cosm = cosmm1*cosphi - sinmm1*sinphi
372 sinm = cosmm1*sinphi + sinmm1*cosphi
373 end do
374
375 if (mi < 0) then
376 phase = sinm
377 dphase = mabs*cosm
378 else
379 phase = cosm
380 dphase = (-mabs)*sinm
381 end if
382
383 pmm = m_one
384 fac = m_one
385
386 if (mabs > m_zero) then
387 do i = 1, mabs
388 pmm = (-pmm)*fac*xysize
389 fac = fac + m_two
390 end do
391 end if
392
393 if (li == mabs) then
394 plgndr = pmm
395 dplg = (-li)*rz*pmm/(xysize**2)
396 else
397 pmmp1 = rz*(2*mabs + 1)*pmm
398 if (li == mabs + 1) then
399 plgndr = pmmp1
400 dplg = -((li*rz*pmmp1 - (mabs + li)*pmm)/(xysize**2))
401 else
402 do l = mabs + 2, li
403 pll = (rz*(2*l - 1)*pmmp1 - (l + mabs - 1)*pmm)/(l - mabs)
404 pmm = pmmp1
405 pmmp1 = pll
406 end do
407 plgndr = pll
408 dplg = -((li*rz*pll - (l + mabs - 1)*pmm)/(xysize**2))
409 end if
410 end if
411
412 ilm0 = li*li + li
413 cmi = c(ilm0 + mi)
414 ylm = cmi*plgndr*phase
415
416 end subroutine ylmr_real
417
418
419 ! ---------------------------------------------------------
429 subroutine weights(N, M, cc, side)
430 integer, intent(in) :: n, m
431 real(real64), intent(out) :: cc(0:,0:,0:)
432 integer, optional, intent(in) :: side
433
434 integer :: i, j, k, mn, side_
435 real(real64) :: c1, c2, c3, c4, c5, xi
436 real(real64), allocatable :: x(:)
437
438 push_sub(weights)
439
440 safe_allocate(x(0:m))
441
442 if (present(side)) then
443 side_ = side
444 else
445 side_ = 0
446 end if
447
448 select case (side_)
449 case (-1)
450 ! grid-points for left-side finite-difference formulas on an equi-spaced grid
451 mn = m
452 x(:) = (/(-i,i=0,mn)/)
453 case (+1)
454 ! grid-points for right-side finite-difference formulas on an equi-spaced grid
455 mn = m
456 x(:) = (/(i,i=0,mn)/)
457 case default
458 ! grid-points for centered finite-difference formulas on an equi-spaced grid
459 mn = m/2
460 x(:) = (/0,(-i,i,i=1,mn)/)
461 end select
462
463
464
465 xi = m_zero ! point at which the approx. are to be accurate
466
467 cc = m_zero
468 cc(0,0,0) = m_one
469
470 c1 = m_one
471 c4 = x(0) - xi
472
473 do j = 1, m
474 mn = min(j,n)
475 c2 = m_one
476 c5 = c4
477 c4 = x(j) - xi
478
479 do k = 0, j - 1
480 c3 = x(j) - x(k)
481 c2 = c2*c3
482
483 if (j <= n) cc(k, j - 1, j) = m_zero
484 cc(k, j, 0) = c4*cc(k, j - 1, 0)/c3
485
486 do i = 1, mn
487 cc(k, j, i) = (c4*cc(k, j - 1, i) - i*cc(k, j - 1, i - 1))/c3
488 end do
489
490 cc(j, j, 0) = -c1*c5*cc(j - 1, j - 1, 0) / c2
491 end do
492
493 do i = 1, mn
494 cc(j, j, i) = c1*(i*cc(j - 1, j - 1, i - 1) - c5*cc(j - 1, j - 1, i))/c2
495 end do
496
497 c1 = c2
498 end do
499
500 safe_deallocate_a(x)
501 pop_sub(weights)
502 end subroutine weights
503
504 ! ---------------------------------------------------------
505 real(real64) pure function ddelta(i, j)
506 integer, intent(in) :: i
507 integer, intent(in) :: j
508
509 ! no push_sub in pure function
510
511 if (i == j) then
512 ddelta = m_one
513 else
514 ddelta = m_zero
515 end if
516
517 end function ddelta
518
519
520 ! ---------------------------------------------------------
523 subroutine make_idx_set(n, out, length, in)
524 integer, intent(in) :: n
525 integer, allocatable, intent(out) :: out(:)
526 integer, intent(out) :: length
527 integer, optional, intent(in) :: in(:)
528
529 integer :: i
530
531 push_sub(make_idx_set)
532
533 if (present(in)) then
534 length = ubound(in, 1)
535 safe_allocate(out(1:length))
536 out = in
537 else
538 length = n
539 safe_allocate(out(1:length))
540 do i = 1, length
541 out(i) = i
542 end do
543 end if
544
545 pop_sub(make_idx_set)
546 end subroutine make_idx_set
547
548
549 ! ---------------------------------------------------------
552 logical function member(n, a)
553 integer, intent(in) :: n
554 integer, intent(in) :: a(:)
555
556 integer :: i
557
558 ! no push_sub, called too frequently
559
560 member = .false.
561
562 do i = 1, ubound(a, 1)
563 if (a(i) == n) then
564 member = .true.
565 exit
566 end if
567 end do
568
569 end function member
570
571
572 ! ---------------------------------------------------------
573 subroutine interpolation_coefficients(nn, xa, xx, cc)
574 integer, intent(in) :: nn
575 real(real64), intent(in) :: xa(:)
576 real(real64), intent(in) :: xx
577 real(real64), intent(out) :: cc(:)
578
579 integer :: ii, kk
580
581 ! no push_sub, called too frequently
582
583 do ii = 1, nn
584 cc(ii) = m_one
585 do kk = 1, nn
586 if (kk == ii) cycle
587 cc(ii) = cc(ii)*(xx - xa(kk))/(xa(ii) - xa(kk))
588 end do
589 end do
590
591 end subroutine interpolation_coefficients
592
593
594 ! ---------------------------------------------------------
596 logical pure function even(n)
597 integer, intent(in) :: n
599 even = (mod(n, 2) == 0)
600
601 end function even
602
603
604 ! ---------------------------------------------------------
606 logical pure function odd(n)
607 integer, intent(in) :: n
608
609 odd = .not. even(n)
610
611 end function odd
612
613 ! ---------------------------------------------------------
617 real(real64), intent(in) :: x(:)
618 real(real64), intent(out) :: u(:)
619
620 integer :: n, k, j
621 real(real64) :: sumx2
622 real(real64), allocatable :: xx(:)
623
625
626 n = size(x)
627 assert(n>1)
628 assert(size(u) == n-1)
629
630 ! These lines make the code less machine-dependent.
631 safe_allocate(xx(1:n))
632 do j = 1, n
633 if (abs(x(j))<1.0e-8_real64) then
634 xx(j) = m_zero
635 else
636 xx(j) = x(j)
637 end if
638 end do
639
640 u = m_zero
641 do k = 1, n-1
642 sumx2 = m_zero
643 do j = k+1, n
644 sumx2 = sumx2 + xx(j)**2
645 end do
646 if (abs(sumx2) <= m_epsilon .and. abs(xx(k)) <= m_epsilon) exit
647 if (k < n - 1) then
648 u(k) = atan2(sqrt(sumx2), xx(k))
649 else
650 u(n-1) = atan2(xx(n), xx(n-1))
651 end if
652 end do
653
654
656 end subroutine cartesian2hyperspherical
657 ! ---------------------------------------------------------
658
659
660 ! ---------------------------------------------------------
663 subroutine hyperspherical2cartesian(u, x)
664 real(real64), intent(in) :: u(:)
665 real(real64), intent(out) :: x(:)
667 integer :: n, j, k
668
670
671 n = size(x)
672 assert(n>1)
673 assert(size(u) == n-1)
674
675 if (n == 2) then
676 x(1) = cos(u(1))
677 x(2) = sin(u(1))
678 elseif (n == 3) then
679 x(1) = cos(u(1))
680 x(2) = sin(u(1))*cos(u(2))
681 x(3) = sin(u(1))*sin(u(2))
682 else
683 x(1) = cos(u(1))
684 x(2) = sin(u(1))*cos(u(2))
685 x(3) = sin(u(1))*sin(u(2))*cos(u(3))
686 do j = 4, n - 1
687 x(j) = m_one
688 do k = 1, j - 1
689 x(j) = x(j) * sin(u(k))
690 end do
691 x(j) = x(j) * cos(u(j))
692 end do
693 x(n) = m_one
694 do k = 1, n - 2
695 x(n) = x(n) * sin(u(k))
696 end do
697 x(n) = x(n) * sin(u(n-1))
698 end if
701 end subroutine hyperspherical2cartesian
702 ! ---------------------------------------------------------
703
704
705 ! ---------------------------------------------------------
709 subroutine hypersphere_grad_matrix(grad_matrix, r, x)
710 real(real64), intent(out) :: grad_matrix(:,:)
711 real(real64), intent(in) :: r
712 real(real64), intent(in) :: x(:)
713
714 integer :: n, l, m
715
717
718 n = size(x)+1 ! the dimension of the matrix is (n-1)x(n)
719
720 ! --- l=1 ---
721 grad_matrix = m_one
722 grad_matrix(1,1) = -r*sin(x(1))
723 do m = 2, n-1
724 grad_matrix(m,1) = m_zero
725 end do
726
727 ! --- l=2..(n-1) ---
728 do l=2, n-1
729 do m=1, n-1
730 if (m == l) then
731 grad_matrix(m,l) = -r*grad_matrix(m,l)*product(sin(x(1:l)))
732 elseif (m < l) then
733 grad_matrix(m,l) = grad_matrix(m,l)*r*cos(x(m))*cos(x(l))*product(sin(x(1:m-1)))*product(sin(x(m+1:l-1)))
734 else
735 grad_matrix(m,l) = m_zero
736 end if
737 end do
738 end do
739
740 ! --- l=n ---
741 do m=1, n-1
742 grad_matrix(m,n) = r*cos(x(m))*grad_matrix(m,n)*product(sin(x(1:m-1)))*product(sin(x(m+1:n-1)))
743 end do
744
746 end subroutine hypersphere_grad_matrix
747
748 ! ---------------------------------------------------------
749 integer(int64) pure function pad88(size, blk)
750 integer(int64), intent(in) :: size
751 integer(int64), intent(in) :: blk
752
753 integer(int64) :: mm
754
755 mm = mod(size, blk)
756 if (mm == 0) then
757 pad88 = size
758 else
759 pad88 = size + blk - mm
760 end if
761 end function pad88
762
763 ! ---------------------------------------------------------
764 integer(int64) pure function pad48(size, blk)
765 integer, intent(in) :: size
766 integer(int64), intent(in) :: blk
767
768 pad48 = pad88(int(size, int64), blk)
769 end function pad48
770
771 ! ---------------------------------------------------------
772 integer(int64) pure function pad8(size, blk)
773 integer(int64), intent(in) :: size
774 integer, intent(in) :: blk
775
776 pad8 = pad88(size, int(blk, int64))
777 end function pad8
778
779 ! ---------------------------------------------------------
780 integer pure function pad4(size, blk)
781 integer, intent(in) :: size
782 integer, intent(in) :: blk
783
784 pad4 = int(pad88(int(size, int64), int(blk, int64)), int32)
785 end function pad4
786
787 ! ---------------------------------------------------------
788
793 integer pure function pad_pow2(size)
794 integer, intent(in) :: size
795
796 integer :: mm, mm2
797
798 mm = size
799 pad_pow2 = 1
800
801 ! loop below never terminates otherwise! just pick 1 as value.
802 if (size == 0) return
803
804 ! first we divide by two all the times we can, so we catch exactly
805 ! the case when size is already a power of two
806 do
807 mm2 = mm/2
808 if (mm - 2*mm2 /= 0) exit
810 mm = mm2
811 end do
812
813 ! the rest is handled by log2
814 if (mm /= 1) then
815 pad_pow2 = pad_pow2*2**(ceiling(log(real(mm, real64) )/log(2.0_8)))
816 end if
817
818 end function pad_pow2
819
820 ! -------------------------------------------------------
821
822 real(real64) pure function dlog2(xx)
823 real(real64), intent(in) :: xx
824
825 dlog2 = log(xx)/log(m_two)
826 end function dlog2
827
828 ! -------------------------------------------------------
829
830 integer pure function ilog2(xx)
831 integer, intent(in) :: xx
832
833 ilog2 = nint(log2(real(xx, real64) ))
834 end function ilog2
835
836 ! -------------------------------------------------------
837
838 integer(int64) pure function llog2(xx)
839 integer(int64), intent(in) :: xx
840
841 llog2 = nint(log2(real(xx, real64) ), kind=int64)
842 end function llog2
843
844 ! -------------------------------------------------------
846 complex(real64) pure function exponential(z)
847 complex(real64), intent(in) :: z
848
849 exponential = exp(z)
850 end function exponential
851
852 ! -------------------------------------------------------
858 complex(real64) pure function phi1(z)
859 complex(real64), intent(in) :: z
860
861 real(real64), parameter :: cut = 0.002_real64
862
863 if (abs(z) < cut) then
864 phi1 = m_one + m_half*z*(m_one + z/m_three*(m_one + m_fourth*z*(m_one + z/m_five)))
865 else
866 phi1 = (exp(z) - m_z1)/z
867 end if
868 end function phi1
869
870 ! -------------------------------------------------------
877 complex(real64) pure function phi2(z)
878 complex(real64), intent(in) :: z
879
880 real(real64), parameter :: cut = 0.002_real64
881
882 if (abs(z) < cut) then
883 phi2 = m_half + z/6.0_real64*(m_one + m_fourth*z*(m_one + z/m_five*(m_one + z/6.0_real64)))
884 else
885 phi2 = (exp(z) - z - m_z1)/z**2
886 end if
887 end function phi2
888
889 ! -------------------------------------------------------
891 real(real64) pure function square_root(x)
892 real(real64), intent(in) :: x
893
894 square_root = sqrt(x)
895 end function square_root
896
897 ! -------------------------------------------------------
898
899 logical function is_prime(n)
900 integer, intent(in) :: n
901
902 integer :: i, root
903
904 push_sub(is_prime)
905
906 if (n < 1) then
907 message(1) = "Internal error: is_prime does not take negative numbers."
908 call messages_fatal(1)
909 end if
910 if (n == 1) then
911 is_prime = .false.
912 pop_sub(is_prime)
913 return
914 end if
916 root = nint(sqrt(real(n, real64) ))
917 do i = 2, root
918 if (mod(n,i) == 0) then
919 is_prime = .false.
920 pop_sub(is_prime)
921 return
922 end if
923 end do
924
925 is_prime = .true.
926 pop_sub(is_prime)
927 end function is_prime
928
929 ! ---------------------------------------------------------
934 subroutine generate_rotation_matrix(R, ff, tt)
935 real(real64), intent(out) :: r(:,:)
936 real(real64), intent(in) :: ff(:)
937 real(real64), intent(in) :: tt(:)
938
939 integer :: dim, i, j
940 real(real64) :: th, uv, uu, vv, ft
941 real(real64), allocatable :: axis(:), u(:), v(:), f(:), t(:), p(:)
942
944
945 dim = size(ff,1)
946
947 assert((dim < 3) .or. (dim > 2))
948 assert(size(tt,1) == dim)
949 assert((size(r,1) == dim) .and. (size(r,2) == dim))
950
951 safe_allocate(u(1:dim))
952 safe_allocate(v(1:dim))
953 safe_allocate(f(1:dim))
954 safe_allocate(t(1:dim))
955
956
957 !normalize
958 f = ff / norm2(ff)
959 t = tt / norm2(tt)
960
961 ft = dot_product(f,t)
962
963 if (abs(ft) < m_one) then
964 select case (dim)
965 case (2)
966 th = acos(ft)
967 r(1,1) = cos(th)
968 r(1,2) = -sin(th)
969
970 r(2,1) = sin(th)
971 r(2,2) = cos(th)
972
973 case (3)
974 if (.false.) then
975 !Old implementation
976 safe_allocate(axis(1:dim))
977 th = acos(ft)
978
979 u = f / dot_product(f,f)
980 v = t /dot_product(t,t)
981
982 axis = dcross_product(u,v)
983 axis = axis / norm2(axis)
985 r(1,1) = cos(th) + axis(1)**2 * (1 - cos(th))
986 r(1,2) = axis(1)*axis(2)*(1-cos(th)) + axis(3)*sin(th)
987 r(1,3) = axis(1)*axis(3)*(1-cos(th)) - axis(2)*sin(th)
988
989 r(2,1) = axis(2)*axis(1)*(1-cos(th)) - axis(3)*sin(th)
990 r(2,2) = cos(th) + axis(2)**2 * (1 - cos(th))
991 r(2,3) = axis(2)*axis(3)*(1-cos(th)) + axis(1)*sin(th)
993 r(3,1) = axis(3)*axis(1)*(1-cos(th)) + axis(2)*sin(th)
994 r(3,2) = axis(3)*axis(2)*(1-cos(th)) - axis(1)*sin(th)
995 r(3,3) = cos(th) + axis(3)**2 * (1 - cos(th))
996
997 safe_deallocate_a(axis)
998 end if
999
1000 if (.true.) then
1001 !Naive implementation
1002 th = acos(ft)
1003 u = dcross_product(f,t)
1004 u = u / norm2(u)
1005
1006 r(1,1) = u(1)**2 + (1-u(1)**2)*cos(th)
1007 r(1,2) = u(1)*u(2)*(1-cos(th)) - u(3)*sin(th)
1008 r(1,3) = u(1)*u(3) + u(2)*sin(th)
1009
1010 r(2,1) = u(1)*u(2)*(1-cos(th)) + u(3)*sin(th)
1011 r(2,2) = u(2)**2 + (1-u(2)**2)*cos(th)
1012 r(2,3) = u(2)*u(3)*(1-cos(th)) - u(1)*sin(th)
1013
1014 r(3,1) = u(1)*u(3)*(1-cos(th)) - u(2)*sin(th)
1015 r(3,2) = u(2)*u(3)*(1-cos(th)) + u(1)*sin(th)
1016 r(3,3) = u(3)**2 + (1-u(3)**2)*cos(th)
1017 end if
1018
1019 if (.false.) then
1020 !Fast
1021 safe_allocate(p(1:dim))
1022
1023 if (abs(f(1)) <= abs(f(2)) .and. abs(f(1)) < abs(f(3))) then
1024 p = (/m_one, m_zero, m_zero/)
1025 else if (abs(f(2)) < abs(f(1)) .and. abs(f(2)) <= abs(f(3))) then
1026 p = (/m_zero, m_one, m_zero/)
1027 else if (abs(f(3)) <= abs(f(1)) .and. abs(f(3)) < abs(f(2))) then
1028 p = (/m_zero, m_zero, m_one/)
1029 end if
1030
1031 u = p - f
1032 v = p - t
1033
1034 uu = dot_product(u,u)
1035 vv = dot_product(v,v)
1036 uv = dot_product(u,v)
1037
1038 do i = 1,3
1039 do j = 1,3
1040
1041 r(i,j) = ddelta(i,j) - m_two * u(i)*u(j)/uu - m_two * v(i)*v(j)/vv &
1042 + 4_real64*uv * v(i)*u(j) /(uu*vv)
1043
1044 end do
1045 end do
1046
1047 safe_deallocate_a(p)
1048 end if
1049
1050
1051 end select
1052
1053 else
1054
1055 r = m_zero
1056 do i = 1,dim
1057 r(i,i) = m_one
1058 end do
1059
1060 end if
1061
1062
1063 safe_deallocate_a(u)
1064 safe_deallocate_a(v)
1065
1067 end subroutine generate_rotation_matrix
1068
1069 ! ---------------------------------------------------------
1078 subroutine numder_ridders(x, h, res, err, f)
1079 implicit none
1080 real(real64), intent(in) :: x, h
1081 real(real64), intent(out) :: res, err
1082 interface
1083 subroutine f(x, fx)
1084 use, intrinsic :: iso_fortran_env
1085 implicit none
1086 real(real64), intent(in) :: x
1087 real(real64), intent(inout) :: fx
1088 end subroutine f
1089 end interface
1090
1091 real(real64) :: con = 1.4_real64, big = 1.0e30_real64, safe = m_two
1092 integer :: ntab = 20, i, j
1093 real(real64) :: errt, fac, hh, fx1, fx2
1094 real(real64), allocatable :: a(:, :)
1095
1096 push_sub(numder_ridders)
1097
1098 if (abs(h) <= m_epsilon) then
1099 message(1) = "h must be nonzero in numder_ridders"
1100 call messages_fatal(1)
1101 end if
1102
1103 safe_allocate(a(1:ntab, 1:ntab))
1104
1105 hh = h
1106 call f(x+hh, fx1)
1107 call f(x-hh, fx2)
1108 a(1,1) = (fx1-fx2) / (m_two*hh)
1109 err = big
1110 do i = 2, ntab
1111 hh = hh / con
1112 call f(x+hh, fx1)
1113 call f(x-hh, fx2)
1114 a(1,i) = (fx1-fx2) / (m_two*hh)
1115 fac = con**2
1116 do j = 2, i
1117 a(j,i) = (a(j-1,i)*fac-a(j-1,i-1)) / (fac-m_one)
1118 fac = con**2*fac
1119 errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
1120 if (errt .le. err) then
1121 err = errt
1122 res = a(j,i)
1123 end if
1124 end do
1125 if (abs(a(i,i)-a(i-1,i-1)) .ge. safe*err) return
1126 end do
1127
1128 safe_deallocate_a(a)
1129 pop_sub(numder_ridders)
1130 end subroutine numder_ridders
1131
1132 ! ---------------------------------------------------------
1133 pure function dzcross_product(a, b) result(c)
1134 real(real64), intent(in) :: a(:)
1135 complex(real64), intent(in) :: b(:)
1136
1137 complex(real64) :: c(1:3)
1138
1139 c(1) = a(2)*b(3) - a(3)*b(2)
1140 c(2) = a(3)*b(1) - a(1)*b(3)
1141 c(3) = a(1)*b(2) - a(2)*b(1)
1142
1143 end function dzcross_product
1144
1145 ! ---------------------------------------------------------
1146 pure function zdcross_product(a, b) result(c)
1147 complex(real64), intent(in) :: a(:)
1148 real(real64), intent(in) :: b(:)
1149
1150 complex(real64) :: c(1:3)
1151
1152 c(1) = a(2)*b(3) - a(3)*b(2)
1153 c(2) = a(3)*b(1) - a(1)*b(3)
1154 c(3) = a(1)*b(2) - a(2)*b(1)
1155
1156 end function zdcross_product
1157
1158
1159!*****************************************************************************80
1160!
1161!! LM_POLYNOMIAL evaluates Laguerre polynomials Lm(n,m,x).
1162!
1163! First terms:
1164!
1165! M = 0
1166!
1167! Lm(0,0,X) = 1
1168! Lm(1,0,X) = -X + 1
1169! Lm(2,0,X) = X^2 - 4 X + 2
1170! Lm(3,0,X) = -X^3 + 9 X^2 - 18 X + 6
1171! Lm(4,0,X) = X^4 - 16 X^3 + 72 X^2 - 96 X + 24
1172! Lm(5,0,X) = -X^5 + 25 X^4 - 200 X^3 + 600 X^2 - 600 x + 120
1173! Lm(6,0,X) = X^6 - 36 X^5 + 450 X^4 - 2400 X^3 + 5400 X^2 - 4320 X + 720
1174!
1175! M = 1
1176!
1177! Lm(0,1,X) = 0
1178! Lm(1,1,X) = -1,
1179! Lm(2,1,X) = 2 X - 4,
1180! Lm(3,1,X) = -3 X^2 + 18 X - 18,
1181! Lm(4,1,X) = 4 X^3 - 48 X^2 + 144 X - 96
1182!
1183! M = 2
1184!
1185! Lm(0,2,X) = 0
1186! Lm(1,2,X) = 0,
1187! Lm(2,2,X) = 2,
1188! Lm(3,2,X) = -6 X + 18,
1189! Lm(4,2,X) = 12 X^2 - 96 X + 144
1190!
1191! M = 3
1192!
1193! Lm(0,3,X) = 0
1194! Lm(1,3,X) = 0,
1195! Lm(2,3,X) = 0,
1196! Lm(3,3,X) = -6,
1197! Lm(4,3,X) = 24 X - 96
1198!
1199! M = 4
1200!
1201! Lm(0,4,X) = 0
1202! Lm(1,4,X) = 0
1203! Lm(2,4,X) = 0
1204! Lm(3,4,X) = 0
1205! Lm(4,4,X) = 24
1206!
1207! Recursion:
1208!
1209! Lm(0,M,X) = 1
1210! Lm(1,M,X) = (M+1-X)
1211!
1212! if 2 <= N:
1213!
1214! Lm(N,M,X) = ( (M+2*N-1-X) * Lm(N-1,M,X)
1215! + (1-M-N) * Lm(N-2,M,X) ) / N
1216!
1217! Special values:
1218!
1219! For M = 0, the associated Laguerre polynomials Lm(N,M,X) are equal
1220! to the Laguerre polynomials L(N,X).
1221!
1222! Licensing:
1223!
1224! This code is distributed under the GNU LGPL license.
1225!
1226! Modified:
1227!
1228! 08 February 2003
1229!
1230! Author:
1231!
1232! John Burkardt
1233!
1234! Reference:
1235!
1236! Milton Abramowitz, Irene Stegun,
1237! Handbook of Mathematical Functions,
1238! National Bureau of Standards, 1964,
1239! ISBN: 0-486-61272-4,
1240! LC: QA47.A34.
1241!
1242! Parameters:
1243!
1244! Input, integer ( kind = 4 ) MM, the number of evaluation points.
1245!
1246! Input, integer ( kind = 4 ) N, the highest order polynomial to compute.
1247! Note that polynomials 0 through N will be computed.
1248!
1249! Input, integer ( kind = 4 ) M, the parameter. M must be nonnegative.
1250!
1251! Input, real ( kind = rk ) X(MM), the evaluation points.
1252!
1253! Output, real ( kind = rk ) CX(MM,0:N), the associated Laguerre polynomials
1254! of degrees 0 through N evaluated at the evaluation points.
1255!
1256! Taken from
1257! https://people.sc.fsu.edu/~jburkardt/f_src/laguerre_polynomial/laguerre_polynomial.html
1258! and adapted to Octopus by N. Tancogne-Dejean
1259 subroutine generalized_laguerre_polynomial ( np, nn, mm, xx, cx )
1260 integer, intent(in) :: np
1261 integer, intent(in) :: nn, mm
1262 real(real64), intent(in) :: xx(np)
1263 real(real64), intent(out) :: cx(np)
1264
1265 integer ii
1266 real(real64), allocatable :: cx_tmp(:,:)
1267
1268 assert(mm >= 0)
1269
1270 safe_allocate(cx_tmp(1:np, 0:nn))
1271
1272 if (nn < 0) then
1273 cx = m_zero
1274 return
1275 end if
1276
1277 if (nn == 0) then
1278 cx(1:np) = m_one
1279 return
1280 end if
1281
1282 cx_tmp(1:np,0) = m_one
1283 cx_tmp(1:np,1) = real(mm + 1, real64) - xx(1:np)
1284
1285 do ii = 2, nn
1286 cx_tmp(1:np, ii) = &
1287 (( real(mm + 2*ii - 1, real64) - xx(1:np)) * cx_tmp(1:np,ii-1) &
1288 + real(-mm -ii + 1, real64) * cx_tmp(1:np,ii-2)) / real(ii, real64)
1289 end do
1290
1291 cx(1:np) = cx_tmp(1:np, nn)
1292
1293 safe_deallocate_a(cx_tmp)
1294
1295 end subroutine generalized_laguerre_polynomial
1296
1297 subroutine dupper_triangular_to_hermitian(nn, aa)
1298 integer, intent(in) :: nn
1299 real(real64), intent(inout) :: aa(:, :)
1300
1301 integer :: ii, jj
1302
1303 do jj = 1, nn
1304 do ii = 1, jj-1
1305 aa(jj, ii) = aa(ii, jj)
1306 end do
1307 end do
1308 end subroutine dupper_triangular_to_hermitian
1309
1310 subroutine zupper_triangular_to_hermitian(nn, aa)
1311 integer, intent(in) :: nn
1312 complex(real64), intent(inout) :: aa(:, :)
1313
1314 integer :: ii, jj
1315
1316 do jj = 1, nn
1317 do ii = 1, jj-1
1318 aa(jj, ii) = conjg(aa(ii, jj))
1319 end do
1320 aa(jj, jj) = real(aa(jj, jj), real64)
1321 end do
1322 end subroutine zupper_triangular_to_hermitian
1323
1324 subroutine dlower_triangular_to_hermitian(nn, aa)
1325 integer, intent(in) :: nn
1326 real(real64), intent(inout) :: aa(:, :)
1327
1328 integer :: ii, jj
1329
1330 do jj = 1, nn
1331 do ii = 1, jj-1
1332 aa(ii, jj) = aa(jj, ii)
1333 end do
1334 end do
1335 end subroutine dlower_triangular_to_hermitian
1336
1337 subroutine zlower_triangular_to_hermitian(nn, aa)
1338 integer, intent(in) :: nn
1339 complex(real64), intent(inout) :: aa(:, :)
1340
1341 integer :: ii, jj
1342
1343 do jj = 1, nn
1344 do ii = 1, jj-1
1345 aa(ii, jj) = conjg(aa(jj, ii))
1346 end do
1347 aa(jj, jj) = real(aa(jj, jj), real64)
1348 end do
1349 end subroutine zlower_triangular_to_hermitian
1350
1351
1352 subroutine dsymmetrize_matrix(nn, aa)
1353 integer, intent(in) :: nn
1354 real(real64), intent(inout) :: aa(:, :)
1355
1356 integer :: ii, jj
1357
1358 do ii = 1, nn
1359 do jj = ii+1, nn
1360 aa(ii, jj) = m_half * (aa(ii, jj) + aa(jj, ii))
1361 aa(jj, ii) = aa(ii, jj)
1362 end do
1363 end do
1364 end subroutine dsymmetrize_matrix
1365
1366 subroutine dzero_small_elements_matrix(nn, aa, tol)
1367 integer, intent(in) :: nn
1368 real(real64), intent(inout) :: aa(:, :)
1369 real(real64) :: tol
1370
1371 where(abs(aa)<tol)
1372 aa = m_zero
1373 end where
1374 end subroutine dzero_small_elements_matrix
1375
1376 ! ---------------------------------------------------------
1378 logical pure function is_diagonal(dim, matrix, tol)
1379 integer, intent(in) :: dim
1380 real(real64), intent(in) :: matrix(:, :)
1381 real(real64), intent(in) :: tol
1382
1383 integer :: ii, jj
1384 real(real64) :: max_diag
1385
1386 max_diag = m_zero
1387 do ii = 1, dim
1388 max_diag = max(abs(matrix(ii, ii)), max_diag)
1389 end do
1391 is_diagonal = .true.
1392 do ii = 1, dim
1393 do jj = 1, dim
1394 if (ii == jj) cycle
1395 if (abs(matrix(ii, jj)) > tol*max_diag) is_diagonal = .false.
1396 end do
1397 end do
1398
1399 end function is_diagonal
1400
1401#include "undef.F90"
1402#include "complex.F90"
1403#include "math_inc.F90"
1404
1405#include "undef.F90"
1406#include "real.F90"
1407#include "math_inc.F90"
1408
1409end module math_oct_m
1410
1411!! Local Variables:
1412!! mode: f90
1413!! coding: utf-8
1414!! End:
This is the common interface to a simple-minded polynomical interpolation procedure (simple use of th...
Definition: math.F90:187
double acos(double __x) __attribute__((__nothrow__
double log(double __x) __attribute__((__nothrow__
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
double sqrt(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
double atan2(double __y, double __x) __attribute__((__nothrow__
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_four
Definition: global.F90:192
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
complex(real64), parameter, public m_z0
Definition: global.F90:198
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_three
Definition: global.F90:191
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public dsymmetrize_matrix(nn, aa)
Definition: math.F90:1446
complex(real64) pure function, public phi2(z)
Compute phi2(z) = (phi1(z)-1)/z = (exp(z) - z - 1)/z^2.
Definition: math.F90:971
integer pure function ilog2(xx)
Definition: math.F90:924
logical pure function, public even(n)
Returns if n is even.
Definition: math.F90:690
pure real(real64) function, dimension(dim, dim) ddiagonal_matrix(dim, diag)
Currently only returns a matrix whose diagonal elements are all the same. Note that the integer versi...
Definition: math.F90:1756
subroutine, public generalized_laguerre_polynomial(np, nn, mm, xx, cx)
Definition: math.F90:1353
complex(real64) pure function, public exponential(z)
Wrapper for exponential.
Definition: math.F90:940
subroutine dinterpolate_0(xa, ya, x, y)
Definition: math.F90:1834
subroutine, public ylmr_cmplx(xx, li, mi, ylm)
Computes spherical harmonics ylm at position (x, y, z)
Definition: math.F90:302
subroutine zinterpolate_0(xa, ya, x, y)
Definition: math.F90:1643
subroutine dinterpolate_1(xa, ya, x, y)
Definition: math.F90:1806
real(real64) pure function dlog2(xx)
Definition: math.F90:916
integer(int64) pure function pad88(size, blk)
Definition: math.F90:843
real(real64) pure function, public square_root(x)
Wrapper for sqrt.
Definition: math.F90:985
logical function, public member(n, a)
Considers a(1:ubound(a, 1)) as an integer set and checks if n is a member of it.
Definition: math.F90:646
integer pure function pad4(size, blk)
Definition: math.F90:874
subroutine dlower_triangular_to_hermitian(nn, aa)
Definition: math.F90:1418
pure complex(real64) function, dimension(1:3), public dzcross_product(a, b)
Definition: math.F90:1227
subroutine, public cartesian2hyperspherical(x, u)
Performs a transformation of an n-dimensional vector from Cartesian coordinates to hyperspherical coo...
Definition: math.F90:710
integer pure function, public pad_pow2(size)
create array size, which is padded to powers of 2
Definition: math.F90:887
subroutine zinterpolate_1(xa, ya, x, y)
Definition: math.F90:1615
pure real(real64) function, dimension(1:3), public dcross_product(a, b)
Definition: math.F90:1859
integer(int64) pure function pad48(size, blk)
Definition: math.F90:858
subroutine, public ylmr_real(xx, li, mi, ylm)
This is a Numerical Recipes-based subroutine computes real spherical harmonics ylm at position (x,...
Definition: math.F90:373
subroutine dupper_triangular_to_hermitian(nn, aa)
Definition: math.F90:1391
subroutine, public interpolation_coefficients(nn, xa, xx, cc)
Definition: math.F90:667
subroutine, public generate_rotation_matrix(R, ff, tt)
Generates a rotation matrix R to rotate a vector f to t.
Definition: math.F90:1028
elemental logical function zis_close_scalar(x, y, rtol, atol)
Same as dis_close_scalar for complex numbers.
Definition: math.F90:235
logical pure function, public odd(n)
Returns if n is odd.
Definition: math.F90:700
real(real64) pure function, public ddelta(i, j)
Definition: math.F90:599
subroutine, public make_idx_set(n, out, length, in)
Construct out(1:length) = (/1, ..., n/) if in is not present, out(1:length) = in otherwise.
Definition: math.F90:617
logical function, public is_prime(n)
Definition: math.F90:993
complex(real64) pure function, public phi1(z)
Compute phi1(z) = (exp(z)-1)/z.
Definition: math.F90:952
pure complex(real64) function, dimension(1:3), public zcross_product(a, b)
Definition: math.F90:1668
subroutine, public weights(N, M, cc, side)
Compute the weights for finite-difference calculations:
Definition: math.F90:523
subroutine zinterpolate_2(xa, ya, x, y)
Definition: math.F90:1580
pure complex(real64) function, dimension(dim, dim) zdiagonal_matrix(dim, diag)
Currently only returns a matrix whose diagonal elements are all the same. Note that the integer versi...
Definition: math.F90:1565
subroutine zlower_triangular_to_hermitian(nn, aa)
Definition: math.F90:1431
logical pure function, public is_diagonal(dim, matrix, tol)
Returns true is the matrix of size dim x dim is diagonal, given a relative tolerance.
Definition: math.F90:1472
subroutine zupper_triangular_to_hermitian(nn, aa)
Definition: math.F90:1404
pure integer function, dimension(dim, dim) idiagonal_matrix(dim, diag)
Currently only returns a matrix whose diagonal elements are all the same. Note that the real and comp...
Definition: math.F90:251
integer(int64) pure function pad8(size, blk)
Definition: math.F90:866
subroutine dinterpolate_2(xa, ya, x, y)
Definition: math.F90:1771
subroutine, public numder_ridders(x, h, res, err, f)
Numerical derivative (Ridder`s algorithm).
Definition: math.F90:1172
integer(int64) pure function llog2(xx)
Definition: math.F90:932
recursive real(real64) function, public hermite(n, x)
Definition: math.F90:266
elemental logical function dis_close_scalar(x, y, rtol, atol)
Are and equal within a tolerance.
Definition: math.F90:222
pure complex(real64) function, dimension(1:3), public zdcross_product(a, b)
Definition: math.F90:1240
subroutine, public hyperspherical2cartesian(u, x)
Performs the inverse transformation of cartesian2hyperspherical.
Definition: math.F90:757
subroutine, public dzero_small_elements_matrix(nn, aa, tol)
Definition: math.F90:1460
subroutine, public hypersphere_grad_matrix(grad_matrix, r, x)
Gives the hyperspherical gradient matrix, which contains the derivatives of the Cartesian coordinates...
Definition: math.F90:803
recursive integer function, public factorial(n)
Definition: math.F90:286
int true(void)