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 Full-mesh scalar/spinor Bloch currents from CP2K image operators and analytic AO moments.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE qs_kubo_bloch
13 : USE cell_types, ONLY: cell_type
14 : USE cp_control_types, ONLY: dft_control_type
15 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
16 : dbcsr_p_type
17 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
18 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
19 : USE input_constants, ONLY: do_method_gapw,&
20 : do_method_gpw
21 : USE kinds, ONLY: dp,&
22 : int_8
23 : USE kpoint_k_r_trafo_simple, ONLY: replicate_rs_matrices,&
24 : rs_to_kp
25 : USE kpoint_methods, ONLY: kpoint_init_cell_index
26 : USE kpoint_types, ONLY: kpoint_create,&
27 : kpoint_release,&
28 : kpoint_type
29 : USE kubo_bloch, ONLY: bloch_current,&
30 : bloch_current_from_frame
31 : USE message_passing, ONLY: mp_para_env_type
32 : USE qs_environment_types, ONLY: get_qs_env,&
33 : qs_env_release,&
34 : qs_environment_type
35 : USE qs_finite_ao, ONLY: finite_ao_soc
36 : USE qs_gamma2kp, ONLY: create_kp_from_gamma
37 : USE qs_mo_types, ONLY: mo_set_type
38 : USE qs_moments, ONLY: build_local_moment_matrix_rs_img
39 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
40 : USE qs_property_symmetry, ONLY: property_symmetry_frame,&
41 : property_symmetry_plan
42 : #include "./base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 : PRIVATE
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kubo_bloch'
47 : PUBLIC :: setup_bloch_transport
48 : CONTAINS
49 : ! **************************************************************************************************
50 : !> \brief Build a full Gamma-centered property mesh, independent of SCF reduction.
51 : !> \param qs_env Converged GPW/GAPW environment
52 : !> \param grid Full property mesh; nonperiodic directions must be one
53 : !> \param do_soc Add complete GTH SOC after a restricted SCF
54 : !> \param projection Cartesian projection onto periodic transport directions
55 : !> \param max_memory_mb Safety bound on replicated array memory per rank, in MiB
56 : !> \param energies Complete spectrum, flattened (spin,k) channels
57 : !> \param current Anti-Hermitian -i*velocity per channel
58 : !> \param occupations Maximum occupations including normalized k weights
59 : !> \param symmetry Reconstruct symmetry-equivalent property eigenframes
60 : !> \param backend Geometric symmetry backend
61 : !> \param eps_symmetry Geometric tolerance
62 : !> \param eps_operators Target frame and current validation tolerance
63 : !> \note This dense reference stores image operators on every rank; eigensystems are distributed by k.
64 : ! **************************************************************************************************
65 14 : SUBROUTINE setup_bloch_transport(qs_env, grid, do_soc, projection, max_memory_mb, energies, current, occupations, &
66 : symmetry, backend, eps_symmetry, eps_operators)
67 : TYPE(qs_environment_type), POINTER :: qs_env
68 : INTEGER, INTENT(IN) :: grid(3)
69 : LOGICAL, INTENT(IN) :: do_soc
70 : REAL(dp), INTENT(IN) :: projection(3, 3)
71 : INTEGER, INTENT(IN) :: max_memory_mb
72 : REAL(dp), ALLOCATABLE, INTENT(OUT) :: energies(:, :)
73 : COMPLEX(dp), ALLOCATABLE, INTENT(OUT) :: current(:, :, :, :)
74 : REAL(dp), ALLOCATABLE, INTENT(OUT) :: occupations(:)
75 : LOGICAL, INTENT(IN) :: symmetry
76 : INTEGER, INTENT(IN) :: backend
77 : REAL(dp), INTENT(IN) :: eps_symmetry, eps_operators
78 :
79 14 : COMPLEX(dp), ALLOCATABLE :: a(:, :), cur(:, :, :), d(:, :, :), &
80 14 : dh(:, :, :), dsk(:, :, :), h(:, :), &
81 14 : s(:, :), v(:, :, :), source(:, :), TARGET(:, :), &
82 14 : source_current(:, :, :), expected(:, :, :)
83 : INTEGER :: axis, channel, failed, idir, ik, n, nao, &
84 : ni, nk, nspin, spin, status, ir, nr, iw, &
85 : trial, ntry, operation_try, alternatives
86 14 : INTEGER, ALLOCATABLE :: representative(:), operation(:)
87 : INTEGER(int_8) :: bytes
88 : LOGICAL :: do_kpoints, antiunitary
89 : REAL(dp) :: k(3), rotation(3, 3), residuals(2), errors(3), &
90 : candidate_errors(3), best_errors(3)
91 14 : REAL(dp), ALLOCATABLE :: mesh(:, :), representative_k(:, :), source_energies(:)
92 14 : REAL(dp), ALLOCATABLE :: ds(:, :, :, :), hs(:, :, :, :), &
93 14 : ss(:, :, :, :), vs(:, :, :, :)
94 : TYPE(cell_type), POINTER :: cell
95 14 : TYPE(dbcsr_p_type), POINTER :: dr(:, :), hr(:, :), sr(:, :), vr(:, :)
96 : TYPE(dft_control_type), POINTER :: dft_control
97 : TYPE(kpoint_type), POINTER :: images, scf, property_points
98 14 : TYPE(mo_set_type), POINTER :: mos(:)
99 : TYPE(mp_para_env_type), POINTER :: para_env
100 14 : TYPE(neighbor_list_set_p_type), POINTER :: all_pairs(:)
101 : TYPE(qs_environment_type), POINTER :: qp
102 :
103 14 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints, dft_control=dft_control)
104 14 : SELECT CASE (dft_control%qs_control%method_id)
105 : CASE (do_method_gpw, do_method_gapw)
106 : CASE DEFAULT
107 14 : CPABORT('Kubo BLOCH currently supports GPW/GAPW only.')
108 : END SELECT
109 14 : NULLIFY (qp, images, dr, vr, property_points)
110 14 : IF (do_kpoints) THEN
111 12 : qp => qs_env
112 : ELSE
113 2 : CALL create_kp_from_gamma(qs_env, qp)
114 : END IF
115 : CALL get_qs_env(qp, cell=cell, para_env=para_env, matrix_ks_kp=hr, matrix_s_kp=sr, &
116 14 : sab_all=all_pairs, kpoints=scf, mos=mos)
117 14 : CALL dbcsr_get_info(sr(1, 1)%matrix, nfullrows_total=nao)
118 14 : n = nao
119 14 : IF (do_soc) n = 2*nao
120 14 : nspin = SIZE(hr, 1)
121 56 : nk = PRODUCT(grid)
122 14 : CALL kpoint_create(images)
123 14 : CALL kpoint_init_cell_index(images, all_pairs, para_env, ni)
124 : bytes = 8_int_8*nao*nao*ni*(nspin + 4 + MERGE(3, 0, do_soc)) + &
125 14 : 16_int_8*n*n*(3_int_8*nspin*nk + 48_int_8)
126 14 : IF (max_memory_mb < 1 .OR. bytes > INT(max_memory_mb, int_8)*1024*1024) THEN
127 0 : CPABORT('Kubo BLOCH arrays exceed MAX_MEMORY_MB per MPI rank.')
128 : END IF
129 2130 : ALLOCATE (dr(3, ni))
130 14 : CALL build_local_moment_matrix_rs_img(qp, dr)
131 196 : ALLOCATE (hs(nspin, nao, nao, ni), ss(1, nao, nao, ni), ds(3, nao, nao, ni), SOURCE=0.0_dp)
132 14 : CALL replicate_rs_matrices(hr, scf, hs, images%cell_to_index)
133 14 : CALL replicate_rs_matrices(sr, scf, ss, images%cell_to_index)
134 14 : CALL replicate_rs_matrices(dr, images, ds, images%cell_to_index)
135 14 : CALL dbcsr_deallocate_matrix_set(dr)
136 14 : IF (do_soc) THEN
137 10 : CALL finite_ao_soc(qp, vr)
138 40 : ALLOCATE (vs(3, nao, nao, ni), SOURCE=0.0_dp)
139 10 : CALL replicate_rs_matrices(vr, scf, vs, images%cell_to_index)
140 10 : CALL dbcsr_deallocate_matrix_set(vr)
141 : END IF
142 154 : ALLOCATE (energies(n, nspin*nk), current(n, n, 3, nspin*nk), occupations(nspin*nk))
143 0 : ALLOCATE (h(n, n), s(n, n), dh(n, n, 3), dsk(n, n, 3), d(n, n, 3), cur(n, n, 3), &
144 308 : a(nao, nao), v(nao, nao, 3))
145 : ALLOCATE (mesh(3, nk), source(n, n), TARGET(n, n), source_energies(n), &
146 182 : source_current(n, n, 3), expected(n, n, 3))
147 108 : DO ik = 1, nk
148 : mesh(:, ik) = REAL([MOD(ik - 1, grid(1)), MOD((ik - 1)/grid(1), grid(2)), &
149 390 : (ik - 1)/(grid(1)*grid(2))], dp)/REAL(grid, dp)
150 : END DO
151 14 : IF (symmetry) THEN
152 : CALL property_symmetry_plan(qp, mesh, grid, projection, backend, eps_symmetry, &
153 8 : property_points, representative, operation)
154 8 : nr = property_points%nkp
155 24 : ALLOCATE (representative_k(3, nr))
156 120 : representative_k(:, :) = property_points%xkp
157 : ELSE
158 6 : nr = nk
159 12 : ALLOCATE (representative_k(3, nr))
160 142 : representative_k(:, :) = mesh
161 24 : ALLOCATE (representative(nk), operation(nk))
162 74 : representative(:) = [(ik, ik=1, nk)]
163 6 : operation = 0
164 : END IF
165 14 : energies = 0.0_dp
166 14 : current = 0.0_dp
167 14 : failed = 0
168 14 : errors = 0.0_dp
169 14 : alternatives = 0
170 28 : DO spin = 1, nspin
171 108 : DO ik = 1, nk
172 94 : channel = ik + (spin - 1)*nk
173 94 : occupations(channel) = mos(spin)%maxocc/REAL(nk, dp)
174 108 : IF (do_soc) occupations(channel) = 1.0_dp/REAL(nk, dp)
175 : END DO
176 90 : DO ir = 1, nr
177 62 : IF (MOD(ir - 1, para_env%num_pe) /= para_env%mepos) CYCLE
178 124 : k = representative_k(:, ir)
179 31 : CALL operators_at_k()
180 31 : CALL bloch_current(h, s, dh, dsk, d, source_energies, source_current, status, source)
181 31 : IF (status /= 0) THEN
182 0 : failed = MAX(failed, ABS(status))
183 0 : CYCLE
184 : END IF
185 289 : DO ik = 1, nk
186 213 : IF (representative(ik) /= ir) CYCLE
187 47 : channel = ik + (spin - 1)*nk
188 1477 : energies(:, channel) = source_energies
189 172802 : cur(:, :, :) = source_current
190 47 : IF (operation(ik) /= 0) THEN
191 64 : k = mesh(:, ik)
192 16 : CALL operators_at_k()
193 16 : ntry = property_points%kp_sym(ir)%kpoint_sym%nwred
194 64 : best_errors = HUGE(1.0_dp)
195 16 : DO trial = 0, ntry
196 16 : operation_try = operation(ik)
197 16 : IF (trial > 0) operation_try = trial
198 0 : IF (trial > 0 .AND. operation_try == operation(ik)) CYCLE
199 : CALL property_symmetry_frame(property_points, ir, operation_try, k, source, do_soc, TARGET, &
200 16 : rotation, antiunitary, status)
201 16 : IF (status /= 0) CYCLE
202 208 : IF (MAXVAL(ABS(MATMUL(rotation, MATMUL(projection, TRANSPOSE(rotation))) - projection)) &
203 1264 : > 1.e-9_dp) CYCLE
204 : CALL bloch_current_from_frame(h, s, dh, dsk, d, source_energies, TARGET, cur, &
205 16 : eps_operators, residuals, status)
206 96 : candidate_errors = [residuals, 0.0_dp]
207 16 : IF (status == 0) THEN
208 16 : expected = 0.0_dp
209 64 : DO axis = 1, 3
210 208 : DO idir = 1, 3
211 192 : IF (antiunitary) THEN
212 : expected(:, :, axis) = expected(:, :, axis) + &
213 12654 : rotation(axis, idir)*CONJG(source_current(:, :, idir))
214 : ELSE
215 254718 : expected(:, :, axis) = expected(:, :, axis) + rotation(axis, idir)*source_current(:, :, idir)
216 : END IF
217 : END DO
218 : END DO
219 178280 : candidate_errors(3) = MAXVAL(ABS(cur - expected))/MAX(1.0_dp, MAXVAL(ABS(cur)))
220 : END IF
221 112 : IF (MAXVAL(candidate_errors) < MAXVAL(best_errors)) best_errors = candidate_errors
222 64 : IF (status /= 0 .OR. MAXVAL(candidate_errors) > eps_operators) CYCLE
223 16 : best_errors = candidate_errors
224 16 : IF (trial > 0) alternatives = alternatives + 1
225 48 : EXIT
226 : END DO
227 64 : errors = MAX(errors, best_errors)
228 16 : IF (trial > ntry) THEN
229 0 : failed = MAX(failed, 5)
230 0 : CYCLE
231 : END IF
232 : END IF
233 250 : DO axis = 1, 3
234 777 : DO idir = 1, 3
235 518406 : current(:, :, axis, channel) = current(:, :, axis, channel) + projection(axis, idir)*cur(:, :, idir)
236 : END DO
237 : END DO
238 : END DO
239 : END DO
240 : END DO
241 14 : CALL para_env%max(failed)
242 14 : CALL para_env%max(errors)
243 14 : CALL para_env%sum(alternatives)
244 14 : iw = cp_logger_get_default_io_unit()
245 14 : IF (symmetry .AND. iw > 0) THEN
246 4 : WRITE (iw, '(T2,A,I0,A,I0)') 'KUBO_TRANSPORT| Property symmetry representatives: ', nr, ' / ', nk
247 4 : WRITE (iw, '(T2,A,3ES14.5)') 'KUBO_TRANSPORT| Symmetry errors S/H/current: ', errors
248 4 : WRITE (iw, '(T2,A,I0)') 'KUBO_TRANSPORT| Alternative validated orbit operations: ', alternatives
249 : END IF
250 14 : IF (failed /= 0) THEN
251 0 : IF (symmetry) THEN
252 0 : CPABORT('Bloch Kubo property symmetry failed target metric/Hamiltonian/current covariance.')
253 : ELSE
254 0 : CPABORT('Bloch Kubo failed metric, Hermiticity or AO-connection validation.')
255 : END IF
256 : END IF
257 14 : CALL para_env%sum(energies)
258 14 : CALL para_env%sum(current)
259 14 : CALL kpoint_release(images)
260 14 : IF (ASSOCIATED(property_points)) CALL kpoint_release(property_points)
261 58 : IF (.NOT. do_kpoints) THEN
262 2 : CALL qs_env_release(qp)
263 2 : DEALLOCATE (qp)
264 : END IF
265 : CONTAINS
266 : ! **************************************************************************************************
267 : !> \brief Fourier-transform the frozen image operators at the current (k,spin).
268 : ! **************************************************************************************************
269 47 : SUBROUTINE operators_at_k()
270 : INTEGER :: component, direction
271 :
272 47 : CALL rs_to_kp(hs(spin, :, :, :), a, images%index_to_cell, k)
273 47 : CALL spin_diagonal(a, h)
274 47 : IF (do_soc) THEN
275 124 : DO component = 1, 3
276 124 : CALL rs_to_kp(vs(component, :, :, :), v(:, :, component), images%index_to_cell, k)
277 : END DO
278 31 : CALL add_soc(h, v)
279 : END IF
280 47 : CALL rs_to_kp(ss(1, :, :, :), a, images%index_to_cell, k)
281 47 : CALL spin_diagonal(a, s)
282 188 : DO component = 1, 3
283 141 : CALL rs_to_kp(hs(spin, :, :, :), a, images%index_to_cell, k, component, cell%hmat)
284 141 : CALL spin_diagonal(a, dh(:, :, component))
285 141 : IF (do_soc) THEN
286 372 : DO direction = 1, 3
287 372 : CALL rs_to_kp(vs(direction, :, :, :), v(:, :, direction), images%index_to_cell, k, component, cell%hmat)
288 : END DO
289 93 : CALL add_soc(dh(:, :, component), v)
290 : END IF
291 141 : CALL rs_to_kp(ss(1, :, :, :), a, images%index_to_cell, k, component, cell%hmat)
292 141 : CALL spin_diagonal(a, dsk(:, :, component))
293 141 : CALL rs_to_kp(ds(component, :, :, :), a, images%index_to_cell, k)
294 188 : CALL spin_diagonal(a, d(:, :, component))
295 : END DO
296 47 : END SUBROUTINE operators_at_k
297 : END SUBROUTINE setup_bloch_transport
298 : ! **************************************************************************************************
299 : !> \brief Duplicate a scalar operator in the physical spin basis.
300 : !> \param a Scalar operator
301 : !> \param out Scalar or doubled spinor operator
302 : ! **************************************************************************************************
303 517 : SUBROUTINE spin_diagonal(a, out)
304 : COMPLEX(dp), INTENT(IN) :: a(:, :)
305 : COMPLEX(dp), INTENT(OUT) :: out(:, :)
306 :
307 : INTEGER :: n
308 :
309 517 : n = SIZE(a, 1)
310 633435 : out = 0.0_dp
311 186131 : out(:n, :n) = a
312 154099 : IF (SIZE(out, 1) == 2*n) out(n + 1:, n + 1:) = a
313 517 : END SUBROUTINE spin_diagonal
314 : ! **************************************************************************************************
315 : !> \brief Add i*sigma times the real-spherical-harmonic SOC components, including their derivatives.
316 : !> \param h Spinor Hamiltonian or derivative
317 : !> \param v Fourier-transformed real antisymmetric SOC components (without the factor i)
318 : ! **************************************************************************************************
319 124 : SUBROUTINE add_soc(h, v)
320 : COMPLEX(dp), INTENT(INOUT) :: h(:, :)
321 : COMPLEX(dp), INTENT(IN) :: v(:, :, :)
322 :
323 : COMPLEX(dp), PARAMETER :: iu = (0.0_dp, 1.0_dp)
324 :
325 : INTEGER :: n
326 :
327 124 : n = SIZE(v, 1)
328 55972 : h(:n, :n) = h(:n, :n) + iu*v(:, :, 3)
329 55972 : h(n + 1:, n + 1:) = h(n + 1:, n + 1:) - iu*v(:, :, 3)
330 55972 : h(:n, n + 1:) = h(:n, n + 1:) + iu*v(:, :, 1) - v(:, :, 2)
331 55972 : h(n + 1:, :n) = h(n + 1:, :n) + iu*v(:, :, 1) + v(:, :, 2)
332 124 : END SUBROUTINE add_soc
333 : END MODULE qs_kubo_bloch
|
