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 Finite projected-Hamiltonian current in a nonorthogonal, possibly spinor AO basis.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE kubo_projected
13 : USE ieee_arithmetic, ONLY: ieee_is_finite
14 : #if defined(__HAS_IEEE_EXCEPTIONS)
15 : USE ieee_exceptions, ONLY: IEEE_ALL, ieee_get_halting_mode, ieee_set_halting_mode
16 : #endif
17 : USE kinds, ONLY: dp
18 : IMPLICIT NONE
19 : PRIVATE
20 : PUBLIC :: projected_current, projected_eigensystem
21 : CONTAINS
22 : ! **************************************************************************************************
23 : !> \brief Compute C^dagger (H S^-1 X - X S^-1 H) C without forming S^-1.
24 : !> \param h Complete scalar/spinor covariant Hamiltonian, not an occupied-state projection
25 : !> \param s Positive definite overlap, with explicit spin blocks when applicable
26 : !> \param positions Three analytic covariant Cartesian first moments
27 : !> \param energies Eigenvalues of H C = S C E
28 : !> \param current Anti-Hermitian commutators in the eigenbasis; velocity is i*current in atomic units
29 : !> \param status 0 success; -1 invalid matrices; -2 ill-conditioned metric; positive LAPACK failure
30 : !> \note This is the commutator of projected operators. It is not the projection of the
31 : !> continuum commutator in an incomplete basis, nor a periodic Berry velocity.
32 : ! **************************************************************************************************
33 14 : SUBROUTINE projected_current(h, s, positions, energies, current, status)
34 : COMPLEX(dp), INTENT(IN) :: h(:, :), s(:, :), positions(:, :, :)
35 : REAL(dp), INTENT(OUT) :: energies(:)
36 : COMPLEX(dp), INTENT(OUT) :: current(:, :, :)
37 : INTEGER, INTENT(OUT) :: status
38 :
39 14 : COMPLEX(dp), ALLOCATABLE :: c(:, :), moment(:, :)
40 : INTEGER :: axis, i, j, n
41 :
42 14 : status = -1
43 106 : energies = 0.0_dp
44 3668 : current = 0.0_dp
45 14 : n = SIZE(h, 1)
46 42 : IF (n < 1 .OR. SIZE(h, 2) /= n .OR. ANY(SHAPE(s) /= [n, n])) RETURN
47 56 : IF (SIZE(energies) /= n .OR. ANY(SHAPE(current) /= [n, n, 3])) RETURN
48 56 : IF (ANY(SHAPE(positions) /= [n, n, 3])) RETURN
49 14 : IF (.NOT. valid_hermitian(h) .OR. .NOT. valid_hermitian(s)) RETURN
50 50 : DO axis = 1, 3
51 50 : IF (.NOT. valid_hermitian(positions(:, :, axis))) RETURN
52 : END DO
53 72 : ALLOCATE (c(n, n), moment(n, n))
54 12 : CALL projected_eigensystem(h, s, energies, c, status)
55 12 : IF (status /= 0) RETURN
56 10 : status = -1
57 40 : DO axis = 1, 3
58 507216 : moment(:, :) = MATMUL(CONJG(TRANSPOSE(c)), MATMUL(positions(:, :, axis), c))
59 268 : DO j = 1, n
60 3402 : DO i = 1, n
61 3372 : current(i, j, axis) = (energies(i) - energies(j))*moment(i, j)
62 : END DO
63 : END DO
64 : END DO
65 3412 : IF (.NOT. ALL(ieee_is_finite(REAL(current, dp)))) RETURN
66 3412 : IF (.NOT. ALL(ieee_is_finite(AIMAG(current)))) RETURN
67 10 : status = 0
68 14 : END SUBROUTINE projected_current
69 : ! **************************************************************************************************
70 : !> \brief Complete generalized eigensystem with a positive, well-conditioned AO metric.
71 : !> \param h Hermitian covariant Hamiltonian
72 : !> \param s Hermitian positive definite overlap
73 : !> \param energies Complete spectrum
74 : !> \param c Metric-orthonormal eigenvectors
75 : !> \param status 0 success; -1 invalid matrices; -2 metric failure; positive LAPACK failure
76 : ! **************************************************************************************************
77 27817 : SUBROUTINE projected_eigensystem(h, s, energies, c, status)
78 : COMPLEX(dp), INTENT(IN) :: h(:, :), s(:, :)
79 : REAL(dp), INTENT(OUT) :: energies(:)
80 : COMPLEX(dp), INTENT(OUT) :: c(:, :)
81 : INTEGER, INTENT(OUT) :: status
82 :
83 27817 : COMPLEX(dp), ALLOCATABLE :: b(:, :), work(:)
84 27817 : REAL(dp), ALLOCATABLE :: rwork(:)
85 : INTEGER :: n, info
86 : #if defined(__HAS_IEEE_EXCEPTIONS)
87 : LOGICAL :: halt(SIZE(IEEE_ALL))
88 : #endif
89 :
90 27817 : status = -1
91 84255 : energies = 0.0_dp
92 223471 : c = 0.0_dp
93 27817 : n = SIZE(h, 1)
94 83451 : IF (n < 1 .OR. SIZE(h, 2) /= n .OR. ANY(SHAPE(s) /= [n, n])) RETURN
95 83451 : IF (SIZE(energies) /= n .OR. ANY(SHAPE(c) /= [n, n])) RETURN
96 27817 : IF (.NOT. valid_hermitian(h) .OR. .NOT. valid_hermitian(s)) RETURN
97 222536 : ALLOCATE (b(n, n), work(MAX(1, 2*n)), rwork(MAX(1, 3*n)))
98 223471 : b(:, :) = s
99 : #if defined(__HAS_IEEE_EXCEPTIONS)
100 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
101 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
102 : #endif
103 27817 : CALL zheev('N', 'U', n, b, n, energies, work, SIZE(work), rwork, info)
104 : #if defined(__HAS_IEEE_EXCEPTIONS)
105 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
106 : #endif
107 27817 : status = ABS(info) + 1
108 84255 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(energies))) RETURN
109 27817 : status = -2
110 84255 : IF (MAXVAL(energies) <= 0.0_dp) RETURN
111 140693 : IF (MINVAL(energies) <= 1.e-10_dp*MAXVAL(energies)) RETURN
112 223429 : c = h
113 223429 : b(:, :) = s
114 : #if defined(__HAS_IEEE_EXCEPTIONS)
115 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
116 : #endif
117 27815 : CALL zhegv(1, 'V', 'U', n, c, n, b, n, energies, work, SIZE(work), rwork, info)
118 : #if defined(__HAS_IEEE_EXCEPTIONS)
119 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
120 : #endif
121 27815 : status = ABS(info) + 1
122 84245 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(energies))) RETURN
123 223429 : IF (.NOT. ALL(ieee_is_finite(REAL(c, dp)))) RETURN
124 223429 : IF (.NOT. ALL(ieee_is_finite(AIMAG(c)))) RETURN
125 27815 : status = 0
126 27817 : END SUBROUTINE projected_eigensystem
127 : ! **************************************************************************************************
128 : !> \brief Reject nonfinite or non-Hermitian data instead of discarding a triangle.
129 : !> \param a Covariant matrix
130 : !> \return Whether both components are finite and Hermitian within tolerance
131 : ! **************************************************************************************************
132 55700 : LOGICAL FUNCTION valid_hermitian(a) RESULT(valid)
133 : COMPLEX(dp), INTENT(IN) :: a(:, :)
134 :
135 905896 : valid = ALL(ieee_is_finite(REAL(a, dp))) .AND. ALL(ieee_is_finite(AIMAG(a)))
136 55700 : IF (.NOT. valid) RETURN
137 905896 : valid = MAXVAL(ABS(a - CONJG(TRANSPOSE(a)))) <= 1.e-10_dp*MAX(TINY(1.0_dp), MAXVAL(ABS(a)))
138 55700 : END FUNCTION valid_hermitian
139 30 : END MODULE kubo_projected
|
