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 27829 : 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 27829 : COMPLEX(dp), ALLOCATABLE :: b(:, :), work(:)
84 27829 : REAL(dp), ALLOCATABLE :: rwork(:)
85 27829 : INTEGER, ALLOCATABLE :: iwork(:)
86 : INTEGER :: n, info
87 : #if defined(__HAS_IEEE_EXCEPTIONS)
88 : LOGICAL :: halt(SIZE(IEEE_ALL))
89 : #endif
90 :
91 27829 : status = -1
92 84515 : energies = 0.0_dp
93 229539 : c = 0.0_dp
94 27829 : n = SIZE(h, 1)
95 83487 : IF (n < 1 .OR. SIZE(h, 2) /= n .OR. ANY(SHAPE(s) /= [n, n])) RETURN
96 83487 : IF (SIZE(energies) /= n .OR. ANY(SHAPE(c) /= [n, n])) RETURN
97 27829 : IF (.NOT. valid_hermitian(h) .OR. .NOT. valid_hermitian(s)) RETURN
98 278290 : ALLOCATE (b(n, n), work(MAX(1, 2*n + n*n)), rwork(MAX(1, 1 + 5*n + 2*n*n)), iwork(MAX(1, 3 + 5*n)))
99 229539 : b(:, :) = s
100 : #if defined(__HAS_IEEE_EXCEPTIONS)
101 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
102 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
103 : #endif
104 27829 : CALL zheevd('N', 'U', n, b, n, energies, work, SIZE(work), rwork, SIZE(rwork), iwork, SIZE(iwork), info)
105 : #if defined(__HAS_IEEE_EXCEPTIONS)
106 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
107 : #endif
108 27829 : status = ABS(info) + 1
109 84515 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(energies))) RETURN
110 27829 : status = -2
111 84515 : IF (MAXVAL(energies) <= 0.0_dp) RETURN
112 141201 : IF (MINVAL(energies) <= 1.e-10_dp*MAXVAL(energies)) RETURN
113 229497 : c = h
114 229497 : b(:, :) = s
115 : #if defined(__HAS_IEEE_EXCEPTIONS)
116 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
117 : #endif
118 27827 : CALL zhegvd(1, 'V', 'U', n, c, n, b, n, energies, work, SIZE(work), rwork, SIZE(rwork), iwork, SIZE(iwork), info)
119 : #if defined(__HAS_IEEE_EXCEPTIONS)
120 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
121 : #endif
122 27827 : status = ABS(info) + 1
123 84505 : IF (info /= 0 .OR. .NOT. ALL(ieee_is_finite(energies))) RETURN
124 229497 : IF (.NOT. ALL(ieee_is_finite(REAL(c, dp)))) RETURN
125 229497 : IF (.NOT. ALL(ieee_is_finite(AIMAG(c)))) RETURN
126 27827 : status = 0
127 27829 : END SUBROUTINE projected_eigensystem
128 : ! **************************************************************************************************
129 : !> \brief Reject nonfinite or non-Hermitian data instead of discarding a triangle.
130 : !> \param a Covariant matrix
131 : !> \return Whether both components are finite and Hermitian within tolerance
132 : ! **************************************************************************************************
133 55724 : LOGICAL FUNCTION valid_hermitian(a) RESULT(valid)
134 : COMPLEX(dp), INTENT(IN) :: a(:, :)
135 :
136 930168 : valid = ALL(ieee_is_finite(REAL(a, dp))) .AND. ALL(ieee_is_finite(AIMAG(a)))
137 55724 : IF (.NOT. valid) RETURN
138 930168 : valid = MAXVAL(ABS(a - CONJG(TRANSPOSE(a)))) <= 1.e-10_dp*MAX(TINY(1.0_dp), MAXVAL(ABS(a)))
139 55724 : END FUNCTION valid_hermitian
140 30 : END MODULE kubo_projected
|
