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 Metric-aware Bloch velocities without eigenvector derivatives or band-gap divisions.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE kubo_bloch
13 : USE ieee_arithmetic, ONLY: ieee_is_finite
14 : USE kinds, ONLY: dp
15 : USE kubo_projected, ONLY: projected_eigensystem
16 :
17 : IMPLICIT NONE
18 : PRIVATE
19 : PUBLIC :: bloch_current, bloch_current_from_frame, fermi_slope
20 : CONTAINS
21 : ! **************************************************************************************************
22 : !> \brief Evaluate v_nm = C^H(dH - E_n*dS)C + i*(E_n-E_m)*C^H*D*C.
23 : !> \param h Full covariant Hamiltonian, including SOC when applicable
24 : !> \param s Overlap in the same Bloch convention
25 : !> \param dh Cartesian k derivatives of h, in Ha*bohr
26 : !> \param ds Cartesian k derivatives of s, in bohr
27 : !> \param d Cell-gauge image first moments, satisfying D-D^H=-i*dS, in bohr
28 : !> \param energies Generalized spectrum
29 : !> \param current Anti-Hermitian current -i*v in the eigenbasis, including diagonal velocities
30 : !> \param status 0 success; -1 invalid data; -2 metric failure; -3 connection failure; positive LAPACK failure
31 : !> \param coefficients Optional complete eigenframe for symmetry reconstruction
32 : !> \note For a general k-dependent basis D transforms as B^H*D*B-i*(dB)^H*S*B.
33 : !> The connection terms are essential even when the spectrum is accurately reproduced.
34 : ! **************************************************************************************************
35 27805 : SUBROUTINE bloch_current(h, s, dh, ds, d, energies, current, status, coefficients)
36 : COMPLEX(dp), INTENT(IN) :: h(:, :), s(:, :), dh(:, :, :), &
37 : ds(:, :, :), d(:, :, :)
38 : REAL(dp), INTENT(OUT) :: energies(:)
39 : COMPLEX(dp), INTENT(OUT) :: current(:, :, :)
40 : INTEGER, INTENT(OUT) :: status
41 : COMPLEX(dp), INTENT(OUT), OPTIONAL :: coefficients(:, :)
42 :
43 27805 : COMPLEX(dp), ALLOCATABLE :: c(:, :)
44 : INTEGER :: n
45 : REAL(dp) :: residuals(2)
46 :
47 27805 : n = SIZE(h, 1)
48 27805 : status = -1
49 84159 : energies = 0.0_dp
50 694690 : current = 0.0_dp
51 27805 : IF (PRESENT(coefficients)) THEN
52 28087 : coefficients = 0.0_dp
53 183 : IF (ANY(SHAPE(coefficients) /= [n, n])) RETURN
54 : END IF
55 111220 : ALLOCATE (c(n, n))
56 27805 : CALL projected_eigensystem(h, s, energies, c, status)
57 27805 : IF (status /= 0) RETURN
58 27805 : CALL bloch_current_from_frame(h, s, dh, ds, d, energies, c, current, 1.e-8_dp, residuals, status)
59 55831 : IF (status == 0 .AND. PRESENT(coefficients)) coefficients = c
60 27805 : END SUBROUTINE bloch_current
61 : ! **************************************************************************************************
62 : !> \brief Evaluate the current in a validated reconstructed eigenframe, without rediagonalization.
63 : !> \param h Target Hamiltonian
64 : !> \param s Target AO metric
65 : !> \param dh Cartesian Hamiltonian derivatives
66 : !> \param ds Cartesian overlap derivatives
67 : !> \param d Cell-gauge image first moments
68 : !> \param energies Source spectrum, including complete degenerate blocks
69 : !> \param c Candidate target eigenframe
70 : !> \param current Target current in exactly the supplied gauge
71 : !> \param tolerance Maximum orthonormality and relative projected-eigenproblem errors
72 : !> \param residuals Metric error and energy-scaled projected-Hamiltonian error
73 : !> \param status 0 success; -1 invalid data; -3 connection failure; -4 invalid eigenframe
74 : ! **************************************************************************************************
75 27911 : SUBROUTINE bloch_current_from_frame(h, s, dh, ds, d, energies, c, current, tolerance, residuals, status)
76 : COMPLEX(dp), INTENT(IN) :: h(:, :), s(:, :), dh(:, :, :), &
77 : ds(:, :, :), d(:, :, :)
78 : REAL(dp), INTENT(IN) :: energies(:)
79 : COMPLEX(dp), INTENT(IN) :: c(:, :)
80 : COMPLEX(dp), INTENT(OUT) :: current(:, :, :)
81 : REAL(dp), INTENT(IN) :: tolerance
82 : REAL(dp), INTENT(OUT) :: residuals(2)
83 : INTEGER, INTENT(OUT) :: status
84 :
85 : COMPLEX(dp), PARAMETER :: iu = (0.0_dp, 1.0_dp)
86 :
87 27911 : COMPLEX(dp), ALLOCATABLE :: dc(:, :), hc(:, :), sc(:, :), v(:, :)
88 : INTEGER :: axis, i, j, n
89 : REAL(dp) :: scale
90 :
91 27911 : status = -1
92 785810 : current = 0.0_dp
93 83733 : residuals = HUGE(1.0_dp)
94 27911 : n = SIZE(h, 1)
95 139555 : IF (n < 1 .OR. ANY(SHAPE(h) /= [n, n]) .OR. ANY(SHAPE(s) /= [n, n])) RETURN
96 111644 : IF (SIZE(energies) /= n .OR. ANY(SHAPE(current) /= [n, n, 3])) RETURN
97 195377 : IF (ANY(SHAPE(dh) /= [n, n, 3]) .OR. ANY(SHAPE(ds) /= [n, n, 3])) RETURN
98 111644 : IF (ANY(SHAPE(d) /= [n, n, 3])) RETURN
99 83733 : IF (ANY(SHAPE(c) /= [n, n])) RETURN
100 27911 : IF (.NOT. ieee_is_finite(tolerance) .OR. tolerance <= 0.0_dp .OR. tolerance >= 1.0_dp) RETURN
101 85069 : IF (.NOT. ALL(ieee_is_finite(energies))) RETURN
102 505266 : IF (.NOT. ALL(ieee_is_finite(REAL(c, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(c)))) RETURN
103 505266 : IF (.NOT. ALL(ieee_is_finite(REAL(h, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(h)))) RETURN
104 505266 : IF (.NOT. ALL(ieee_is_finite(REAL(s, dp))) .OR. .NOT. ALL(ieee_is_finite(AIMAG(s)))) RETURN
105 505266 : IF (MAXVAL(ABS(h - CONJG(TRANSPOSE(h)))) > 1.e-9_dp*MAX(1.0_dp, MAXVAL(ABS(h)))) RETURN
106 505266 : IF (MAXVAL(ABS(s - CONJG(TRANSPOSE(s)))) > 1.e-9_dp*MAX(1.0_dp, MAXVAL(ABS(s)))) RETURN
107 27911 : IF (.NOT. finite(dh) .OR. .NOT. finite(ds) .OR. .NOT. finite(d)) RETURN
108 111554 : DO axis = 1, 3
109 757479 : scale = MAX(TINY(1.0_dp), MAXVAL(ABS(dh(:, :, axis))))
110 757479 : IF (MAXVAL(ABS(dh(:, :, axis) - CONJG(TRANSPOSE(dh(:, :, axis))))) > 1.e-9_dp*scale) RETURN
111 757479 : scale = MAX(TINY(1.0_dp), MAXVAL(ABS(ds(:, :, axis))))
112 757479 : IF (MAXVAL(ABS(ds(:, :, axis) - CONJG(TRANSPOSE(ds(:, :, axis))))) > 1.e-9_dp*scale) RETURN
113 83673 : status = -3
114 1431285 : scale = MAX(1.0_dp, MAXVAL(ABS(ds(:, :, axis))), MAXVAL(ABS(d(:, :, axis))))
115 757479 : IF (MAXVAL(ABS(d(:, :, axis) - CONJG(TRANSPOSE(d(:, :, axis))) + iu*ds(:, :, axis))) &
116 : > 1.e-9_dp*scale) RETURN
117 111524 : status = -1
118 : END DO
119 278810 : ALLOCATE (hc(n, n), sc(n, n), dc(n, n), v(n, n))
120 2080856 : sc(:, :) = MATMUL(CONJG(TRANSPOSE(c)), MATMUL(s, c))
121 2080856 : hc(:, :) = MATMUL(CONJG(TRANSPOSE(c)), MATMUL(h, c))
122 84979 : DO i = 1, n
123 57098 : sc(i, i) = sc(i, i) - 1.0_dp
124 84979 : hc(i, i) = hc(i, i) - energies(i)
125 : END DO
126 617706 : residuals = [MAXVAL(ABS(sc)), MAXVAL(ABS(hc))/MAX(1.0_dp, MAXVAL(ABS(energies)))]
127 27881 : status = -4
128 167286 : IF (.NOT. ALL(ieee_is_finite(residuals)) .OR. MAXVAL(residuals) > tolerance) RETURN
129 27821 : status = -1
130 111284 : DO axis = 1, 3
131 7798320 : hc(:, :) = MATMUL(CONJG(TRANSPOSE(c)), MATMUL(dh(:, :, axis), c))
132 7798320 : sc(:, :) = MATMUL(CONJG(TRANSPOSE(c)), MATMUL(ds(:, :, axis), c))
133 7798320 : dc(:, :) = MATMUL(CONJG(TRANSPOSE(c)), MATMUL(d(:, :, axis), c))
134 254397 : DO j = 1, n
135 756009 : DO i = 1, n
136 672546 : v(i, j) = hc(i, j) - energies(i)*sc(i, j) + iu*(energies(i) - energies(j))*dc(i, j)
137 : END DO
138 : END DO
139 783830 : current(:, :, axis) = -iu*v
140 : END DO
141 27821 : IF (.NOT. finite(current)) RETURN
142 27821 : status = 0
143 27911 : END SUBROUTINE bloch_current_from_frame
144 : ! **************************************************************************************************
145 : !> \brief Finite check on all components, including connection matrices which are not Hermitian.
146 : !> \param a Matrix triplet
147 : !> \return Whether all entries are finite
148 : ! **************************************************************************************************
149 111554 : LOGICAL FUNCTION finite(a)
150 : COMPLEX(dp), INTENT(IN) :: a(:, :, :)
151 :
152 6282520 : finite = ALL(ieee_is_finite(REAL(a, dp))) .AND. ALL(ieee_is_finite(AIMAG(a)))
153 111554 : END FUNCTION finite
154 : ! **************************************************************************************************
155 : !> \brief Stable Fermi divided difference including the intraband and degenerate limits.
156 : !> \param ea First energy
157 : !> \param eb Second energy
158 : !> \param mu Chemical potential
159 : !> \param kt Positive temperature in energy units
160 : !> \return (f(ea)-f(eb))/(ea-eb) for unit occupation, nonpositive
161 : ! **************************************************************************************************
162 1102708 : PURE REAL(dp) FUNCTION fermi_slope(ea, eb, mu, kt) RESULT(slope)
163 : REAL(dp), INTENT(IN) :: ea, eb, mu, kt
164 : REAL(dp) :: delta, f
165 :
166 1102708 : delta = ea - eb
167 1102708 : IF (ABS(delta) < 1.e-4_dp*kt) THEN
168 72356 : f = occupation((ea + eb)*0.5_dp)
169 72356 : slope = -f*(1.0_dp - f)/kt*(1.0_dp + (delta/kt)**2*(1.0_dp - 6*f + 6*f*f)/24.0_dp)
170 : ELSE
171 1030352 : slope = (occupation(ea) - occupation(eb))/delta
172 : END IF
173 : CONTAINS
174 : ! **************************************************************************************************
175 : !> \brief Stable logistic occupation.
176 : !> \param e Energy
177 : !> \return Unit-occupation Fermi function
178 : ! **************************************************************************************************
179 2133060 : PURE REAL(dp) FUNCTION occupation(e) RESULT(f)
180 : REAL(dp), INTENT(IN) :: e
181 :
182 : REAL(dp) :: x
183 :
184 2133060 : x = EXP(-ABS((e - mu)/kt))
185 2133060 : f = x/(1 + x)
186 2133060 : IF (e < mu) f = 1.0_dp - f
187 2133060 : END FUNCTION occupation
188 : END FUNCTION fermi_slope
189 306151 : END MODULE kubo_bloch
|
