(git:6344619)
Loading...
Searching...
No Matches
kubo_bloch.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
20CONTAINS
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 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 COMPLEX(dp), ALLOCATABLE :: c(:, :)
44 INTEGER :: n
45 REAL(dp) :: residuals(2)
46
47 n = SIZE(h, 1)
48 status = -1
49 energies = 0.0_dp
50 current = 0.0_dp
51 IF (PRESENT(coefficients)) THEN
52 coefficients = 0.0_dp
53 IF (any(shape(coefficients) /= [n, n])) RETURN
54 END IF
55 ALLOCATE (c(n, n))
56 CALL projected_eigensystem(h, s, energies, c, status)
57 IF (status /= 0) RETURN
58 CALL bloch_current_from_frame(h, s, dh, ds, d, energies, c, current, 1.e-8_dp, residuals, status)
59 IF (status == 0 .AND. PRESENT(coefficients)) coefficients = c
60 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 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 COMPLEX(dp), ALLOCATABLE :: dc(:, :), hc(:, :), sc(:, :), v(:, :)
88 INTEGER :: axis, i, j, n
89 REAL(dp) :: scale
90
91 status = -1
92 current = 0.0_dp
93 residuals = huge(1.0_dp)
94 n = SIZE(h, 1)
95 IF (n < 1 .OR. any(shape(h) /= [n, n]) .OR. any(shape(s) /= [n, n])) RETURN
96 IF (SIZE(energies) /= n .OR. any(shape(current) /= [n, n, 3])) RETURN
97 IF (any(shape(dh) /= [n, n, 3]) .OR. any(shape(ds) /= [n, n, 3])) RETURN
98 IF (any(shape(d) /= [n, n, 3])) RETURN
99 IF (any(shape(c) /= [n, n])) RETURN
100 IF (.NOT. ieee_is_finite(tolerance) .OR. tolerance <= 0.0_dp .OR. tolerance >= 1.0_dp) RETURN
101 IF (.NOT. all(ieee_is_finite(energies))) RETURN
102 IF (.NOT. all(ieee_is_finite(real(c, dp))) .OR. .NOT. all(ieee_is_finite(aimag(c)))) RETURN
103 IF (.NOT. all(ieee_is_finite(real(h, dp))) .OR. .NOT. all(ieee_is_finite(aimag(h)))) RETURN
104 IF (.NOT. all(ieee_is_finite(real(s, dp))) .OR. .NOT. all(ieee_is_finite(aimag(s)))) RETURN
105 IF (maxval(abs(h - conjg(transpose(h)))) > 1.e-9_dp*max(1.0_dp, maxval(abs(h)))) RETURN
106 IF (maxval(abs(s - conjg(transpose(s)))) > 1.e-9_dp*max(1.0_dp, maxval(abs(s)))) RETURN
107 IF (.NOT. finite(dh) .OR. .NOT. finite(ds) .OR. .NOT. finite(d)) RETURN
108 DO axis = 1, 3
109 scale = max(tiny(1.0_dp), maxval(abs(dh(:, :, axis))))
110 IF (maxval(abs(dh(:, :, axis) - conjg(transpose(dh(:, :, axis))))) > 1.e-9_dp*scale) RETURN
111 scale = max(tiny(1.0_dp), maxval(abs(ds(:, :, axis))))
112 IF (maxval(abs(ds(:, :, axis) - conjg(transpose(ds(:, :, axis))))) > 1.e-9_dp*scale) RETURN
113 status = -3
114 scale = max(1.0_dp, maxval(abs(ds(:, :, axis))), maxval(abs(d(:, :, axis))))
115 IF (maxval(abs(d(:, :, axis) - conjg(transpose(d(:, :, axis))) + iu*ds(:, :, axis))) &
116 > 1.e-9_dp*scale) RETURN
117 status = -1
118 END DO
119 ALLOCATE (hc(n, n), sc(n, n), dc(n, n), v(n, n))
120 sc(:, :) = matmul(conjg(transpose(c)), matmul(s, c))
121 hc(:, :) = matmul(conjg(transpose(c)), matmul(h, c))
122 DO i = 1, n
123 sc(i, i) = sc(i, i) - 1.0_dp
124 hc(i, i) = hc(i, i) - energies(i)
125 END DO
126 residuals = [maxval(abs(sc)), maxval(abs(hc))/max(1.0_dp, maxval(abs(energies)))]
127 status = -4
128 IF (.NOT. all(ieee_is_finite(residuals)) .OR. maxval(residuals) > tolerance) RETURN
129 status = -1
130 DO axis = 1, 3
131 hc(:, :) = matmul(conjg(transpose(c)), matmul(dh(:, :, axis), c))
132 sc(:, :) = matmul(conjg(transpose(c)), matmul(ds(:, :, axis), c))
133 dc(:, :) = matmul(conjg(transpose(c)), matmul(d(:, :, axis), c))
134 DO j = 1, n
135 DO i = 1, n
136 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 current(:, :, axis) = -iu*v
140 END DO
141 IF (.NOT. finite(current)) RETURN
142 status = 0
143 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 LOGICAL FUNCTION finite(a)
150 COMPLEX(dp), INTENT(IN) :: a(:, :, :)
151
152 finite = all(ieee_is_finite(real(a, dp))) .AND. all(ieee_is_finite(aimag(a)))
153 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 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 delta = ea - eb
167 IF (abs(delta) < 1.e-4_dp*kt) THEN
168 f = occupation((ea + eb)*0.5_dp)
169 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 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 PURE REAL(dp) FUNCTION occupation(e) RESULT(f)
180 REAL(dp), INTENT(IN) :: e
181
182 REAL(dp) :: x
183
184 x = exp(-abs((e - mu)/kt))
185 f = x/(1 + x)
186 IF (e < mu) f = 1.0_dp - f
187 END FUNCTION occupation
188 END FUNCTION fermi_slope
189END MODULE kubo_bloch
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kubo_bloch
Metric-aware Bloch velocities without eigenvector derivatives or band-gap divisions.
Definition kubo_bloch.F:12
kubo_bloch::bloch_current
subroutine, public bloch_current(h, s, dh, ds, d, energies, current, status, coefficients)
Evaluate v_nm = C^H(dH - E_n*dS)C + i*(E_n-E_m)*C^H*D*C.
Definition kubo_bloch.F:36
kubo_bloch::bloch_current_from_frame
subroutine, public bloch_current_from_frame(h, s, dh, ds, d, energies, c, current, tolerance, residuals, status)
Evaluate the current in a validated reconstructed eigenframe, without rediagonalization.
Definition kubo_bloch.F:76
kubo_bloch::fermi_slope
pure real(dp) function, public fermi_slope(ea, eb, mu, kt)
Stable Fermi divided difference including the intraband and degenerate limits.
Definition kubo_bloch.F:163
kubo_projected
Finite projected-Hamiltonian current in a nonorthogonal, possibly spinor AO basis.
Definition kubo_projected.F:12
kubo_projected::projected_eigensystem
subroutine, public projected_eigensystem(h, s, energies, c, status)
Complete generalized eigensystem with a positive, well-conditioned AO metric.
Definition kubo_projected.F:78