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-volume Kubo-Greenwood transport from converged Quickstep matrices.
10 : ! **************************************************************************************************
11 : MODULE qs_kubo_transport
12 : USE bibliography, ONLY: EsteveParedes2023,&
13 : KuhneHeskeProdan2020,&
14 : cite_reference
15 : USE cell_types, ONLY: cell_type,&
16 : pbc
17 : USE cp_blacs_env, ONLY: cp_blacs_env_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
19 : dbcsr_p_type,&
20 : dbcsr_type
21 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
22 : dbcsr_deallocate_matrix_set
23 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
24 : cp_fm_struct_release,&
25 : cp_fm_struct_type
26 : USE cp_fm_types, ONLY: cp_fm_create,&
27 : cp_fm_get_submatrix,&
28 : cp_fm_release,&
29 : cp_fm_type
30 : USE cp_log_handling, ONLY: cp_get_default_logger,&
31 : cp_logger_get_default_io_unit,&
32 : cp_logger_type
33 : USE ieee_arithmetic, ONLY: ieee_is_finite
34 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 : section_vals_type,&
36 : section_vals_val_get
37 : USE kinds, ONLY: dp
38 : USE kubo_bloch, ONLY: fermi_slope
39 : USE kubo_hall, ONLY: hall_response
40 : USE kubo_projected, ONLY: projected_current
41 : USE mathconstants, ONLY: pi
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE particle_types, ONLY: particle_type
44 : USE physcon, ONLY: a_bohr,&
45 : angstrom,&
46 : e_charge,&
47 : evolt,&
48 : h_bar,&
49 : kelvin
50 : USE qs_environment_types, ONLY: get_qs_env,&
51 : qs_environment_type
52 : USE qs_finite_ao, ONLY: finite_ao_positions,&
53 : finite_ao_soc
54 : USE qs_kubo_bloch, ONLY: setup_bloch_transport
55 : USE qs_mo_types, ONLY: get_mo_set,&
56 : mo_set_type
57 : USE spectral_localizer_dbcsr, ONLY: localizer_ao_operators,&
58 : localizer_dbcsr_dense,&
59 : localizer_dbcsr_release,&
60 : localizer_matrix_type
61 : USE string_utilities, ONLY: uppercase
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 :
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kubo_transport'
69 :
70 : PUBLIC :: qs_scf_post_kubo_transport
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief Share full finite AO/SOC operators with the localizer, retaining complex matrix elements.
76 : !> \param qs_env Converged finite Quickstep environment
77 : !> \param spin Scalar spin channel, one for SOC
78 : !> \param do_soc Include GTH SOC
79 : !> \param energies Complete scalar/spinor eigenvalues
80 : !> \param current Three projected commutators in the eigenbasis
81 : ! **************************************************************************************************
82 8 : SUBROUTINE setup_projected_transport(qs_env, spin, do_soc, energies, current)
83 : TYPE(qs_environment_type), POINTER :: qs_env
84 : INTEGER, INTENT(IN) :: spin
85 : LOGICAL, INTENT(IN) :: do_soc
86 : REAL(dp), INTENT(OUT) :: energies(:)
87 : COMPLEX(dp), INTENT(OUT) :: current(:, :, :)
88 :
89 16 : COMPLEX(dp), ALLOCATABLE :: component(:, :), hd(:, :), sd(:, :), &
90 : xd(:, :, :)
91 : INTEGER :: axis, n, status
92 16 : TYPE(dbcsr_p_type), POINTER :: h(:), positions(:), s(:), soc(:, :)
93 40 : TYPE(localizer_matrix_type) :: metric, operators(4)
94 : TYPE(mp_para_env_type), POINTER :: para_env
95 :
96 8 : NULLIFY (positions, soc)
97 8 : CALL get_qs_env(qs_env, matrix_ks=h, matrix_s=s, para_env=para_env)
98 8 : CALL finite_ao_positions(qs_env, positions)
99 8 : IF (do_soc) THEN
100 2 : CALL finite_ao_soc(qs_env, soc)
101 2 : CALL localizer_ao_operators(h(spin)%matrix, s(1)%matrix, positions, operators, metric, soc)
102 : ELSE
103 6 : CALL localizer_ao_operators(h(spin)%matrix, s(1)%matrix, positions, operators, metric)
104 : END IF
105 8 : CALL localizer_dbcsr_dense(operators(1), para_env, hd)
106 8 : CALL localizer_dbcsr_dense(metric, para_env, sd)
107 8 : n = SIZE(hd, 1)
108 40 : ALLOCATE (xd(n, n, 3))
109 32 : DO axis = 1, 3
110 24 : CALL localizer_dbcsr_dense(operators(axis + 1), para_env, component)
111 6216 : xd(:, :, axis) = component
112 32 : DEALLOCATE (component)
113 : END DO
114 8 : status = 0
115 8 : IF (para_env%is_source()) CALL projected_current(hd, sd, xd, energies, current, status)
116 8 : CALL para_env%bcast(status, para_env%source)
117 8 : IF (status /= 0) THEN
118 0 : CPABORT("Kubo projected AO matrices or eigensystem invalid.")
119 : END IF
120 232 : CALL para_env%bcast(energies, para_env%source)
121 12440 : CALL para_env%bcast(current, para_env%source)
122 8 : CALL dbcsr_deallocate_matrix_set(positions)
123 8 : IF (do_soc) CALL dbcsr_deallocate_matrix_set(soc)
124 40 : DO axis = 1, 4
125 40 : CALL localizer_dbcsr_release(operators(axis))
126 : END DO
127 8 : CALL localizer_dbcsr_release(metric)
128 24 : END SUBROUTINE setup_projected_transport
129 :
130 : ! **************************************************************************************************
131 : !> \brief Compute and print the finite-volume Kubo-Greenwood conductivity tensor.
132 : !> \param qs_env Quickstep environment
133 : ! **************************************************************************************************
134 14227 : SUBROUTINE qs_scf_post_kubo_transport(qs_env)
135 : TYPE(qs_environment_type), POINTER :: qs_env
136 :
137 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_kubo_transport'
138 :
139 : CHARACTER(LEN=32) :: density_unit, measure_unit, method, &
140 : periodic_label, sigma_iso_unit, &
141 : sigma_tensor_unit
142 : CHARACTER(LEN=512) :: header
143 : COMPLEX(KIND=dp), ALLOCATABLE, &
144 14227 : DIMENSION(:, :, :, :) :: dHH
145 : INTEGER :: axis, current_operator, hall_status, handle, imu, ispin, max_ao, max_memory_mb, &
146 : nao, natom, nelectron_total, neutral_grid, nk, nmu, nperiodic, nspins, nstate, &
147 : output_unit, symmetry_backend, transport_ndim
148 14227 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_atom
149 14227 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
150 14227 : INTEGER, POINTER :: grid(:)
151 : LOGICAL :: allow_mo_reuse, do_hall, do_kpoints, &
152 : do_soc, kubo_active, &
153 : neutral_mu_explicit, ot_active, &
154 : property_symmetry
155 14227 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: reused_mos
156 14227 : REAL(dp), ALLOCATABLE :: scalar_current(:, :, :)
157 : REAL(KIND=dp) :: density_factor, dissipation, emax, emin, eps_operators, eps_symmetry, &
158 : iso_factor, measure, measure_m, mu, mu0, mu_step, ne, nh, sigma_iso, temperature, &
159 : tensor_factor
160 14227 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: maxocc
161 14227 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eig, s_dense
162 : REAL(KIND=dp), DIMENSION(3, 3) :: hall, hall_part, projection, sigma, &
163 : sigma_out
164 14227 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy_range
165 : TYPE(cell_type), POINTER :: cell
166 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
167 : TYPE(cp_logger_type), POINTER :: logger
168 14227 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
169 14227 : TYPE(dbcsr_p_type), POINTER :: images(:, :)
170 14227 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
171 : TYPE(mp_para_env_type), POINTER :: para_env
172 14227 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
173 : TYPE(section_vals_type), POINTER :: input_section, kubo_section, ot_section
174 :
175 14227 : CALL timeset(routineN, handle)
176 :
177 14227 : NULLIFY (blacs_env, cell, energy_range, input_section, kubo_section, logger, matrix_ks, matrix_s, &
178 14227 : mos, ot_section, para_env, particle_set, row_blk_sizes)
179 :
180 14227 : CALL get_qs_env(qs_env, input=input_section)
181 14227 : kubo_section => section_vals_get_subs_vals(input_section, "PROPERTIES%KUBO_TRANSPORT")
182 14227 : CALL section_vals_val_get(kubo_section, "_SECTION_PARAMETERS_", l_val=kubo_active)
183 14227 : IF (.NOT. kubo_active) THEN
184 14197 : CALL timestop(handle)
185 14197 : RETURN
186 : END IF
187 30 : CALL cite_reference(KuhneHeskeProdan2020)
188 :
189 30 : CALL section_vals_val_get(kubo_section, "METHOD", c_val=method)
190 30 : CALL uppercase(method)
191 30 : SELECT CASE (TRIM(method))
192 : CASE ("DIAGONALIZATION")
193 : CASE ("TD", "CHEBYSHEV")
194 0 : CPABORT("KUBO_TRANSPORT METHOD "//TRIM(method)//" is reserved but not implemented yet.")
195 : CASE DEFAULT
196 30 : CPABORT("Unknown KUBO_TRANSPORT METHOD: "//TRIM(method))
197 : END SELECT
198 30 : ot_section => section_vals_get_subs_vals(input_section, "DFT%SCF%OT")
199 30 : CALL section_vals_val_get(ot_section, "_SECTION_PARAMETERS_", l_val=ot_active)
200 30 : allow_mo_reuse = .NOT. ot_active
201 :
202 30 : CALL section_vals_val_get(kubo_section, "CURRENT_OPERATOR", i_val=current_operator)
203 30 : CALL section_vals_val_get(kubo_section, "SOC", l_val=do_soc)
204 30 : CALL section_vals_val_get(kubo_section, "HALL_RESPONSE", l_val=do_hall)
205 30 : IF (do_hall .AND. current_operator == 0) THEN
206 0 : CPABORT("Kubo HALL_RESPONSE requires CURRENT_OPERATOR PROJECTED_AO or BLOCH.")
207 : END IF
208 30 : CALL section_vals_val_get(kubo_section, "SYMMETRY", l_val=property_symmetry)
209 30 : IF (property_symmetry .AND. current_operator /= 2) THEN
210 0 : CPABORT('Kubo property SYMMETRY requires BLOCH.')
211 : END IF
212 30 : CALL get_qs_env(qs_env, do_kpoints=do_kpoints, cell=cell)
213 30 : IF (do_kpoints .AND. current_operator /= 2) THEN
214 0 : CPABORT("Kubo explicit k-points require CURRENT_OPERATOR BLOCH.")
215 : END IF
216 30 : IF (do_soc .AND. current_operator == 0) THEN
217 0 : CPABORT("Kubo SOC requires CURRENT_OPERATOR PROJECTED_AO or BLOCH.")
218 : END IF
219 66 : IF (current_operator == 1 .AND. ANY(cell%perd /= 0)) THEN
220 0 : CPABORT("Kubo PROJECTED_AO requires PERIODIC NONE.")
221 : END IF
222 : CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, mos=mos, natom=natom, &
223 30 : nelectron_total=nelectron_total, para_env=para_env, particle_set=particle_set)
224 30 : CALL get_qs_env(qs_env, matrix_ks_kp=images)
225 30 : matrix_ks => images(:, 1)
226 30 : CALL get_qs_env(qs_env, matrix_s_kp=images)
227 30 : matrix_s => images(:, 1)
228 30 : CPASSERT(ASSOCIATED(matrix_s))
229 30 : CPASSERT(ASSOCIATED(matrix_ks))
230 30 : CPASSERT(ASSOCIATED(cell))
231 30 : CPASSERT(ASSOCIATED(particle_set))
232 :
233 30 : CALL section_vals_val_get(kubo_section, "TEMPERATURE", r_val=temperature)
234 30 : CALL section_vals_val_get(kubo_section, "DISSIPATION", r_val=dissipation)
235 30 : CALL section_vals_val_get(kubo_section, "ENERGY_RANGE", r_vals=energy_range)
236 30 : CALL section_vals_val_get(kubo_section, "NEUTRAL_MU", explicit=neutral_mu_explicit)
237 30 : IF (neutral_mu_explicit) CALL section_vals_val_get(kubo_section, "NEUTRAL_MU", r_val=mu0)
238 30 : CALL section_vals_val_get(kubo_section, "N_MU", i_val=nmu)
239 30 : CALL section_vals_val_get(kubo_section, "NEUTRAL_GRID", i_val=neutral_grid)
240 :
241 30 : IF (temperature <= 0.0_dp) CPABORT("KUBO_TRANSPORT TEMPERATURE must be positive.")
242 30 : IF (dissipation <= 0.0_dp) CPABORT("KUBO_TRANSPORT DISSIPATION must be positive.")
243 30 : IF (nmu < 1) CPABORT("KUBO_TRANSPORT N_MU must be at least one.")
244 30 : IF (neutral_grid < 10) CPABORT("KUBO_TRANSPORT NEUTRAL_GRID must be at least ten.")
245 :
246 30 : nspins = SIZE(matrix_ks)
247 30 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao, row_blk_size=row_blk_sizes)
248 30 : nstate = nao
249 30 : IF (current_operator /= 0) THEN
250 20 : CALL section_vals_val_get(kubo_section, "MAX_AO", i_val=max_ao)
251 20 : IF (max_ao < 1 .OR. nao > max_ao) CPABORT("Kubo AO current exceeds MAX_AO.")
252 20 : IF (.NOT. qs_env%scf_convergence_available .OR. .NOT. qs_env%scf_converged) THEN
253 0 : CPABORT("Kubo AO current requires a converged SCF.")
254 : END IF
255 : END IF
256 30 : IF (do_soc) THEN
257 12 : IF (nspins /= 1) CPABORT("Kubo SOC requires a restricted SCF.")
258 12 : nstate = 2*nao
259 : END IF
260 30 : CALL build_ao_atom_map(row_blk_sizes, natom, ao_atom)
261 30 : IF (current_operator == 0) CALL dbcsr_to_dense(matrix_s(1)%matrix, blacs_env, para_env, s_dense)
262 : CALL setup_transport_geometry(cell, transport_ndim, nperiodic, periodic_label, projection, &
263 30 : measure, measure_unit)
264 30 : CALL transport_unit_labels(transport_ndim, density_unit, sigma_iso_unit, sigma_tensor_unit)
265 :
266 30 : IF (current_operator == 2) THEN
267 14 : CALL cite_reference(EsteveParedes2023)
268 14 : CALL section_vals_val_get(kubo_section, "MP_GRID", i_vals=grid)
269 14 : CALL section_vals_val_get(kubo_section, "MAX_MEMORY_MB", i_val=max_memory_mb)
270 14 : CALL section_vals_val_get(kubo_section, "SYMMETRY_BACKEND", i_val=symmetry_backend)
271 14 : CALL section_vals_val_get(kubo_section, "EPS_SYMMETRY", r_val=eps_symmetry)
272 14 : CALL section_vals_val_get(kubo_section, "EPS_SYMMETRY_OPERATORS", r_val=eps_operators)
273 14 : IF (.NOT. ieee_is_finite(eps_symmetry) .OR. eps_symmetry <= 0.0_dp .OR. eps_symmetry >= 1.0_dp) THEN
274 0 : CPABORT('Invalid Kubo EPS_SYMMETRY.')
275 : END IF
276 14 : IF (.NOT. ieee_is_finite(eps_operators) .OR. eps_operators <= 0.0_dp .OR. eps_operators >= 1.0_dp) THEN
277 0 : CPABORT('Invalid Kubo EPS_SYMMETRY_OPERATORS.')
278 : END IF
279 56 : IF (COUNT(cell%perd /= 0) == 0) CPABORT("Kubo BLOCH requires a periodic cell.")
280 98 : IF (ANY(grid < 1) .OR. ANY(grid /= 1 .AND. cell%perd == 0)) CPABORT("Invalid Kubo MP_GRID.")
281 : nk = 1
282 56 : DO axis = 1, 3
283 42 : IF (grid(axis) > 1000000/nk) CPABORT("Kubo MP_GRID exceeds one million points.")
284 56 : nk = nk*grid(axis)
285 : END DO
286 : CALL setup_bloch_transport(qs_env, grid, do_soc, projection, max_memory_mb, eig, dHH, maxocc, &
287 14 : property_symmetry, symmetry_backend, eps_symmetry, eps_operators)
288 56 : ALLOCATE (reused_mos(SIZE(eig, 2)), SOURCE=.FALSE.)
289 : ELSE
290 208 : ALLOCATE (eig(nstate, nspins), dHH(nstate, nstate, 3, nspins), maxocc(nspins), reused_mos(nspins))
291 34 : DO ispin = 1, nspins
292 18 : maxocc(ispin) = mos(ispin)%maxocc
293 18 : IF (do_soc) maxocc(ispin) = 1.0_dp
294 34 : IF (current_operator == 1) THEN
295 8 : CALL setup_projected_transport(qs_env, ispin, do_soc, eig(:, ispin), dHH(:, :, :, ispin))
296 8 : reused_mos(ispin) = .FALSE.
297 : ELSE
298 50 : ALLOCATE (scalar_current(nao, nao, 3))
299 : CALL setup_spin_transport(matrix_ks(ispin)%matrix, s_dense, blacs_env, para_env, &
300 : particle_set, cell, projection, ao_atom, mos(ispin), allow_mo_reuse, &
301 10 : eig(:, ispin), scalar_current, reused_mos(ispin))
302 3340 : dHH(:, :, :, ispin) = CMPLX(scalar_current, 0.0_dp, dp)
303 10 : DEALLOCATE (scalar_current)
304 : END IF
305 : END DO
306 : END IF
307 :
308 3214 : emin = MINVAL(eig)
309 3214 : emax = MAXVAL(eig)
310 30 : IF (.NOT. neutral_mu_explicit) THEN
311 : CALL find_neutral_mu(eig, maxocc, temperature, REAL(nelectron_total, KIND=dp), neutral_grid, &
312 0 : emin, emax, mu0)
313 : END IF
314 :
315 30 : IF (energy_range(1) < energy_range(2)) THEN
316 30 : emin = energy_range(1)
317 30 : emax = energy_range(2)
318 : END IF
319 :
320 30 : measure_m = measure*a_bohr**transport_ndim
321 30 : density_factor = 1.0_dp/measure_m
322 30 : CALL transport_output_factors(transport_ndim, iso_factor, tensor_factor)
323 :
324 30 : logger => cp_get_default_logger()
325 30 : output_unit = cp_logger_get_default_io_unit(logger)
326 :
327 30 : IF (output_unit > 0) THEN
328 15 : WRITE (output_unit, '(/,T2,A)') "Kubo-Greenwood finite-volume transport"
329 15 : WRITE (output_unit, '(T3,A,T38,A)') "Method:", TRIM(method)
330 15 : IF (current_operator == 1) THEN
331 3 : WRITE (output_unit, '(T3,A)') "Current: metric-aware projected AO commutator, not continuum velocity."
332 3 : WRITE (output_unit, '(T3,A,L1)') "GTH SOC in the complete AO spinor space: ", do_soc
333 12 : ELSE IF (current_operator == 2) THEN
334 7 : WRITE (output_unit, '(T3,A)') "Current: Bloch H/S derivatives and analytic AO connection."
335 28 : WRITE (output_unit, '(T3,A,3I6)') "Full property mesh: ", grid
336 7 : WRITE (output_unit, '(T3,A,L1)') "GTH SOC in the complete AO spinor space: ", do_soc
337 : END IF
338 15 : IF (current_operator /= 0) THEN
339 10 : IF (do_hall) THEN
340 6 : WRITE (output_unit, '(T3,A)') "Dissipative and antisymmetric DC charge response; scalar relaxation rate."
341 6 : WRITE (output_unit, '(T3,A)') "Tensor convention: first index current, second index electric field."
342 : ELSE
343 4 : WRITE (output_unit, '(T3,A)') "Dissipative symmetric conductivity only; no Hall response."
344 : END IF
345 : END IF
346 15 : WRITE (output_unit, '(T3,A,T38,A)') "Cell periodicity:", TRIM(periodic_label)
347 15 : IF (nperiodic == 0) THEN
348 6 : WRITE (output_unit, '(T3,A,T38,A)') "Transport normalization:", "3D finite box"
349 : ELSE
350 9 : WRITE (output_unit, '(T3,A,T38,I12)') "Transport dimensionality:", transport_ndim
351 : END IF
352 16 : IF (ALL(reused_mos)) THEN
353 1 : WRITE (output_unit, '(T3,A,T38,A)') "Eigensystem source:", "CP2K canonical MOs"
354 69 : ELSE IF (ANY(reused_mos)) THEN
355 0 : WRITE (output_unit, '(T3,A,T38,A)') "Eigensystem source:", "mixed canonical/internal"
356 : ELSE
357 14 : WRITE (output_unit, '(T3,A,T38,A)') "Eigensystem source:", "internal diagonalization"
358 : END IF
359 15 : WRITE (output_unit, '(T3,A,T38,I12)') "Number of atoms:", natom
360 15 : WRITE (output_unit, '(T3,A,T38,I12)') "Number of AO functions:", nao
361 15 : WRITE (output_unit, '(T3,A,T38,F12.4)') "Temperature [K]:", temperature*kelvin
362 15 : WRITE (output_unit, '(T3,A,T38,F12.4)') "Dissipation [K]:", dissipation*kelvin
363 15 : WRITE (output_unit, '(T3,A,T38,F12.6)') "Dissipation [eV]:", dissipation*evolt
364 15 : WRITE (output_unit, '(T3,A,T38,F12.6)') "Neutral chemical potential [eV]:", mu0*evolt
365 15 : WRITE (output_unit, '(T3,A,T38,ES12.4,1X,A)') "Transport measure:", &
366 30 : measure*angstrom**transport_ndim, TRIM(measure_unit)
367 : header = "mu[Ha] mu-mu0[eV] Ne["//TRIM(density_unit)//"] "// &
368 : "Nh["//TRIM(density_unit)//"] sigma_iso["//TRIM(sigma_iso_unit)//"] "// &
369 : "sigma_xx["//TRIM(sigma_tensor_unit)//"] sigma_xy["// &
370 : TRIM(sigma_tensor_unit)//"] sigma_xz["//TRIM(sigma_tensor_unit)// &
371 : "] sigma_yx["//TRIM(sigma_tensor_unit)//"] sigma_yy["// &
372 : TRIM(sigma_tensor_unit)//"] sigma_yz["//TRIM(sigma_tensor_unit)// &
373 : "] sigma_zx["//TRIM(sigma_tensor_unit)//"] sigma_zy["// &
374 15 : TRIM(sigma_tensor_unit)//"] sigma_zz["//TRIM(sigma_tensor_unit)//"]"
375 15 : WRITE (output_unit, '(/,T3,A)') TRIM(header)
376 : END IF
377 :
378 30 : IF (nmu == 1) THEN
379 : mu_step = 0.0_dp
380 : ELSE
381 0 : mu_step = (emax - emin)/REAL(nmu - 1, KIND=dp)
382 : END IF
383 :
384 60 : DO imu = 1, nmu
385 30 : mu = emin + REAL(imu - 1, KIND=dp)*mu_step
386 30 : CALL electron_hole_density(eig, maxocc, mu, mu0, temperature, ne, nh)
387 : CALL conductivity_at_mu(eig, dHH, maxocc, mu, temperature, dissipation, measure, &
388 44 : transport_ndim, sigma, current_operator == 2 .OR. do_hall)
389 30 : IF (do_hall) THEN
390 12 : hall = 0.0_dp
391 76 : DO ispin = 1, SIZE(eig, 2)
392 : CALL hall_response(eig(:, ispin), dHH(:, :, :, ispin), mu, temperature, dissipation, &
393 64 : hall_part, hall_status)
394 64 : IF (hall_status /= 0) THEN
395 0 : CPABORT("Invalid Kubo Hall current, spectrum, or relaxation parameters.")
396 : END IF
397 844 : hall = hall + maxocc(ispin)*hall_part
398 : END DO
399 156 : hall = (e_charge**2/h_bar)*hall*angstrom**(2 - transport_ndim)/measure
400 156 : sigma = sigma + hall
401 : END IF
402 90 : IF (output_unit > 0) THEN
403 : sigma_iso = (sigma(1, 1) + sigma(2, 2) + sigma(3, 3))/ &
404 15 : REAL(transport_ndim, KIND=dp)*iso_factor
405 195 : sigma_out(:, :) = sigma(:, :)*tensor_factor
406 15 : WRITE (output_unit, '(T3,F12.6,F14.6,2ES16.7,10ES17.8)') mu, (mu - mu0)*evolt, &
407 15 : ne*density_factor, nh*density_factor, &
408 15 : sigma_iso, &
409 15 : sigma_out(1, 1), sigma_out(1, 2), sigma_out(1, 3), &
410 15 : sigma_out(2, 1), sigma_out(2, 2), sigma_out(2, 3), &
411 30 : sigma_out(3, 1), sigma_out(3, 2), sigma_out(3, 3)
412 2 : SELECT CASE (transport_ndim)
413 : CASE (1)
414 2 : WRITE (output_unit, '(T3,A,1X,ES17.8)') "KUBO_TRANSPORT| sigma_iso[S*m]", sigma_iso
415 : CASE (2)
416 3 : WRITE (output_unit, '(T3,A,1X,ES17.8)') "KUBO_TRANSPORT| sigma_iso[S]", sigma_iso
417 : CASE DEFAULT
418 15 : WRITE (output_unit, '(T3,A,1X,ES17.8)') "KUBO_TRANSPORT| sigma_iso[S/cm]", sigma_iso
419 : END SELECT
420 15 : IF (do_hall) THEN
421 6 : WRITE (output_unit, '(T3,A,1X,A,3ES17.8)') "KUBO_TRANSPORT| hall_xy_yz_zx", &
422 6 : "["//TRIM(sigma_tensor_unit)//"]", &
423 12 : tensor_factor*hall(1, 2), tensor_factor*hall(2, 3), tensor_factor*hall(3, 1)
424 6 : WRITE (output_unit, '(T3,A,1X,ES17.8)') "KUBO_TRANSPORT| hall_max", &
425 84 : tensor_factor*MAXVAL(ABS(hall))
426 : END IF
427 : END IF
428 : END DO
429 :
430 30 : DEALLOCATE (ao_atom, dHH, eig, maxocc, reused_mos)
431 30 : IF (ALLOCATED(s_dense)) DEALLOCATE (s_dense)
432 :
433 30 : CALL timestop(handle)
434 :
435 28454 : END SUBROUTINE qs_scf_post_kubo_transport
436 :
437 : ! **************************************************************************************************
438 : !> \brief Build a global AO -> atom map from DBCSR atomic block sizes.
439 : !> \param row_blk_sizes ...
440 : !> \param natom ...
441 : !> \param ao_atom ...
442 : ! **************************************************************************************************
443 30 : SUBROUTINE build_ao_atom_map(row_blk_sizes, natom, ao_atom)
444 : INTEGER, DIMENSION(:), INTENT(IN) :: row_blk_sizes
445 : INTEGER, INTENT(IN) :: natom
446 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ao_atom
447 :
448 : INTEGER :: iao, iatom, nao, u
449 :
450 78 : nao = SUM(row_blk_sizes(1:natom))
451 90 : ALLOCATE (ao_atom(nao))
452 30 : iao = 0
453 78 : DO iatom = 1, natom
454 478 : DO u = 1, row_blk_sizes(iatom)
455 400 : iao = iao + 1
456 448 : ao_atom(iao) = iatom
457 : END DO
458 : END DO
459 :
460 30 : END SUBROUTINE build_ao_atom_map
461 :
462 : ! **************************************************************************************************
463 : !> \brief Copy a DBCSR matrix to a replicated dense matrix.
464 : !> \param matrix ...
465 : !> \param blacs_env ...
466 : !> \param para_env ...
467 : !> \param dense ...
468 : ! **************************************************************************************************
469 20 : SUBROUTINE dbcsr_to_dense(matrix, blacs_env, para_env, dense)
470 : TYPE(dbcsr_type), INTENT(IN) :: matrix
471 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
472 : TYPE(mp_para_env_type), POINTER :: para_env
473 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
474 : INTENT(OUT) :: dense
475 :
476 : INTEGER :: ncol, nrow
477 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
478 : TYPE(cp_fm_type) :: fm
479 :
480 20 : NULLIFY (fm_struct)
481 :
482 20 : CALL dbcsr_get_info(matrix, nfullrows_total=nrow, nfullcols_total=ncol)
483 80 : ALLOCATE (dense(nrow, ncol))
484 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
485 20 : nrow_global=nrow, ncol_global=ncol)
486 20 : CALL cp_fm_create(fm, fm_struct)
487 20 : CALL copy_dbcsr_to_fm(matrix, fm)
488 20 : CALL cp_fm_get_submatrix(fm, dense, n_rows=nrow, n_cols=ncol)
489 20 : CALL cp_fm_release(fm)
490 20 : CALL cp_fm_struct_release(fm_struct)
491 :
492 40 : END SUBROUTINE dbcsr_to_dense
493 :
494 : ! **************************************************************************************************
495 : !> \brief Set up the periodic transport subspace and its normalization measure.
496 : !> \param cell ...
497 : !> \param transport_ndim ...
498 : !> \param nperiodic ...
499 : !> \param periodic_label ...
500 : !> \param projection ...
501 : !> \param measure ...
502 : !> \param measure_unit ...
503 : ! **************************************************************************************************
504 30 : SUBROUTINE setup_transport_geometry(cell, transport_ndim, nperiodic, periodic_label, projection, &
505 : measure, measure_unit)
506 : TYPE(cell_type), POINTER :: cell
507 : INTEGER, INTENT(OUT) :: transport_ndim, nperiodic
508 : CHARACTER(LEN=*), INTENT(OUT) :: periodic_label
509 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: projection
510 : REAL(KIND=dp), INTENT(OUT) :: measure
511 : CHARACTER(LEN=*), INTENT(OUT) :: measure_unit
512 :
513 : INTEGER :: idir, jdir, kdir
514 : REAL(KIND=dp) :: detg, invg11, invg12, invg22, norm2
515 : REAL(KIND=dp), DIMENSION(3) :: v1, v2
516 : REAL(KIND=dp), DIMENSION(3, 2) :: basis2
517 :
518 0 : CPASSERT(ASSOCIATED(cell))
519 :
520 120 : nperiodic = COUNT(cell%perd(1:3) == 1)
521 30 : CALL transport_periodicity_label(cell%perd, periodic_label)
522 30 : projection = 0.0_dp
523 :
524 42 : SELECT CASE (nperiodic)
525 : CASE (0)
526 12 : transport_ndim = 3
527 12 : measure = cell%deth
528 48 : DO idir = 1, 3
529 48 : projection(idir, idir) = 1.0_dp
530 : END DO
531 12 : measure_unit = "Angstrom^3"
532 : CASE (1)
533 4 : transport_ndim = 1
534 4 : kdir = 0
535 16 : DO idir = 1, 3
536 16 : IF (cell%perd(idir) == 1) kdir = idir
537 : END DO
538 16 : v1(:) = cell%hmat(:, kdir)
539 16 : norm2 = DOT_PRODUCT(v1, v1)
540 4 : IF (norm2 <= 0.0_dp) CPABORT("KUBO_TRANSPORT periodic cell vector has zero length.")
541 4 : measure = SQRT(norm2)
542 16 : DO jdir = 1, 3
543 52 : DO idir = 1, 3
544 48 : projection(idir, jdir) = v1(idir)*v1(jdir)/norm2
545 : END DO
546 : END DO
547 4 : measure_unit = "Angstrom"
548 : CASE (2)
549 6 : transport_ndim = 2
550 6 : kdir = 0
551 24 : DO idir = 1, 3
552 24 : IF (cell%perd(idir) == 1) THEN
553 12 : kdir = kdir + 1
554 48 : basis2(:, kdir) = cell%hmat(:, idir)
555 : END IF
556 : END DO
557 24 : v1(:) = basis2(:, 1)
558 24 : v2(:) = basis2(:, 2)
559 24 : norm2 = DOT_PRODUCT(v1, v1)
560 24 : invg22 = DOT_PRODUCT(v2, v2)
561 24 : invg12 = DOT_PRODUCT(v1, v2)
562 6 : detg = norm2*invg22 - invg12*invg12
563 6 : IF (detg <= 0.0_dp) CPABORT("KUBO_TRANSPORT periodic cell vectors are linearly dependent.")
564 6 : measure = SQRT(detg)
565 6 : invg11 = invg22/detg
566 6 : invg22 = norm2/detg
567 6 : invg12 = -invg12/detg
568 24 : DO jdir = 1, 3
569 78 : DO idir = 1, 3
570 : projection(idir, jdir) = basis2(idir, 1)*invg11*basis2(jdir, 1) + &
571 : basis2(idir, 1)*invg12*basis2(jdir, 2) + &
572 : basis2(idir, 2)*invg12*basis2(jdir, 1) + &
573 72 : basis2(idir, 2)*invg22*basis2(jdir, 2)
574 : END DO
575 : END DO
576 6 : measure_unit = "Angstrom^2"
577 : CASE DEFAULT
578 8 : transport_ndim = 3
579 8 : measure = cell%deth
580 32 : DO idir = 1, 3
581 32 : projection(idir, idir) = 1.0_dp
582 : END DO
583 38 : measure_unit = "Angstrom^3"
584 : END SELECT
585 :
586 30 : IF (measure <= 0.0_dp) CPABORT("KUBO_TRANSPORT transport normalization measure is not positive.")
587 :
588 30 : END SUBROUTINE setup_transport_geometry
589 :
590 : ! **************************************************************************************************
591 : !> \brief Human-readable periodicity label.
592 : !> \param perd ...
593 : !> \param label ...
594 : ! **************************************************************************************************
595 30 : SUBROUTINE transport_periodicity_label(perd, label)
596 : INTEGER, DIMENSION(3), INTENT(IN) :: perd
597 : CHARACTER(LEN=*), INTENT(OUT) :: label
598 :
599 30 : label = ""
600 30 : IF (perd(1) == 1) label = TRIM(label)//"X"
601 30 : IF (perd(2) == 1) label = TRIM(label)//"Y"
602 30 : IF (perd(3) == 1) label = TRIM(label)//"Z"
603 30 : IF (LEN_TRIM(label) == 0) label = "NONE"
604 :
605 30 : END SUBROUTINE transport_periodicity_label
606 :
607 : ! **************************************************************************************************
608 : !> \brief Output unit labels for a transport dimensionality.
609 : !> \param transport_ndim ...
610 : !> \param density_unit ...
611 : !> \param sigma_iso_unit ...
612 : !> \param sigma_tensor_unit ...
613 : ! **************************************************************************************************
614 30 : SUBROUTINE transport_unit_labels(transport_ndim, density_unit, sigma_iso_unit, sigma_tensor_unit)
615 : INTEGER, INTENT(IN) :: transport_ndim
616 : CHARACTER(LEN=*), INTENT(OUT) :: density_unit, sigma_iso_unit, &
617 : sigma_tensor_unit
618 :
619 34 : SELECT CASE (transport_ndim)
620 : CASE (1)
621 4 : density_unit = "m^-1"
622 4 : sigma_iso_unit = "S*m"
623 4 : sigma_tensor_unit = "S*m"
624 : CASE (2)
625 6 : density_unit = "m^-2"
626 6 : sigma_iso_unit = "S"
627 6 : sigma_tensor_unit = "S"
628 : CASE DEFAULT
629 20 : density_unit = "m^-3"
630 20 : sigma_iso_unit = "S/cm"
631 20 : sigma_tensor_unit = "S/m"
632 : END SELECT
633 :
634 30 : END SUBROUTINE transport_unit_labels
635 :
636 : ! **************************************************************************************************
637 : !> \brief Output conversion factors from internal dimensional transport units.
638 : !> \param transport_ndim ...
639 : !> \param iso_factor ...
640 : !> \param tensor_factor ...
641 : ! **************************************************************************************************
642 30 : SUBROUTINE transport_output_factors(transport_ndim, iso_factor, tensor_factor)
643 : INTEGER, INTENT(IN) :: transport_ndim
644 : REAL(KIND=dp), INTENT(OUT) :: iso_factor, tensor_factor
645 :
646 34 : SELECT CASE (transport_ndim)
647 : CASE (1)
648 4 : iso_factor = 1.0E-10_dp
649 4 : tensor_factor = 1.0E-10_dp
650 : CASE (2)
651 6 : iso_factor = 1.0_dp
652 6 : tensor_factor = 1.0_dp
653 : CASE DEFAULT
654 20 : iso_factor = 1.0E8_dp
655 30 : tensor_factor = 1.0E10_dp
656 : END SELECT
657 :
658 30 : END SUBROUTINE transport_output_factors
659 :
660 : ! **************************************************************************************************
661 : !> \brief Prepare eigenvalues and the transformed [X,H] matrices for one spin channel.
662 : !> \param ks_matrix ...
663 : !> \param s_dense ...
664 : !> \param blacs_env ...
665 : !> \param para_env ...
666 : !> \param particle_set ...
667 : !> \param cell ...
668 : !> \param projection ...
669 : !> \param ao_atom ...
670 : !> \param mo_set ...
671 : !> \param allow_mo_reuse ...
672 : !> \param eig ...
673 : !> \param dHH ...
674 : !> \param reused_mos ...
675 : ! **************************************************************************************************
676 10 : SUBROUTINE setup_spin_transport(ks_matrix, s_dense, blacs_env, para_env, particle_set, cell, &
677 10 : projection, ao_atom, mo_set, allow_mo_reuse, eig, dHH, reused_mos)
678 : TYPE(dbcsr_type), INTENT(IN) :: ks_matrix
679 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: s_dense
680 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
681 : TYPE(mp_para_env_type), POINTER :: para_env
682 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
683 : TYPE(cell_type), POINTER :: cell
684 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: projection
685 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_atom
686 : TYPE(mo_set_type), INTENT(IN) :: mo_set
687 : LOGICAL, INTENT(IN) :: allow_mo_reuse
688 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eig
689 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: dHH
690 : LOGICAL, INTENT(OUT) :: reused_mos
691 :
692 : INTEGER :: idir, mo_nao, nao, nmo
693 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_dense, gmat, h_dense, horth, op, &
694 10 : op_orth, uvec, work
695 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
696 : TYPE(cp_fm_type), POINTER :: mo_coeff
697 :
698 10 : NULLIFY (mo_coeff, mo_eigenvalues)
699 10 : nao = SIZE(s_dense, 1)
700 10 : reused_mos = .FALSE.
701 :
702 10 : CALL dbcsr_to_dense(ks_matrix, blacs_env, para_env, h_dense)
703 60 : ALLOCATE (op(nao, nao), work(nao, nao))
704 :
705 10 : IF (allow_mo_reuse) THEN
706 : CALL get_mo_set(mo_set=mo_set, nao=mo_nao, nmo=nmo, eigenvalues=mo_eigenvalues, &
707 2 : mo_coeff=mo_coeff)
708 : reused_mos = mo_nao == nao .AND. nmo >= nao .AND. ASSOCIATED(mo_eigenvalues) .AND. &
709 2 : ASSOCIATED(mo_coeff)
710 2 : IF (reused_mos) reused_mos = SIZE(mo_eigenvalues) >= nao
711 : END IF
712 :
713 10 : IF (reused_mos) THEN
714 6 : ALLOCATE (coeff_dense(nao, nao))
715 22 : eig(:) = mo_eigenvalues(1:nao)
716 2 : CALL cp_fm_get_submatrix(mo_coeff, coeff_dense, n_rows=nao, n_cols=nao)
717 8 : DO idir = 1, 3
718 6 : CALL build_commutator_kernel(h_dense, particle_set, cell, projection, ao_atom, idir, op)
719 8 : CALL transform_to_eigenbasis(op, coeff_dense, dHH(:, :, idir), work)
720 : END DO
721 2 : DEALLOCATE (coeff_dense)
722 : ELSE
723 72 : ALLOCATE (gmat(nao, nao), horth(nao, nao), op_orth(nao, nao), uvec(nao, nao))
724 8 : CALL inverse_sqrt_symmetric(s_dense, gmat)
725 :
726 16888 : work(:, :) = MATMUL(h_dense, gmat)
727 8888 : horth(:, :) = MATMUL(gmat, work)
728 8 : CALL symmetrize(horth)
729 :
730 888 : uvec(:, :) = horth
731 8 : CALL diagonalize_symmetric(uvec, eig)
732 :
733 32 : DO idir = 1, 3
734 24 : CALL build_commutator_kernel(h_dense, particle_set, cell, projection, ao_atom, idir, op)
735 26664 : work(:, :) = MATMUL(op, gmat)
736 26664 : op_orth(:, :) = MATMUL(gmat, work)
737 32 : CALL transform_to_eigenbasis(op_orth, uvec, dHH(:, :, idir), work)
738 : END DO
739 8 : DEALLOCATE (gmat, horth, op_orth, uvec)
740 : END IF
741 :
742 10 : DEALLOCATE (h_dense, op, work)
743 :
744 20 : END SUBROUTINE setup_spin_transport
745 :
746 : ! **************************************************************************************************
747 : !> \brief Symmetric inverse square root via dense diagonalization.
748 : !> \param matrix ...
749 : !> \param inverse_sqrt ...
750 : ! **************************************************************************************************
751 8 : SUBROUTINE inverse_sqrt_symmetric(matrix, inverse_sqrt)
752 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
753 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: inverse_sqrt
754 :
755 : INTEGER :: i, info, lwork, n
756 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig, work
757 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: scaled, vectors
758 :
759 8 : n = SIZE(matrix, 1)
760 8 : lwork = MAX(1, 8*n)
761 80 : ALLOCATE (eig(n), scaled(n, n), vectors(n, n), work(lwork))
762 :
763 888 : vectors(:, :) = matrix
764 8 : CALL dsyev("V", "L", n, vectors, n, eig, work, lwork, info)
765 8 : IF (info /= 0) CPABORT("Diagonalization of the overlap matrix failed.")
766 88 : IF (MINVAL(eig) <= 0.0_dp) CPABORT("Overlap matrix is not positive definite.")
767 :
768 888 : scaled(:, :) = vectors
769 88 : DO i = 1, n
770 888 : scaled(:, i) = scaled(:, i)/SQRT(eig(i))
771 : END DO
772 8 : CALL dgemm("N", "T", n, n, n, 1.0_dp, scaled, n, vectors, n, 0.0_dp, inverse_sqrt, n)
773 :
774 8 : DEALLOCATE (eig, scaled, vectors, work)
775 :
776 8 : END SUBROUTINE inverse_sqrt_symmetric
777 :
778 : ! **************************************************************************************************
779 : !> \brief Dense symmetric diagonalization.
780 : !> \param matrix ...
781 : !> \param eig ...
782 : ! **************************************************************************************************
783 8 : SUBROUTINE diagonalize_symmetric(matrix, eig)
784 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix
785 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eig
786 :
787 : INTEGER :: info, lwork, n
788 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
789 :
790 8 : n = SIZE(matrix, 1)
791 8 : lwork = MAX(1, 8*n)
792 24 : ALLOCATE (work(lwork))
793 8 : CALL dsyev("V", "L", n, matrix, n, eig, work, lwork, info)
794 8 : IF (info /= 0) CPABORT("Diagonalization of the orthogonalized KS matrix failed.")
795 8 : DEALLOCATE (work)
796 :
797 8 : END SUBROUTINE diagonalize_symmetric
798 :
799 : ! **************************************************************************************************
800 : !> \brief Enforce exact symmetry after dense products.
801 : !> \param matrix ...
802 : ! **************************************************************************************************
803 8 : SUBROUTINE symmetrize(matrix)
804 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix
805 :
806 : INTEGER :: i, j, n
807 : REAL(KIND=dp) :: value
808 :
809 8 : n = SIZE(matrix, 1)
810 88 : DO j = 1, n
811 448 : DO i = j + 1, n
812 360 : value = 0.5_dp*(matrix(i, j) + matrix(j, i))
813 360 : matrix(i, j) = value
814 440 : matrix(j, i) = value
815 : END DO
816 : END DO
817 :
818 8 : END SUBROUTINE symmetrize
819 :
820 : ! **************************************************************************************************
821 : !> \brief Build element-wise kernel (P*(r_col-r_row))_idir * matrix(row,col).
822 : !> \param matrix ...
823 : !> \param particle_set ...
824 : !> \param cell ...
825 : !> \param projection ...
826 : !> \param ao_atom ...
827 : !> \param idir ...
828 : !> \param op ...
829 : ! **************************************************************************************************
830 30 : SUBROUTINE build_commutator_kernel(matrix, particle_set, cell, projection, ao_atom, idir, op)
831 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
832 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
833 : TYPE(cell_type), POINTER :: cell
834 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: projection
835 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_atom
836 : INTEGER, INTENT(IN) :: idir
837 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: op
838 :
839 : INTEGER :: i, j, n
840 : REAL(KIND=dp), DIMENSION(3) :: dr, projected_dr
841 :
842 30 : n = SIZE(matrix, 1)
843 330 : DO j = 1, n
844 3330 : DO i = 1, n
845 12000 : dr(:) = pbc(particle_set(ao_atom(j))%r(:) - particle_set(ao_atom(i))%r(:), cell)
846 39000 : projected_dr(:) = MATMUL(projection, dr)
847 3300 : op(i, j) = projected_dr(idir)*matrix(i, j)
848 : END DO
849 : END DO
850 :
851 30 : END SUBROUTINE build_commutator_kernel
852 :
853 : ! **************************************************************************************************
854 : !> \brief Transform an AO matrix to the eigenvector basis.
855 : !> \param op ...
856 : !> \param uvec ...
857 : !> \param transformed ...
858 : !> \param work ...
859 : ! **************************************************************************************************
860 30 : SUBROUTINE transform_to_eigenbasis(op, uvec, transformed, work)
861 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: op, uvec
862 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: transformed
863 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: work
864 :
865 : INTEGER :: n
866 :
867 30 : n = SIZE(op, 1)
868 30 : CALL dgemm("N", "N", n, n, n, 1.0_dp, op, n, uvec, n, 0.0_dp, work, n)
869 30 : CALL dgemm("T", "N", n, n, n, 1.0_dp, uvec, n, work, n, 0.0_dp, transformed, n)
870 :
871 30 : END SUBROUTINE transform_to_eigenbasis
872 :
873 : ! **************************************************************************************************
874 : !> \brief Fermi occupation with overflow protection.
875 : !> \param eig ...
876 : !> \param mu ...
877 : !> \param kT ...
878 : !> \param maxocc ...
879 : !> \return ...
880 : ! **************************************************************************************************
881 6144 : PURE FUNCTION fermi_occ(eig, mu, kT, maxocc) RESULT(occ)
882 : REAL(KIND=dp), INTENT(IN) :: eig, mu, kT, maxocc
883 : REAL(KIND=dp) :: occ
884 :
885 : REAL(KIND=dp) :: x
886 :
887 6144 : x = (eig - mu)/kT
888 6144 : IF (x > 80.0_dp) THEN
889 : occ = 0.0_dp
890 1824 : ELSE IF (x < -80.0_dp) THEN
891 1036 : occ = maxocc
892 : ELSE
893 788 : occ = maxocc/(1.0_dp + EXP(x))
894 : END IF
895 :
896 6144 : END FUNCTION fermi_occ
897 :
898 : ! **************************************************************************************************
899 : !> \brief Fermi vacancy below the neutral point.
900 : !> \param eig ...
901 : !> \param mu ...
902 : !> \param kT ...
903 : !> \param maxocc ...
904 : !> \return ...
905 : ! **************************************************************************************************
906 912 : PURE FUNCTION fermi_hole(eig, mu, kT, maxocc) RESULT(hole)
907 : REAL(KIND=dp), INTENT(IN) :: eig, mu, kT, maxocc
908 : REAL(KIND=dp) :: hole
909 :
910 912 : hole = maxocc - fermi_occ(eig, mu, kT, maxocc)
911 :
912 912 : END FUNCTION fermi_hole
913 :
914 : ! **************************************************************************************************
915 : !> \brief Total electron count at a chemical potential.
916 : !> \param eig ...
917 : !> \param maxocc ...
918 : !> \param mu ...
919 : !> \param kT ...
920 : !> \return ...
921 : ! **************************************************************************************************
922 0 : PURE FUNCTION electron_count(eig, maxocc, mu, kT) RESULT(count)
923 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
924 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
925 : REAL(KIND=dp), INTENT(IN) :: mu, kT
926 : REAL(KIND=dp) :: count
927 :
928 : INTEGER :: ispin, n
929 :
930 0 : count = 0.0_dp
931 0 : DO ispin = 1, SIZE(eig, 2)
932 0 : DO n = 1, SIZE(eig, 1)
933 0 : count = count + fermi_occ(eig(n, ispin), mu, kT, maxocc(ispin))
934 : END DO
935 : END DO
936 :
937 0 : END FUNCTION electron_count
938 :
939 : ! **************************************************************************************************
940 : !> \brief Determine the neutral chemical potential using the legacy two-step definition.
941 : !> \param eig ...
942 : !> \param maxocc ...
943 : !> \param kT ...
944 : !> \param neutral_electrons ...
945 : !> \param neutral_grid ...
946 : !> \param emin ...
947 : !> \param emax ...
948 : !> \param mu0 ...
949 : ! **************************************************************************************************
950 0 : SUBROUTINE find_neutral_mu(eig, maxocc, kT, neutral_electrons, neutral_grid, emin, emax, mu0)
951 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
952 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
953 : REAL(KIND=dp), INTENT(IN) :: kT, neutral_electrons
954 : INTEGER, INTENT(IN) :: neutral_grid
955 : REAL(KIND=dp), INTENT(IN) :: emin, emax
956 : REAL(KIND=dp), INTENT(OUT) :: mu0
957 :
958 : INTEGER :: i
959 : REAL(KIND=dp) :: best, count, mu, mup, ne, nh, &
960 : previous_count
961 :
962 0 : mup = 0.5_dp*(emin + emax)
963 0 : previous_count = electron_count(eig, maxocc, emin, kT)
964 0 : DO i = 1, neutral_grid
965 0 : mu = emin + (emax - emin)*REAL(i, KIND=dp)/REAL(neutral_grid, KIND=dp)
966 0 : count = electron_count(eig, maxocc, mu, kT)
967 0 : IF ((count - neutral_electrons)*(previous_count - neutral_electrons) <= 0.0_dp) mup = mu
968 0 : previous_count = count
969 : END DO
970 :
971 0 : mu0 = mup
972 0 : best = HUGE(1.0_dp)
973 0 : DO i = 1, neutral_grid
974 0 : mu = emin + (emax - emin)*REAL(i, KIND=dp)/REAL(neutral_grid, KIND=dp)
975 0 : CALL electron_hole_density(eig, maxocc, mu, mup, kT, ne, nh)
976 0 : IF (ABS(nh - ne) <= best) THEN
977 0 : mu0 = mu
978 0 : best = ABS(nh - ne)
979 : END IF
980 : END DO
981 :
982 0 : END SUBROUTINE find_neutral_mu
983 :
984 : ! **************************************************************************************************
985 : !> \brief Electron and hole counts relative to a neutral point.
986 : !> \param eig ...
987 : !> \param maxocc ...
988 : !> \param mu ...
989 : !> \param mu0 ...
990 : !> \param kT ...
991 : !> \param ne ...
992 : !> \param nh ...
993 : ! **************************************************************************************************
994 30 : SUBROUTINE electron_hole_density(eig, maxocc, mu, mu0, kT, ne, nh)
995 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
996 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
997 : REAL(KIND=dp), INTENT(IN) :: mu, mu0, kT
998 : REAL(KIND=dp), INTENT(OUT) :: ne, nh
999 :
1000 : INTEGER :: ispin, n
1001 :
1002 30 : ne = 0.0_dp
1003 30 : nh = 0.0_dp
1004 142 : DO ispin = 1, SIZE(eig, 2)
1005 3214 : DO n = 1, SIZE(eig, 1)
1006 3072 : IF (eig(n, ispin) <= mu0) nh = nh + fermi_hole(eig(n, ispin), mu, kT, maxocc(ispin))
1007 3184 : IF (eig(n, ispin) >= mu0) ne = ne + fermi_occ(eig(n, ispin), mu, kT, maxocc(ispin))
1008 : END DO
1009 : END DO
1010 :
1011 30 : END SUBROUTINE electron_hole_density
1012 :
1013 : ! **************************************************************************************************
1014 : !> \brief Conductivity tensor at one chemical potential.
1015 : !> \param eig ...
1016 : !> \param dHH ...
1017 : !> \param maxocc ...
1018 : !> \param mu ...
1019 : !> \param kT ...
1020 : !> \param dissipation ...
1021 : !> \param measure ...
1022 : !> \param transport_ndim ...
1023 : !> \param sigma ...
1024 : !> \param bloch Use the stable Fermi limit for Bloch currents or the complete DC tensor
1025 : ! **************************************************************************************************
1026 30 : SUBROUTINE conductivity_at_mu(eig, dHH, maxocc, mu, kT, dissipation, measure, transport_ndim, sigma, bloch)
1027 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
1028 : COMPLEX(KIND=dp), DIMENSION(:, :, :, :), &
1029 : INTENT(IN) :: dHH
1030 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
1031 : REAL(KIND=dp), INTENT(IN) :: mu, kT, dissipation, measure
1032 : INTEGER, INTENT(IN) :: transport_ndim
1033 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: sigma
1034 : LOGICAL, INTENT(IN) :: bloch
1035 :
1036 : INTEGER :: idir, ispin, jdir, m, n, nao
1037 : REAL(KIND=dp) :: delta, eps, factor, g0
1038 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occ
1039 :
1040 30 : nao = SIZE(eig, 1)
1041 30 : eps = 1.0E-12_dp
1042 30 : g0 = e_charge**2/(pi*h_bar)
1043 30 : sigma = 0.0_dp
1044 :
1045 90 : ALLOCATE (occ(nao))
1046 :
1047 142 : DO ispin = 1, SIZE(eig, 2)
1048 3184 : DO n = 1, nao
1049 3184 : occ(n) = fermi_occ(eig(n, ispin), mu, kT, maxocc(ispin))
1050 : END DO
1051 :
1052 478 : DO idir = 1, 3
1053 1456 : DO jdir = 1, 3
1054 28992 : DO n = 1, nao
1055 1065168 : DO m = 1, nao
1056 1036512 : delta = eig(n, ispin) - eig(m, ispin)
1057 1036512 : IF (bloch) THEN
1058 : factor = maxocc(ispin)*fermi_slope(eig(n, ispin), eig(m, ispin), mu, kT)* &
1059 1022112 : dissipation/(dissipation**2 + delta**2)
1060 : ELSE
1061 14400 : factor = (occ(n) - occ(m))/(delta + eps)*dissipation/(dissipation**2 + delta**2)
1062 : END IF
1063 : sigma(jdir, idir) = sigma(jdir, idir) + &
1064 1064160 : REAL(dHH(m, n, jdir, ispin)*dHH(n, m, idir, ispin), dp)*factor
1065 : END DO
1066 : END DO
1067 : END DO
1068 : END DO
1069 : END DO
1070 :
1071 390 : sigma = pi*g0*sigma*angstrom**(2 - transport_ndim)/measure
1072 :
1073 30 : DEALLOCATE (occ)
1074 :
1075 30 : END SUBROUTINE conductivity_at_mu
1076 :
1077 : END MODULE qs_kubo_transport
|
