(git:58e3e09)
mathlib.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Collection of simple mathematical functions and subroutines
10 !> \par History
11 !> FUNCTION angle updated and FUNCTION dihedral angle added; cleaned
12 !> (13.03.2004,MK)
13 !> \author MK (15.11.1998)
14 ! **************************************************************************************************
15 MODULE mathlib
16 
17  USE kinds, ONLY: default_string_length,&
18  dp
19  USE mathconstants, ONLY: euler,&
20  fac,&
21  pi
22 #include "../base/base_uses.f90"
23 
24  IMPLICIT NONE
25 
26  PRIVATE
27 
28  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mathlib'
29  REAL(KIND=dp), PARAMETER :: eps_geo = 1.0e-6_dp
30 
31  ! Public subroutines
32 
33  PUBLIC :: build_rotmat, &
34  jacobi, &
35  diamat_all, &
36  invmat, &
37  invmat_symm, &
38  invert_matrix, &
42  unit_matrix, diag, &
43  erfc_cutoff, &
45 
46  ! Public functions
47 
48  PUBLIC :: angle, &
49  binomial, &
50  binomial_gen, &
51  multinomial, &
52  det_3x3, &
54  gcd, &
55  inv_3x3, &
56  lcm, &
58  pswitch, &
59  rotate_vector, &
62  get_diag, &
63  set_diag
64 
65  INTERFACE det_3x3
66  MODULE PROCEDURE det_3x3_1, det_3x3_2
67  END INTERFACE
68 
69  INTERFACE invert_matrix
70  MODULE PROCEDURE invert_matrix_d, invert_matrix_z
71  END INTERFACE
72 
73  INTERFACE set_diag
74  MODULE PROCEDURE set_diag_scalar_d, set_diag_scalar_z
75  END INTERFACE
76 
77  INTERFACE swap
78  MODULE PROCEDURE swap_scalar, swap_vector
79  END INTERFACE
80 
81  INTERFACE unit_matrix
82  MODULE PROCEDURE unit_matrix_d, unit_matrix_z
83  END INTERFACE
84 
85 CONTAINS
86 
87 ! **************************************************************************************************
88 !> \brief Polynomial (5th degree) switching function
89 !> f(a) = 1 .... f(b) = 0 with f'(a) = f"(a) = f'(b) = f"(b) = 0
90 !> \param x ...
91 !> \param a ...
92 !> \param b ...
93 !> \param order ...
94 !> \return =0 : f(x)
95 !> \return =1 : f'(x)
96 !> \return =2 : f"(x)
97 ! **************************************************************************************************
98  FUNCTION pswitch(x, a, b, order) RESULT(fx)
99  REAL(kind=dp) :: x, a, b
100  INTEGER :: order
101  REAL(kind=dp) :: fx
102 
103  REAL(kind=dp) :: u, u2, u3
104 
105  cpassert(b > a)
106  IF (x < a .OR. x > b) THEN
107  ! outside switching intervall
108  IF (order > 0) THEN
109  ! derivatives are 0
110  fx = 0.0_dp
111  ELSE
112  IF (x < a) THEN
113  ! x < a => f(x) = 1
114  fx = 1.0_dp
115  ELSE
116  ! x > b => f(x) = 0
117  fx = 0.0_dp
118  END IF
119  END IF
120  ELSE
121  ! renormalized coordinate
122  u = (x - a)/(b - a)
123  SELECT CASE (order)
124  CASE (0)
125  u2 = u*u
126  u3 = u2*u
127  fx = 1._dp - 10._dp*u3 + 15._dp*u2*u2 - 6._dp*u2*u3
128  CASE (1)
129  u2 = u*u
130  fx = -30._dp*u2 + 60._dp*u*u2 - 30._dp*u2*u2
131  fx = fx/(b - a)
132  CASE (2)
133  u2 = u*u
134  fx = -60._dp*u + 180._dp*u2 - 120._dp*u*u2
135  fx = fx/(b - a)**2
136  CASE DEFAULT
137  cpabort('order not defined')
138  END SELECT
139  END IF
140 
141  END FUNCTION pswitch
142 
143 ! **************************************************************************************************
144 !> \brief determines if a value is not normal (e.g. for Inf and Nan)
145 !> based on IO to work also under optimization.
146 !> \param a input value
147 !> \return TRUE for NaN and Inf
148 ! **************************************************************************************************
149  LOGICAL FUNCTION abnormal_value(a)
150  REAL(kind=dp) :: a
151 
152  CHARACTER(LEN=32) :: buffer
153 
154  abnormal_value = .false.
155  ! the function should work when compiled with -ffast-math and similar
156  ! unfortunately, that option asserts that all numbers are normals,
157  ! which the compiler uses to optimize the function to .FALSE. if based on the IEEE module
158  ! therefore, pass this to the Fortran runtime/printf, if things are NaN or Inf, error out.
159  WRITE (buffer, *) a
160  IF (index(buffer, "N") .NE. 0 .OR. index(buffer, "n") .NE. 0) abnormal_value = .true.
161 
162  END FUNCTION abnormal_value
163 
164 ! **************************************************************************************************
165 !> \brief Calculation of the angle between the vectors a and b.
166 !> The angle is returned in radians.
167 !> \param a ...
168 !> \param b ...
169 !> \return ...
170 !> \date 14.10.1998
171 !> \author MK
172 !> \version 1.0
173 ! **************************************************************************************************
174  PURE FUNCTION angle(a, b) RESULT(angle_ab)
175 
176  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: a, b
177  REAL(kind=dp) :: angle_ab
178 
179  REAL(kind=dp) :: length_of_a, length_of_b
180  REAL(kind=dp), DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
181 
182  length_of_a = sqrt(dot_product(a, a))
183  length_of_b = sqrt(dot_product(b, b))
184 
185  IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
186  a_norm(:) = a(:)/length_of_a
187  b_norm(:) = b(:)/length_of_b
188  angle_ab = acos(min(max(dot_product(a_norm, b_norm), -1.0_dp), 1.0_dp))
189  ELSE
190  angle_ab = 0.0_dp
191  END IF
192 
193  END FUNCTION angle
194 
195 ! **************************************************************************************************
196 !> \brief The binomial coefficient n over k for 0 <= k <= n is calculated,
197 !> otherwise zero is returned.
198 !> \param n ...
199 !> \param k ...
200 !> \return ...
201 !> \date 08.03.1999
202 !> \author MK
203 !> \version 1.0
204 ! **************************************************************************************************
205  ELEMENTAL FUNCTION binomial(n, k) RESULT(n_over_k)
206 
207  INTEGER, INTENT(IN) :: n, k
208  REAL(kind=dp) :: n_over_k
209 
210  IF ((k >= 0) .AND. (k <= n)) THEN
211  n_over_k = fac(n)/(fac(n - k)*fac(k))
212  ELSE
213  n_over_k = 0.0_dp
214  END IF
215 
216  END FUNCTION binomial
217 
218 ! **************************************************************************************************
219 !> \brief The generalized binomial coefficient z over k for 0 <= k <= n is calculated.
220 !> (z) z*(z-1)*...*(z-k+2)*(z-k+1)
221 !> ( ) = ---------------------------
222 !> (k) k!
223 !> \param z ...
224 !> \param k ...
225 !> \return ...
226 !> \date 11.11.2019
227 !> \author FS
228 !> \version 1.0
229 ! **************************************************************************************************
230  ELEMENTAL FUNCTION binomial_gen(z, k) RESULT(z_over_k)
231 
232  REAL(kind=dp), INTENT(IN) :: z
233  INTEGER, INTENT(IN) :: k
234  REAL(kind=dp) :: z_over_k
235 
236  INTEGER :: i
237 
238  IF (k >= 0) THEN
239  z_over_k = 1.0_dp
240  DO i = 1, k
241  z_over_k = z_over_k*(z - i + 1)/real(i, dp)
242  END DO
243  ELSE
244  z_over_k = 0.0_dp
245  END IF
246 
247  END FUNCTION binomial_gen
248 
249 ! **************************************************************************************************
250 !> \brief Calculates the multinomial coefficients
251 !> \param n ...
252 !> \param k ...
253 !> \return ...
254 !> \author Ole Schuett
255 ! **************************************************************************************************
256  PURE FUNCTION multinomial(n, k) RESULT(res)
257 
258  INTEGER, INTENT(IN) :: n
259  INTEGER, DIMENSION(:), INTENT(IN) :: k
260  REAL(kind=dp) :: res
261 
262  INTEGER :: i
263  REAL(kind=dp) :: denom
264 
265  IF (all(k >= 0) .AND. sum(k) == n) THEN
266  denom = 1.0_dp
267  DO i = 1, SIZE(k)
268  denom = denom*fac(k(i))
269  END DO
270  res = fac(n)/denom
271  ELSE
272  res = 0.0_dp
273  END IF
274 
275  END FUNCTION multinomial
276 
277 ! **************************************************************************************************
278 !> \brief The rotation matrix rotmat which rotates a vector about a
279 !> rotation axis defined by the vector a is build up.
280 !> The rotation angle is phi (radians).
281 !> \param phi ...
282 !> \param a ...
283 !> \param rotmat ...
284 !> \date 16.10.1998
285 !> \author MK
286 !> \version 1.0
287 ! **************************************************************************************************
288  PURE SUBROUTINE build_rotmat(phi, a, rotmat)
289  REAL(kind=dp), INTENT(IN) :: phi
290  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
291  REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT) :: rotmat
292 
293  REAL(kind=dp) :: cosp, cost, length_of_a, sinp
294  REAL(kind=dp), DIMENSION(3) :: d
295 
296  length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
297  ! Check the length of the vector a
298  IF (length_of_a > eps_geo) THEN
299 
300  d(:) = a(:)/length_of_a
301 
302  cosp = cos(phi)
303  sinp = sin(phi)
304  cost = 1.0_dp - cosp
305 
306  rotmat(1, 1) = d(1)*d(1)*cost + cosp
307  rotmat(1, 2) = d(1)*d(2)*cost - d(3)*sinp
308  rotmat(1, 3) = d(1)*d(3)*cost + d(2)*sinp
309  rotmat(2, 1) = d(2)*d(1)*cost + d(3)*sinp
310  rotmat(2, 2) = d(2)*d(2)*cost + cosp
311  rotmat(2, 3) = d(2)*d(3)*cost - d(1)*sinp
312  rotmat(3, 1) = d(3)*d(1)*cost - d(2)*sinp
313  rotmat(3, 2) = d(3)*d(2)*cost + d(1)*sinp
314  rotmat(3, 3) = d(3)*d(3)*cost + cosp
315  ELSE
316  CALL unit_matrix(rotmat)
317  END IF
318 
319  END SUBROUTINE build_rotmat
320 
321 ! **************************************************************************************************
322 !> \brief Returns the determinante of the 3x3 matrix a.
323 !> \param a ...
324 !> \return ...
325 !> \date 13.03.2004
326 !> \author MK
327 !> \version 1.0
328 ! **************************************************************************************************
329  PURE FUNCTION det_3x3_1(a) RESULT(det_a)
330  REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: a
331  REAL(kind=dp) :: det_a
332 
333  det_a = a(1, 1)*(a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2)) + &
334  a(1, 2)*(a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3)) + &
335  a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
336 
337  END FUNCTION det_3x3_1
338 
339 ! **************************************************************************************************
340 !> \brief Returns the determinante of the 3x3 matrix a given by its columns.
341 !> \param a1 ...
342 !> \param a2 ...
343 !> \param a3 ...
344 !> \return ...
345 !> \date 13.03.2004
346 !> \author MK
347 !> \version 1.0
348 ! **************************************************************************************************
349  PURE FUNCTION det_3x3_2(a1, a2, a3) RESULT(det_a)
350  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3
351  REAL(kind=dp) :: det_a
352 
353  det_a = a1(1)*(a2(2)*a3(3) - a3(2)*a2(3)) + &
354  a2(1)*(a3(2)*a1(3) - a1(2)*a3(3)) + &
355  a3(1)*(a1(2)*a2(3) - a2(2)*a1(3))
356 
357  END FUNCTION det_3x3_2
358 
359 ! **************************************************************************************************
360 !> \brief Diagonalize the symmetric n by n matrix a using the LAPACK
361 !> library. Only the upper triangle of matrix a is used.
362 !> Externals (LAPACK 3.0)
363 !> \param a ...
364 !> \param eigval ...
365 !> \param dac ...
366 !> \date 29.03.1999
367 !> \par Variables
368 !> - a : Symmetric matrix to be diagonalized (input; upper triangle) ->
369 !> - eigenvectors of the matrix a (output).
370 !> - dac : If true, then the divide-and-conquer algorithm is applied.
371 !> - eigval : Eigenvalues of the matrix a (output).
372 !> \author MK
373 !> \version 1.0
374 ! **************************************************************************************************
375  SUBROUTINE diamat_all(a, eigval, dac)
376  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
377  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigval
378  LOGICAL, INTENT(IN), OPTIONAL :: dac
379 
380  CHARACTER(len=*), PARAMETER :: routinen = 'diamat_all'
381 
382  INTEGER :: handle, info, liwork, lwork, n, nb
383  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
384  INTEGER, EXTERNAL :: ilaenv
385  LOGICAL :: divide_and_conquer
386  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
387 
388  EXTERNAL dsyev, dsyevd
389 
390  CALL timeset(routinen, handle)
391 
392  ! Get the size of the matrix a
393  n = SIZE(a, 1)
394 
395  ! Check the size of matrix a
396  IF (SIZE(a, 2) /= n) THEN
397  cpabort("Check the size of matrix a (parameter #1)")
398  END IF
399 
400  ! Check the size of vector eigval
401  IF (SIZE(eigval) /= n) THEN
402  cpabort("The dimension of vector eigval is too small")
403  END IF
404 
405  ! Check, if the divide-and-conquer algorithm is requested
406 
407  IF (PRESENT(dac)) THEN
408  divide_and_conquer = dac
409  ELSE
410  divide_and_conquer = .false.
411  END IF
412 
413  ! Get the optimal work storage size
414 
415  IF (divide_and_conquer) THEN
416  lwork = 2*n**2 + 6*n + 1
417  liwork = 5*n + 3
418  ELSE
419  nb = ilaenv(1, "DSYTRD", "U", n, -1, -1, -1)
420  lwork = (nb + 2)*n
421  END IF
422 
423  ! Allocate work storage
424 
425  ALLOCATE (work(lwork))
426  IF (divide_and_conquer) THEN
427  ALLOCATE (iwork(liwork))
428  END IF
429 
430  ! Diagonalize the matrix a
431 
432  info = 0
433  IF (divide_and_conquer) THEN
434  CALL dsyevd("V", "U", n, a, n, eigval, work, lwork, iwork, liwork, info)
435  ELSE
436  CALL dsyev("V", "U", n, a, n, eigval, work, lwork, info)
437  END IF
438 
439  IF (info /= 0) THEN
440  IF (divide_and_conquer) THEN
441  cpabort("The matrix diagonalization with dsyevd failed")
442  ELSE
443  cpabort("The matrix diagonalization with dsyev failed")
444  END IF
445  END IF
446 
447  ! Release work storage
448 
449  DEALLOCATE (work)
450 
451  IF (divide_and_conquer) THEN
452  DEALLOCATE (iwork)
453  END IF
454 
455  CALL timestop(handle)
456 
457  END SUBROUTINE diamat_all
458 
459 ! **************************************************************************************************
460 !> \brief Returns the dihedral angle, i.e. the angle between the planes
461 !> defined by the vectors (-ab,bc) and (cd,-bc).
462 !> The dihedral angle is returned in radians.
463 !> \param ab ...
464 !> \param bc ...
465 !> \param cd ...
466 !> \return ...
467 !> \date 13.03.2004
468 !> \author MK
469 !> \version 1.0
470 ! **************************************************************************************************
471  PURE FUNCTION dihedral_angle(ab, bc, cd) RESULT(dihedral_angle_abcd)
472  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: ab, bc, cd
473  REAL(kind=dp) :: dihedral_angle_abcd
474 
475  REAL(kind=dp) :: det_abcd
476  REAL(kind=dp), DIMENSION(3) :: abc, bcd
477 
478  abc = vector_product(bc, -ab)
479  bcd = vector_product(cd, -bc)
480  ! Calculate the normal vectors of the planes
481  ! defined by the points a,b,c and b,c,d
482 
483  det_abcd = det_3x3(abc, bcd, -bc)
484  dihedral_angle_abcd = sign(1.0_dp, det_abcd)*angle(abc, bcd)
485 
486  END FUNCTION dihedral_angle
487 
488 ! **************************************************************************************************
489 !> \brief Return the diagonal elements of matrix a as a vector.
490 !> \param a ...
491 !> \return ...
492 !> \date 20.11.1998
493 !> \author MK
494 !> \version 1.0
495 ! **************************************************************************************************
496  PURE FUNCTION get_diag(a) RESULT(a_diag)
497  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a
498  REAL(kind=dp), &
499  DIMENSION(MIN(SIZE(a, 1), SIZE(a, 2))) :: a_diag
500 
501  INTEGER :: i, n
502 
503  n = min(SIZE(a, 1), SIZE(a, 2))
504 
505  DO i = 1, n
506  a_diag(i) = a(i, i)
507  END DO
508 
509  END FUNCTION get_diag
510 
511 ! **************************************************************************************************
512 !> \brief Returns the inverse of the 3 x 3 matrix a.
513 !> \param a ...
514 !> \return ...
515 !> \date 13.03.2004
516 !> \author MK
517 !> \version 1.0
518 ! **************************************************************************************************
519  PURE FUNCTION inv_3x3(a) RESULT(a_inv)
520 
521  REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: a
522  REAL(kind=dp), DIMENSION(3, 3) :: a_inv
523 
524  REAL(kind=dp) :: det_a
525 
526  det_a = 1.0_dp/det_3x3(a)
527 
528  a_inv(1, 1) = (a(2, 2)*a(3, 3) - a(3, 2)*a(2, 3))*det_a
529  a_inv(2, 1) = (a(2, 3)*a(3, 1) - a(3, 3)*a(2, 1))*det_a
530  a_inv(3, 1) = (a(2, 1)*a(3, 2) - a(3, 1)*a(2, 2))*det_a
531 
532  a_inv(1, 2) = (a(1, 3)*a(3, 2) - a(3, 3)*a(1, 2))*det_a
533  a_inv(2, 2) = (a(1, 1)*a(3, 3) - a(3, 1)*a(1, 3))*det_a
534  a_inv(3, 2) = (a(1, 2)*a(3, 1) - a(3, 2)*a(1, 1))*det_a
535 
536  a_inv(1, 3) = (a(1, 2)*a(2, 3) - a(2, 2)*a(1, 3))*det_a
537  a_inv(2, 3) = (a(1, 3)*a(2, 1) - a(2, 3)*a(1, 1))*det_a
538  a_inv(3, 3) = (a(1, 1)*a(2, 2) - a(2, 1)*a(1, 2))*det_a
539 
540  END FUNCTION inv_3x3
541 
542 ! **************************************************************************************************
543 !> \brief returns inverse of matrix using the lapack routines DGETRF and DGETRI
544 !> \param a ...
545 !> \param info ...
546 ! **************************************************************************************************
547  SUBROUTINE invmat(a, info)
548  REAL(kind=dp), INTENT(INOUT) :: a(:, :)
549  INTEGER, INTENT(OUT) :: info
550 
551  INTEGER :: lwork, n
552  INTEGER, ALLOCATABLE :: ipiv(:)
553  REAL(kind=dp), ALLOCATABLE :: work(:)
554 
555  n = SIZE(a, 1)
556  lwork = 20*n
557  ALLOCATE (ipiv(n))
558  ALLOCATE (work(lwork))
559  ipiv = 0
560  work = 0._dp
561  info = 0
562  CALL dgetrf(n, n, a, n, ipiv, info)
563  IF (info == 0) THEN
564  CALL dgetri(n, a, n, ipiv, work, lwork, info)
565  END IF
566  DEALLOCATE (ipiv, work)
567  END SUBROUTINE invmat
568 
569 ! **************************************************************************************************
570 !> \brief returns inverse of real symmetric, positive definite matrix
571 !> \param a matrix
572 !> \param cholesky_triangle if cholesky decomposition of a was already done
573 !> using dpotrf, indicating if the upper or lower triangle of a is
574 !> stored. If not given, cholesky decomposition of a will be done
575 !> before inversion.
576 !> \author Dorothea Golze [02.2015]
577 ! **************************************************************************************************
578  SUBROUTINE invmat_symm(a, cholesky_triangle)
579  REAL(kind=dp), INTENT(INOUT) :: a(:, :)
580  CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: cholesky_triangle
581 
582  CHARACTER(LEN=*), PARAMETER :: routinen = 'invmat_symm'
583 
584  CHARACTER(LEN=1) :: my_triangle
585  INTEGER :: handle, i, info, n
586 
587  CALL timeset(routinen, handle)
588 
589  n = SIZE(a, 1)
590  info = 0
591 
592  IF (PRESENT(cholesky_triangle)) THEN
593  my_triangle = cholesky_triangle
594  ELSE
595  my_triangle = "U"
596  END IF
597 
598  ! do cholesky decomposition
599  IF (.NOT. PRESENT(cholesky_triangle)) THEN
600  CALL dpotrf(my_triangle, n, a, n, info)
601  IF (info /= 0) THEN
602  cpabort("DPOTRF failed")
603  END IF
604  END IF
605 
606  ! do inversion using the cholesky decomposition
607  CALL dpotri(my_triangle, n, a, n, info)
608  IF (info /= 0) THEN
609  cpabort("Matrix inversion failed")
610  END IF
611 
612  IF (my_triangle == "U") THEN
613  DO i = 1, n - 1
614  a(i + 1:n, i) = a(i, i + 1:n)
615  END DO
616  ELSE
617  DO i = 1, n - 1
618  a(i, i + 1:n) = a(i + 1:n, i)
619  END DO
620  END IF
621 
622  CALL timestop(handle)
623 
624  END SUBROUTINE invmat_symm
625 
626 ! **************************************************************************************************
627 !> \brief Compute the inverse of the n by n real matrix a using the LAPACK
628 !> library
629 !> \param a ...
630 !> \param a_inverse ...
631 !> \param eval_error ...
632 !> \param option ...
633 !> \param improve ...
634 !> \date 23.03.1999
635 !> \par Variables
636 !> - a : Real matrix to be inverted (input).
637 !> - a_inverse: Inverse of the matrix a (output).
638 !> - a_lu : LU factorization of matrix a.
639 !> - a_norm : Norm of matrix a.
640 !> - error : Estimated error of the inversion.
641 !> - r_cond : Reciprocal condition number of the matrix a.
642 !> - trans : "N" => invert a
643 !> - "T" => invert transpose(a)
644 !> \author MK
645 !> \version 1.0
646 !> \note NB add improve argument, used to disable call to dgerfs
647 ! **************************************************************************************************
648  SUBROUTINE invert_matrix_d(a, a_inverse, eval_error, option, improve)
649  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a
650  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: a_inverse
651  REAL(kind=dp), INTENT(OUT) :: eval_error
652  CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: option
653  LOGICAL, INTENT(IN), OPTIONAL :: improve
654 
655  CHARACTER(LEN=1) :: norm, trans
656  CHARACTER(LEN=default_string_length) :: message
657  INTEGER :: info, iter, n
658  INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv, iwork
659  LOGICAL :: do_improve
660  REAL(kind=dp) :: a_norm, old_eval_error, r_cond
661  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, work
662  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
663  REAL(kind=dp), EXTERNAL :: dlange
664 
665  EXTERNAL dgecon, dgerfs, dgetrf, dgetrs
666 
667  ! Check for optional parameter
668  IF (PRESENT(option)) THEN
669  trans = option
670  ELSE
671  trans = "N"
672  END IF
673 
674  IF (PRESENT(improve)) THEN
675  do_improve = improve
676  ELSE
677  do_improve = .true.
678  END IF
679 
680  ! Get the dimension of matrix a
681  n = SIZE(a, 1)
682 
683  ! Check array dimensions
684  IF (n == 0) THEN
685  cpabort("Matrix to be inverted of zero size")
686  END IF
687 
688  IF (n /= SIZE(a, 2)) THEN
689  cpabort("Check the array bounds of parameter #1")
690  END IF
691 
692  IF ((n /= SIZE(a_inverse, 1)) .OR. &
693  (n /= SIZE(a_inverse, 2))) THEN
694  cpabort("Check the array bounds of parameter #2")
695  END IF
696 
697  ! Allocate work storage
698  ALLOCATE (a_lu(n, n))
699 
700  ALLOCATE (b(n, n))
701 
702  ALLOCATE (berr(n))
703 
704  ALLOCATE (ferr(n))
705 
706  ALLOCATE (ipiv(n))
707 
708  ALLOCATE (iwork(n))
709 
710  ALLOCATE (work(4*n))
711 
712  a_lu(1:n, 1:n) = a(1:n, 1:n)
713 
714  ! Compute the LU factorization of the matrix a
715  CALL dgetrf(n, n, a_lu, n, ipiv, info)
716 
717  IF (info /= 0) THEN
718  cpabort("The LU factorization in dgetrf failed")
719  END IF
720 
721  ! Compute the norm of the matrix a
722 
723  IF (trans == "N") THEN
724  norm = '1'
725  ELSE
726  norm = 'I'
727  END IF
728 
729  a_norm = dlange(norm, n, n, a, n, work)
730 
731  ! Compute the reciprocal of the condition number of a
732 
733  CALL dgecon(norm, n, a_lu, n, a_norm, r_cond, work, iwork, info)
734 
735  IF (info /= 0) THEN
736  cpabort("The computation of the condition number in dgecon failed")
737  END IF
738 
739  IF (r_cond < epsilon(0.0_dp)) THEN
740  WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
741  CALL cp_abort(__location__, &
742  "Bad condition number "//trim(message)//" (smaller than the machine "// &
743  "working precision)")
744  END IF
745 
746  ! Solve a system of linear equations using the LU factorization computed by dgetrf
747 
748  CALL unit_matrix(a_inverse)
749 
750  CALL dgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
751 
752  IF (info /= 0) THEN
753  cpabort("Solving the system of linear equations in dgetrs failed")
754  END IF
755 
756  ! Improve the computed solution iteratively
757  CALL unit_matrix(b) ! Initialize right-hand sides
758 
759  eval_error = 0.0_dp
760 
761  IF (do_improve) THEN
762  DO iter = 1, 10
763 
764  CALL dgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
765  work, iwork, info)
766 
767  IF (info /= 0) THEN
768  cpabort("Improving the computed solution in dgerfs failed")
769  END IF
770 
771  old_eval_error = eval_error
772  eval_error = maxval(ferr)
773 
774  IF (abs(eval_error - old_eval_error) <= epsilon(1.0_dp)) EXIT
775 
776  END DO
777  END IF
778 
779  ! Release work storage
780  DEALLOCATE (work)
781  DEALLOCATE (iwork)
782  DEALLOCATE (ipiv)
783  DEALLOCATE (ferr)
784  DEALLOCATE (berr)
785  DEALLOCATE (b)
786  DEALLOCATE (a_lu)
787 
788  END SUBROUTINE invert_matrix_d
789 
790 ! **************************************************************************************************
791 !> \brief Compute the inverse of the n by n complex matrix a using the LAPACK
792 !> library
793 !> \param a ...
794 !> \param a_inverse ...
795 !> \param eval_error ...
796 !> \param option ...
797 !> \date 08.06.2009
798 !> \par Variables
799 !> - a : Complex matrix to be inverted (input).
800 !> - a_inverse: Inverse of the matrix a (output).
801 !> - a_lu : LU factorization of matrix a.
802 !> - a_norm : Norm of matrix a.
803 !> - error : Estimated error of the inversion.
804 !> - r_cond : Reciprocal condition number of the matrix a.
805 !> - trans : "N" => invert a
806 !> - "T" => invert transpose(a)
807 !> \author MK
808 !> \version 1.0
809 ! **************************************************************************************************
810  SUBROUTINE invert_matrix_z(a, a_inverse, eval_error, option)
811  COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: a
812  COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: a_inverse
813  REAL(kind=dp), INTENT(OUT) :: eval_error
814  CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: option
815 
816  CHARACTER(LEN=1) :: norm, trans
817  CHARACTER(LEN=default_string_length) :: message
818  COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
819  COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
820  INTEGER :: info, iter, n
821  INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
822  REAL(kind=dp) :: a_norm, old_eval_error, r_cond
823  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, rwork
824  REAL(kind=dp), EXTERNAL :: zlange
825 
826  EXTERNAL zgecon, zgerfs, zgetrf, zgetrs
827 
828  ! Check for optional parameter
829  IF (PRESENT(option)) THEN
830  trans = option
831  ELSE
832  trans = "N"
833  END IF
834 
835  ! Get the dimension of matrix a
836  n = SIZE(a, 1)
837 
838  ! Check array dimensions
839  IF (n == 0) THEN
840  cpabort("Matrix to be inverted of zero size")
841  END IF
842 
843  IF (n /= SIZE(a, 2)) THEN
844  cpabort("Check the array bounds of parameter #1")
845  END IF
846 
847  IF ((n /= SIZE(a_inverse, 1)) .OR. &
848  (n /= SIZE(a_inverse, 2))) THEN
849  cpabort("Check the array bounds of parameter #2")
850  END IF
851 
852  ! Allocate work storage
853  ALLOCATE (a_lu(n, n))
854 
855  ALLOCATE (b(n, n))
856 
857  ALLOCATE (berr(n))
858 
859  ALLOCATE (ferr(n))
860 
861  ALLOCATE (ipiv(n))
862 
863  ALLOCATE (rwork(2*n))
864 
865  ALLOCATE (work(2*n))
866 
867  a_lu(1:n, 1:n) = a(1:n, 1:n)
868 
869  ! Compute the LU factorization of the matrix a
870  CALL zgetrf(n, n, a_lu, n, ipiv, info)
871 
872  IF (info /= 0) THEN
873  cpabort("The LU factorization in dgetrf failed")
874  END IF
875 
876  ! Compute the norm of the matrix a
877 
878  IF (trans == "N") THEN
879  norm = '1'
880  ELSE
881  norm = 'I'
882  END IF
883 
884  a_norm = zlange(norm, n, n, a, n, work)
885 
886  ! Compute the reciprocal of the condition number of a
887 
888  CALL zgecon(norm, n, a_lu, n, a_norm, r_cond, work, rwork, info)
889 
890  IF (info /= 0) THEN
891  cpabort("The computation of the condition number in dgecon failed")
892  END IF
893 
894  IF (r_cond < epsilon(0.0_dp)) THEN
895  WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
896  CALL cp_abort(__location__, &
897  "Bad condition number "//trim(message)//" (smaller than the machine "// &
898  "working precision)")
899  END IF
900 
901  ! Solve a system of linear equations using the LU factorization computed by dgetrf
902 
903  CALL unit_matrix(a_inverse)
904 
905  CALL zgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
906 
907  IF (info /= 0) THEN
908  cpabort("Solving the system of linear equations in dgetrs failed")
909  END IF
910 
911  ! Improve the computed solution iteratively
912  CALL unit_matrix(b) ! Initialize right-hand sides
913 
914  eval_error = 0.0_dp
915 
916  DO iter = 1, 10
917 
918  CALL zgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
919  work, rwork, info)
920 
921  IF (info /= 0) THEN
922  cpabort("Improving the computed solution in dgerfs failed")
923  END IF
924 
925  old_eval_error = eval_error
926  eval_error = maxval(ferr)
927 
928  IF (abs(eval_error - old_eval_error) <= epsilon(1.0_dp)) EXIT
929 
930  END DO
931 
932  ! Release work storage
933  DEALLOCATE (work)
934  DEALLOCATE (rwork)
935  DEALLOCATE (ipiv)
936  DEALLOCATE (ferr)
937  DEALLOCATE (berr)
938  DEALLOCATE (b)
939  DEALLOCATE (a_lu)
940 
941  END SUBROUTINE invert_matrix_z
942 
943 ! **************************************************************************************************
944 !> \brief returns the pseudoinverse of a real, square matrix using singular
945 !> value decomposition
946 !> \param a matrix a
947 !> \param a_pinverse pseudoinverse of matrix a
948 !> \param rskip parameter for setting small singular values to zero
949 !> \param determinant determinant of matrix a (optional output)
950 !> \param sval array holding singular values of matrix a (optional output)
951 !> \author Dorothea Golze [02.2015]
952 ! **************************************************************************************************
953  SUBROUTINE get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
954 
955  REAL(kind=dp), DIMENSION(:, :) :: a, a_pinverse
956  REAL(kind=dp), INTENT(IN) :: rskip
957  REAL(kind=dp), INTENT(OUT), OPTIONAL :: determinant
958  REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
959  OPTIONAL, POINTER :: sval
960 
961  CHARACTER(LEN=*), PARAMETER :: routinen = 'get_pseudo_inverse_svd'
962 
963  INTEGER :: handle, i, info, lwork, n
964  INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
965  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sig, work
966  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sig_plus, temp_mat, u, vt
967 
968  CALL timeset(routinen, handle)
969 
970  n = SIZE(a, 1)
971  ALLOCATE (u(n, n), vt(n, n), sig(n), sig_plus(n, n), iwork(8*n), work(1), temp_mat(n, n))
972  u(:, :) = 0.0_dp
973  vt(:, :) = 0.0_dp
974  sig(:) = 0.0_dp
975  sig_plus = 0.0_dp
976  work = 0.0_dp
977  iwork = 0
978  IF (PRESENT(determinant)) determinant = 1.0_dp
979 
980  ! work size query
981  lwork = -1
982  CALL dgesdd('A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
983  lwork, iwork(1), info)
984 
985  IF (info /= 0) THEN
986  cpabort("ERROR in DGESDD: Could not retrieve work array sizes")
987  END IF
988  lwork = int(work(1))
989  DEALLOCATE (work)
990  ALLOCATE (work(lwork))
991 
992  ! do SVD
993  CALL dgesdd('A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
994  lwork, iwork(1), info)
995 
996  IF (info /= 0) THEN
997  cpabort("SVD failed")
998  END IF
999 
1000  IF (PRESENT(sval)) THEN
1001  cpassert(.NOT. ASSOCIATED(sval))
1002  ALLOCATE (sval(n))
1003  sval(:) = sig
1004  END IF
1005 
1006  ! set singular values that are too small to zero
1007  DO i = 1, n
1008  IF (sig(i) > rskip*maxval(sig)) THEN
1009  IF (PRESENT(determinant)) &
1010  determinant = determinant*sig(i)
1011  sig_plus(i, i) = 1._dp/sig(i)
1012  ELSE
1013  sig_plus(i, i) = 0.0_dp
1014  END IF
1015  END DO
1016 
1017  ! build pseudoinverse: V*sig_plus*UT
1018  CALL dgemm("N", "T", n, n, n, 1._dp, sig_plus, n, u, n, 0._dp, temp_mat, n)
1019  CALL dgemm("T", "N", n, n, n, 1._dp, vt, n, temp_mat, n, 0._dp, a_pinverse, n)
1020 
1021  DEALLOCATE (u, vt, sig, iwork, work, sig_plus, temp_mat)
1022 
1023  CALL timestop(handle)
1024 
1025  END SUBROUTINE get_pseudo_inverse_svd
1026 
1027 ! **************************************************************************************************
1028 !> \brief returns the pseudoinverse of a real, symmetric and positive definite
1029 !> matrix using diagonalization.
1030 !> \param a matrix a
1031 !> \param a_pinverse pseudoinverse of matrix a
1032 !> \param rskip parameter for setting small eigenvalues to zero
1033 !> \author Dorothea Golze [02.2015]
1034 ! **************************************************************************************************
1035  SUBROUTINE get_pseudo_inverse_diag(a, a_pinverse, rskip)
1036 
1037  REAL(kind=dp), DIMENSION(:, :) :: a, a_pinverse
1038  REAL(kind=dp), INTENT(IN) :: rskip
1039 
1040  CHARACTER(LEN=*), PARAMETER :: routinen = 'get_pseudo_inverse_diag'
1041 
1042  INTEGER :: handle, i, info, lwork, n
1043  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig, work
1044  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dinv, temp_mat
1045 
1046  CALL timeset(routinen, handle)
1047 
1048  info = 0
1049  n = SIZE(a, 1)
1050  ALLOCATE (dinv(n, n), eig(n), work(1), temp_mat(n, n))
1051  dinv(:, :) = 0.0_dp
1052  eig(:) = 0.0_dp
1053  work(:) = 0.0_dp
1054  temp_mat = 0.0_dp
1055 
1056  ! work size query
1057  lwork = -1
1058  CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1059  IF (info /= 0) THEN
1060  cpabort("ERROR in DSYEV: Could not retrieve work array sizes")
1061  END IF
1062  lwork = int(work(1))
1063  DEALLOCATE (work)
1064  ALLOCATE (work(lwork))
1065  work = 0.0_dp
1066 
1067  ! get eigenvalues and eigenvectors
1068  CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1069 
1070  IF (info /= 0) THEN
1071  cpabort("Matrix diagonalization failed")
1072  END IF
1073 
1074  ! set eigenvalues that are too small to zero
1075  DO i = 1, n
1076  IF (eig(i) > rskip*maxval(eig)) THEN
1077  dinv(i, i) = 1.0_dp/eig(i)
1078  ELSE
1079  dinv(i, i) = 0._dp
1080  END IF
1081  END DO
1082 
1083  ! build pseudoinverse: U*dinv*UT
1084  CALL dgemm("N", "T", n, n, n, 1._dp, dinv, n, a, n, 0._dp, temp_mat, n)
1085  CALL dgemm("N", "N", n, n, n, 1._dp, a, n, temp_mat, n, 0._dp, a_pinverse, n)
1086 
1087  DEALLOCATE (eig, work, dinv, temp_mat)
1088 
1089  CALL timestop(handle)
1090 
1091  END SUBROUTINE get_pseudo_inverse_diag
1092 
1093 ! **************************************************************************************************
1094 !> \brief Reflection of the vector a through a mirror plane defined by the
1095 !> normal vector b. The reflected vector a is stored in a_mirror.
1096 !> \param a ...
1097 !> \param b ...
1098 !> \return ...
1099 !> \date 16.10.1998
1100 !> \author MK
1101 !> \version 1.0
1102 ! **************************************************************************************************
1103  PURE FUNCTION reflect_vector(a, b) RESULT(a_mirror)
1104 
1105  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a, b
1106  REAL(kind=dp), DIMENSION(3) :: a_mirror
1107 
1108  REAL(kind=dp) :: length_of_b, scapro
1109  REAL(kind=dp), DIMENSION(3) :: d
1110 
1111  length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1112 
1113  IF (length_of_b > eps_geo) THEN
1114 
1115  d(:) = b(:)/length_of_b
1116 
1117  ! Calculate the mirror image a_mirror of the vector a
1118  scapro = a(1)*d(1) + a(2)*d(2) + a(3)*d(3)
1119 
1120  a_mirror(:) = a(:) - 2.0_dp*scapro*d(:)
1121 
1122  ELSE
1123 
1124  a_mirror(:) = 0.0_dp
1125 
1126  END IF
1127 
1128  END FUNCTION reflect_vector
1129 
1130 ! **************************************************************************************************
1131 !> \brief Rotation of the vector a about an rotation axis defined by the
1132 !> vector b. The rotation angle is phi (radians). The rotated vector
1133 !> a is stored in a_rot.
1134 !> \param a ...
1135 !> \param phi ...
1136 !> \param b ...
1137 !> \return ...
1138 !> \date 16.10.1998
1139 !> \author MK
1140 !> \version 1.0
1141 ! **************************************************************************************************
1142  PURE FUNCTION rotate_vector(a, phi, b) RESULT(a_rot)
1143  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1144  REAL(kind=dp), INTENT(IN) :: phi
1145  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: b
1146  REAL(kind=dp), DIMENSION(3) :: a_rot
1147 
1148  REAL(kind=dp) :: length_of_b
1149  REAL(kind=dp), DIMENSION(3, 3) :: rotmat
1150 
1151  length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1152  IF (length_of_b > eps_geo) THEN
1153 
1154  ! Build up the rotation matrix rotmat
1155  CALL build_rotmat(phi, b, rotmat)
1156 
1157  ! Rotate the vector a by phi about the axis defined by vector b
1158  a_rot(:) = matmul(rotmat, a)
1159 
1160  ELSE
1161 
1162  a_rot(:) = 0.0_dp
1163 
1164  END IF
1165 
1166  END FUNCTION rotate_vector
1167 
1168 ! **************************************************************************************************
1169 !> \brief Set the diagonal elements of matrix a to b.
1170 !> \param a ...
1171 !> \param b ...
1172 !> \date 20.11.1998
1173 !> \author MK
1174 !> \version 1.0
1175 ! **************************************************************************************************
1176  PURE SUBROUTINE set_diag_scalar_d(a, b)
1177  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1178  REAL(kind=dp), INTENT(IN) :: b
1179 
1180  INTEGER :: i, n
1181 
1182  n = min(SIZE(a, 1), SIZE(a, 2))
1183  DO i = 1, n
1184  a(i, i) = b
1185  END DO
1186 
1187  END SUBROUTINE set_diag_scalar_d
1188 
1189 ! **************************************************************************************************
1190 !> \brief ...
1191 !> \param a ...
1192 !> \param b ...
1193 ! **************************************************************************************************
1194  PURE SUBROUTINE set_diag_scalar_z(a, b)
1195  COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1196  COMPLEX(KIND=dp), INTENT(IN) :: b
1197 
1198  INTEGER :: i, n
1199 
1200  n = min(SIZE(a, 1), SIZE(a, 2))
1201  DO i = 1, n
1202  a(i, i) = b
1203  END DO
1204 
1205  END SUBROUTINE set_diag_scalar_z
1206 
1207 ! **************************************************************************************************
1208 !> \brief Symmetrize the matrix a.
1209 !> \param a ...
1210 !> \param option ...
1211 !> \date 16.10.1998
1212 !> \author MK
1213 !> \version 1.0
1214 ! **************************************************************************************************
1215  SUBROUTINE symmetrize_matrix(a, option)
1216  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1217  CHARACTER(LEN=*), INTENT(IN) :: option
1218 
1219  INTEGER :: i, n
1220 
1221  n = min(SIZE(a, 1), SIZE(a, 2))
1222 
1223  IF (option == "lower_to_upper") THEN
1224  DO i = 1, n - 1
1225  a(i, i + 1:n) = a(i + 1:n, i)
1226  END DO
1227  ELSE IF (option == "upper_to_lower") THEN
1228  DO i = 1, n - 1
1229  a(i + 1:n, i) = a(i, i + 1:n)
1230  END DO
1231  ELSE IF (option == "anti_lower_to_upper") THEN
1232  DO i = 1, n - 1
1233  a(i, i + 1:n) = -a(i + 1:n, i)
1234  END DO
1235  ELSE IF (option == "anti_upper_to_lower") THEN
1236  DO i = 1, n - 1
1237  a(i + 1:n, i) = -a(i, i + 1:n)
1238  END DO
1239  ELSE
1240  cpabort("Invalid option <"//trim(option)//"> was specified for parameter #2")
1241  END IF
1242 
1243  END SUBROUTINE symmetrize_matrix
1244 
1245 ! **************************************************************************************************
1246 !> \brief Set the matrix a to be a unit matrix.
1247 !> \param a ...
1248 !> \date 16.10.1998
1249 !> \author MK
1250 !> \version 1.0
1251 ! **************************************************************************************************
1252  PURE SUBROUTINE unit_matrix_d(a)
1253  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1254 
1255  a(:, :) = 0.0_dp
1256  CALL set_diag(a, 1.0_dp)
1257 
1258  END SUBROUTINE unit_matrix_d
1259 
1260 ! **************************************************************************************************
1261 !> \brief ...
1262 !> \param a ...
1263 ! **************************************************************************************************
1264  PURE SUBROUTINE unit_matrix_z(a)
1265  COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1266 
1267  a(:, :) = (0.0_dp, 0.0_dp)
1268  CALL set_diag(a, (1.0_dp, 0.0_dp))
1269 
1270  END SUBROUTINE unit_matrix_z
1271 
1272 ! **************************************************************************************************
1273 !> \brief Calculation of the vector product c = a x b.
1274 !> \param a ...
1275 !> \param b ...
1276 !> \return ...
1277 !> \date 16.10.1998
1278 !> \author MK
1279 !> \version 1.0
1280 ! **************************************************************************************************
1281  PURE FUNCTION vector_product(a, b) RESULT(c)
1282  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a, b
1283  REAL(kind=dp), DIMENSION(3) :: c
1284 
1285  c(1) = a(2)*b(3) - a(3)*b(2)
1286  c(2) = a(3)*b(1) - a(1)*b(3)
1287  c(3) = a(1)*b(2) - a(2)*b(1)
1288 
1289  END FUNCTION vector_product
1290 
1291 ! **************************************************************************************************
1292 !> \brief computes the greatest common divisor of two number
1293 !> \param a ...
1294 !> \param b ...
1295 !> \return ...
1296 !> \author Joost VandeVondele
1297 ! **************************************************************************************************
1298  ELEMENTAL FUNCTION gcd(a, b)
1299  INTEGER, INTENT(IN) :: a, b
1300  INTEGER :: gcd
1301 
1302  INTEGER :: aa, ab, l, rem, s
1303 
1304  aa = abs(a)
1305  ab = abs(b)
1306  IF (aa < ab) THEN
1307  s = aa
1308  l = ab
1309  ELSE
1310  s = ab
1311  l = aa
1312  END IF
1313  IF (s .NE. 0) THEN
1314  DO
1315  rem = mod(l, s)
1316  IF (rem == 0) EXIT
1317  l = s
1318  s = rem
1319  END DO
1320  gcd = s
1321  ELSE
1322  gcd = l
1323  END IF
1324  END FUNCTION gcd
1325 
1326 ! **************************************************************************************************
1327 !> \brief computes the least common multiplier of two numbers
1328 !> \param a ...
1329 !> \param b ...
1330 !> \return ...
1331 !> \author Joost VandeVondele
1332 ! **************************************************************************************************
1333  ELEMENTAL FUNCTION lcm(a, b)
1334  INTEGER, INTENT(IN) :: a, b
1335  INTEGER :: lcm
1336 
1337  INTEGER :: tmp
1338 
1339  tmp = gcd(a, b)
1340  IF (tmp == 0) THEN
1341  lcm = 0
1342  ELSE
1343  ! could still overflow if the true lcm is larger than maxint
1344  lcm = abs((a/tmp)*b)
1345  END IF
1346  END FUNCTION lcm
1347 
1348 ! **************************************************************************************************
1349 !> \brief computes the exponential integral
1350 !> Ei(x) = Int(exp(-x*t)/t,t=1..infinity) x>0
1351 !> \param x ...
1352 !> \return ...
1353 !> \author JGH (adapted from Numerical recipies)
1354 ! **************************************************************************************************
1355  FUNCTION ei(x)
1356  REAL(dp) :: x, ei
1357 
1358  INTEGER, PARAMETER :: maxit = 100
1359  REAL(dp), PARAMETER :: eps = epsilon(0.0_dp), &
1360  fpmin = tiny(0.0_dp)
1361 
1362  INTEGER :: k
1363  REAL(dp) :: fact, prev, sum1, term
1364 
1365  IF (x <= 0._dp) THEN
1366  cpabort("Invalid argument")
1367  END IF
1368 
1369  IF (x < fpmin) THEN
1370  ei = log(x) + euler
1371  ELSE IF (x <= -log(eps)) THEN
1372  sum1 = 0._dp
1373  fact = 1._dp
1374  DO k = 1, maxit
1375  fact = fact*x/real(k, dp)
1376  term = fact/real(k, dp)
1377  sum1 = sum1 + term
1378  IF (term < eps*sum1) EXIT
1379  END DO
1380  ei = sum1 + log(x) + euler
1381  ELSE
1382  sum1 = 0._dp
1383  term = 1._dp
1384  DO k = 1, maxit
1385  prev = term
1386  term = term*real(k, dp)/x
1387  IF (term < eps) EXIT
1388  IF (term < prev) THEN
1389  sum1 = sum1 + term
1390  ELSE
1391  sum1 = sum1 - prev
1392  EXIT
1393  END IF
1394  END DO
1395  ei = exp(x)*(1._dp + sum1)/x
1396  END IF
1397 
1398  END FUNCTION ei
1399 
1400 ! **************************************************************************************************
1401 !> \brief computes the exponential integral
1402 !> En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,1,..
1403 !> Note: Ei(-x) = -E1(x)
1404 !> \param n ...
1405 !> \param x ...
1406 !> \return ...
1407 !> \par History
1408 !> 05.2007 Created
1409 !> \author Manuel Guidon (adapted from Numerical recipies)
1410 ! **************************************************************************************************
1411  ELEMENTAL IMPURE FUNCTION expint(n, x)
1412  INTEGER, INTENT(IN) :: n
1413  REAL(dp), INTENT(IN) :: x
1414  REAL(dp) :: expint
1415 
1416  INTEGER, PARAMETER :: maxit = 100
1417  REAL(dp), PARAMETER :: eps = 6.e-14_dp, euler = 0.5772156649015328606065120_dp, &
1418  fpmin = tiny(0.0_dp)
1419 
1420  INTEGER :: i, ii, nm1
1421  REAL(dp) :: a, b, c, d, del, fact, h, psi
1422 
1423  nm1 = n - 1
1424 
1425  IF (n .LT. 0 .OR. x .LT. 0.0_dp .OR. (x .EQ. 0.0_dp .AND. (n .EQ. 0 .OR. n .EQ. 1))) THEN
1426  cpabort("Invalid argument")
1427  ELSE IF (n .EQ. 0) THEN !Special case.
1428  expint = exp(-x)/x
1429  ELSE IF (x .EQ. 0.0_dp) THEN !Another special case.
1430  expint = 1.0_dp/nm1
1431  ELSE IF (x .GT. 1.0_dp) THEN !Lentz's algorithm (5.2).
1432  b = x + n
1433  c = 1.0_dp/fpmin
1434  d = 1.0_dp/b
1435  h = d
1436  DO i = 1, maxit
1437  a = -i*(nm1 + i)
1438  b = b + 2.0_dp
1439  d = 1.0_dp/(a*d + b)
1440  c = b + a/c
1441  del = c*d
1442  h = h*del
1443  IF (abs(del - 1.0_dp) .LT. eps) THEN
1444  expint = h*exp(-x)
1445  RETURN
1446  END IF
1447  END DO
1448  cpabort("continued fraction failed in expint")
1449  ELSE !Evaluate series.
1450  IF (nm1 .NE. 0) THEN !Set first term.
1451  expint = 1.0_dp/nm1
1452  ELSE
1453  expint = -log(x) - euler
1454  END IF
1455  fact = 1.0_dp
1456  DO i = 1, maxit
1457  fact = -fact*x/i
1458  IF (i .NE. nm1) THEN
1459  del = -fact/(i - nm1)
1460  ELSE
1461  psi = -euler !Compute I(n).
1462  DO ii = 1, nm1
1463  psi = psi + 1.0_dp/ii
1464  END DO
1465  del = fact*(-log(x) + psi)
1466  END IF
1467  expint = expint + del
1468  IF (abs(del) .LT. abs(expint)*eps) RETURN
1469  END DO
1470  cpabort("series failed in expint")
1471  END IF
1472  RETURN
1473  END FUNCTION expint
1474 
1475 ! **************************************************************************************************
1476 !> \brief Jacobi matrix diagonalization. The eigenvalues are returned in
1477 !> vector d and the eigenvectors are returned in matrix v in ascending
1478 !> order.
1479 !>
1480 !> \param a ...
1481 !> \param d ...
1482 !> \param v ...
1483 !> \par History
1484 !> - Creation (20.11.98, Matthias Krack)
1485 ! **************************************************************************************************
1486  SUBROUTINE jacobi(a, d, v)
1487  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1488  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: d
1489  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: v
1490 
1491  INTEGER :: n
1492 
1493  n = SIZE(d(:))
1494 
1495  ! Diagonalize matrix a
1496  CALL diag(n, a, d, v)
1497 
1498  ! Sort eigenvalues and eigenvector in ascending order
1499  CALL eigsrt(n, d, v)
1500 
1501  END SUBROUTINE jacobi
1502 
1503 ! **************************************************************************************************
1504 !> \brief Diagonalize matrix a. The eigenvalues are returned in vector d
1505 !> and the eigenvectors are returned in matrix v.
1506 !>
1507 !> \param n matrix/vector extent (problem size)
1508 !> \param a matrix to be diagonalised
1509 !> \param d vector of eigenvalues
1510 !> \param v matrix of eigenvectors
1511 !> \par History
1512 !> - Creation (20.11.98, Matthias Krack)
1513 ! **************************************************************************************************
1514  SUBROUTINE diag(n, a, d, v)
1515 
1516  INTEGER, INTENT(IN) :: n
1517  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1518  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: d
1519  REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: v
1520 
1521  REAL(kind=dp), PARAMETER :: a_eps = 1.0e-10_dp, d_eps = 1.0e-3_dp
1522 
1523  INTEGER :: i, ip, iq
1524  REAL(kind=dp) :: a_max, apq, c, d_min, dip, diq, g, h, s, &
1525  t, tau, theta, tresh
1526  REAL(kind=dp), DIMENSION(n) :: b, z
1527 
1528  a_max = 0.0_dp
1529  DO ip = 1, n - 1
1530  a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1531  b(ip) = a(ip, ip) ! get_diag(a)
1532  END DO
1533  b(n) = a(n, n)
1534 
1535  CALL unit_matrix(v)
1536 
1537  ! Go for 50 iterations
1538  DO i = 1, 50
1539  d = b
1540  d_min = max(d_eps, minval(abs(b)))
1541  IF (a_max < a_eps*d_min) RETURN
1542  tresh = merge(a_max, 0.0_dp, (i < 4))
1543  z = 0.0_dp
1544  DO ip = 1, n - 1
1545  DO iq = ip + 1, n
1546  dip = d(ip)
1547  diq = d(iq)
1548  apq = a(ip, iq)
1549  g = 100.0_dp*abs(apq)
1550  IF (tresh < abs(apq)) THEN
1551  h = diq - dip
1552  IF ((abs(h) + g) .NE. abs(h)) THEN
1553  theta = 0.5_dp*h/apq
1554  t = 1.0_dp/(abs(theta) + sqrt(1.0_dp + theta**2))
1555  IF (theta < 0.0_dp) t = -t
1556  ELSE
1557  t = apq/h
1558  END IF
1559  c = 1.0_dp/sqrt(1.0_dp + t**2)
1560  s = t*c
1561  tau = s/(1.0_dp + c)
1562  h = t*apq
1563  z(ip) = z(ip) - h
1564  z(iq) = z(iq) + h
1565  d(ip) = dip - h
1566  d(iq) = diq + h
1567  a(ip, iq) = 0.0_dp
1568  CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
1569  CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
1570  CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
1571  CALL jrotate(v(:, ip), v(:, iq), s, tau)
1572  ELSE IF ((4 < i) .AND. &
1573  ((abs(dip) + g) == abs(dip)) .AND. &
1574  ((abs(diq) + g) == abs(diq))) THEN
1575  a(ip, iq) = 0.0_dp
1576  END IF
1577  END DO
1578  END DO
1579  b = b + z
1580  a_max = 0.0_dp
1581  DO ip = 1, n - 1
1582  a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1583  END DO
1584  END DO
1585  WRITE (*, '(/,T2,A,/)') 'Too many iterations in jacobi'
1586 
1587  END SUBROUTINE diag
1588 
1589 ! **************************************************************************************************
1590 !> \brief Perform a Jacobi rotation of the vectors a and b.
1591 !>
1592 !> \param a ...
1593 !> \param b ...
1594 !> \param ss ...
1595 !> \param tt ...
1596 !> \par History
1597 !> - Creation (20.11.98, Matthias Krack)
1598 ! **************************************************************************************************
1599  PURE SUBROUTINE jrotate(a, b, ss, tt)
1600  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1601  REAL(kind=dp), INTENT(IN) :: ss, tt
1602 
1603  REAL(kind=dp) :: u, v
1604 
1605  u = 1.0_dp - ss*tt
1606  v = ss/u
1607 
1608  a = a*u - b*ss
1609  b = b*(u + ss*v) + a*v
1610 
1611  END SUBROUTINE jrotate
1612 
1613 ! **************************************************************************************************
1614 !> \brief Sort the values in vector d in ascending order and swap the
1615 !> corresponding columns of matrix v.
1616 !>
1617 !> \param n ...
1618 !> \param d ...
1619 !> \param v ...
1620 !> \par History
1621 !> - Creation (20.11.98, Matthias Krack)
1622 ! **************************************************************************************************
1623  SUBROUTINE eigsrt(n, d, v)
1624  INTEGER, INTENT(IN) :: n
1625  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: d
1626  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: v
1627 
1628  INTEGER :: i, j
1629 
1630  DO i = 1, n - 1
1631  j = sum(minloc(d(i:n))) + i - 1
1632  IF (j /= i) THEN
1633  CALL swap(d(i), d(j))
1634  CALL swap(v(:, i), v(:, j))
1635  END IF
1636  END DO
1637 
1638  END SUBROUTINE eigsrt
1639 
1640 ! **************************************************************************
1641 !> \brief Swap two scalars
1642 !>
1643 !> \param a ...
1644 !> \param b ...
1645 !> \par History
1646 !> - Creation (20.11.98, Matthias Krack)
1647 ! **************************************************************************************************
1648  ELEMENTAL SUBROUTINE swap_scalar(a, b)
1649  REAL(kind=dp), INTENT(INOUT) :: a, b
1650 
1651  REAL(kind=dp) :: c
1652 
1653  c = a
1654  a = b
1655  b = c
1656 
1657  END SUBROUTINE swap_scalar
1658 
1659 ! **************************************************************************
1660 !> \brief Swap two vectors
1661 !>
1662 !> \param a ...
1663 !> \param b ...
1664 !> \par History
1665 !> - Creation (20.11.98, Matthias Krack)
1666 ! **************************************************************************************************
1667  SUBROUTINE swap_vector(a, b)
1668  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1669 
1670  INTEGER :: i, n
1671  REAL(kind=dp) :: c
1672 
1673  n = SIZE(a)
1674 
1675  IF (n /= SIZE(b)) THEN
1676  cpabort("Check the array bounds of the parameters")
1677  END IF
1678 
1679  DO i = 1, n
1680  c = a(i)
1681  a(i) = b(i)
1682  b(i) = c
1683  END DO
1684 
1685  END SUBROUTINE swap_vector
1686 
1687 ! **************************************************************************************************
1688 !> \brief - compute a truncation radius for the shortrange operator
1689 !> \param eps target accuracy!> \param omg screening parameter
1690 !> \param omg ...
1691 !> \param r_cutoff cutoff radius
1692 !> \par History
1693 !> 10.2012 created [Hossein Banihashemian]
1694 !> 05.2019 moved here from hfx_types (A. Bussy)
1695 !> \author Hossein Banihashemian
1696 ! **************************************************************************************************
1697  SUBROUTINE erfc_cutoff(eps, omg, r_cutoff)
1698  IMPLICIT NONE
1699 
1700  REAL(dp), INTENT(in) :: eps, omg
1701  REAL(dp), INTENT(out) :: r_cutoff
1702 
1703  CHARACTER(LEN=*), PARAMETER :: routinen = 'erfc_cutoff'
1704 
1705  REAL(dp), PARAMETER :: abstol = 1e-10_dp, soltol = 1e-16_dp
1706  REAL(dp) :: r0, f0, fprime0, delta_r
1707  INTEGER :: iter, handle
1708  INTEGER, PARAMETER :: itermax = 1000
1709 
1710  CALL timeset(routinen, handle)
1711 
1712  ! initial guess assuming that we are in the asymptotic regime of the erf, and the solution is about 10.
1713  r0 = sqrt(-log(eps*omg*10**2))/omg
1714  CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1715 
1716  DO iter = 1, itermax
1717  delta_r = f0/fprime0
1718  r0 = r0 - delta_r
1719  CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1720  IF (abs(delta_r) .LT. abstol .OR. abs(f0) .LT. soltol) EXIT
1721  END DO
1722  cpassert(iter <= itermax)
1723  r_cutoff = r0
1724 
1725  CALL timestop(handle)
1726  CONTAINS
1727 ! **************************************************************************************************
1728 !> \brief ...
1729 !> \param r ...
1730 !> \param eps ...
1731 !> \param omega ...
1732 !> \param fn ...
1733 !> \param df ...
1734 ! **************************************************************************************************
1735  ELEMENTAL SUBROUTINE eval_transc_func(r, eps, omega, fn, df)
1736  REAL(dp), INTENT(in) :: r, eps, omega
1737  REAL(dp), INTENT(out) :: fn, df
1738 
1739  fn = erfc(omega*r) - r*eps
1740  df = -omega*2*exp(-omega**2*r**2)/sqrt(pi) - eps
1741  END SUBROUTINE eval_transc_func
1742  END SUBROUTINE erfc_cutoff
1743 
1744 ! **************************************************************************************************
1745 !> \brief Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd
1746 !> \param matrix ...
1747 !> \param eigenvectors ...
1748 !> \param eigenvalues ...
1749 !> \author A. Bussy
1750 ! **************************************************************************************************
1751  SUBROUTINE complex_diag(matrix, eigenvectors, eigenvalues)
1752 
1753  COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix, eigenvectors
1754  REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1755 
1756  COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1757  INTEGER :: info, liwork, lrwork, lwork, n
1758  INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
1759  REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: rwork
1760 
1761  eigenvectors(:, :) = matrix(:, :)
1762 
1763  n = SIZE(matrix, 1)
1764  ALLOCATE (iwork(1), rwork(1), work(1))
1765 
1766  ! work space query
1767  lwork = -1
1768  lrwork = -1
1769  liwork = -1
1770 
1771  CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1772 
1773  lwork = ceiling(real(work(1), kind=dp))
1774  lrwork = ceiling(real(rwork(1), kind=dp))
1775  liwork = iwork(1)
1776 
1777  DEALLOCATE (iwork, rwork, work)
1778 
1779  ALLOCATE (iwork(liwork))
1780  iwork(:) = 0
1781  ALLOCATE (rwork(lrwork))
1782  rwork(:) = 0.0_dp
1783  ALLOCATE (work(lwork))
1784  work(:) = cmplx(0.0_dp, 0.0_dp, kind=dp)
1785 
1786  !diagonalization proper
1787  CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1788 
1789  IF (info /= 0) &
1790  cpabort("Diagonalisation of a complex matrix failed")
1791  DEALLOCATE (iwork, rwork, work)
1792 
1793  END SUBROUTINE complex_diag
1794 
1795 END MODULE mathlib
real(kind=dp) function det_3x3(a)
...
Definition: dumpdcd.F:1091
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public euler
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition: mathconstants.F:37
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
subroutine, public get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
returns the pseudoinverse of a real, square matrix using singular value decomposition
Definition: mathlib.F:954
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
Definition: mathlib.F:1487
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
Definition: mathlib.F:206
elemental integer function, public lcm(a, b)
computes the least common multiplier of two numbers
Definition: mathlib.F:1334
subroutine, public get_pseudo_inverse_diag(a, a_pinverse, rskip)
returns the pseudoinverse of a real, symmetric and positive definite matrix using diagonalization.
Definition: mathlib.F:1036
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
Definition: mathlib.F:579
elemental real(kind=dp) function, public binomial_gen(z, k)
The generalized binomial coefficient z over k for 0 <= k <= n is calculated. (z) z*(z-1)*....
Definition: mathlib.F:231
pure real(kind=dp) function, dimension(min(size(a, 1), size(a, 2))), public get_diag(a)
Return the diagonal elements of matrix a as a vector.
Definition: mathlib.F:497
pure real(kind=dp) function, dimension(3), public reflect_vector(a, b)
Reflection of the vector a through a mirror plane defined by the normal vector b. The reflected vecto...
Definition: mathlib.F:1104
pure real(kind=dp) function, public angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Definition: mathlib.F:175
pure subroutine, public build_rotmat(phi, a, rotmat)
The rotation matrix rotmat which rotates a vector about a rotation axis defined by the vector a is bu...
Definition: mathlib.F:289
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition: mathlib.F:1299
subroutine, public symmetrize_matrix(a, option)
Symmetrize the matrix a.
Definition: mathlib.F:1216
subroutine, public complex_diag(matrix, eigenvectors, eigenvalues)
Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd.
Definition: mathlib.F:1752
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition: mathlib.F:520
subroutine, public invmat(a, info)
returns inverse of matrix using the lapack routines DGETRF and DGETRI
Definition: mathlib.F:548
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition: mathlib.F:376
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
Definition: mathlib.F:1515
pure real(kind=dp) function, dimension(3), public rotate_vector(a, phi, b)
Rotation of the vector a about an rotation axis defined by the vector b. The rotation angle is phi (r...
Definition: mathlib.F:1143
subroutine, public erfc_cutoff(eps, omg, r_cutoff)
compute a truncation radius for the shortrange operator
Definition: mathlib.F:1698
logical function, public abnormal_value(a)
determines if a value is not normal (e.g. for Inf and Nan) based on IO to work also under optimizatio...
Definition: mathlib.F:150
elemental impure real(dp) function, public expint(n, x)
computes the exponential integral En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,...
Definition: mathlib.F:1412
pure real(kind=dp) function, public dihedral_angle(ab, bc, cd)
Returns the dihedral angle, i.e. the angle between the planes defined by the vectors (-ab,...
Definition: mathlib.F:472
real(kind=dp) function, public pswitch(x, a, b, order)
Polynomial (5th degree) switching function f(a) = 1 .... f(b) = 0 with f'(a) = f"(a) = f'(b) = f"(b) ...
Definition: mathlib.F:99
pure real(kind=dp) function, dimension(3), public vector_product(a, b)
Calculation of the vector product c = a x b.
Definition: mathlib.F:1282
pure real(kind=dp) function, public multinomial(n, k)
Calculates the multinomial coefficients.
Definition: mathlib.F:257