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: KuhneHeskeProdan2020,&
13 : cite_reference
14 : USE cell_types, ONLY: cell_type,&
15 : pbc
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
18 : dbcsr_p_type,&
19 : dbcsr_type
20 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
21 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
22 : cp_fm_struct_release,&
23 : cp_fm_struct_type
24 : USE cp_fm_types, ONLY: cp_fm_create,&
25 : cp_fm_get_submatrix,&
26 : cp_fm_release,&
27 : cp_fm_type
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_get_default_io_unit,&
30 : cp_logger_type
31 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
32 : section_vals_type,&
33 : section_vals_val_get
34 : USE kinds, ONLY: dp
35 : USE mathconstants, ONLY: pi
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE particle_types, ONLY: particle_type
38 : USE physcon, ONLY: a_bohr,&
39 : angstrom,&
40 : e_charge,&
41 : evolt,&
42 : h_bar,&
43 : kelvin
44 : USE qs_environment_types, ONLY: get_qs_env,&
45 : qs_environment_type
46 : USE qs_mo_types, ONLY: get_mo_set,&
47 : mo_set_type
48 : USE string_utilities, ONLY: uppercase
49 : #include "./base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 :
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kubo_transport'
56 :
57 : PUBLIC :: qs_scf_post_kubo_transport
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief Compute and print the finite-volume Kubo-Greenwood conductivity tensor.
63 : !> \param qs_env Quickstep environment
64 : ! **************************************************************************************************
65 12667 : SUBROUTINE qs_scf_post_kubo_transport(qs_env)
66 : TYPE(qs_environment_type), POINTER :: qs_env
67 :
68 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_kubo_transport'
69 :
70 : CHARACTER(LEN=32) :: density_unit, measure_unit, method, &
71 : periodic_label, sigma_iso_unit, &
72 : sigma_tensor_unit
73 : CHARACTER(LEN=512) :: header
74 : INTEGER :: handle, imu, ispin, nao, natom, &
75 : nelectron_total, neutral_grid, nmu, &
76 : nperiodic, nspins, output_unit, &
77 : transport_ndim
78 12667 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_atom
79 12667 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
80 : LOGICAL :: allow_mo_reuse, do_kpoints, kubo_active, &
81 : neutral_mu_explicit, ot_active
82 12667 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: reused_mos
83 : REAL(KIND=dp) :: density_factor, dissipation, emax, emin, iso_factor, measure, measure_m, &
84 : mu, mu0, mu_step, ne, nh, sigma_iso, temperature, tensor_factor
85 12667 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: maxocc
86 12667 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eig, s_dense
87 12667 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dHH
88 : REAL(KIND=dp), DIMENSION(3, 3) :: projection, sigma, sigma_out
89 12667 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy_range
90 : TYPE(cell_type), POINTER :: cell
91 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
92 : TYPE(cp_logger_type), POINTER :: logger
93 12667 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
94 12667 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
95 : TYPE(mp_para_env_type), POINTER :: para_env
96 12667 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
97 : TYPE(section_vals_type), POINTER :: input_section, kubo_section, ot_section
98 :
99 12667 : CALL timeset(routineN, handle)
100 :
101 12667 : NULLIFY (blacs_env, cell, energy_range, input_section, kubo_section, logger, matrix_ks, matrix_s, &
102 12667 : mos, ot_section, para_env, particle_set, row_blk_sizes)
103 :
104 12667 : CALL get_qs_env(qs_env, input=input_section)
105 12667 : kubo_section => section_vals_get_subs_vals(input_section, "PROPERTIES%KUBO_TRANSPORT")
106 12667 : CALL section_vals_val_get(kubo_section, "_SECTION_PARAMETERS_", l_val=kubo_active)
107 12667 : IF (.NOT. kubo_active) THEN
108 12657 : CALL timestop(handle)
109 12657 : RETURN
110 : END IF
111 10 : CALL cite_reference(KuhneHeskeProdan2020)
112 :
113 10 : CALL section_vals_val_get(kubo_section, "METHOD", c_val=method)
114 10 : CALL uppercase(method)
115 10 : SELECT CASE (TRIM(method))
116 : CASE ("DIAGONALIZATION")
117 : CASE ("TD", "CHEBYSHEV")
118 0 : CPABORT("KUBO_TRANSPORT METHOD "//TRIM(method)//" is reserved but not implemented yet.")
119 : CASE DEFAULT
120 10 : CPABORT("Unknown KUBO_TRANSPORT METHOD: "//TRIM(method))
121 : END SELECT
122 10 : ot_section => section_vals_get_subs_vals(input_section, "DFT%SCF%OT")
123 10 : CALL section_vals_val_get(ot_section, "_SECTION_PARAMETERS_", l_val=ot_active)
124 10 : allow_mo_reuse = .NOT. ot_active
125 :
126 : CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, do_kpoints=do_kpoints, &
127 : matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, natom=natom, &
128 10 : nelectron_total=nelectron_total, para_env=para_env, particle_set=particle_set)
129 :
130 10 : IF (do_kpoints) CPABORT("KUBO_TRANSPORT currently expects a finite Gamma-point supercell.")
131 10 : CPASSERT(ASSOCIATED(matrix_s))
132 10 : CPASSERT(ASSOCIATED(matrix_ks))
133 10 : CPASSERT(ASSOCIATED(cell))
134 10 : CPASSERT(ASSOCIATED(particle_set))
135 :
136 10 : CALL section_vals_val_get(kubo_section, "TEMPERATURE", r_val=temperature)
137 10 : CALL section_vals_val_get(kubo_section, "DISSIPATION", r_val=dissipation)
138 10 : CALL section_vals_val_get(kubo_section, "ENERGY_RANGE", r_vals=energy_range)
139 10 : CALL section_vals_val_get(kubo_section, "NEUTRAL_MU", explicit=neutral_mu_explicit)
140 10 : IF (neutral_mu_explicit) CALL section_vals_val_get(kubo_section, "NEUTRAL_MU", r_val=mu0)
141 10 : CALL section_vals_val_get(kubo_section, "N_MU", i_val=nmu)
142 10 : CALL section_vals_val_get(kubo_section, "NEUTRAL_GRID", i_val=neutral_grid)
143 :
144 10 : IF (temperature <= 0.0_dp) CPABORT("KUBO_TRANSPORT TEMPERATURE must be positive.")
145 10 : IF (dissipation <= 0.0_dp) CPABORT("KUBO_TRANSPORT DISSIPATION must be positive.")
146 10 : IF (nmu < 1) CPABORT("KUBO_TRANSPORT N_MU must be at least one.")
147 10 : IF (neutral_grid < 10) CPABORT("KUBO_TRANSPORT NEUTRAL_GRID must be at least ten.")
148 :
149 10 : nspins = SIZE(matrix_ks)
150 10 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao, row_blk_size=row_blk_sizes)
151 10 : CALL build_ao_atom_map(row_blk_sizes, natom, ao_atom)
152 10 : CALL dbcsr_to_dense(matrix_s(1)%matrix, blacs_env, para_env, s_dense)
153 : CALL setup_transport_geometry(cell, transport_ndim, nperiodic, periodic_label, projection, &
154 10 : measure, measure_unit)
155 10 : CALL transport_unit_labels(transport_ndim, density_unit, sigma_iso_unit, sigma_tensor_unit)
156 :
157 130 : ALLOCATE (eig(nao, nspins), dHH(nao, nao, 3, nspins), maxocc(nspins), reused_mos(nspins))
158 :
159 20 : DO ispin = 1, nspins
160 10 : maxocc(ispin) = mos(ispin)%maxocc
161 : CALL setup_spin_transport(matrix_ks(ispin)%matrix, s_dense, blacs_env, para_env, &
162 : particle_set, cell, projection, ao_atom, mos(ispin), allow_mo_reuse, &
163 20 : eig(:, ispin), dHH(:, :, :, ispin), reused_mos(ispin))
164 : END DO
165 :
166 120 : emin = MINVAL(eig)
167 120 : emax = MAXVAL(eig)
168 10 : IF (.NOT. neutral_mu_explicit) THEN
169 : CALL find_neutral_mu(eig, maxocc, temperature, REAL(nelectron_total, KIND=dp), neutral_grid, &
170 0 : emin, emax, mu0)
171 : END IF
172 :
173 10 : IF (energy_range(1) < energy_range(2)) THEN
174 10 : emin = energy_range(1)
175 10 : emax = energy_range(2)
176 : END IF
177 :
178 10 : measure_m = measure*a_bohr**transport_ndim
179 10 : density_factor = 1.0_dp/measure_m
180 10 : CALL transport_output_factors(transport_ndim, iso_factor, tensor_factor)
181 :
182 10 : logger => cp_get_default_logger()
183 10 : output_unit = cp_logger_get_default_io_unit(logger)
184 :
185 10 : IF (output_unit > 0) THEN
186 5 : WRITE (output_unit, '(/,T2,A)') "Kubo-Greenwood finite-volume transport"
187 5 : WRITE (output_unit, '(T3,A,T38,A)') "Method:", TRIM(method)
188 5 : WRITE (output_unit, '(T3,A,T38,A)') "Cell periodicity:", TRIM(periodic_label)
189 5 : IF (nperiodic == 0) THEN
190 3 : WRITE (output_unit, '(T3,A,T38,A)') "Transport normalization:", "3D finite box"
191 : ELSE
192 2 : WRITE (output_unit, '(T3,A,T38,I12)') "Transport dimensionality:", transport_ndim
193 : END IF
194 6 : IF (ALL(reused_mos)) THEN
195 1 : WRITE (output_unit, '(T3,A,T38,A)') "Eigensystem source:", "CP2K canonical MOs"
196 8 : ELSEIF (ANY(reused_mos)) THEN
197 0 : WRITE (output_unit, '(T3,A,T38,A)') "Eigensystem source:", "mixed canonical/internal"
198 : ELSE
199 4 : WRITE (output_unit, '(T3,A,T38,A)') "Eigensystem source:", "internal diagonalization"
200 : END IF
201 5 : WRITE (output_unit, '(T3,A,T38,I12)') "Number of atoms:", natom
202 5 : WRITE (output_unit, '(T3,A,T38,I12)') "Number of AO functions:", nao
203 5 : WRITE (output_unit, '(T3,A,T38,F12.4)') "Temperature [K]:", temperature*kelvin
204 5 : WRITE (output_unit, '(T3,A,T38,F12.4)') "Dissipation [K]:", dissipation*kelvin
205 5 : WRITE (output_unit, '(T3,A,T38,F12.6)') "Dissipation [eV]:", dissipation*evolt
206 5 : WRITE (output_unit, '(T3,A,T38,F12.6)') "Neutral chemical potential [eV]:", mu0*evolt
207 5 : WRITE (output_unit, '(T3,A,T38,ES12.4,1X,A)') "Transport measure:", &
208 10 : measure*angstrom**transport_ndim, TRIM(measure_unit)
209 : header = "mu[Ha] mu-mu0[eV] Ne["//TRIM(density_unit)//"] "// &
210 : "Nh["//TRIM(density_unit)//"] sigma_iso["//TRIM(sigma_iso_unit)//"] "// &
211 : "sigma_xx["//TRIM(sigma_tensor_unit)//"] sigma_xy["// &
212 : TRIM(sigma_tensor_unit)//"] sigma_xz["//TRIM(sigma_tensor_unit)// &
213 : "] sigma_yx["//TRIM(sigma_tensor_unit)//"] sigma_yy["// &
214 : TRIM(sigma_tensor_unit)//"] sigma_yz["//TRIM(sigma_tensor_unit)// &
215 : "] sigma_zx["//TRIM(sigma_tensor_unit)//"] sigma_zy["// &
216 5 : TRIM(sigma_tensor_unit)//"] sigma_zz["//TRIM(sigma_tensor_unit)//"]"
217 5 : WRITE (output_unit, '(/,T3,A)') TRIM(header)
218 : END IF
219 :
220 10 : IF (nmu == 1) THEN
221 : mu_step = 0.0_dp
222 : ELSE
223 0 : mu_step = (emax - emin)/REAL(nmu - 1, KIND=dp)
224 : END IF
225 :
226 20 : DO imu = 1, nmu
227 10 : mu = emin + REAL(imu - 1, KIND=dp)*mu_step
228 10 : CALL electron_hole_density(eig, maxocc, mu, mu0, temperature, ne, nh)
229 : CALL conductivity_at_mu(eig, dHH, maxocc, mu, temperature, dissipation, measure, &
230 10 : transport_ndim, sigma)
231 30 : IF (output_unit > 0) THEN
232 : sigma_iso = (sigma(1, 1) + sigma(2, 2) + sigma(3, 3))/ &
233 5 : REAL(transport_ndim, KIND=dp)*iso_factor
234 65 : sigma_out(:, :) = sigma(:, :)*tensor_factor
235 5 : WRITE (output_unit, '(T3,F12.6,F14.6,2ES16.7,10ES17.8)') mu, (mu - mu0)*evolt, &
236 5 : ne*density_factor, nh*density_factor, &
237 5 : sigma_iso, &
238 5 : sigma_out(1, 1), sigma_out(1, 2), sigma_out(1, 3), &
239 5 : sigma_out(2, 1), sigma_out(2, 2), sigma_out(2, 3), &
240 10 : sigma_out(3, 1), sigma_out(3, 2), sigma_out(3, 3)
241 1 : SELECT CASE (transport_ndim)
242 : CASE (1)
243 1 : WRITE (output_unit, '(T3,A,1X,ES17.8)') "KUBO_TRANSPORT| sigma_iso[S*m]", sigma_iso
244 : CASE (2)
245 1 : WRITE (output_unit, '(T3,A,1X,ES17.8)') "KUBO_TRANSPORT| sigma_iso[S]", sigma_iso
246 : CASE DEFAULT
247 5 : WRITE (output_unit, '(T3,A,1X,ES17.8)') "KUBO_TRANSPORT| sigma_iso[S/cm]", sigma_iso
248 : END SELECT
249 : END IF
250 : END DO
251 :
252 10 : DEALLOCATE (ao_atom, dHH, eig, maxocc, reused_mos, s_dense)
253 :
254 10 : CALL timestop(handle)
255 :
256 38071 : END SUBROUTINE qs_scf_post_kubo_transport
257 :
258 : ! **************************************************************************************************
259 : !> \brief Build a global AO -> atom map from DBCSR atomic block sizes.
260 : !> \param row_blk_sizes ...
261 : !> \param natom ...
262 : !> \param ao_atom ...
263 : ! **************************************************************************************************
264 10 : SUBROUTINE build_ao_atom_map(row_blk_sizes, natom, ao_atom)
265 : INTEGER, DIMENSION(:), INTENT(IN) :: row_blk_sizes
266 : INTEGER, INTENT(IN) :: natom
267 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ao_atom
268 :
269 : INTEGER :: iao, iatom, nao, u
270 :
271 30 : nao = SUM(row_blk_sizes(1:natom))
272 30 : ALLOCATE (ao_atom(nao))
273 10 : iao = 0
274 30 : DO iatom = 1, natom
275 130 : DO u = 1, row_blk_sizes(iatom)
276 100 : iao = iao + 1
277 120 : ao_atom(iao) = iatom
278 : END DO
279 : END DO
280 :
281 10 : END SUBROUTINE build_ao_atom_map
282 :
283 : ! **************************************************************************************************
284 : !> \brief Copy a DBCSR matrix to a replicated dense matrix.
285 : !> \param matrix ...
286 : !> \param blacs_env ...
287 : !> \param para_env ...
288 : !> \param dense ...
289 : ! **************************************************************************************************
290 20 : SUBROUTINE dbcsr_to_dense(matrix, blacs_env, para_env, dense)
291 : TYPE(dbcsr_type), INTENT(IN) :: matrix
292 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
293 : TYPE(mp_para_env_type), POINTER :: para_env
294 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
295 : INTENT(OUT) :: dense
296 :
297 : INTEGER :: ncol, nrow
298 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
299 : TYPE(cp_fm_type) :: fm
300 :
301 20 : NULLIFY (fm_struct)
302 :
303 20 : CALL dbcsr_get_info(matrix, nfullrows_total=nrow, nfullcols_total=ncol)
304 80 : ALLOCATE (dense(nrow, ncol))
305 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
306 20 : nrow_global=nrow, ncol_global=ncol)
307 20 : CALL cp_fm_create(fm, fm_struct)
308 20 : CALL copy_dbcsr_to_fm(matrix, fm)
309 20 : CALL cp_fm_get_submatrix(fm, dense, n_rows=nrow, n_cols=ncol)
310 20 : CALL cp_fm_release(fm)
311 20 : CALL cp_fm_struct_release(fm_struct)
312 :
313 40 : END SUBROUTINE dbcsr_to_dense
314 :
315 : ! **************************************************************************************************
316 : !> \brief Set up the periodic transport subspace and its normalization measure.
317 : !> \param cell ...
318 : !> \param transport_ndim ...
319 : !> \param nperiodic ...
320 : !> \param periodic_label ...
321 : !> \param projection ...
322 : !> \param measure ...
323 : !> \param measure_unit ...
324 : ! **************************************************************************************************
325 10 : SUBROUTINE setup_transport_geometry(cell, transport_ndim, nperiodic, periodic_label, projection, &
326 : measure, measure_unit)
327 : TYPE(cell_type), POINTER :: cell
328 : INTEGER, INTENT(OUT) :: transport_ndim, nperiodic
329 : CHARACTER(LEN=*), INTENT(OUT) :: periodic_label
330 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: projection
331 : REAL(KIND=dp), INTENT(OUT) :: measure
332 : CHARACTER(LEN=*), INTENT(OUT) :: measure_unit
333 :
334 : INTEGER :: idir, jdir, kdir
335 : REAL(KIND=dp) :: detg, invg11, invg12, invg22, norm2
336 : REAL(KIND=dp), DIMENSION(3) :: v1, v2
337 : REAL(KIND=dp), DIMENSION(3, 2) :: basis2
338 :
339 0 : CPASSERT(ASSOCIATED(cell))
340 :
341 40 : nperiodic = COUNT(cell%perd(1:3) == 1)
342 10 : CALL transport_periodicity_label(cell%perd, periodic_label)
343 10 : projection = 0.0_dp
344 :
345 16 : SELECT CASE (nperiodic)
346 : CASE (0)
347 6 : transport_ndim = 3
348 6 : measure = cell%deth
349 24 : DO idir = 1, 3
350 24 : projection(idir, idir) = 1.0_dp
351 : END DO
352 6 : measure_unit = "Angstrom^3"
353 : CASE (1)
354 2 : transport_ndim = 1
355 2 : kdir = 0
356 8 : DO idir = 1, 3
357 8 : IF (cell%perd(idir) == 1) kdir = idir
358 : END DO
359 8 : v1(:) = cell%hmat(:, kdir)
360 8 : norm2 = DOT_PRODUCT(v1, v1)
361 2 : IF (norm2 <= 0.0_dp) CPABORT("KUBO_TRANSPORT periodic cell vector has zero length.")
362 2 : measure = SQRT(norm2)
363 8 : DO jdir = 1, 3
364 26 : DO idir = 1, 3
365 24 : projection(idir, jdir) = v1(idir)*v1(jdir)/norm2
366 : END DO
367 : END DO
368 2 : measure_unit = "Angstrom"
369 : CASE (2)
370 2 : transport_ndim = 2
371 2 : kdir = 0
372 8 : DO idir = 1, 3
373 8 : IF (cell%perd(idir) == 1) THEN
374 4 : kdir = kdir + 1
375 16 : basis2(:, kdir) = cell%hmat(:, idir)
376 : END IF
377 : END DO
378 8 : v1(:) = basis2(:, 1)
379 8 : v2(:) = basis2(:, 2)
380 8 : norm2 = DOT_PRODUCT(v1, v1)
381 8 : invg22 = DOT_PRODUCT(v2, v2)
382 8 : invg12 = DOT_PRODUCT(v1, v2)
383 2 : detg = norm2*invg22 - invg12*invg12
384 2 : IF (detg <= 0.0_dp) CPABORT("KUBO_TRANSPORT periodic cell vectors are linearly dependent.")
385 2 : measure = SQRT(detg)
386 2 : invg11 = invg22/detg
387 2 : invg22 = norm2/detg
388 2 : invg12 = -invg12/detg
389 8 : DO jdir = 1, 3
390 26 : DO idir = 1, 3
391 : projection(idir, jdir) = basis2(idir, 1)*invg11*basis2(jdir, 1) + &
392 : basis2(idir, 1)*invg12*basis2(jdir, 2) + &
393 : basis2(idir, 2)*invg12*basis2(jdir, 1) + &
394 24 : basis2(idir, 2)*invg22*basis2(jdir, 2)
395 : END DO
396 : END DO
397 2 : measure_unit = "Angstrom^2"
398 : CASE DEFAULT
399 0 : transport_ndim = 3
400 0 : measure = cell%deth
401 0 : DO idir = 1, 3
402 0 : projection(idir, idir) = 1.0_dp
403 : END DO
404 10 : measure_unit = "Angstrom^3"
405 : END SELECT
406 :
407 10 : IF (measure <= 0.0_dp) CPABORT("KUBO_TRANSPORT transport normalization measure is not positive.")
408 :
409 10 : END SUBROUTINE setup_transport_geometry
410 :
411 : ! **************************************************************************************************
412 : !> \brief Human-readable periodicity label.
413 : !> \param perd ...
414 : !> \param label ...
415 : ! **************************************************************************************************
416 10 : SUBROUTINE transport_periodicity_label(perd, label)
417 : INTEGER, DIMENSION(3), INTENT(IN) :: perd
418 : CHARACTER(LEN=*), INTENT(OUT) :: label
419 :
420 10 : label = ""
421 10 : IF (perd(1) == 1) label = TRIM(label)//"X"
422 10 : IF (perd(2) == 1) label = TRIM(label)//"Y"
423 10 : IF (perd(3) == 1) label = TRIM(label)//"Z"
424 10 : IF (LEN_TRIM(label) == 0) label = "NONE"
425 :
426 10 : END SUBROUTINE transport_periodicity_label
427 :
428 : ! **************************************************************************************************
429 : !> \brief Output unit labels for a transport dimensionality.
430 : !> \param transport_ndim ...
431 : !> \param density_unit ...
432 : !> \param sigma_iso_unit ...
433 : !> \param sigma_tensor_unit ...
434 : ! **************************************************************************************************
435 10 : SUBROUTINE transport_unit_labels(transport_ndim, density_unit, sigma_iso_unit, sigma_tensor_unit)
436 : INTEGER, INTENT(IN) :: transport_ndim
437 : CHARACTER(LEN=*), INTENT(OUT) :: density_unit, sigma_iso_unit, &
438 : sigma_tensor_unit
439 :
440 12 : SELECT CASE (transport_ndim)
441 : CASE (1)
442 2 : density_unit = "m^-1"
443 2 : sigma_iso_unit = "S*m"
444 2 : sigma_tensor_unit = "S*m"
445 : CASE (2)
446 2 : density_unit = "m^-2"
447 2 : sigma_iso_unit = "S"
448 2 : sigma_tensor_unit = "S"
449 : CASE DEFAULT
450 6 : density_unit = "m^-3"
451 6 : sigma_iso_unit = "S/cm"
452 6 : sigma_tensor_unit = "S/m"
453 : END SELECT
454 :
455 10 : END SUBROUTINE transport_unit_labels
456 :
457 : ! **************************************************************************************************
458 : !> \brief Output conversion factors from internal dimensional transport units.
459 : !> \param transport_ndim ...
460 : !> \param iso_factor ...
461 : !> \param tensor_factor ...
462 : ! **************************************************************************************************
463 10 : SUBROUTINE transport_output_factors(transport_ndim, iso_factor, tensor_factor)
464 : INTEGER, INTENT(IN) :: transport_ndim
465 : REAL(KIND=dp), INTENT(OUT) :: iso_factor, tensor_factor
466 :
467 12 : SELECT CASE (transport_ndim)
468 : CASE (1)
469 2 : iso_factor = 1.0E-10_dp
470 2 : tensor_factor = 1.0E-10_dp
471 : CASE (2)
472 2 : iso_factor = 1.0_dp
473 2 : tensor_factor = 1.0_dp
474 : CASE DEFAULT
475 6 : iso_factor = 1.0E8_dp
476 10 : tensor_factor = 1.0E10_dp
477 : END SELECT
478 :
479 10 : END SUBROUTINE transport_output_factors
480 :
481 : ! **************************************************************************************************
482 : !> \brief Prepare eigenvalues and the transformed [X,H] matrices for one spin channel.
483 : !> \param ks_matrix ...
484 : !> \param s_dense ...
485 : !> \param blacs_env ...
486 : !> \param para_env ...
487 : !> \param particle_set ...
488 : !> \param cell ...
489 : !> \param projection ...
490 : !> \param ao_atom ...
491 : !> \param mo_set ...
492 : !> \param allow_mo_reuse ...
493 : !> \param eig ...
494 : !> \param dHH ...
495 : !> \param reused_mos ...
496 : ! **************************************************************************************************
497 10 : SUBROUTINE setup_spin_transport(ks_matrix, s_dense, blacs_env, para_env, particle_set, cell, &
498 10 : projection, ao_atom, mo_set, allow_mo_reuse, eig, dHH, reused_mos)
499 : TYPE(dbcsr_type), INTENT(IN) :: ks_matrix
500 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: s_dense
501 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
502 : TYPE(mp_para_env_type), POINTER :: para_env
503 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
504 : TYPE(cell_type), POINTER :: cell
505 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: projection
506 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_atom
507 : TYPE(mo_set_type), INTENT(IN) :: mo_set
508 : LOGICAL, INTENT(IN) :: allow_mo_reuse
509 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eig
510 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: dHH
511 : LOGICAL, INTENT(OUT) :: reused_mos
512 :
513 : INTEGER :: idir, mo_nao, nao, nmo
514 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_dense, gmat, h_dense, horth, op, &
515 10 : op_orth, uvec, work
516 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
517 : TYPE(cp_fm_type), POINTER :: mo_coeff
518 :
519 10 : NULLIFY (mo_coeff, mo_eigenvalues)
520 10 : nao = SIZE(s_dense, 1)
521 10 : reused_mos = .FALSE.
522 :
523 10 : CALL dbcsr_to_dense(ks_matrix, blacs_env, para_env, h_dense)
524 60 : ALLOCATE (op(nao, nao), work(nao, nao))
525 :
526 10 : IF (allow_mo_reuse) THEN
527 : CALL get_mo_set(mo_set=mo_set, nao=mo_nao, nmo=nmo, eigenvalues=mo_eigenvalues, &
528 2 : mo_coeff=mo_coeff)
529 : reused_mos = mo_nao == nao .AND. nmo >= nao .AND. ASSOCIATED(mo_eigenvalues) .AND. &
530 2 : ASSOCIATED(mo_coeff)
531 2 : IF (reused_mos) reused_mos = SIZE(mo_eigenvalues) >= nao
532 : END IF
533 :
534 10 : IF (reused_mos) THEN
535 6 : ALLOCATE (coeff_dense(nao, nao))
536 22 : eig(:) = mo_eigenvalues(1:nao)
537 2 : CALL cp_fm_get_submatrix(mo_coeff, coeff_dense, n_rows=nao, n_cols=nao)
538 8 : DO idir = 1, 3
539 6 : CALL build_commutator_kernel(h_dense, particle_set, cell, projection, ao_atom, idir, op)
540 8 : CALL transform_to_eigenbasis(op, coeff_dense, dHH(:, :, idir), work)
541 : END DO
542 2 : DEALLOCATE (coeff_dense)
543 : ELSE
544 72 : ALLOCATE (gmat(nao, nao), horth(nao, nao), op_orth(nao, nao), uvec(nao, nao))
545 8 : CALL inverse_sqrt_symmetric(s_dense, gmat)
546 :
547 16888 : work(:, :) = MATMUL(h_dense, gmat)
548 8888 : horth(:, :) = MATMUL(gmat, work)
549 8 : CALL symmetrize(horth)
550 :
551 888 : uvec(:, :) = horth
552 8 : CALL diagonalize_symmetric(uvec, eig)
553 :
554 32 : DO idir = 1, 3
555 24 : CALL build_commutator_kernel(h_dense, particle_set, cell, projection, ao_atom, idir, op)
556 26664 : work(:, :) = MATMUL(op, gmat)
557 26664 : op_orth(:, :) = MATMUL(gmat, work)
558 32 : CALL transform_to_eigenbasis(op_orth, uvec, dHH(:, :, idir), work)
559 : END DO
560 8 : DEALLOCATE (gmat, horth, op_orth, uvec)
561 : END IF
562 :
563 10 : DEALLOCATE (h_dense, op, work)
564 :
565 20 : END SUBROUTINE setup_spin_transport
566 :
567 : ! **************************************************************************************************
568 : !> \brief Symmetric inverse square root via dense diagonalization.
569 : !> \param matrix ...
570 : !> \param inverse_sqrt ...
571 : ! **************************************************************************************************
572 8 : SUBROUTINE inverse_sqrt_symmetric(matrix, inverse_sqrt)
573 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
574 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: inverse_sqrt
575 :
576 : INTEGER :: i, info, lwork, n
577 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig, work
578 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: scaled, vectors
579 :
580 8 : n = SIZE(matrix, 1)
581 8 : lwork = MAX(1, 8*n)
582 80 : ALLOCATE (eig(n), scaled(n, n), vectors(n, n), work(lwork))
583 :
584 888 : vectors(:, :) = matrix
585 8 : CALL dsyev("V", "L", n, vectors, n, eig, work, lwork, info)
586 8 : IF (info /= 0) CPABORT("Diagonalization of the overlap matrix failed.")
587 88 : IF (MINVAL(eig) <= 0.0_dp) CPABORT("Overlap matrix is not positive definite.")
588 :
589 888 : scaled(:, :) = vectors
590 88 : DO i = 1, n
591 888 : scaled(:, i) = scaled(:, i)/SQRT(eig(i))
592 : END DO
593 8 : CALL dgemm("N", "T", n, n, n, 1.0_dp, scaled, n, vectors, n, 0.0_dp, inverse_sqrt, n)
594 :
595 8 : DEALLOCATE (eig, scaled, vectors, work)
596 :
597 8 : END SUBROUTINE inverse_sqrt_symmetric
598 :
599 : ! **************************************************************************************************
600 : !> \brief Dense symmetric diagonalization.
601 : !> \param matrix ...
602 : !> \param eig ...
603 : ! **************************************************************************************************
604 8 : SUBROUTINE diagonalize_symmetric(matrix, eig)
605 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix
606 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eig
607 :
608 : INTEGER :: info, lwork, n
609 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
610 :
611 8 : n = SIZE(matrix, 1)
612 8 : lwork = MAX(1, 8*n)
613 24 : ALLOCATE (work(lwork))
614 8 : CALL dsyev("V", "L", n, matrix, n, eig, work, lwork, info)
615 8 : IF (info /= 0) CPABORT("Diagonalization of the orthogonalized KS matrix failed.")
616 8 : DEALLOCATE (work)
617 :
618 8 : END SUBROUTINE diagonalize_symmetric
619 :
620 : ! **************************************************************************************************
621 : !> \brief Enforce exact symmetry after dense products.
622 : !> \param matrix ...
623 : ! **************************************************************************************************
624 8 : SUBROUTINE symmetrize(matrix)
625 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix
626 :
627 : INTEGER :: i, j, n
628 : REAL(KIND=dp) :: value
629 :
630 8 : n = SIZE(matrix, 1)
631 88 : DO j = 1, n
632 448 : DO i = j + 1, n
633 360 : value = 0.5_dp*(matrix(i, j) + matrix(j, i))
634 360 : matrix(i, j) = value
635 440 : matrix(j, i) = value
636 : END DO
637 : END DO
638 :
639 8 : END SUBROUTINE symmetrize
640 :
641 : ! **************************************************************************************************
642 : !> \brief Build element-wise kernel (P*(r_col-r_row))_idir * matrix(row,col).
643 : !> \param matrix ...
644 : !> \param particle_set ...
645 : !> \param cell ...
646 : !> \param projection ...
647 : !> \param ao_atom ...
648 : !> \param idir ...
649 : !> \param op ...
650 : ! **************************************************************************************************
651 30 : SUBROUTINE build_commutator_kernel(matrix, particle_set, cell, projection, ao_atom, idir, op)
652 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
653 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
654 : TYPE(cell_type), POINTER :: cell
655 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: projection
656 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_atom
657 : INTEGER, INTENT(IN) :: idir
658 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: op
659 :
660 : INTEGER :: i, j, n
661 : REAL(KIND=dp), DIMENSION(3) :: dr, projected_dr
662 :
663 30 : n = SIZE(matrix, 1)
664 330 : DO j = 1, n
665 3330 : DO i = 1, n
666 12000 : dr(:) = pbc(particle_set(ao_atom(j))%r(:) - particle_set(ao_atom(i))%r(:), cell)
667 39000 : projected_dr(:) = MATMUL(projection, dr)
668 3300 : op(i, j) = projected_dr(idir)*matrix(i, j)
669 : END DO
670 : END DO
671 :
672 30 : END SUBROUTINE build_commutator_kernel
673 :
674 : ! **************************************************************************************************
675 : !> \brief Transform an AO matrix to the eigenvector basis.
676 : !> \param op ...
677 : !> \param uvec ...
678 : !> \param transformed ...
679 : !> \param work ...
680 : ! **************************************************************************************************
681 30 : SUBROUTINE transform_to_eigenbasis(op, uvec, transformed, work)
682 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: op, uvec
683 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: transformed
684 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: work
685 :
686 : INTEGER :: n
687 :
688 30 : n = SIZE(op, 1)
689 30 : CALL dgemm("N", "N", n, n, n, 1.0_dp, op, n, uvec, n, 0.0_dp, work, n)
690 30 : CALL dgemm("T", "N", n, n, n, 1.0_dp, uvec, n, work, n, 0.0_dp, transformed, n)
691 :
692 30 : END SUBROUTINE transform_to_eigenbasis
693 :
694 : ! **************************************************************************************************
695 : !> \brief Fermi occupation with overflow protection.
696 : !> \param eig ...
697 : !> \param mu ...
698 : !> \param kT ...
699 : !> \param maxocc ...
700 : !> \return ...
701 : ! **************************************************************************************************
702 200 : PURE FUNCTION fermi_occ(eig, mu, kT, maxocc) RESULT(occ)
703 : REAL(KIND=dp), INTENT(IN) :: eig, mu, kT, maxocc
704 : REAL(KIND=dp) :: occ
705 :
706 : REAL(KIND=dp) :: x
707 :
708 200 : x = (eig - mu)/kT
709 200 : IF (x > 80.0_dp) THEN
710 : occ = 0.0_dp
711 20 : ELSEIF (x < -80.0_dp) THEN
712 0 : occ = maxocc
713 : ELSE
714 20 : occ = maxocc/(1.0_dp + EXP(x))
715 : END IF
716 :
717 200 : END FUNCTION fermi_occ
718 :
719 : ! **************************************************************************************************
720 : !> \brief Fermi vacancy below the neutral point.
721 : !> \param eig ...
722 : !> \param mu ...
723 : !> \param kT ...
724 : !> \param maxocc ...
725 : !> \return ...
726 : ! **************************************************************************************************
727 10 : PURE FUNCTION fermi_hole(eig, mu, kT, maxocc) RESULT(hole)
728 : REAL(KIND=dp), INTENT(IN) :: eig, mu, kT, maxocc
729 : REAL(KIND=dp) :: hole
730 :
731 10 : hole = maxocc - fermi_occ(eig, mu, kT, maxocc)
732 :
733 10 : END FUNCTION fermi_hole
734 :
735 : ! **************************************************************************************************
736 : !> \brief Total electron count at a chemical potential.
737 : !> \param eig ...
738 : !> \param maxocc ...
739 : !> \param mu ...
740 : !> \param kT ...
741 : !> \return ...
742 : ! **************************************************************************************************
743 0 : PURE FUNCTION electron_count(eig, maxocc, mu, kT) RESULT(count)
744 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
745 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
746 : REAL(KIND=dp), INTENT(IN) :: mu, kT
747 : REAL(KIND=dp) :: count
748 :
749 : INTEGER :: ispin, n
750 :
751 0 : count = 0.0_dp
752 0 : DO ispin = 1, SIZE(eig, 2)
753 0 : DO n = 1, SIZE(eig, 1)
754 0 : count = count + fermi_occ(eig(n, ispin), mu, kT, maxocc(ispin))
755 : END DO
756 : END DO
757 :
758 0 : END FUNCTION electron_count
759 :
760 : ! **************************************************************************************************
761 : !> \brief Determine the neutral chemical potential using the legacy two-step definition.
762 : !> \param eig ...
763 : !> \param maxocc ...
764 : !> \param kT ...
765 : !> \param neutral_electrons ...
766 : !> \param neutral_grid ...
767 : !> \param emin ...
768 : !> \param emax ...
769 : !> \param mu0 ...
770 : ! **************************************************************************************************
771 0 : SUBROUTINE find_neutral_mu(eig, maxocc, kT, neutral_electrons, neutral_grid, emin, emax, mu0)
772 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
773 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
774 : REAL(KIND=dp), INTENT(IN) :: kT, neutral_electrons
775 : INTEGER, INTENT(IN) :: neutral_grid
776 : REAL(KIND=dp), INTENT(IN) :: emin, emax
777 : REAL(KIND=dp), INTENT(OUT) :: mu0
778 :
779 : INTEGER :: i
780 : REAL(KIND=dp) :: best, count, mu, mup, ne, nh, &
781 : previous_count
782 :
783 0 : mup = 0.5_dp*(emin + emax)
784 0 : previous_count = electron_count(eig, maxocc, emin, kT)
785 0 : DO i = 1, neutral_grid
786 0 : mu = emin + (emax - emin)*REAL(i, KIND=dp)/REAL(neutral_grid, KIND=dp)
787 0 : count = electron_count(eig, maxocc, mu, kT)
788 0 : IF ((count - neutral_electrons)*(previous_count - neutral_electrons) <= 0.0_dp) mup = mu
789 0 : previous_count = count
790 : END DO
791 :
792 0 : mu0 = mup
793 0 : best = HUGE(1.0_dp)
794 0 : DO i = 1, neutral_grid
795 0 : mu = emin + (emax - emin)*REAL(i, KIND=dp)/REAL(neutral_grid, KIND=dp)
796 0 : CALL electron_hole_density(eig, maxocc, mu, mup, kT, ne, nh)
797 0 : IF (ABS(nh - ne) <= best) THEN
798 0 : mu0 = mu
799 0 : best = ABS(nh - ne)
800 : END IF
801 : END DO
802 :
803 0 : END SUBROUTINE find_neutral_mu
804 :
805 : ! **************************************************************************************************
806 : !> \brief Electron and hole counts relative to a neutral point.
807 : !> \param eig ...
808 : !> \param maxocc ...
809 : !> \param mu ...
810 : !> \param mu0 ...
811 : !> \param kT ...
812 : !> \param ne ...
813 : !> \param nh ...
814 : ! **************************************************************************************************
815 10 : SUBROUTINE electron_hole_density(eig, maxocc, mu, mu0, kT, ne, nh)
816 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
817 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
818 : REAL(KIND=dp), INTENT(IN) :: mu, mu0, kT
819 : REAL(KIND=dp), INTENT(OUT) :: ne, nh
820 :
821 : INTEGER :: ispin, n
822 :
823 10 : ne = 0.0_dp
824 10 : nh = 0.0_dp
825 20 : DO ispin = 1, SIZE(eig, 2)
826 120 : DO n = 1, SIZE(eig, 1)
827 100 : IF (eig(n, ispin) <= mu0) nh = nh + fermi_hole(eig(n, ispin), mu, kT, maxocc(ispin))
828 110 : IF (eig(n, ispin) >= mu0) ne = ne + fermi_occ(eig(n, ispin), mu, kT, maxocc(ispin))
829 : END DO
830 : END DO
831 :
832 10 : END SUBROUTINE electron_hole_density
833 :
834 : ! **************************************************************************************************
835 : !> \brief Conductivity tensor at one chemical potential.
836 : !> \param eig ...
837 : !> \param dHH ...
838 : !> \param maxocc ...
839 : !> \param mu ...
840 : !> \param kT ...
841 : !> \param dissipation ...
842 : !> \param measure ...
843 : !> \param transport_ndim ...
844 : !> \param sigma ...
845 : ! **************************************************************************************************
846 10 : SUBROUTINE conductivity_at_mu(eig, dHH, maxocc, mu, kT, dissipation, measure, transport_ndim, sigma)
847 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
848 : REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dHH
849 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: maxocc
850 : REAL(KIND=dp), INTENT(IN) :: mu, kT, dissipation, measure
851 : INTEGER, INTENT(IN) :: transport_ndim
852 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: sigma
853 :
854 : INTEGER :: idir, ispin, jdir, m, n, nao
855 : REAL(KIND=dp) :: delta, eps, factor, g0
856 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occ
857 :
858 10 : nao = SIZE(eig, 1)
859 10 : eps = 1.0E-12_dp
860 10 : g0 = e_charge**2/(pi*h_bar)
861 10 : sigma = 0.0_dp
862 :
863 30 : ALLOCATE (occ(nao))
864 :
865 20 : DO ispin = 1, SIZE(eig, 2)
866 110 : DO n = 1, nao
867 110 : occ(n) = fermi_occ(eig(n, ispin), mu, kT, maxocc(ispin))
868 : END DO
869 :
870 50 : DO idir = 1, 3
871 130 : DO jdir = 1, 3
872 1020 : DO n = 1, nao
873 9990 : DO m = 1, nao
874 9000 : delta = eig(n, ispin) - eig(m, ispin)
875 9000 : factor = (occ(n) - occ(m))/(delta + eps)*dissipation/(dissipation**2 + delta**2)
876 : sigma(jdir, idir) = sigma(jdir, idir) + &
877 9900 : dHH(m, n, jdir, ispin)*dHH(n, m, idir, ispin)*factor
878 : END DO
879 : END DO
880 : END DO
881 : END DO
882 : END DO
883 :
884 130 : sigma = pi*g0*sigma*angstrom**(2 - transport_ndim)/measure
885 :
886 10 : DEALLOCATE (occ)
887 :
888 10 : END SUBROUTINE conductivity_at_mu
889 :
890 : END MODULE qs_kubo_transport
|
