Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : z_one,&
23 : z_zero
24 : #include "../base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mathlib'
31 : REAL(KIND=dp), PARAMETER :: eps_geo = 1.0E-6_dp
32 :
33 : ! Public subroutines
34 :
35 : PUBLIC :: build_rotmat, &
36 : jacobi, &
37 : diamat_all, &
38 : invmat, &
39 : invmat_symm, &
40 : invert_matrix, &
41 : get_pseudo_inverse_svd, &
42 : get_pseudo_inverse_diag, &
43 : symmetrize_matrix, &
44 : unit_matrix, diag, &
45 : erfc_cutoff, &
46 : diag_antisym, &
47 : diag_complex, &
48 : gemm_square, &
49 : geeig_right
50 :
51 : ! Public functions
52 :
53 : PUBLIC :: angle, &
54 : binomial, &
55 : binomial_gen, &
56 : multinomial, &
57 : det_3x3, &
58 : dihedral_angle, &
59 : gcd, &
60 : inv_3x3, &
61 : lcm, &
62 : vector_product, &
63 : pswitch, &
64 : rotate_vector, &
65 : reflect_vector, &
66 : expint, abnormal_value, &
67 : get_diag, &
68 : set_diag
69 :
70 : INTERFACE det_3x3
71 : MODULE PROCEDURE det_3x3_1, det_3x3_2
72 : END INTERFACE
73 :
74 : INTERFACE invert_matrix
75 : MODULE PROCEDURE invert_matrix_d, invert_matrix_z
76 : END INTERFACE
77 :
78 : INTERFACE set_diag
79 : MODULE PROCEDURE set_diag_scalar_d, set_diag_scalar_z
80 : END INTERFACE
81 :
82 : INTERFACE swap
83 : MODULE PROCEDURE swap_scalar, swap_vector
84 : END INTERFACE
85 :
86 : INTERFACE unit_matrix
87 : MODULE PROCEDURE unit_matrix_d, unit_matrix_z
88 : END INTERFACE
89 :
90 : INTERFACE gemm_square
91 : MODULE PROCEDURE zgemm_square_2, zgemm_square_3, dgemm_square_2, dgemm_square_3
92 : END INTERFACE
93 :
94 : CONTAINS
95 :
96 : ! **************************************************************************************************
97 : !> \brief Polynomial (5th degree) switching function
98 : !> f(a) = 1 .... f(b) = 0 with f'(a) = f"(a) = f'(b) = f"(b) = 0
99 : !> \param x ...
100 : !> \param a ...
101 : !> \param b ...
102 : !> \param order ...
103 : !> \return =0 : f(x)
104 : !> \return =1 : f'(x)
105 : !> \return =2 : f"(x)
106 : ! **************************************************************************************************
107 16030 : FUNCTION pswitch(x, a, b, order) RESULT(fx)
108 : REAL(KIND=dp) :: x, a, b
109 : INTEGER :: order
110 : REAL(KIND=dp) :: fx
111 :
112 : REAL(KIND=dp) :: u, u2, u3
113 :
114 16030 : CPASSERT(b > a)
115 16030 : IF (x < a .OR. x > b) THEN
116 : ! outside switching intervall
117 14990 : IF (order > 0) THEN
118 : ! derivatives are 0
119 : fx = 0.0_dp
120 : ELSE
121 7495 : IF (x < a) THEN
122 : ! x < a => f(x) = 1
123 : fx = 1.0_dp
124 : ELSE
125 : ! x > b => f(x) = 0
126 7288 : fx = 0.0_dp
127 : END IF
128 : END IF
129 : ELSE
130 : ! renormalized coordinate
131 1040 : u = (x - a)/(b - a)
132 1560 : SELECT CASE (order)
133 : CASE (0)
134 520 : u2 = u*u
135 520 : u3 = u2*u
136 520 : fx = 1._dp - 10._dp*u3 + 15._dp*u2*u2 - 6._dp*u2*u3
137 : CASE (1)
138 520 : u2 = u*u
139 520 : fx = -30._dp*u2 + 60._dp*u*u2 - 30._dp*u2*u2
140 520 : fx = fx/(b - a)
141 : CASE (2)
142 0 : u2 = u*u
143 0 : fx = -60._dp*u + 180._dp*u2 - 120._dp*u*u2
144 0 : fx = fx/(b - a)**2
145 : CASE DEFAULT
146 1040 : CPABORT('order not defined')
147 : END SELECT
148 : END IF
149 :
150 16030 : END FUNCTION pswitch
151 :
152 : ! **************************************************************************************************
153 : !> \brief determines if a value is not normal (e.g. for Inf and Nan)
154 : !> based on IO to work also under optimization.
155 : !> \param a input value
156 : !> \return TRUE for NaN and Inf
157 : ! **************************************************************************************************
158 347436 : LOGICAL FUNCTION abnormal_value(a)
159 : REAL(KIND=dp) :: a
160 :
161 : CHARACTER(LEN=32) :: buffer
162 :
163 347436 : abnormal_value = .FALSE.
164 : ! the function should work when compiled with -ffast-math and similar
165 : ! unfortunately, that option asserts that all numbers are normals,
166 : ! which the compiler uses to optimize the function to .FALSE. if based on the IEEE module
167 : ! therefore, pass this to the Fortran runtime/printf, if things are NaN or Inf, error out.
168 347436 : WRITE (buffer, *) a
169 347436 : IF (INDEX(buffer, "N") /= 0 .OR. INDEX(buffer, "n") /= 0) abnormal_value = .TRUE.
170 :
171 347436 : END FUNCTION abnormal_value
172 :
173 : ! **************************************************************************************************
174 : !> \brief Calculation of the angle between the vectors a and b.
175 : !> The angle is returned in radians.
176 : !> \param a ...
177 : !> \param b ...
178 : !> \return ...
179 : !> \date 14.10.1998
180 : !> \author MK
181 : !> \version 1.0
182 : ! **************************************************************************************************
183 625925 : PURE FUNCTION angle(a, b) RESULT(angle_ab)
184 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, b
185 : REAL(KIND=dp) :: angle_ab
186 :
187 : REAL(KIND=dp) :: length_of_a, length_of_b
188 625925 : REAL(KIND=dp), DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
189 :
190 2503700 : length_of_a = SQRT(DOT_PRODUCT(a, a))
191 2503700 : length_of_b = SQRT(DOT_PRODUCT(b, b))
192 :
193 625925 : IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
194 2503700 : a_norm(:) = a(:)/length_of_a
195 2503700 : b_norm(:) = b(:)/length_of_b
196 2503700 : angle_ab = ACOS(MIN(MAX(DOT_PRODUCT(a_norm, b_norm), -1.0_dp), 1.0_dp))
197 : ELSE
198 : angle_ab = 0.0_dp
199 : END IF
200 :
201 625925 : END FUNCTION angle
202 :
203 : ! **************************************************************************************************
204 : !> \brief The binomial coefficient n over k for 0 <= k <= n is calculated,
205 : !> otherwise zero is returned.
206 : !> \param n ...
207 : !> \param k ...
208 : !> \return ...
209 : !> \date 08.03.1999
210 : !> \author MK
211 : !> \version 1.0
212 : ! **************************************************************************************************
213 343928383 : ELEMENTAL FUNCTION binomial(n, k) RESULT(n_over_k)
214 : INTEGER, INTENT(IN) :: n, k
215 : REAL(KIND=dp) :: n_over_k
216 :
217 343928383 : IF ((k >= 0) .AND. (k <= n)) THEN
218 341136369 : n_over_k = fac(n)/(fac(n - k)*fac(k))
219 : ELSE
220 : n_over_k = 0.0_dp
221 : END IF
222 :
223 343928383 : END FUNCTION binomial
224 :
225 : ! **************************************************************************************************
226 : !> \brief The generalized binomial coefficient z over k for 0 <= k <= n is calculated.
227 : !> (z) z*(z-1)*...*(z-k+2)*(z-k+1)
228 : !> ( ) = ---------------------------
229 : !> (k) k!
230 : !> \param z ...
231 : !> \param k ...
232 : !> \return ...
233 : !> \date 11.11.2019
234 : !> \author FS
235 : !> \version 1.0
236 : ! **************************************************************************************************
237 171640 : ELEMENTAL FUNCTION binomial_gen(z, k) RESULT(z_over_k)
238 : REAL(KIND=dp), INTENT(IN) :: z
239 : INTEGER, INTENT(IN) :: k
240 : REAL(KIND=dp) :: z_over_k
241 :
242 : INTEGER :: i
243 :
244 171640 : IF (k >= 0) THEN
245 : z_over_k = 1.0_dp
246 257460 : DO i = 1, k
247 257460 : z_over_k = z_over_k*(z - i + 1)/REAL(i, dp)
248 : END DO
249 : ELSE
250 : z_over_k = 0.0_dp
251 : END IF
252 :
253 171640 : END FUNCTION binomial_gen
254 :
255 : ! **************************************************************************************************
256 : !> \brief Calculates the multinomial coefficients
257 : !> \param n ...
258 : !> \param k ...
259 : !> \return ...
260 : !> \author Ole Schuett
261 : ! **************************************************************************************************
262 8982 : PURE FUNCTION multinomial(n, k) RESULT(res)
263 : INTEGER, INTENT(IN) :: n
264 : INTEGER, DIMENSION(:), INTENT(IN) :: k
265 : REAL(KIND=dp) :: res
266 :
267 : INTEGER :: i
268 : REAL(KIND=dp) :: denom
269 :
270 71856 : IF (ALL(k >= 0) .AND. SUM(k) == n) THEN
271 5280 : denom = 1.0_dp
272 21120 : DO i = 1, SIZE(k)
273 21120 : denom = denom*fac(k(i))
274 : END DO
275 5280 : res = fac(n)/denom
276 : ELSE
277 : res = 0.0_dp
278 : END IF
279 :
280 8982 : END FUNCTION multinomial
281 :
282 : ! **************************************************************************************************
283 : !> \brief The rotation matrix rotmat which rotates a vector about a
284 : !> rotation axis defined by the vector a is build up.
285 : !> The rotation angle is phi (radians).
286 : !> \param phi ...
287 : !> \param a ...
288 : !> \param rotmat ...
289 : !> \date 16.10.1998
290 : !> \author MK
291 : !> \version 1.0
292 : ! **************************************************************************************************
293 15759 : PURE SUBROUTINE build_rotmat(phi, a, rotmat)
294 : REAL(KIND=dp), INTENT(IN) :: phi
295 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a
296 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: rotmat
297 :
298 : REAL(KIND=dp) :: cosp, cost, length_of_a, sinp
299 : REAL(KIND=dp), DIMENSION(3) :: d
300 :
301 15759 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
302 : ! Check the length of the vector a
303 15759 : IF (length_of_a > eps_geo) THEN
304 :
305 63036 : d(:) = a(:)/length_of_a
306 :
307 15759 : cosp = COS(phi)
308 15759 : sinp = SIN(phi)
309 15759 : cost = 1.0_dp - cosp
310 :
311 15759 : rotmat(1, 1) = d(1)*d(1)*cost + cosp
312 15759 : rotmat(1, 2) = d(1)*d(2)*cost - d(3)*sinp
313 15759 : rotmat(1, 3) = d(1)*d(3)*cost + d(2)*sinp
314 15759 : rotmat(2, 1) = d(2)*d(1)*cost + d(3)*sinp
315 15759 : rotmat(2, 2) = d(2)*d(2)*cost + cosp
316 15759 : rotmat(2, 3) = d(2)*d(3)*cost - d(1)*sinp
317 15759 : rotmat(3, 1) = d(3)*d(1)*cost - d(2)*sinp
318 15759 : rotmat(3, 2) = d(3)*d(2)*cost + d(1)*sinp
319 15759 : rotmat(3, 3) = d(3)*d(3)*cost + cosp
320 : ELSE
321 0 : CALL unit_matrix(rotmat)
322 : END IF
323 :
324 15759 : END SUBROUTINE build_rotmat
325 :
326 : ! **************************************************************************************************
327 : !> \brief Returns the determinante of the 3x3 matrix a.
328 : !> \param a ...
329 : !> \return ...
330 : !> \date 13.03.2004
331 : !> \author MK
332 : !> \version 1.0
333 : ! **************************************************************************************************
334 16463106 : PURE FUNCTION det_3x3_1(a) RESULT(det_a)
335 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
336 : REAL(KIND=dp) :: det_a
337 :
338 : det_a = a(1, 1)*(a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2)) + &
339 : a(1, 2)*(a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3)) + &
340 16463106 : a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
341 :
342 16463106 : END FUNCTION det_3x3_1
343 :
344 : ! **************************************************************************************************
345 : !> \brief Returns the determinante of the 3x3 matrix a given by its columns.
346 : !> \param a1 ...
347 : !> \param a2 ...
348 : !> \param a3 ...
349 : !> \return ...
350 : !> \date 13.03.2004
351 : !> \author MK
352 : !> \version 1.0
353 : ! **************************************************************************************************
354 6270 : PURE FUNCTION det_3x3_2(a1, a2, a3) RESULT(det_a)
355 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3
356 : REAL(KIND=dp) :: det_a
357 :
358 : det_a = a1(1)*(a2(2)*a3(3) - a3(2)*a2(3)) + &
359 : a2(1)*(a3(2)*a1(3) - a1(2)*a3(3)) + &
360 6270 : a3(1)*(a1(2)*a2(3) - a2(2)*a1(3))
361 :
362 6270 : END FUNCTION det_3x3_2
363 :
364 : ! **************************************************************************************************
365 : !> \brief Diagonalize the symmetric n by n matrix a using the LAPACK
366 : !> library. Only the upper triangle of matrix a is used.
367 : !> Externals (LAPACK 3.0)
368 : !> \param a ...
369 : !> \param eigval ...
370 : !> \param dac ...
371 : !> \date 29.03.1999
372 : !> \par Variables
373 : !> - a : Symmetric matrix to be diagonalized (input; upper triangle) ->
374 : !> - eigenvectors of the matrix a (output).
375 : !> - dac : If true, then the divide-and-conquer algorithm is applied.
376 : !> - eigval : Eigenvalues of the matrix a (output).
377 : !> \author MK
378 : !> \version 1.0
379 : ! **************************************************************************************************
380 74805 : SUBROUTINE diamat_all(a, eigval, dac)
381 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
382 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigval
383 : LOGICAL, INTENT(IN), OPTIONAL :: dac
384 :
385 : CHARACTER(len=*), PARAMETER :: routineN = 'diamat_all'
386 :
387 : INTEGER :: handle, info, liwork, lwork, n, nb
388 74805 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
389 : INTEGER, EXTERNAL :: ilaenv
390 : LOGICAL :: divide_and_conquer
391 74805 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
392 :
393 : EXTERNAL dsyev, dsyevd
394 :
395 74805 : CALL timeset(routineN, handle)
396 :
397 : ! Get the size of the matrix a
398 74805 : n = SIZE(a, 1)
399 :
400 : ! Check the size of matrix a
401 74805 : IF (SIZE(a, 2) /= n) THEN
402 0 : CPABORT("Check the size of matrix a (parameter #1)")
403 : END IF
404 :
405 : ! Check the size of vector eigval
406 74805 : IF (SIZE(eigval) /= n) THEN
407 0 : CPABORT("The dimension of vector eigval is too small")
408 : END IF
409 :
410 : ! Check, if the divide-and-conquer algorithm is requested
411 :
412 74805 : IF (PRESENT(dac)) THEN
413 205 : divide_and_conquer = dac
414 : ELSE
415 : divide_and_conquer = .FALSE.
416 : END IF
417 :
418 : ! Get the optimal work storage size
419 :
420 205 : IF (divide_and_conquer) THEN
421 205 : lwork = 2*n**2 + 6*n + 1
422 205 : liwork = 5*n + 3
423 : ELSE
424 74600 : nb = ilaenv(1, "DSYTRD", "U", n, -1, -1, -1)
425 74600 : lwork = (nb + 2)*n
426 : END IF
427 :
428 : ! Allocate work storage
429 :
430 224415 : ALLOCATE (work(lwork))
431 74805 : IF (divide_and_conquer) THEN
432 615 : ALLOCATE (iwork(liwork))
433 : END IF
434 :
435 : ! Diagonalize the matrix a
436 :
437 74805 : info = 0
438 : IF (divide_and_conquer) THEN
439 205 : CALL dsyevd("V", "U", n, a, n, eigval, work, lwork, iwork, liwork, info)
440 : ELSE
441 74990 : CALL dsyev("V", "U", n, a, n, eigval, work, lwork, info)
442 : END IF
443 :
444 74805 : IF (info /= 0) THEN
445 0 : IF (divide_and_conquer) THEN
446 0 : CPABORT("The matrix diagonalization with dsyevd failed")
447 : ELSE
448 0 : CPABORT("The matrix diagonalization with dsyev failed")
449 : END IF
450 : END IF
451 :
452 : ! Release work storage
453 74805 : DEALLOCATE (work)
454 :
455 74805 : IF (divide_and_conquer) THEN
456 205 : DEALLOCATE (iwork)
457 : END IF
458 :
459 74805 : CALL timestop(handle)
460 :
461 149610 : END SUBROUTINE diamat_all
462 :
463 : ! **************************************************************************************************
464 : !> \brief Returns the dihedral angle, i.e. the angle between the planes
465 : !> defined by the vectors (-ab,bc) and (cd,-bc).
466 : !> The dihedral angle is returned in radians.
467 : !> \param ab ...
468 : !> \param bc ...
469 : !> \param cd ...
470 : !> \return ...
471 : !> \date 13.03.2004
472 : !> \author MK
473 : !> \version 1.0
474 : ! **************************************************************************************************
475 116 : PURE FUNCTION dihedral_angle(ab, bc, cd) RESULT(dihedral_angle_abcd)
476 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ab, bc, cd
477 : REAL(KIND=dp) :: dihedral_angle_abcd
478 :
479 : REAL(KIND=dp) :: det_abcd
480 : REAL(KIND=dp), DIMENSION(3) :: abc, bcd
481 :
482 464 : abc = vector_product(bc, -ab)
483 464 : bcd = vector_product(cd, -bc)
484 : ! Calculate the normal vectors of the planes
485 : ! defined by the points a,b,c and b,c,d
486 :
487 464 : det_abcd = det_3x3(abc, bcd, -bc)
488 116 : dihedral_angle_abcd = SIGN(1.0_dp, det_abcd)*angle(abc, bcd)
489 :
490 116 : END FUNCTION dihedral_angle
491 :
492 : ! **************************************************************************************************
493 : !> \brief Return the diagonal elements of matrix a as a vector.
494 : !> \param a ...
495 : !> \return ...
496 : !> \date 20.11.1998
497 : !> \author MK
498 : !> \version 1.0
499 : ! **************************************************************************************************
500 55926 : PURE FUNCTION get_diag(a) RESULT(a_diag)
501 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: a
502 : REAL(KIND=dp), &
503 : DIMENSION(MIN(SIZE(a, 1), SIZE(a, 2))) :: a_diag
504 :
505 : INTEGER :: i, n
506 :
507 55926 : n = MIN(SIZE(a, 1), SIZE(a, 2))
508 :
509 5641638 : DO i = 1, n
510 5641638 : a_diag(i) = a(i, i)
511 : END DO
512 :
513 55926 : END FUNCTION get_diag
514 :
515 : ! **************************************************************************************************
516 : !> \brief Returns the inverse of the 3 x 3 matrix a.
517 : !> \param a ...
518 : !> \return ...
519 : !> \date 13.03.2004
520 : !> \author MK
521 : !> \version 1.0
522 : ! **************************************************************************************************
523 15704959 : PURE FUNCTION inv_3x3(a) RESULT(a_inv)
524 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
525 : REAL(KIND=dp), DIMENSION(3, 3) :: a_inv
526 :
527 : REAL(KIND=dp) :: det_a
528 :
529 15704959 : det_a = 1.0_dp/det_3x3(a)
530 :
531 15704959 : a_inv(1, 1) = (a(2, 2)*a(3, 3) - a(3, 2)*a(2, 3))*det_a
532 15704959 : a_inv(2, 1) = (a(2, 3)*a(3, 1) - a(3, 3)*a(2, 1))*det_a
533 15704959 : a_inv(3, 1) = (a(2, 1)*a(3, 2) - a(3, 1)*a(2, 2))*det_a
534 :
535 15704959 : a_inv(1, 2) = (a(1, 3)*a(3, 2) - a(3, 3)*a(1, 2))*det_a
536 15704959 : a_inv(2, 2) = (a(1, 1)*a(3, 3) - a(3, 1)*a(1, 3))*det_a
537 15704959 : a_inv(3, 2) = (a(1, 2)*a(3, 1) - a(3, 2)*a(1, 1))*det_a
538 :
539 15704959 : a_inv(1, 3) = (a(1, 2)*a(2, 3) - a(2, 2)*a(1, 3))*det_a
540 15704959 : a_inv(2, 3) = (a(1, 3)*a(2, 1) - a(2, 3)*a(1, 1))*det_a
541 15704959 : a_inv(3, 3) = (a(1, 1)*a(2, 2) - a(2, 1)*a(1, 2))*det_a
542 :
543 15704959 : END FUNCTION inv_3x3
544 :
545 : ! **************************************************************************************************
546 : !> \brief returns inverse of matrix using the lapack routines DGETRF and DGETRI
547 : !> \param a ...
548 : !> \param info ...
549 : ! **************************************************************************************************
550 6190 : SUBROUTINE invmat(a, info)
551 : REAL(KIND=dp), INTENT(INOUT) :: a(:, :)
552 : INTEGER, INTENT(OUT) :: info
553 :
554 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invmat'
555 :
556 : INTEGER :: handle, lwork, n
557 6190 : INTEGER, ALLOCATABLE :: ipiv(:)
558 6190 : REAL(KIND=dp), ALLOCATABLE :: work(:)
559 :
560 6190 : CALL timeset(routineN, handle)
561 :
562 6190 : n = SIZE(a, 1)
563 6190 : lwork = 20*n
564 18570 : ALLOCATE (ipiv(n))
565 18570 : ALLOCATE (work(lwork))
566 59582 : ipiv = 0
567 1074030 : work = 0._dp
568 6190 : info = 0
569 437640 : CALL dgetrf(n, n, a, n, ipiv, info)
570 6190 : IF (info == 0) THEN
571 437592 : CALL dgetri(n, a, n, ipiv, work, lwork, info)
572 : END IF
573 6190 : DEALLOCATE (ipiv, work)
574 :
575 6190 : CALL timestop(handle)
576 :
577 6190 : END SUBROUTINE invmat
578 :
579 : ! **************************************************************************************************
580 : !> \brief returns inverse of real symmetric, positive definite matrix
581 : !> \param a matrix
582 : !> \param potrf if cholesky decomposition of a was already done using dpotrf.
583 : !> If not given, cholesky decomposition of a will be done before inversion.
584 : !> \param uplo indicating if the upper or lower triangle of a is stored.
585 : !> \author Dorothea Golze [02.2015]
586 : ! **************************************************************************************************
587 2430 : SUBROUTINE invmat_symm(a, potrf, uplo)
588 : REAL(KIND=dp), INTENT(INOUT) :: a(:, :)
589 : LOGICAL, INTENT(IN), OPTIONAL :: potrf
590 : CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: uplo
591 :
592 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invmat_symm'
593 :
594 : CHARACTER(LEN=1) :: myuplo
595 : INTEGER :: handle, info, n
596 : LOGICAL :: do_potrf
597 :
598 2430 : CALL timeset(routineN, handle)
599 :
600 2430 : do_potrf = .TRUE.
601 2430 : IF (PRESENT(potrf)) do_potrf = potrf
602 :
603 2430 : myuplo = 'U'
604 2430 : IF (PRESENT(uplo)) myuplo = uplo
605 :
606 2430 : n = SIZE(a, 1)
607 2430 : info = 0
608 :
609 : ! do cholesky decomposition
610 2430 : IF (do_potrf) THEN
611 3413 : CALL dpotrf(myuplo, n, a, n, info)
612 2213 : IF (info /= 0) CPABORT("DPOTRF failed")
613 : END IF
614 :
615 : ! do inversion using the cholesky decomposition
616 3630 : CALL dpotri(myuplo, n, a, n, info)
617 2430 : IF (info /= 0) CPABORT("Matrix inversion failed")
618 :
619 : ! complete the matrix
620 2430 : IF ((myuplo == "U") .OR. (myuplo == "u")) THEN
621 2430 : CALL symmetrize_matrix(a, "upper_to_lower")
622 : ELSE
623 0 : CALL symmetrize_matrix(a, "lower_to_upper")
624 : END IF
625 :
626 2430 : CALL timestop(handle)
627 :
628 2430 : END SUBROUTINE invmat_symm
629 :
630 : ! **************************************************************************************************
631 : !> \brief Compute the inverse of the n by n real matrix a using the LAPACK
632 : !> library
633 : !> \param a ...
634 : !> \param a_inverse ...
635 : !> \param eval_error ...
636 : !> \param option ...
637 : !> \param improve ...
638 : !> \date 23.03.1999
639 : !> \par Variables
640 : !> - a : Real matrix to be inverted (input).
641 : !> - a_inverse: Inverse of the matrix a (output).
642 : !> - a_lu : LU factorization of matrix a.
643 : !> - a_norm : Norm of matrix a.
644 : !> - error : Estimated error of the inversion.
645 : !> - r_cond : Reciprocal condition number of the matrix a.
646 : !> - trans : "N" => invert a
647 : !> - "T" => invert transpose(a)
648 : !> \author MK
649 : !> \version 1.0
650 : !> \note NB add improve argument, used to disable call to dgerfs
651 : ! **************************************************************************************************
652 2258 : SUBROUTINE invert_matrix_d(a, a_inverse, eval_error, option, improve)
653 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: a
654 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: a_inverse
655 : REAL(KIND=dp), INTENT(OUT) :: eval_error
656 : CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: option
657 : LOGICAL, INTENT(IN), OPTIONAL :: improve
658 :
659 : CHARACTER(LEN=1) :: norm, trans
660 : CHARACTER(LEN=default_string_length) :: message
661 : INTEGER :: info, iter, n
662 2258 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv, iwork
663 : LOGICAL :: do_improve
664 : REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
665 2258 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, work
666 2258 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
667 : REAL(KIND=dp), EXTERNAL :: dlange
668 :
669 : EXTERNAL dgecon, dgerfs, dgetrf, dgetrs
670 :
671 : ! Check for optional parameter
672 2258 : IF (PRESENT(option)) THEN
673 498 : trans = option
674 : ELSE
675 1760 : trans = "N"
676 : END IF
677 :
678 2258 : IF (PRESENT(improve)) THEN
679 536 : do_improve = improve
680 : ELSE
681 : do_improve = .TRUE.
682 : END IF
683 :
684 : ! Get the dimension of matrix a
685 2258 : n = SIZE(a, 1)
686 :
687 : ! Check array dimensions
688 2258 : IF (n == 0) THEN
689 0 : CPABORT("Matrix to be inverted of zero size")
690 : END IF
691 :
692 2258 : IF (n /= SIZE(a, 2)) THEN
693 0 : CPABORT("Check the array bounds of parameter #1")
694 : END IF
695 :
696 2258 : IF ((n /= SIZE(a_inverse, 1)) .OR. &
697 : (n /= SIZE(a_inverse, 2))) THEN
698 0 : CPABORT("Check the array bounds of parameter #2")
699 : END IF
700 :
701 : ! Allocate work storage
702 9032 : ALLOCATE (a_lu(n, n))
703 6774 : ALLOCATE (b(n, n))
704 6774 : ALLOCATE (berr(n))
705 4516 : ALLOCATE (ferr(n))
706 6774 : ALLOCATE (ipiv(n))
707 4516 : ALLOCATE (iwork(n))
708 6774 : ALLOCATE (work(4*n))
709 :
710 2461906 : a_lu(1:n, 1:n) = a(1:n, 1:n)
711 :
712 : ! Compute the LU factorization of the matrix a
713 2258 : CALL dgetrf(n, n, a_lu, n, ipiv, info)
714 :
715 2258 : IF (info /= 0) THEN
716 0 : CPABORT("The LU factorization in dgetrf failed")
717 : END IF
718 :
719 : ! Compute the norm of the matrix a
720 :
721 2258 : IF (trans == "N") THEN
722 2008 : norm = '1'
723 : ELSE
724 250 : norm = 'I'
725 : END IF
726 :
727 2258 : a_norm = dlange(norm, n, n, a, n, work)
728 :
729 : ! Compute the reciprocal of the condition number of a
730 :
731 2258 : CALL dgecon(norm, n, a_lu, n, a_norm, r_cond, work, iwork, info)
732 :
733 2258 : IF (info /= 0) THEN
734 0 : CPABORT("The computation of the condition number in dgecon failed")
735 : END IF
736 :
737 2258 : IF (r_cond < EPSILON(0.0_dp)) THEN
738 0 : WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
739 : CALL cp_abort(__LOCATION__, &
740 : "Bad condition number "//TRIM(message)//" (smaller than the machine "// &
741 0 : "working precision)")
742 : END IF
743 :
744 : ! Solve a system of linear equations using the LU factorization computed by dgetrf
745 :
746 2258 : CALL unit_matrix(a_inverse)
747 :
748 2258 : CALL dgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
749 :
750 2258 : IF (info /= 0) THEN
751 0 : CPABORT("Solving the system of linear equations in dgetrs failed")
752 : END IF
753 :
754 : ! Improve the computed solution iteratively
755 2258 : CALL unit_matrix(b) ! Initialize right-hand sides
756 :
757 2258 : eval_error = 0.0_dp
758 :
759 2258 : IF (do_improve) THEN
760 6490 : DO iter = 1, 10
761 :
762 : CALL dgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
763 6252 : work, iwork, info)
764 :
765 6252 : IF (info /= 0) THEN
766 0 : CPABORT("Improving the computed solution in dgerfs failed")
767 : END IF
768 :
769 6252 : old_eval_error = eval_error
770 171918 : eval_error = MAXVAL(ferr)
771 :
772 6490 : IF (ABS(eval_error - old_eval_error) <= EPSILON(1.0_dp)) EXIT
773 :
774 : END DO
775 : END IF
776 :
777 : ! Release work storage
778 2258 : DEALLOCATE (work)
779 2258 : DEALLOCATE (iwork)
780 2258 : DEALLOCATE (ipiv)
781 2258 : DEALLOCATE (ferr)
782 2258 : DEALLOCATE (berr)
783 2258 : DEALLOCATE (b)
784 2258 : DEALLOCATE (a_lu)
785 :
786 2258 : END SUBROUTINE invert_matrix_d
787 :
788 : ! **************************************************************************************************
789 : !> \brief Compute the inverse of the n by n complex matrix a using the LAPACK
790 : !> library
791 : !> \param a ...
792 : !> \param a_inverse ...
793 : !> \param eval_error ...
794 : !> \param option ...
795 : !> \date 08.06.2009
796 : !> \par Variables
797 : !> - a : Complex matrix to be inverted (input).
798 : !> - a_inverse: Inverse of the matrix a (output).
799 : !> - a_lu : LU factorization of matrix a.
800 : !> - a_norm : Norm of matrix a.
801 : !> - error : Estimated error of the inversion.
802 : !> - r_cond : Reciprocal condition number of the matrix a.
803 : !> - trans : "N" => invert a
804 : !> - "T" => invert transpose(a)
805 : !> \author MK
806 : !> \version 1.0
807 : ! **************************************************************************************************
808 0 : SUBROUTINE invert_matrix_z(a, a_inverse, eval_error, option)
809 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: a
810 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: a_inverse
811 : REAL(KIND=dp), INTENT(OUT) :: eval_error
812 : CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: option
813 :
814 : CHARACTER(LEN=1) :: norm, trans
815 : CHARACTER(LEN=default_string_length) :: message
816 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
817 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
818 : INTEGER :: info, iter, n
819 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
820 : REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
821 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, rwork
822 : REAL(KIND=dp), EXTERNAL :: zlange
823 :
824 : EXTERNAL zgecon, zgerfs, zgetrf, zgetrs
825 :
826 : ! Check for optional parameter
827 0 : IF (PRESENT(option)) THEN
828 0 : trans = option
829 : ELSE
830 0 : trans = "N"
831 : END IF
832 :
833 : ! Get the dimension of matrix a
834 0 : n = SIZE(a, 1)
835 :
836 : ! Check array dimensions
837 0 : IF (n == 0) THEN
838 0 : CPABORT("Matrix to be inverted of zero size")
839 : END IF
840 :
841 0 : IF (n /= SIZE(a, 2)) THEN
842 0 : CPABORT("Check the array bounds of parameter #1")
843 : END IF
844 :
845 0 : IF ((n /= SIZE(a_inverse, 1)) .OR. &
846 : (n /= SIZE(a_inverse, 2))) THEN
847 0 : CPABORT("Check the array bounds of parameter #2")
848 : END IF
849 :
850 : ! Allocate work storage
851 0 : ALLOCATE (a_lu(n, n))
852 0 : ALLOCATE (b(n, n))
853 0 : ALLOCATE (berr(n))
854 0 : ALLOCATE (ferr(n))
855 0 : ALLOCATE (ipiv(n))
856 0 : ALLOCATE (rwork(2*n))
857 0 : ALLOCATE (work(2*n))
858 :
859 0 : a_lu(1:n, 1:n) = a(1:n, 1:n)
860 :
861 : ! Compute the LU factorization of the matrix a
862 0 : CALL zgetrf(n, n, a_lu, n, ipiv, info)
863 :
864 0 : IF (info /= 0) THEN
865 0 : CPABORT("The LU factorization in dgetrf failed")
866 : END IF
867 :
868 : ! Compute the norm of the matrix a
869 :
870 0 : IF (trans == "N") THEN
871 0 : norm = '1'
872 : ELSE
873 0 : norm = 'I'
874 : END IF
875 :
876 0 : a_norm = zlange(norm, n, n, a, n, work)
877 :
878 : ! Compute the reciprocal of the condition number of a
879 :
880 0 : CALL zgecon(norm, n, a_lu, n, a_norm, r_cond, work, rwork, info)
881 :
882 0 : IF (info /= 0) THEN
883 0 : CPABORT("The computation of the condition number in dgecon failed")
884 : END IF
885 :
886 0 : IF (r_cond < EPSILON(0.0_dp)) THEN
887 0 : WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
888 : CALL cp_abort(__LOCATION__, &
889 : "Bad condition number "//TRIM(message)//" (smaller than the machine "// &
890 0 : "working precision)")
891 : END IF
892 :
893 : ! Solve a system of linear equations using the LU factorization computed by dgetrf
894 :
895 0 : CALL unit_matrix(a_inverse)
896 :
897 0 : CALL zgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
898 :
899 0 : IF (info /= 0) THEN
900 0 : CPABORT("Solving the system of linear equations in dgetrs failed")
901 : END IF
902 :
903 : ! Improve the computed solution iteratively
904 0 : CALL unit_matrix(b) ! Initialize right-hand sides
905 :
906 0 : eval_error = 0.0_dp
907 :
908 0 : DO iter = 1, 10
909 :
910 : CALL zgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
911 0 : work, rwork, info)
912 :
913 0 : IF (info /= 0) THEN
914 0 : CPABORT("Improving the computed solution in dgerfs failed")
915 : END IF
916 :
917 0 : old_eval_error = eval_error
918 0 : eval_error = MAXVAL(ferr)
919 :
920 0 : IF (ABS(eval_error - old_eval_error) <= EPSILON(1.0_dp)) EXIT
921 :
922 : END DO
923 :
924 : ! Release work storage
925 0 : DEALLOCATE (work)
926 0 : DEALLOCATE (rwork)
927 0 : DEALLOCATE (ipiv)
928 0 : DEALLOCATE (ferr)
929 0 : DEALLOCATE (berr)
930 0 : DEALLOCATE (b)
931 0 : DEALLOCATE (a_lu)
932 :
933 0 : END SUBROUTINE invert_matrix_z
934 :
935 : ! **************************************************************************************************
936 : !> \brief returns the pseudoinverse of a real, square matrix using singular
937 : !> value decomposition
938 : !> \param a matrix a
939 : !> \param a_pinverse pseudoinverse of matrix a
940 : !> \param rskip parameter for setting small singular values to zero
941 : !> \param determinant determinant of matrix a (optional output)
942 : !> \param sval array holding singular values of matrix a (optional output)
943 : !> \author Dorothea Golze [02.2015]
944 : ! **************************************************************************************************
945 535 : SUBROUTINE get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
946 : REAL(KIND=dp), DIMENSION(:, :) :: a, a_pinverse
947 : REAL(KIND=dp), INTENT(IN) :: rskip
948 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: determinant
949 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
950 : OPTIONAL, POINTER :: sval
951 :
952 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_pseudo_inverse_svd'
953 :
954 : INTEGER :: handle, i, info, lwork, n
955 535 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
956 535 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sig, work
957 535 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sig_plus, temp_mat, u, vt
958 :
959 535 : CALL timeset(routineN, handle)
960 :
961 535 : n = SIZE(a, 1)
962 7490 : ALLOCATE (u(n, n), vt(n, n), sig(n), sig_plus(n, n), iwork(8*n), work(1), temp_mat(n, n))
963 48679 : u(:, :) = 0.0_dp
964 48679 : vt(:, :) = 0.0_dp
965 2821 : sig(:) = 0.0_dp
966 48679 : sig_plus = 0.0_dp
967 1070 : work = 0.0_dp
968 18823 : iwork = 0
969 535 : IF (PRESENT(determinant)) determinant = 1.0_dp
970 :
971 : ! work size query
972 535 : lwork = -1
973 : CALL dgesdd('A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
974 535 : lwork, iwork(1), info)
975 :
976 535 : IF (info /= 0) THEN
977 0 : CPABORT("ERROR in DGESDD: Could not retrieve work array sizes")
978 : END IF
979 535 : lwork = INT(work(1))
980 535 : DEALLOCATE (work)
981 1605 : ALLOCATE (work(lwork))
982 :
983 : ! do SVD
984 : CALL dgesdd('A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
985 535 : lwork, iwork(1), info)
986 :
987 535 : IF (info /= 0) THEN
988 0 : CPABORT("SVD failed")
989 : END IF
990 :
991 535 : IF (PRESENT(sval)) THEN
992 0 : CPASSERT(.NOT. ASSOCIATED(sval))
993 0 : ALLOCATE (sval(n))
994 0 : sval(:) = sig
995 : END IF
996 :
997 : ! set singular values that are too small to zero
998 2821 : DO i = 1, n
999 50965 : IF (sig(i) > rskip*MAXVAL(sig)) THEN
1000 2250 : IF (PRESENT(determinant)) &
1001 0 : determinant = determinant*sig(i)
1002 2250 : sig_plus(i, i) = 1._dp/sig(i)
1003 : ELSE
1004 36 : sig_plus(i, i) = 0.0_dp
1005 : END IF
1006 : END DO
1007 :
1008 : ! build pseudoinverse: V*sig_plus*UT
1009 535 : CALL dgemm("N", "T", n, n, n, 1._dp, sig_plus, n, u, n, 0._dp, temp_mat, n)
1010 535 : CALL dgemm("T", "N", n, n, n, 1._dp, vt, n, temp_mat, n, 0._dp, a_pinverse, n)
1011 :
1012 535 : DEALLOCATE (u, vt, sig, iwork, work, sig_plus, temp_mat)
1013 :
1014 535 : CALL timestop(handle)
1015 :
1016 535 : END SUBROUTINE get_pseudo_inverse_svd
1017 :
1018 : ! **************************************************************************************************
1019 : !> \brief returns the pseudoinverse of a real, symmetric and positive definite
1020 : !> matrix using diagonalization.
1021 : !> \param a matrix a
1022 : !> \param a_pinverse pseudoinverse of matrix a
1023 : !> \param rskip parameter for setting small eigenvalues to zero
1024 : !> \author Dorothea Golze [02.2015]
1025 : ! **************************************************************************************************
1026 1161 : SUBROUTINE get_pseudo_inverse_diag(a, a_pinverse, rskip)
1027 : REAL(KIND=dp), DIMENSION(:, :) :: a, a_pinverse
1028 : REAL(KIND=dp), INTENT(IN) :: rskip
1029 :
1030 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_pseudo_inverse_diag'
1031 :
1032 : INTEGER :: handle, i, info, lwork, n
1033 1161 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig, work
1034 1161 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dinv, temp_mat
1035 :
1036 1161 : CALL timeset(routineN, handle)
1037 :
1038 1161 : info = 0
1039 1161 : n = SIZE(a, 1)
1040 9288 : ALLOCATE (dinv(n, n), eig(n), work(1), temp_mat(n, n))
1041 186519 : dinv(:, :) = 0.0_dp
1042 15075 : eig(:) = 0.0_dp
1043 2322 : work(:) = 0.0_dp
1044 186519 : temp_mat = 0.0_dp
1045 :
1046 : ! work size query
1047 1161 : lwork = -1
1048 1161 : CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1049 1161 : IF (info /= 0) THEN
1050 0 : CPABORT("ERROR in DSYEV: Could not retrieve work array sizes")
1051 : END IF
1052 1161 : lwork = INT(work(1))
1053 1161 : DEALLOCATE (work)
1054 3483 : ALLOCATE (work(lwork))
1055 474237 : work = 0.0_dp
1056 :
1057 : ! get eigenvalues and eigenvectors
1058 1161 : CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1059 :
1060 1161 : IF (info /= 0) THEN
1061 0 : CPABORT("Matrix diagonalization failed")
1062 : END IF
1063 :
1064 : ! set eigenvalues that are too small to zero
1065 15075 : DO i = 1, n
1066 200433 : IF (eig(i) > rskip*MAXVAL(eig)) THEN
1067 12454 : dinv(i, i) = 1.0_dp/eig(i)
1068 : ELSE
1069 1460 : dinv(i, i) = 0._dp
1070 : END IF
1071 : END DO
1072 :
1073 : ! build pseudoinverse: U*dinv*UT
1074 1161 : CALL dgemm("N", "T", n, n, n, 1._dp, dinv, n, a, n, 0._dp, temp_mat, n)
1075 1161 : CALL dgemm("N", "N", n, n, n, 1._dp, a, n, temp_mat, n, 0._dp, a_pinverse, n)
1076 :
1077 1161 : DEALLOCATE (eig, work, dinv, temp_mat)
1078 :
1079 1161 : CALL timestop(handle)
1080 :
1081 1161 : END SUBROUTINE get_pseudo_inverse_diag
1082 :
1083 : ! **************************************************************************************************
1084 : !> \brief Reflection of the vector a through a mirror plane defined by the
1085 : !> normal vector b. The reflected vector a is stored in a_mirror.
1086 : !> \param a ...
1087 : !> \param b ...
1088 : !> \return ...
1089 : !> \date 16.10.1998
1090 : !> \author MK
1091 : !> \version 1.0
1092 : ! **************************************************************************************************
1093 19487 : PURE FUNCTION reflect_vector(a, b) RESULT(a_mirror)
1094 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a, b
1095 : REAL(KIND=dp), DIMENSION(3) :: a_mirror
1096 :
1097 : REAL(KIND=dp) :: length_of_b, scapro
1098 : REAL(KIND=dp), DIMENSION(3) :: d
1099 :
1100 19487 : length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1101 :
1102 19487 : IF (length_of_b > eps_geo) THEN
1103 :
1104 77948 : d(:) = b(:)/length_of_b
1105 :
1106 : ! Calculate the mirror image a_mirror of the vector a
1107 19487 : scapro = a(1)*d(1) + a(2)*d(2) + a(3)*d(3)
1108 :
1109 77948 : a_mirror(:) = a(:) - 2.0_dp*scapro*d(:)
1110 :
1111 : ELSE
1112 :
1113 0 : a_mirror(:) = 0.0_dp
1114 :
1115 : END IF
1116 :
1117 19487 : END FUNCTION reflect_vector
1118 :
1119 : ! **************************************************************************************************
1120 : !> \brief Rotation of the vector a about an rotation axis defined by the
1121 : !> vector b. The rotation angle is phi (radians). The rotated vector
1122 : !> a is stored in a_rot.
1123 : !> \param a ...
1124 : !> \param phi ...
1125 : !> \param b ...
1126 : !> \return ...
1127 : !> \date 16.10.1998
1128 : !> \author MK
1129 : !> \version 1.0
1130 : ! **************************************************************************************************
1131 4098 : PURE FUNCTION rotate_vector(a, phi, b) RESULT(a_rot)
1132 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a
1133 : REAL(KIND=dp), INTENT(IN) :: phi
1134 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: b
1135 : REAL(KIND=dp), DIMENSION(3) :: a_rot
1136 :
1137 : REAL(KIND=dp) :: length_of_b
1138 : REAL(KIND=dp), DIMENSION(3, 3) :: rotmat
1139 :
1140 4098 : length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1141 4098 : IF (length_of_b > eps_geo) THEN
1142 :
1143 : ! Build up the rotation matrix rotmat
1144 3964 : CALL build_rotmat(phi, b, rotmat)
1145 :
1146 : ! Rotate the vector a by phi about the axis defined by vector b
1147 63424 : a_rot(:) = MATMUL(rotmat, a)
1148 :
1149 : ELSE
1150 :
1151 536 : a_rot(:) = 0.0_dp
1152 :
1153 : END IF
1154 :
1155 4098 : END FUNCTION rotate_vector
1156 :
1157 : ! **************************************************************************************************
1158 : !> \brief Set the diagonal elements of matrix a to b.
1159 : !> \param a ...
1160 : !> \param b ...
1161 : !> \date 20.11.1998
1162 : !> \author MK
1163 : !> \version 1.0
1164 : ! **************************************************************************************************
1165 37088 : PURE SUBROUTINE set_diag_scalar_d(a, b)
1166 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1167 : REAL(KIND=dp), INTENT(IN) :: b
1168 :
1169 : INTEGER :: i, n
1170 :
1171 37088 : n = MIN(SIZE(a, 1), SIZE(a, 2))
1172 206599 : DO i = 1, n
1173 206599 : a(i, i) = b
1174 : END DO
1175 :
1176 37088 : END SUBROUTINE set_diag_scalar_d
1177 :
1178 : ! **************************************************************************************************
1179 : !> \brief ...
1180 : !> \param a ...
1181 : !> \param b ...
1182 : ! **************************************************************************************************
1183 0 : PURE SUBROUTINE set_diag_scalar_z(a, b)
1184 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1185 : COMPLEX(KIND=dp), INTENT(IN) :: b
1186 :
1187 : INTEGER :: i, n
1188 :
1189 0 : n = MIN(SIZE(a, 1), SIZE(a, 2))
1190 0 : DO i = 1, n
1191 0 : a(i, i) = b
1192 : END DO
1193 :
1194 0 : END SUBROUTINE set_diag_scalar_z
1195 :
1196 : ! **************************************************************************************************
1197 : !> \brief Symmetrize the matrix a.
1198 : !> \param a ...
1199 : !> \param option ...
1200 : !> \date 16.10.1998
1201 : !> \author MK
1202 : !> \version 1.0
1203 : ! **************************************************************************************************
1204 19038 : SUBROUTINE symmetrize_matrix(a, option)
1205 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1206 : CHARACTER(LEN=*), INTENT(IN) :: option
1207 :
1208 : INTEGER :: i, n
1209 :
1210 19038 : n = MIN(SIZE(a, 1), SIZE(a, 2))
1211 :
1212 19038 : IF (option == "lower_to_upper") THEN
1213 0 : DO i = 1, n - 1
1214 0 : a(i, i + 1:n) = a(i + 1:n, i)
1215 : END DO
1216 19038 : ELSE IF (option == "upper_to_lower") THEN
1217 432207 : DO i = 1, n - 1
1218 22848474 : a(i + 1:n, i) = a(i, i + 1:n)
1219 : END DO
1220 36 : ELSE IF (option == "anti_lower_to_upper") THEN
1221 0 : DO i = 1, n - 1
1222 0 : a(i, i + 1:n) = -a(i + 1:n, i)
1223 : END DO
1224 36 : ELSE IF (option == "anti_upper_to_lower") THEN
1225 564 : DO i = 1, n - 1
1226 4716 : a(i + 1:n, i) = -a(i, i + 1:n)
1227 : END DO
1228 : ELSE
1229 0 : CPABORT("Invalid option <"//TRIM(option)//"> was specified for parameter #2")
1230 : END IF
1231 :
1232 19038 : END SUBROUTINE symmetrize_matrix
1233 :
1234 : ! **************************************************************************************************
1235 : !> \brief Set the matrix a to be a unit matrix.
1236 : !> \param a ...
1237 : !> \date 16.10.1998
1238 : !> \author MK
1239 : !> \version 1.0
1240 : ! **************************************************************************************************
1241 37088 : PURE SUBROUTINE unit_matrix_d(a)
1242 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1243 :
1244 6778712 : a(:, :) = 0.0_dp
1245 37088 : CALL set_diag(a, 1.0_dp)
1246 :
1247 37088 : END SUBROUTINE unit_matrix_d
1248 :
1249 : ! **************************************************************************************************
1250 : !> \brief ...
1251 : !> \param a ...
1252 : ! **************************************************************************************************
1253 0 : PURE SUBROUTINE unit_matrix_z(a)
1254 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1255 :
1256 0 : a(:, :) = (0.0_dp, 0.0_dp)
1257 0 : CALL set_diag(a, (1.0_dp, 0.0_dp))
1258 :
1259 0 : END SUBROUTINE unit_matrix_z
1260 :
1261 : ! **************************************************************************************************
1262 : !> \brief Calculation of the vector product c = a x b.
1263 : !> \param a ...
1264 : !> \param b ...
1265 : !> \return ...
1266 : !> \date 16.10.1998
1267 : !> \author MK
1268 : !> \version 1.0
1269 : ! **************************************************************************************************
1270 7659 : PURE FUNCTION vector_product(a, b) RESULT(c)
1271 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a, b
1272 : REAL(KIND=dp), DIMENSION(3) :: c
1273 :
1274 7659 : c(1) = a(2)*b(3) - a(3)*b(2)
1275 7659 : c(2) = a(3)*b(1) - a(1)*b(3)
1276 7659 : c(3) = a(1)*b(2) - a(2)*b(1)
1277 :
1278 7659 : END FUNCTION vector_product
1279 :
1280 : ! **************************************************************************************************
1281 : !> \brief computes the greatest common divisor of two number
1282 : !> \param a ...
1283 : !> \param b ...
1284 : !> \return ...
1285 : !> \author Joost VandeVondele
1286 : ! **************************************************************************************************
1287 988522 : ELEMENTAL FUNCTION gcd(a, b)
1288 : INTEGER, INTENT(IN) :: a, b
1289 : INTEGER :: gcd
1290 :
1291 : INTEGER :: aa, ab, l, rem, s
1292 :
1293 988522 : aa = ABS(a)
1294 988522 : ab = ABS(b)
1295 988522 : IF (aa < ab) THEN
1296 : s = aa
1297 : l = ab
1298 : ELSE
1299 968612 : s = ab
1300 968612 : l = aa
1301 : END IF
1302 988522 : IF (s /= 0) THEN
1303 : DO
1304 988522 : rem = MOD(l, s)
1305 988522 : IF (rem == 0) EXIT
1306 : l = s
1307 988522 : s = rem
1308 : END DO
1309 : GCD = s
1310 : ELSE
1311 : GCD = l
1312 : END IF
1313 988522 : END FUNCTION gcd
1314 :
1315 : ! **************************************************************************************************
1316 : !> \brief computes the least common multiplier of two numbers
1317 : !> \param a ...
1318 : !> \param b ...
1319 : !> \return ...
1320 : !> \author Joost VandeVondele
1321 : ! **************************************************************************************************
1322 474308 : ELEMENTAL FUNCTION lcm(a, b)
1323 : INTEGER, INTENT(IN) :: a, b
1324 : INTEGER :: lcm
1325 :
1326 : INTEGER :: tmp
1327 :
1328 474308 : tmp = gcd(a, b)
1329 474308 : IF (tmp == 0) THEN
1330 : lcm = 0
1331 : ELSE
1332 : ! could still overflow if the true lcm is larger than maxint
1333 474308 : lcm = ABS((a/tmp)*b)
1334 : END IF
1335 474308 : END FUNCTION lcm
1336 :
1337 : ! **************************************************************************************************
1338 : !> \brief computes the exponential integral
1339 : !> Ei(x) = Int(exp(-x*t)/t,t=1..infinity) x>0
1340 : !> \param x ...
1341 : !> \return ...
1342 : !> \author JGH (adapted from Numerical recipies)
1343 : ! **************************************************************************************************
1344 0 : FUNCTION ei(x)
1345 : REAL(dp) :: x, ei
1346 :
1347 : INTEGER, PARAMETER :: maxit = 100
1348 : REAL(dp), PARAMETER :: eps = EPSILON(0.0_dp), &
1349 : fpmin = TINY(0.0_dp)
1350 :
1351 : INTEGER :: k
1352 : REAL(dp) :: fact, prev, sum1, term
1353 :
1354 0 : IF (x <= 0._dp) THEN
1355 0 : CPABORT("Invalid argument")
1356 : END IF
1357 :
1358 0 : IF (x < fpmin) THEN
1359 0 : ei = LOG(x) + euler
1360 0 : ELSE IF (x <= -LOG(EPS)) THEN
1361 : sum1 = 0._dp
1362 : fact = 1._dp
1363 0 : DO k = 1, maxit
1364 0 : fact = fact*x/REAL(k, dp)
1365 0 : term = fact/REAL(k, dp)
1366 0 : sum1 = sum1 + term
1367 0 : IF (term < eps*sum1) EXIT
1368 : END DO
1369 0 : ei = sum1 + LOG(x) + euler
1370 : ELSE
1371 : sum1 = 0._dp
1372 : term = 1._dp
1373 0 : DO k = 1, maxit
1374 0 : prev = term
1375 0 : term = term*REAL(k, dp)/x
1376 0 : IF (term < eps) EXIT
1377 0 : IF (term < prev) THEN
1378 0 : sum1 = sum1 + term
1379 : ELSE
1380 0 : sum1 = sum1 - prev
1381 0 : EXIT
1382 : END IF
1383 : END DO
1384 0 : ei = EXP(x)*(1._dp + sum1)/x
1385 : END IF
1386 :
1387 0 : END FUNCTION ei
1388 :
1389 : ! **************************************************************************************************
1390 : !> \brief computes the exponential integral
1391 : !> En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,1,..
1392 : !> Note: Ei(-x) = -E1(x)
1393 : !> \param n ...
1394 : !> \param x ...
1395 : !> \return ...
1396 : !> \par History
1397 : !> 05.2007 Created
1398 : !> \author Manuel Guidon (adapted from Numerical recipies)
1399 : ! **************************************************************************************************
1400 262050896 : ELEMENTAL IMPURE FUNCTION expint(n, x)
1401 : INTEGER, INTENT(IN) :: n
1402 : REAL(dp), INTENT(IN) :: x
1403 : REAL(dp) :: expint
1404 :
1405 : INTEGER, PARAMETER :: maxit = 100
1406 : REAL(dp), PARAMETER :: eps = 6.e-14_dp, euler = 0.5772156649015328606065120_dp, &
1407 : fpmin = TINY(0.0_dp)
1408 :
1409 : INTEGER :: i, ii, nm1
1410 : REAL(dp) :: a, b, c, d, del, fact, h, psi
1411 :
1412 262050896 : nm1 = n - 1
1413 :
1414 262050896 : IF (n < 0 .OR. x < 0.0_dp .OR. (x == 0.0_dp .AND. (n == 0 .OR. n == 1))) THEN
1415 0 : CPABORT("Invalid argument")
1416 262050896 : ELSE IF (n == 0) THEN !Special case.
1417 0 : expint = EXP(-x)/x
1418 262050896 : ELSE IF (x == 0.0_dp) THEN !Another special case.
1419 0 : expint = 1.0_dp/nm1
1420 262050896 : ELSE IF (x > 1.0_dp) THEN !Lentz's algorithm (5.2).
1421 213921911 : b = x + n
1422 213921911 : c = 1.0_dp/FPMIN
1423 213921911 : d = 1.0_dp/b
1424 213921911 : h = d
1425 6056783603 : DO i = 1, MAXIT
1426 6056783603 : a = -i*(nm1 + i)
1427 6056783603 : b = b + 2.0_dp
1428 6056783603 : d = 1.0_dp/(a*d + b)
1429 6056783603 : c = b + a/c
1430 6056783603 : del = c*d
1431 6056783603 : h = h*del
1432 6056783603 : IF (ABS(del - 1.0_dp) < EPS) THEN
1433 213921911 : expint = h*EXP(-x)
1434 213921911 : RETURN
1435 : END IF
1436 : END DO
1437 0 : CPABORT("continued fraction failed in expint")
1438 : ELSE !Evaluate series.
1439 48128985 : IF (nm1 /= 0) THEN !Set first term.
1440 0 : expint = 1.0_dp/nm1
1441 : ELSE
1442 48128985 : expint = -LOG(x) - euler
1443 : END IF
1444 : fact = 1.0_dp
1445 423129382 : DO i = 1, MAXIT
1446 423129382 : fact = -fact*x/i
1447 423129382 : IF (i /= nm1) THEN
1448 423129382 : del = -fact/(i - nm1)
1449 : ELSE
1450 : psi = -euler !Compute I(n).
1451 0 : DO ii = 1, nm1
1452 0 : psi = psi + 1.0_dp/ii
1453 : END DO
1454 0 : del = fact*(-LOG(x) + psi)
1455 : END IF
1456 423129382 : expint = expint + del
1457 423129382 : IF (ABS(del) < ABS(expint)*EPS) RETURN
1458 : END DO
1459 0 : CPABORT("series failed in expint")
1460 : END IF
1461 :
1462 : END FUNCTION expint
1463 :
1464 : ! **************************************************************************************************
1465 : !> \brief Jacobi matrix diagonalization. The eigenvalues are returned in
1466 : !> vector d and the eigenvectors are returned in matrix v in ascending
1467 : !> order.
1468 : !>
1469 : !> \param a ...
1470 : !> \param d ...
1471 : !> \param v ...
1472 : !> \par History
1473 : !> - Creation (20.11.98, Matthias Krack)
1474 : ! **************************************************************************************************
1475 32514 : SUBROUTINE jacobi(a, d, v)
1476 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1477 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: d
1478 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: v
1479 :
1480 : INTEGER :: n
1481 :
1482 32514 : n = SIZE(d(:))
1483 :
1484 : ! Diagonalize matrix a
1485 32514 : CALL diag(n, a, d, v)
1486 :
1487 : ! Sort eigenvalues and eigenvector in ascending order
1488 32514 : CALL eigsrt(n, d, v)
1489 :
1490 32514 : END SUBROUTINE jacobi
1491 :
1492 : ! **************************************************************************************************
1493 : !> \brief Diagonalize matrix a. The eigenvalues are returned in vector d
1494 : !> and the eigenvectors are returned in matrix v.
1495 : !>
1496 : !> \param n matrix/vector extent (problem size)
1497 : !> \param a matrix to be diagonalised
1498 : !> \param d vector of eigenvalues
1499 : !> \param v matrix of eigenvectors
1500 : !> \par History
1501 : !> - Creation (20.11.98, Matthias Krack)
1502 : ! **************************************************************************************************
1503 32514 : SUBROUTINE diag(n, a, d, v)
1504 : INTEGER, INTENT(IN) :: n
1505 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1506 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: d
1507 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: v
1508 :
1509 : CHARACTER(len=*), PARAMETER :: routineN = 'diag'
1510 : REAL(KIND=dp), PARAMETER :: a_eps = 1.0E-10_dp, d_eps = 1.0E-3_dp
1511 :
1512 : INTEGER :: handle, i, ip, iq
1513 : REAL(KIND=dp) :: a_max, apq, c, d_min, dip, diq, g, h, s, &
1514 : t, tau, theta, tresh
1515 32514 : REAL(KIND=dp), DIMENSION(n) :: b, z
1516 :
1517 32514 : CALL timeset(routineN, handle)
1518 :
1519 32514 : a_max = 0.0_dp
1520 119685 : DO ip = 1, n - 1
1521 965473 : a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip + 1:n))))
1522 119685 : b(ip) = a(ip, ip) ! get_diag(a)
1523 : END DO
1524 32514 : b(n) = a(n, n)
1525 :
1526 32514 : CALL unit_matrix(v)
1527 :
1528 : ! Go for 50 iterations
1529 157925 : DO i = 1, 50
1530 974476 : d = b
1531 1132401 : d_min = MAX(d_eps, MINVAL(ABS(b)))
1532 157925 : IF (a_max < a_eps*d_min) THEN
1533 32514 : CALL timestop(handle)
1534 : RETURN
1535 : END IF
1536 196970 : tresh = MERGE(a_max, 0.0_dp, (i < 4))
1537 822277 : z = 0.0_dp
1538 696866 : DO ip = 1, n - 1
1539 6399744 : DO iq = ip + 1, n
1540 5702878 : dip = d(ip)
1541 5702878 : diq = d(iq)
1542 5702878 : apq = a(ip, iq)
1543 5702878 : g = 100.0_dp*ABS(apq)
1544 6274333 : IF (tresh < ABS(apq)) THEN
1545 3333457 : h = diq - dip
1546 3333457 : IF ((ABS(h) + g) /= ABS(h)) THEN
1547 2522707 : theta = 0.5_dp*h/apq
1548 2522707 : t = 1.0_dp/(ABS(theta) + SQRT(1.0_dp + theta**2))
1549 2522707 : IF (theta < 0.0_dp) t = -t
1550 : ELSE
1551 810750 : t = apq/h
1552 : END IF
1553 3333457 : c = 1.0_dp/SQRT(1.0_dp + t**2)
1554 3333457 : s = t*c
1555 3333457 : tau = s/(1.0_dp + c)
1556 3333457 : h = t*apq
1557 3333457 : z(ip) = z(ip) - h
1558 3333457 : z(iq) = z(iq) + h
1559 3333457 : d(ip) = dip - h
1560 3333457 : d(iq) = diq + h
1561 3333457 : a(ip, iq) = 0.0_dp
1562 3333457 : CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
1563 3333457 : CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
1564 3333457 : CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
1565 3333457 : CALL jrotate(v(:, ip), v(:, iq), s, tau)
1566 : ELSE IF ((4 < i) .AND. &
1567 0 : ((ABS(dip) + g) == ABS(dip)) .AND. &
1568 2369421 : ((ABS(diq) + g) == ABS(diq))) THEN
1569 0 : a(ip, iq) = 0.0_dp
1570 : END IF
1571 : END DO
1572 : END DO
1573 822277 : b = b + z
1574 : a_max = 0.0_dp
1575 696866 : DO ip = 1, n - 1
1576 6971199 : a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip + 1:n))))
1577 : END DO
1578 : END DO
1579 0 : WRITE (*, '(/,T2,A,/)') 'Too many iterations in jacobi'
1580 :
1581 0 : CALL timestop(handle)
1582 :
1583 : END SUBROUTINE diag
1584 :
1585 : ! **************************************************************************************************
1586 : !> \brief Perform a Jacobi rotation of the vectors a and b.
1587 : !>
1588 : !> \param a ...
1589 : !> \param b ...
1590 : !> \param ss ...
1591 : !> \param tt ...
1592 : !> \par History
1593 : !> - Creation (20.11.98, Matthias Krack)
1594 : ! **************************************************************************************************
1595 13333828 : PURE SUBROUTINE jrotate(a, b, ss, tt)
1596 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1597 : REAL(KIND=dp), INTENT(IN) :: ss, tt
1598 :
1599 : REAL(KIND=dp) :: u, v
1600 :
1601 13333828 : u = 1.0_dp - ss*tt
1602 13333828 : v = ss/u
1603 :
1604 228078108 : a = a*u - b*ss
1605 228078108 : b = b*(u + ss*v) + a*v
1606 :
1607 13333828 : END SUBROUTINE jrotate
1608 :
1609 : ! **************************************************************************************************
1610 : !> \brief Sort the values in vector d in ascending order and swap the
1611 : !> corresponding columns of matrix v.
1612 : !>
1613 : !> \param n ...
1614 : !> \param d ...
1615 : !> \param v ...
1616 : !> \par History
1617 : !> - Creation (20.11.98, Matthias Krack)
1618 : ! **************************************************************************************************
1619 32514 : SUBROUTINE eigsrt(n, d, v)
1620 : INTEGER, INTENT(IN) :: n
1621 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: d
1622 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: v
1623 :
1624 : INTEGER :: i, j
1625 :
1626 119685 : DO i = 1, n - 1
1627 965473 : j = SUM(MINLOC(d(i:n))) + i - 1
1628 119685 : IF (j /= i) THEN
1629 44361 : CALL swap(d(i), d(j))
1630 44361 : CALL swap(v(:, i), v(:, j))
1631 : END IF
1632 : END DO
1633 :
1634 32514 : END SUBROUTINE eigsrt
1635 :
1636 : ! **************************************************************************
1637 : !> \brief Swap two scalars
1638 : !>
1639 : !> \param a ...
1640 : !> \param b ...
1641 : !> \par History
1642 : !> - Creation (20.11.98, Matthias Krack)
1643 : ! **************************************************************************************************
1644 44361 : ELEMENTAL SUBROUTINE swap_scalar(a, b)
1645 : REAL(KIND=dp), INTENT(INOUT) :: a, b
1646 :
1647 : REAL(KIND=dp) :: c
1648 :
1649 44361 : c = a
1650 44361 : a = b
1651 44361 : b = c
1652 :
1653 44361 : END SUBROUTINE swap_scalar
1654 :
1655 : ! **************************************************************************
1656 : !> \brief Swap two vectors
1657 : !>
1658 : !> \param a ...
1659 : !> \param b ...
1660 : !> \par History
1661 : !> - Creation (20.11.98, Matthias Krack)
1662 : ! **************************************************************************************************
1663 44361 : SUBROUTINE swap_vector(a, b)
1664 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1665 :
1666 : INTEGER :: i, n
1667 : REAL(KIND=dp) :: c
1668 :
1669 44361 : n = SIZE(a)
1670 :
1671 44361 : IF (n /= SIZE(b)) THEN
1672 0 : CPABORT("Check the array bounds of the parameters")
1673 : END IF
1674 :
1675 960905 : DO i = 1, n
1676 916544 : c = a(i)
1677 916544 : a(i) = b(i)
1678 960905 : b(i) = c
1679 : END DO
1680 :
1681 44361 : END SUBROUTINE swap_vector
1682 :
1683 : ! **************************************************************************************************
1684 : !> \brief - compute a truncation radius for the shortrange operator
1685 : !> \param eps target accuracy!> \param omg screening parameter
1686 : !> \param omg ...
1687 : !> \param r_cutoff cutoff radius
1688 : !> \par History
1689 : !> 10.2012 created [Hossein Banihashemian]
1690 : !> 05.2019 moved here from hfx_types (A. Bussy)
1691 : !> \author Hossein Banihashemian
1692 : ! **************************************************************************************************
1693 66 : SUBROUTINE erfc_cutoff(eps, omg, r_cutoff)
1694 : IMPLICIT NONE
1695 :
1696 : REAL(dp), INTENT(in) :: eps, omg
1697 : REAL(dp), INTENT(out) :: r_cutoff
1698 :
1699 : CHARACTER(LEN=*), PARAMETER :: routineN = 'erfc_cutoff'
1700 :
1701 : REAL(dp), PARAMETER :: abstol = 1E-10_dp, soltol = 1E-16_dp
1702 : REAL(dp) :: r0, f0, fprime0, delta_r
1703 : INTEGER :: iter, handle
1704 : INTEGER, PARAMETER :: iterMAX = 1000
1705 :
1706 66 : CALL timeset(routineN, handle)
1707 :
1708 : ! initial guess assuming that we are in the asymptotic regime of the erf, and the solution is about 10.
1709 66 : r0 = SQRT(-LOG(eps*omg*10**2))/omg
1710 66 : CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1711 :
1712 638 : DO iter = 1, iterMAX
1713 638 : delta_r = f0/fprime0
1714 638 : r0 = r0 - delta_r
1715 638 : CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1716 638 : IF (ABS(delta_r) < abstol .OR. ABS(f0) < soltol) EXIT
1717 : END DO
1718 66 : CPASSERT(iter <= itermax)
1719 66 : r_cutoff = r0
1720 :
1721 66 : CALL timestop(handle)
1722 : CONTAINS
1723 : ! **************************************************************************************************
1724 : !> \brief ...
1725 : !> \param r ...
1726 : !> \param eps ...
1727 : !> \param omega ...
1728 : !> \param fn ...
1729 : !> \param df ...
1730 : ! **************************************************************************************************
1731 704 : ELEMENTAL SUBROUTINE eval_transc_func(r, eps, omega, fn, df)
1732 : REAL(dp), INTENT(in) :: r, eps, omega
1733 : REAL(dp), INTENT(out) :: fn, df
1734 :
1735 : REAL(dp) :: qr
1736 :
1737 704 : qr = omega*r
1738 704 : fn = ERFC(qr) - r*eps
1739 704 : df = -2.0_dp*oorootpi*omega*EXP(-qr**2) - eps
1740 704 : END SUBROUTINE eval_transc_func
1741 : END SUBROUTINE erfc_cutoff
1742 :
1743 : ! **************************************************************************************************
1744 : !> \brief Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd
1745 : !> \param matrix Hermitian matrix is preserved
1746 : !> \param eigenvectors ...
1747 : !> \param eigenvalues ...
1748 : !> \author A. Bussy
1749 : ! **************************************************************************************************
1750 40563 : SUBROUTINE diag_complex(matrix, eigenvectors, eigenvalues)
1751 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
1752 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eigenvectors
1753 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1754 :
1755 : CHARACTER(len=*), PARAMETER :: routineN = 'diag_complex'
1756 :
1757 40563 : COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1758 : INTEGER :: handle, info, liwork, lrwork, lwork, n
1759 40563 : INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
1760 40563 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: rwork
1761 :
1762 40563 : CALL timeset(routineN, handle)
1763 :
1764 40563 : IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) CPABORT("Expected square matrix")
1765 : ! IF (MAXVAL(ABS(matrix - CONJG(TRANSPOSE(matrix)))) > 1e-14_dp) CPABORT("Expected hermitian matrix")
1766 :
1767 40563 : n = SIZE(matrix, 1)
1768 40563 : ALLOCATE (iwork(1), rwork(1), work(1))
1769 :
1770 : ! work space query
1771 40563 : lwork = -1
1772 40563 : lrwork = -1
1773 40563 : liwork = -1
1774 :
1775 40563 : CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1776 :
1777 40563 : lwork = CEILING(REAL(work(1), KIND=dp))
1778 40563 : lrwork = CEILING(rwork(1))
1779 40563 : liwork = iwork(1)
1780 :
1781 40563 : DEALLOCATE (iwork, rwork, work)
1782 283941 : ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
1783 1608109 : eigenvectors(:, :) = matrix(:, :)
1784 :
1785 : ! final diagonalization
1786 40563 : CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1787 :
1788 40563 : DEALLOCATE (iwork, rwork, work)
1789 :
1790 40563 : IF (info /= 0) CPABORT("Diagonalisation of a complex matrix failed")
1791 :
1792 40563 : CALL timestop(handle)
1793 :
1794 40563 : END SUBROUTINE diag_complex
1795 :
1796 : ! **************************************************************************************************
1797 : !> \brief Helper routine for diagonalizing anti symmetric matrices
1798 : !> \param matrix ...
1799 : !> \param evecs ...
1800 : !> \param evals ...
1801 : ! **************************************************************************************************
1802 3539 : SUBROUTINE diag_antisym(matrix, evecs, evals)
1803 : REAL(dp), DIMENSION(:, :) :: matrix
1804 : COMPLEX(dp), DIMENSION(:, :) :: evecs
1805 : COMPLEX(dp), DIMENSION(:) :: evals
1806 :
1807 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: matrix_c
1808 : INTEGER :: n
1809 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1810 :
1811 3539 : IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) CPABORT("Expected square matrix")
1812 : ! IF (MAXVAL(ABS(matrix + TRANSPOSE(matrix))) > 1e-14_dp) CPABORT("Expected anti-symmetric matrix")
1813 :
1814 3539 : n = SIZE(matrix, 1)
1815 21234 : ALLOCATE (matrix_c(n, n), eigenvalues(n))
1816 :
1817 235869 : matrix_c(:, :) = CMPLX(0.0_dp, -matrix, kind=dp)
1818 3539 : CALL diag_complex(matrix_c, evecs, eigenvalues)
1819 27874 : evals = CMPLX(0.0_dp, eigenvalues, kind=dp)
1820 :
1821 3539 : DEALLOCATE (matrix_c, eigenvalues)
1822 3539 : END SUBROUTINE diag_antisym
1823 : ! **************************************************************************************************
1824 : !> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged
1825 : !> \param A_in Input matrix 1
1826 : !> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1827 : !> \param B_in Input matrix 2
1828 : !> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1829 : !> \param C_out Output matrix
1830 : !> \par History
1831 : !> 11.2025 created [Stepan Marek]
1832 : ! **************************************************************************************************
1833 2174 : SUBROUTINE zgemm_square_2(A_in, A_trans, B_in, B_trans, C_out)
1834 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1835 : CHARACTER, INTENT(IN) :: A_trans
1836 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1837 : CHARACTER, INTENT(IN) :: B_trans
1838 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: C_out
1839 :
1840 : CHARACTER(len=*), PARAMETER :: routineN = 'zgemm_square_2'
1841 :
1842 : INTEGER :: handle, n
1843 :
1844 2174 : n = SIZE(A_in, 1)
1845 2174 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1846 2174 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1847 2174 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1848 2174 : IF (n /= SIZE(C_out, 1)) CPABORT("Incompatible (rows) result array 3 (C).")
1849 2174 : IF (n /= SIZE(C_out, 2)) CPABORT("Incompatible (cols) result array 3 (C).")
1850 2174 : IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
1851 : A_trans == 'T' .OR. A_trans == 't' .OR. &
1852 : A_trans == 'C' .OR. A_trans == 'c')) &
1853 0 : CPABORT("Unknown transpose character for array 1 (A).")
1854 2174 : IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
1855 : B_trans == 'T' .OR. B_trans == 't' .OR. &
1856 : B_trans == 'C' .OR. B_trans == 'c')) &
1857 0 : CPABORT("Unknown transpose character for array 2 (B).")
1858 :
1859 2174 : CALL timeset(routineN, handle)
1860 :
1861 2174 : CALL ZGEMM(A_trans, B_trans, n, n, n, z_one, A_in, n, B_in, n, z_zero, C_out, n)
1862 :
1863 2174 : CALL timestop(handle)
1864 :
1865 2174 : END SUBROUTINE zgemm_square_2
1866 : ! **************************************************************************************************
1867 : !> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, real matrices
1868 : !> \param A_in Input matrix 1
1869 : !> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1870 : !> \param B_in Input matrix 2
1871 : !> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1872 : !> \param C_out Output matrix
1873 : !> \par History
1874 : !> 11.2025 created [Stepan Marek]
1875 : ! **************************************************************************************************
1876 2 : SUBROUTINE dgemm_square_2(A_in, A_trans, B_in, B_trans, C_out)
1877 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1878 : CHARACTER, INTENT(IN) :: A_trans
1879 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1880 : CHARACTER, INTENT(IN) :: B_trans
1881 : REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: C_out
1882 :
1883 : CHARACTER(len=*), PARAMETER :: routineN = 'dgemm_square_2'
1884 :
1885 : INTEGER :: handle, n
1886 :
1887 2 : n = SIZE(A_in, 1)
1888 2 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1889 2 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1890 2 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1891 2 : IF (n /= SIZE(C_out, 1)) CPABORT("Incompatible (rows) result array 3 (C).")
1892 2 : IF (n /= SIZE(C_out, 2)) CPABORT("Incompatible (cols) result array 3 (C).")
1893 2 : IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
1894 : A_trans == 'T' .OR. A_trans == 't' .OR. &
1895 : A_trans == 'C' .OR. A_trans == 'c')) &
1896 0 : CPABORT("Unknown transpose character for array 1 (A).")
1897 2 : IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
1898 : B_trans == 'T' .OR. B_trans == 't' .OR. &
1899 : B_trans == 'C' .OR. B_trans == 'c')) &
1900 0 : CPABORT("Unknown transpose character for array 2 (B).")
1901 :
1902 2 : CALL timeset(routineN, handle)
1903 :
1904 2 : CALL DGEMM(A_trans, B_trans, n, n, n, 1.0_dp, A_in, n, B_in, n, 0.0_dp, C_out, n)
1905 :
1906 2 : CALL timestop(handle)
1907 :
1908 2 : END SUBROUTINE dgemm_square_2
1909 : ! **************************************************************************************************
1910 : !> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, for 3 matrices
1911 : !> \param A_in Input matrix 1
1912 : !> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1913 : !> \param B_in Input matrix 2
1914 : !> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1915 : !> \param C_in Input matrix 3
1916 : !> \param C_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 3
1917 : !> \param D_out Output matrix
1918 : !> \par History
1919 : !> 11.2025 created [Stepan Marek]
1920 : ! **************************************************************************************************
1921 29698 : SUBROUTINE zgemm_square_3(A_in, A_trans, B_in, B_trans, C_in, C_trans, D_out)
1922 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1923 : CHARACTER, INTENT(IN) :: A_trans
1924 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1925 : CHARACTER, INTENT(IN) :: B_trans
1926 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: C_in
1927 : CHARACTER, INTENT(IN) :: C_trans
1928 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: D_out
1929 :
1930 : CHARACTER(len=*), PARAMETER :: routineN = 'zgemm_square_3'
1931 :
1932 29698 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1933 : INTEGER :: handle, n
1934 :
1935 29698 : n = SIZE(A_in, 1)
1936 29698 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1937 29698 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1938 29698 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1939 29698 : IF (n /= SIZE(C_in, 1)) CPABORT("Incompatible (rows) array 3 (C).")
1940 29698 : IF (n /= SIZE(C_in, 2)) CPABORT("Non-square array 3 (C).")
1941 29698 : IF (n /= SIZE(D_out, 1)) CPABORT("Incompatible (rows) result array 4 (D).")
1942 29698 : IF (n /= SIZE(D_out, 2)) CPABORT("Incompatible (cols) result array 4 (D).")
1943 29698 : IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
1944 : A_trans == 'T' .OR. A_trans == 't' .OR. &
1945 : A_trans == 'C' .OR. A_trans == 'c')) &
1946 0 : CPABORT("Unknown transpose character for array 1 (A).")
1947 29698 : IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
1948 : B_trans == 'T' .OR. B_trans == 't' .OR. &
1949 : B_trans == 'C' .OR. B_trans == 'c')) &
1950 0 : CPABORT("Unknown transpose character for array 2 (B).")
1951 29698 : IF (.NOT. (C_trans == 'N' .OR. C_trans == 'n' .OR. &
1952 : C_trans == 'T' .OR. C_trans == 't' .OR. &
1953 : C_trans == 'C' .OR. C_trans == 'c')) &
1954 0 : CPABORT("Unknown transpose character for array 3 (C).")
1955 :
1956 29698 : CALL timeset(routineN, handle)
1957 :
1958 1193752 : ALLOCATE (work(n, n), source=z_zero)
1959 :
1960 29698 : CALL ZGEMM(A_trans, B_trans, n, n, n, z_one, A_in, n, B_in, n, z_zero, work, n)
1961 29698 : CALL ZGEMM('N', C_trans, n, n, n, z_one, work, n, C_in, n, z_zero, D_out, n)
1962 :
1963 29698 : DEALLOCATE (work)
1964 :
1965 29698 : CALL timestop(handle)
1966 :
1967 29698 : END SUBROUTINE zgemm_square_3
1968 : ! **************************************************************************************************
1969 : !> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, for 3 matrices
1970 : !> \param A_in Input matrix 1
1971 : !> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1972 : !> \param B_in Input matrix 2
1973 : !> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1974 : !> \param C_in Input matrix 3
1975 : !> \param C_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 3
1976 : !> \param D_out Output matrix
1977 : !> \par History
1978 : !> 11.2025 created [Stepan Marek]
1979 : ! **************************************************************************************************
1980 4 : SUBROUTINE dgemm_square_3(A_in, A_trans, B_in, B_trans, C_in, C_trans, D_out)
1981 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1982 : CHARACTER, INTENT(IN) :: A_trans
1983 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1984 : CHARACTER, INTENT(IN) :: B_trans
1985 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: C_in
1986 : CHARACTER, INTENT(IN) :: C_trans
1987 : REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: D_out
1988 :
1989 : CHARACTER(len=*), PARAMETER :: routineN = 'dgemm_square_3'
1990 :
1991 : INTEGER :: handle, n
1992 4 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1993 :
1994 4 : n = SIZE(A_in, 1)
1995 4 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1996 4 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1997 4 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1998 4 : IF (n /= SIZE(C_in, 1)) CPABORT("Incompatible (rows) array 3 (C).")
1999 4 : IF (n /= SIZE(C_in, 2)) CPABORT("Non-square array 3 (C).")
2000 4 : IF (n /= SIZE(D_out, 1)) CPABORT("Incompatible (rows) result array 4 (D).")
2001 4 : IF (n /= SIZE(D_out, 2)) CPABORT("Incompatible (cols) result array 4 (D).")
2002 4 : IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
2003 : A_trans == 'T' .OR. A_trans == 't' .OR. &
2004 : A_trans == 'C' .OR. A_trans == 'c')) &
2005 0 : CPABORT("Unknown transpose character for array 1 (A).")
2006 4 : IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
2007 : B_trans == 'T' .OR. B_trans == 't' .OR. &
2008 : B_trans == 'C' .OR. B_trans == 'c')) &
2009 0 : CPABORT("Unknown transpose character for array 2 (B).")
2010 4 : IF (.NOT. (C_trans == 'N' .OR. C_trans == 'n' .OR. &
2011 : C_trans == 'T' .OR. C_trans == 't' .OR. &
2012 : C_trans == 'C' .OR. C_trans == 'c')) &
2013 0 : CPABORT("Unknown transpose character for array 3 (C).")
2014 :
2015 4 : CALL timeset(routineN, handle)
2016 :
2017 64 : ALLOCATE (work(n, n), source=0.0_dp)
2018 :
2019 4 : CALL DGEMM(A_trans, B_trans, n, n, n, 1.0_dp, A_in, n, B_in, n, 0.0_dp, work, n)
2020 4 : CALL DGEMM('N', C_trans, n, n, n, 1.0_dp, work, n, C_in, n, 0.0_dp, D_out, n)
2021 :
2022 4 : DEALLOCATE (work)
2023 :
2024 4 : CALL timestop(handle)
2025 :
2026 4 : END SUBROUTINE dgemm_square_3
2027 :
2028 : ! **************************************************************************************************
2029 : !> \brief Solve the generalized eigenvalue equation for complex matrices A*v = B*v*λ
2030 : !> \param A_in ...
2031 : !> \param B_in ...
2032 : !> \param eigenvalues ...
2033 : !> \param eigenvectors ...
2034 : !> \author Shridhar Shanbhag
2035 : ! **************************************************************************************************
2036 126 : SUBROUTINE geeig_right(A_in, B_in, eigenvalues, eigenvectors)
2037 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in, B_in
2038 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
2039 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eigenvectors
2040 :
2041 : CHARACTER(len=*), PARAMETER :: routineN = 'geeig_right'
2042 : COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
2043 : czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2044 :
2045 126 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
2046 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: A, B, work
2047 : INTEGER :: handle, i, icol, irow, lda, ldb, ldc, &
2048 : nao, nc, ncol, nmo, nx
2049 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
2050 :
2051 126 : CALL timeset(routineN, handle)
2052 :
2053 : ! Test sizes
2054 126 : nao = SIZE(A_in, 1)
2055 126 : nmo = SIZE(eigenvalues)
2056 630 : ALLOCATE (evals(nao), cevals(nao))
2057 1008 : ALLOCATE (work(nao, nao), B(nao, nao), A(nao, nao))
2058 3582 : A(:, :) = A_in(:, :)
2059 3582 : B(:, :) = B_in(:, :)
2060 :
2061 : ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
2062 3582 : B = -B
2063 126 : CALL diag_complex(B, work, evals)
2064 702 : evals(:) = -evals(:)
2065 126 : nc = nao
2066 702 : DO i = 1, nao
2067 702 : IF (evals(i) < -1.0_dp) THEN
2068 0 : nc = i - 1
2069 0 : EXIT
2070 : END IF
2071 : END DO
2072 126 : CPASSERT(nc /= 0)
2073 :
2074 126 : IF (nc /= nao) THEN
2075 0 : IF (nc < nmo) THEN
2076 : ! Copy NULL space definition to last vectors of eigenvectors (if needed)
2077 0 : ncol = nmo - nc
2078 0 : CALL zcopy(ncol*nao, work(1, nc + 1), 1, eigenvectors(1, nc + 1), 1)
2079 : END IF
2080 : ! Set NULL space in eigenvector matrix of S to zero
2081 0 : DO icol = nc + 1, nao
2082 0 : DO irow = 1, nao
2083 0 : work(irow, icol) = czero
2084 : END DO
2085 : END DO
2086 : ! Set small eigenvalues to a dummy save value
2087 0 : evals(nc + 1:nao) = 1.0_dp
2088 : END IF
2089 : ! Calculate U*s**(-1/2)
2090 702 : cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
2091 702 : DO i = 1, MIN(SIZE(work, 2), SIZE(cevals))
2092 702 : CALL zscal(MIN(SIZE(work, 2), SIZE(cevals)), cevals(i), work(1, i), 1)
2093 : END DO
2094 : ! Reduce to get U^(-C) * H * U^(-1)
2095 126 : CALL gemm_square(work, 'C', A, 'N', B)
2096 126 : CALL gemm_square(B, 'N', work, 'N', A)
2097 126 : IF (nc /= nao) THEN
2098 : ! set diagonal values to save large value
2099 0 : DO icol = nc + 1, nao
2100 0 : A(icol, icol) = 10000*cone
2101 : END DO
2102 : END IF
2103 : ! Diagonalize
2104 126 : CALL diag_complex(A, B, evals)
2105 702 : eigenvalues(1:nmo) = evals(1:nmo)
2106 126 : nx = MIN(nc, nmo)
2107 : ! Restore vectors C = U^(-1) * C*
2108 126 : lda = SIZE(work, 1)
2109 126 : ldb = SIZE(B, 1)
2110 126 : ldc = SIZE(eigenvectors, 1)
2111 : CALL ZGEMM("N", "N", nao, nx, nc, cone, work, &
2112 126 : lda, B, ldb, czero, eigenvectors, ldc)
2113 :
2114 126 : DEALLOCATE (evals, cevals, work, B, A)
2115 :
2116 126 : CALL timestop(handle)
2117 126 : END SUBROUTINE geeig_right
2118 :
2119 : END MODULE mathlib
|
