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