(git:f56c6e3)
Loading...
Searching...
No Matches
negf_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief NEGF based quantum transport calculations
10! **************************************************************************************************
11MODULE negf_methods
12 USE bibliography, ONLY: bailey2006,&
13 papior2017,&
14 cite_reference
15 USE cp_blacs_env, ONLY: cp_blacs_env_type
16 USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale,&
17 cp_cfm_scale_and_add,&
18 cp_cfm_trace
19 USE cp_cfm_types, ONLY: &
20 copy_cfm_info_type, cp_cfm_cleanup_copy_general, cp_cfm_create, &
21 cp_cfm_finish_copy_general, cp_cfm_get_info, cp_cfm_get_submatrix, cp_cfm_release, &
22 cp_cfm_set_submatrix, cp_cfm_start_copy_general, cp_cfm_to_fm, cp_cfm_type
23 USE cp_control_types, ONLY: dft_control_type
24 USE cp_dbcsr_api, ONLY: dbcsr_copy,&
25 dbcsr_deallocate_matrix,&
26 dbcsr_init_p,&
27 dbcsr_p_type
28 USE cp_dbcsr_contrib, ONLY: dbcsr_dot
29 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
30 USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
31 cp_fm_scale_and_add,&
32 cp_fm_trace
33 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
34 cp_fm_struct_release,&
35 cp_fm_struct_type
36 USE cp_fm_types, ONLY: cp_fm_copy_general,&
37 cp_fm_create,&
38 cp_fm_get_info,&
39 cp_fm_release,&
40 cp_fm_set_all,&
41 cp_fm_to_fm,&
42 cp_fm_type
43 USE cp_log_handling, ONLY: cp_get_default_logger,&
44 cp_logger_get_default_io_unit,&
45 cp_logger_type
46 USE cp_output_handling, ONLY: cp_p_file,&
47 cp_print_key_finished_output,&
48 cp_print_key_should_output,&
49 cp_print_key_unit_nr,&
50 debug_print_level,&
51 high_print_level
52 USE cp_subsys_types, ONLY: cp_subsys_type
53 USE force_env_types, ONLY: force_env_get,&
54 force_env_p_type,&
55 force_env_type
56 USE global_types, ONLY: global_environment_type
57 USE input_constants, ONLY: negfint_method_cc,&
58 negfint_method_simpson
59 USE input_section_types, ONLY: section_vals_get_subs_vals,&
60 section_vals_type,&
61 section_vals_val_get
62 USE kinds, ONLY: default_string_length,&
63 dp
64 USE kpoint_types, ONLY: get_kpoint_info,&
65 kpoint_type
66 USE machine, ONLY: m_walltime
67 USE mathconstants, ONLY: pi,&
68 twopi,&
69 z_one,&
70 z_zero
71 USE message_passing, ONLY: mp_para_env_type
72 USE negf_control_types, ONLY: negf_control_create,&
73 negf_control_release,&
74 negf_control_type,&
75 read_negf_control
76 USE negf_env_types, ONLY: negf_env_create,&
77 negf_env_release,&
78 negf_env_type
79 USE negf_green_cache, ONLY: green_functions_cache_expand,&
80 green_functions_cache_release,&
81 green_functions_cache_reorder,&
82 green_functions_cache_type
83 USE negf_green_methods, ONLY: do_sancho,&
84 negf_contact_broadening_matrix,&
85 negf_contact_self_energy,&
86 negf_retarded_green_function,&
87 sancho_work_matrices_create,&
88 sancho_work_matrices_release,&
89 sancho_work_matrices_type
90 USE negf_integr_cc, ONLY: &
91 cc_interval_full, cc_interval_half, cc_shape_arc, cc_shape_linear, &
92 ccquad_double_number_of_points, ccquad_init, ccquad_reduce_and_append_zdata, &
93 ccquad_refine_integral, ccquad_release, ccquad_type
94 USE negf_integr_simpson, ONLY: simpsonrule_get_next_nodes,&
95 simpsonrule_init,&
96 simpsonrule_refine_integral,&
97 simpsonrule_release,&
98 simpsonrule_type,&
99 sr_shape_arc,&
100 sr_shape_linear
101 USE negf_matrix_utils, ONLY: invert_cell_to_index,&
102 negf_copy_fm_submat_to_dbcsr,&
103 negf_copy_sym_dbcsr_to_fm_submat
104 USE negf_subgroup_types, ONLY: negf_sub_env_create,&
105 negf_sub_env_release,&
106 negf_subgroup_env_type
107 USE parallel_gemm_api, ONLY: parallel_gemm
108 USE physcon, ONLY: e_charge,&
109 evolt,&
110 kelvin,&
111 seconds
112 USE qs_density_mixing_types, ONLY: direct_mixing_nr,&
113 gspace_mixing_nr
114 USE qs_energy_types, ONLY: qs_energy_type
115 USE qs_environment_types, ONLY: get_qs_env,&
116 qs_environment_type
117 USE qs_gspace_mixing, ONLY: gspace_mixing
118 USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
119 USE qs_mixing_utils, ONLY: mixing_allocate,&
120 mixing_init
121 USE qs_rho_methods, ONLY: qs_rho_update_rho
122 USE qs_rho_types, ONLY: qs_rho_get,&
123 qs_rho_type
124 USE qs_scf_methods, ONLY: scf_env_density_mixing
125 USE qs_subsys_types, ONLY: qs_subsys_type
126 USE string_utilities, ONLY: integer_to_string
127#include "./base/base_uses.f90"
128
129 IMPLICIT NONE
130 PRIVATE
131
132 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_methods'
133 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .true.
134
135 PUBLIC :: do_negf
136
137! **************************************************************************************************
138!> \brief Type to accumulate the total number of points used in integration as well as
139!> the final error estimate
140!> \author Sergey Chulkov
141! **************************************************************************************************
142 TYPE integration_status_type
143 INTEGER :: npoints = -1
144 REAL(kind=dp) :: error = -1.0_dp
145 END TYPE integration_status_type
146
147CONTAINS
148
149! **************************************************************************************************
150!> \brief Perform NEGF calculation.
151!> \param force_env Force environment
152!> \par History
153!> * 01.2017 created [Sergey Chulkov]
154! **************************************************************************************************
155 SUBROUTINE do_negf(force_env)
156 TYPE(force_env_type), POINTER :: force_env
157
158 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_negf'
159
160 CHARACTER(len=100) :: sfmt
161 CHARACTER(len=default_string_length) :: contact_id_str
162 INTEGER :: handle, i, icontact, ispin, j, k, &
163 log_unit, n, ncontacts, npoints, &
164 nspins, print_level, print_unit
165 LOGICAL :: debug_output, should_output, &
166 verbose_output
167 REAL(kind=dp) :: energy_max, energy_min
168 REAL(kind=dp), DIMENSION(2) :: current
169 TYPE(cp_blacs_env_type), POINTER :: blacs_env
170 TYPE(cp_logger_type), POINTER :: logger
171 TYPE(cp_subsys_type), POINTER :: cp_subsys
172 TYPE(dft_control_type), POINTER :: dft_control
173 TYPE(force_env_p_type), DIMENSION(:), POINTER :: sub_force_env
174 TYPE(global_environment_type), POINTER :: global_env
175 TYPE(negf_control_type), POINTER :: negf_control
176 TYPE(negf_env_type) :: negf_env
177 TYPE(negf_subgroup_env_type) :: sub_env
178 TYPE(qs_environment_type), POINTER :: qs_env
179 TYPE(section_vals_type), POINTER :: negf_contact_section, &
180 negf_mixing_section, negf_section, &
181 print_section, root_section
182
183 CALL timeset(routinen, handle)
184 logger => cp_get_default_logger()
185 log_unit = cp_logger_get_default_io_unit()
186
187 CALL cite_reference(bailey2006)
188 CALL cite_reference(papior2017)
189
190 NULLIFY (blacs_env, cp_subsys, global_env, qs_env, root_section, sub_force_env)
191 CALL force_env_get(force_env, globenv=global_env, qs_env=qs_env, root_section=root_section, &
192 sub_force_env=sub_force_env, subsys=cp_subsys)
193
194 CALL get_qs_env(qs_env, blacs_env=blacs_env)
195
196 negf_section => section_vals_get_subs_vals(root_section, "NEGF")
197 negf_contact_section => section_vals_get_subs_vals(negf_section, "CONTACT")
198 negf_mixing_section => section_vals_get_subs_vals(negf_section, "MIXING")
199
200 NULLIFY (negf_control)
201 CALL negf_control_create(negf_control)
202 CALL read_negf_control(negf_control, root_section, cp_subsys)
203 CALL get_qs_env(qs_env, dft_control=dft_control)
204
205 ! print unit, if log_unit > 0, otherwise no output
206 log_unit = cp_print_key_unit_nr(logger, negf_section, "PRINT%PROGRAM_RUN_INFO", extension=".Log")
207
208 IF (log_unit > 0) THEN
209 WRITE (log_unit, '(/,T2,79("-"))')
210 WRITE (log_unit, '(T27,A,T62)') "NEGF calculation is started"
211 WRITE (log_unit, '(T2,79("-"))')
212 END IF
213
214 ! print levels, are used if log_unit > 0
215 IF (log_unit > 0) THEN
216 CALL section_vals_val_get(negf_section, "PRINT%PROGRAM_RUN_INFO%PRINT_LEVEL", i_val=print_level)
217 SELECT CASE (print_level)
218 CASE (high_print_level)
219 verbose_output = .true.
220 CASE (debug_print_level)
221 verbose_output = .true.
222 debug_output = .true.
223 CASE DEFAULT
224 verbose_output = .false.
225 debug_output = .false.
226 END SELECT
227 END IF
228
229 IF (log_unit > 0) THEN
230 WRITE (log_unit, "(/,' THE RELEVANT HAMILTONIAN AND OVERLAP MATRICES FROM DFT')")
231 WRITE (log_unit, "( ' ------------------------------------------------------')")
232 END IF
233
234 CALL negf_sub_env_create(sub_env, negf_control, blacs_env, global_env%blacs_grid_layout, global_env%blacs_repeatable)
235 CALL negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
236
237 IF (log_unit > 0) THEN
238 WRITE (log_unit, "(/,' NEGF| The initial Hamiltonian and Overlap matrices are calculated.')")
239 END IF
240
241 IF (log_unit > 0) THEN
242 DO icontact = 1, SIZE(negf_control%contacts)
243 WRITE (log_unit, "(/,' The electrode',I5)") icontact
244 WRITE (log_unit, "( ' ------------------')")
245 WRITE (log_unit, "(' From the force environment:',I16)") negf_control%contacts(icontact)%force_env_index
246 WRITE (log_unit, "(' Number of atoms:',I27)") SIZE(negf_control%contacts(icontact)%atomlist_bulk)
247 IF (verbose_output) WRITE (log_unit, "(' Atoms belonging to a contact (from the entire sysem):')")
248 IF (verbose_output) WRITE (log_unit, "(16I5)") negf_control%contacts(icontact)%atomlist_bulk
249 WRITE (log_unit, "(' Number of atoms in a primary unit cell:',I4)") SIZE(negf_env%contacts(icontact)%atomlist_cell0)
250 IF (verbose_output) WRITE (log_unit, "(' Atoms belonging to a primary unit cell (from the entire sysem):')")
251 IF (verbose_output) WRITE (log_unit, "(16I5)") negf_env%contacts(icontact)%atomlist_cell0
252 IF (verbose_output) WRITE (log_unit, "(' Direction of an electrode: ',I16)") negf_env%contacts(icontact)%direction_axis
253 n = SIZE(negf_env%contacts(icontact)%h_00(1)%local_data, 1)
254 WRITE (sfmt, "('(',i0,'(E15.5))')") n
255 WRITE (log_unit, "(' The number of atomic orbtals:',I14)") n
256 ! print the electrode Hamiltonians for check and debuging
257 IF (debug_output) THEN
258 DO k = 1, dft_control%nspins
259 WRITE (log_unit, "(' The H_00 electrode Hamiltonian for spin',I2)") k
260 DO i = 1, n
261 WRITE (log_unit, sfmt) (negf_env%contacts(icontact)%h_00(k)%local_data(i, j), j=1, n)
262 END DO
263 WRITE (log_unit, "(' The H_01 electrode Hamiltonian for spin',I2)") k
264 DO i = 1, n
265 WRITE (log_unit, sfmt) (negf_env%contacts(icontact)%h_01(k)%local_data(i, j), j=1, n)
266 END DO
267 END DO
268 WRITE (log_unit, "(' The S_00 overlap matrix')")
269 DO i = 1, n
270 WRITE (log_unit, sfmt) (negf_env%contacts(icontact)%s_00%local_data(i, j), j=1, n)
271 END DO
272 WRITE (log_unit, "(' The S_01 overlap matrix')")
273 DO i = 1, n
274 WRITE (log_unit, sfmt) (negf_env%contacts(icontact)%s_01%local_data(i, j), j=1, n)
275 END DO
276 END IF
277 END DO
278 WRITE (log_unit, "(/,' The full scattering region')")
279 WRITE (log_unit, "( ' --------------------------')")
280 WRITE (log_unit, "(' Number of atoms:',I27)") SIZE(negf_control%atomlist_S_screening)
281 IF (verbose_output) WRITE (log_unit, "(' Atoms belonging to a full scattering region:')")
282 IF (verbose_output) WRITE (log_unit, "(16I5)") negf_control%atomlist_S_screening
283 n = SIZE(negf_env%h_s(1)%local_data, 1)
284 WRITE (sfmt, "('(',i0,'(E15.5))')") n
285 WRITE (log_unit, "(' The number of atomic orbtals:',I14)") n
286 ! print the full scattering region Hamiltonians for check and debuging
287 IF (debug_output) THEN
288 DO k = 1, dft_control%nspins
289 WRITE (log_unit, "(' The H_s Hamiltonian for spin',I2)") k
290 DO i = 1, n
291 WRITE (log_unit, sfmt) (negf_env%h_s(k)%local_data(i, j), j=1, n)
292 END DO
293 END DO
294 WRITE (log_unit, "(' The S_s overlap matrix')")
295 DO i = 1, n
296 WRITE (log_unit, sfmt) (negf_env%s_s%local_data(i, j), j=1, n)
297 END DO
298 END IF
299 ! print the electrode Hamiltonians for check and debuging
300 IF (debug_output) THEN
301 WRITE (log_unit, "(/,' Scattering region - electrode contacts')")
302 WRITE (log_unit, "( ' ---------------------------------------')")
303 DO icontact = 1, SIZE(negf_control%contacts)
304 WRITE (log_unit, "(/,' The contact',I5)") icontact
305 WRITE (log_unit, "( ' ----------------')")
306 DO k = 1, dft_control%nspins
307 WRITE (log_unit, "(' The H_sc Hamiltonian for spin',I2)") k
308 DO i = 1, n
309 WRITE (log_unit, sfmt) &
310 (negf_env%h_sc(k, icontact)%local_data(i, j), j=1, SIZE(negf_env%contacts(icontact)%h_00(1)%local_data, 1))
311 END DO
312 END DO
313 WRITE (log_unit, "(' The S_sc overlap matrix')")
314 DO i = 1, n
315 WRITE (log_unit, sfmt) &
316 (negf_env%s_sc(icontact)%local_data(i, j), j=1, SIZE(negf_env%contacts(icontact)%h_00(1)%local_data, 1))
317 END DO
318 END DO
319 END IF
320 WRITE (log_unit, *)
321 END IF
322
323 ! compute contact Fermi level as well as requested properties
324 ! -----------------------------------------------------------
325 ncontacts = SIZE(negf_control%contacts)
326 DO icontact = 1, ncontacts
327 NULLIFY (qs_env)
328 IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
329 CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env)
330 ELSE
331 CALL force_env_get(force_env, qs_env=qs_env)
332 END IF
333 CALL guess_fermi_level(icontact, negf_env, negf_control, sub_env, qs_env, log_unit)
334
335 print_section => section_vals_get_subs_vals(negf_contact_section, "PRINT", i_rep_section=icontact)
336 should_output = btest(cp_print_key_should_output(logger%iter_info, print_section, "DOS"), cp_p_file)
337
338 IF (should_output) THEN
339 CALL section_vals_val_get(print_section, "DOS%FROM_ENERGY", r_val=energy_min)
340 CALL section_vals_val_get(print_section, "DOS%TILL_ENERGY", r_val=energy_max)
341 CALL section_vals_val_get(print_section, "DOS%N_GRIDPOINTS", i_val=npoints)
342
343 CALL integer_to_string(icontact, contact_id_str)
344 print_unit = cp_print_key_unit_nr(logger, print_section, "DOS", &
345 extension=".dos", &
346 middle_name=trim(adjustl(contact_id_str)), &
347 file_status="REPLACE")
348
349 CALL negf_print_dos(print_unit, energy_min, energy_max, npoints, &
350 v_shift=0.0_dp, negf_env=negf_env, negf_control=negf_control, &
351 sub_env=sub_env, base_contact=icontact, just_contact=icontact)
352
353 CALL cp_print_key_finished_output(print_unit, logger, print_section, "DOS")
354 END IF
355 END DO
356
357 ! NEGF procedure
358 ! --------------
359 IF (ncontacts > 1) THEN
360 CALL force_env_get(force_env, qs_env=qs_env)
361
362 ! shift potential
363 ! ---------------
364 CALL shift_potential(negf_env, negf_control, sub_env, qs_env, base_contact=1, log_unit=log_unit)
365
366 ! self-consistent density
367 ! -----------------------
368 CALL converge_density(negf_env, negf_control, sub_env, qs_env, negf_control%v_shift, base_contact=1, log_unit=log_unit)
369
370 ! current
371 ! -------
372 CALL get_qs_env(qs_env, dft_control=dft_control)
373
374 nspins = dft_control%nspins
375
376 cpassert(nspins <= 2)
377 DO ispin = 1, nspins
378 ! compute the electric current flown through a pair of electrodes
379 ! contact_id1 -> extended molecule -> contact_id2.
380 ! Only extended systems with two electrodes are supported at the moment,
381 ! so for the time being the contacts' indices are hardcoded.
382 current(ispin) = negf_compute_current(contact_id1=1, contact_id2=2, &
383 v_shift=negf_control%v_shift, &
384 negf_env=negf_env, &
385 negf_control=negf_control, &
386 sub_env=sub_env, &
387 ispin=ispin, &
388 blacs_env_global=blacs_env)
389 END DO
390
391 IF (log_unit > 0) THEN
392 IF (nspins > 1) THEN
393 WRITE (log_unit, '(/,T2,A,T60,ES20.7E2)') "NEGF| Alpha-spin electric current (A)", current(1)
394 WRITE (log_unit, '(T2,A,T60,ES20.7E2)') "NEGF| Beta-spin electric current (A)", current(2)
395 ELSE
396 WRITE (log_unit, '(/,T2,A,T60,ES20.7E2)') "NEGF| Electric current (A)", 2.0_dp*current(1)
397 END IF
398 END IF
399
400 ! density of states
401 ! -----------------
402 print_section => section_vals_get_subs_vals(negf_section, "PRINT")
403 should_output = btest(cp_print_key_should_output(logger%iter_info, print_section, "DOS"), cp_p_file)
404
405 IF (should_output) THEN
406 CALL section_vals_val_get(print_section, "DOS%FROM_ENERGY", r_val=energy_min)
407 CALL section_vals_val_get(print_section, "DOS%TILL_ENERGY", r_val=energy_max)
408 CALL section_vals_val_get(print_section, "DOS%N_GRIDPOINTS", i_val=npoints)
409
410 CALL integer_to_string(0, contact_id_str)
411 print_unit = cp_print_key_unit_nr(logger, print_section, "DOS", &
412 extension=".dos", &
413 middle_name=trim(adjustl(contact_id_str)), &
414 file_status="REPLACE")
415
416 CALL negf_print_dos(print_unit, energy_min, energy_max, npoints, negf_control%v_shift, &
417 negf_env=negf_env, negf_control=negf_control, &
418 sub_env=sub_env, base_contact=1)
419
420 CALL cp_print_key_finished_output(print_unit, logger, print_section, "DOS")
421 END IF
422
423 ! transmission coefficient
424 ! ------------------------
425 should_output = btest(cp_print_key_should_output(logger%iter_info, print_section, "TRANSMISSION"), cp_p_file)
426
427 IF (should_output) THEN
428 CALL section_vals_val_get(print_section, "TRANSMISSION%FROM_ENERGY", r_val=energy_min)
429 CALL section_vals_val_get(print_section, "TRANSMISSION%TILL_ENERGY", r_val=energy_max)
430 CALL section_vals_val_get(print_section, "TRANSMISSION%N_GRIDPOINTS", i_val=npoints)
431
432 CALL integer_to_string(0, contact_id_str)
433 print_unit = cp_print_key_unit_nr(logger, print_section, "TRANSMISSION", &
434 extension=".transm", &
435 middle_name=trim(adjustl(contact_id_str)), &
436 file_status="REPLACE")
437
438 CALL negf_print_transmission(print_unit, energy_min, energy_max, npoints, negf_control%v_shift, &
439 negf_env=negf_env, negf_control=negf_control, &
440 sub_env=sub_env, contact_id1=1, contact_id2=2)
441
442 CALL cp_print_key_finished_output(print_unit, logger, print_section, "TRANSMISSION")
443 END IF
444
445 END IF
446
447 IF (log_unit > 0) THEN
448 WRITE (log_unit, '(/,T2,79("-"))')
449 WRITE (log_unit, '(T27,A,T62)') "NEGF calculation is finished"
450 WRITE (log_unit, '(T2,79("-"))')
451 END IF
452
453 CALL negf_env_release(negf_env)
454 CALL negf_sub_env_release(sub_env)
455
456 CALL negf_control_release(negf_control)
457 CALL timestop(handle)
458 END SUBROUTINE do_negf
459
460! **************************************************************************************************
461!> \brief Compute the contact's Fermi level.
462!> \param contact_id index of the contact
463!> \param negf_env NEGF environment
464!> \param negf_control NEGF control
465!> \param sub_env NEGF parallel (sub)group environment
466!> \param qs_env QuickStep environment
467!> \param log_unit output unit
468!> \par History
469!> * 10.2017 created [Sergey Chulkov]
470! **************************************************************************************************
471 SUBROUTINE guess_fermi_level(contact_id, negf_env, negf_control, sub_env, qs_env, log_unit)
472 INTEGER, INTENT(in) :: contact_id
473 TYPE(negf_env_type), INTENT(in) :: negf_env
474 TYPE(negf_control_type), POINTER :: negf_control
475 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
476 TYPE(qs_environment_type), POINTER :: qs_env
477 INTEGER, INTENT(in) :: log_unit
478
479 CHARACTER(LEN=*), PARAMETER :: routinen = 'guess_fermi_level'
480 TYPE(cp_fm_type), PARAMETER :: fm_dummy = cp_fm_type()
481
482 CHARACTER(len=default_string_length) :: temperature_str
483 COMPLEX(kind=dp) :: lbound_cpath, lbound_lpath, ubound_lpath
484 INTEGER :: direction_axis_abs, handle, image, &
485 ispin, nao, nimages, nspins, step
486 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
487 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
488 LOGICAL :: do_kpoints
489 REAL(kind=dp) :: delta_au, energy_ubound_minus_fermi, fermi_level_guess, fermi_level_max, &
490 fermi_level_min, nelectrons_guess, nelectrons_max, nelectrons_min, nelectrons_qs_cell0, &
491 nelectrons_qs_cell1, offset_au, rscale, t1, t2, trace
492 TYPE(cp_blacs_env_type), POINTER :: blacs_env_global
493 TYPE(cp_fm_struct_type), POINTER :: fm_struct
494 TYPE(cp_fm_type) :: rho_ao_fm
495 TYPE(cp_fm_type), POINTER :: matrix_s_fm
496 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_qs_kp
497 TYPE(dft_control_type), POINTER :: dft_control
498 TYPE(green_functions_cache_type) :: g_surf_cache
499 TYPE(integration_status_type) :: stats
500 TYPE(kpoint_type), POINTER :: kpoints
501 TYPE(mp_para_env_type), POINTER :: para_env_global
502 TYPE(qs_energy_type), POINTER :: energy
503 TYPE(qs_rho_type), POINTER :: rho_struct
504 TYPE(qs_subsys_type), POINTER :: subsys
505
506 CALL timeset(routinen, handle)
507
508 IF (negf_control%contacts(contact_id)%compute_fermi_level) THEN
509 CALL get_qs_env(qs_env, &
510 blacs_env=blacs_env_global, &
511 dft_control=dft_control, &
512 do_kpoints=do_kpoints, &
513 kpoints=kpoints, &
514 matrix_s_kp=matrix_s_kp, &
515 para_env=para_env_global, &
516 rho=rho_struct, subsys=subsys)
517 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_qs_kp)
518
519 nimages = dft_control%nimages
520 nspins = dft_control%nspins
521 direction_axis_abs = abs(negf_env%contacts(contact_id)%direction_axis)
522
523 cpassert(SIZE(negf_env%contacts(contact_id)%h_00) == nspins)
524
525 IF (sub_env%ngroups > 1) THEN
526 NULLIFY (matrix_s_fm, fm_struct)
527
528 CALL cp_fm_get_info(negf_env%contacts(contact_id)%s_00, nrow_global=nao)
529 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env_global)
530 CALL cp_fm_create(rho_ao_fm, fm_struct)
531
532 ALLOCATE (matrix_s_fm)
533 CALL cp_fm_create(matrix_s_fm, fm_struct)
534 CALL cp_fm_struct_release(fm_struct)
535
536 IF (sub_env%group_distribution(sub_env%mepos_global) == 0) THEN
537 CALL cp_fm_copy_general(negf_env%contacts(contact_id)%s_00, matrix_s_fm, para_env_global)
538 ELSE
539 CALL cp_fm_copy_general(fm_dummy, matrix_s_fm, para_env_global)
540 END IF
541 ELSE
542 matrix_s_fm => negf_env%contacts(contact_id)%s_00
543 CALL cp_fm_get_info(matrix_s_fm, matrix_struct=fm_struct)
544 CALL cp_fm_create(rho_ao_fm, fm_struct)
545 END IF
546
547 IF (do_kpoints) THEN
548 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
549 ELSE
550 ALLOCATE (cell_to_index(0:0, 0:0, 0:0))
551 cell_to_index(0, 0, 0) = 1
552 END IF
553
554 ALLOCATE (index_to_cell(3, nimages))
555 CALL invert_cell_to_index(cell_to_index, nimages, index_to_cell)
556 IF (.NOT. do_kpoints) DEALLOCATE (cell_to_index)
557
558 IF (nspins == 1) THEN
559 ! spin-restricted calculation: number of electrons must be doubled
560 rscale = 2.0_dp
561 ELSE
562 rscale = 1.0_dp
563 END IF
564
565 ! compute the refence number of electrons using the electron density
566 nelectrons_qs_cell0 = 0.0_dp
567 nelectrons_qs_cell1 = 0.0_dp
568 DO image = 1, nimages
569 IF (index_to_cell(direction_axis_abs, image) == 0) THEN
570 DO ispin = 1, nspins
571 CALL dbcsr_dot(rho_ao_qs_kp(ispin, image)%matrix, matrix_s_kp(1, image)%matrix, trace)
572 nelectrons_qs_cell0 = nelectrons_qs_cell0 + trace
573 END DO
574 ELSE IF (abs(index_to_cell(direction_axis_abs, image)) == 1) THEN
575 DO ispin = 1, nspins
576 CALL dbcsr_dot(rho_ao_qs_kp(ispin, image)%matrix, matrix_s_kp(1, image)%matrix, trace)
577 nelectrons_qs_cell1 = nelectrons_qs_cell1 + trace
578 END DO
579 END IF
580 END DO
581 DEALLOCATE (index_to_cell)
582
583 IF (log_unit > 0) THEN
584 WRITE (temperature_str, '(F11.3)') negf_control%contacts(contact_id)%temperature*kelvin
585 WRITE (log_unit, '(/,T2,A,I3)') "COMPUTE FERMI LEVEL OF CONTACT ", contact_id
586 WRITE (log_unit, "( ' ----------------------------------')")
587 WRITE (log_unit, '(A)') " Temperature "//trim(adjustl(temperature_str))//" Kelvin"
588 WRITE (log_unit, '(/,T2,A,T60,F20.10,/)') "Electronic density of the isolated contact unit cell:", &
589 -1.0_dp*(nelectrons_qs_cell0 + nelectrons_qs_cell1)
590 WRITE (log_unit, '(T3,A)') "Step Integration method Time Fermi level Convergence (density)"
591 WRITE (log_unit, '(T3,78("-"))')
592 END IF
593
594 nelectrons_qs_cell0 = 0.0_dp
595 DO ispin = 1, nspins
596 CALL cp_fm_trace(negf_env%contacts(contact_id)%rho_00(ispin), &
597 negf_env%contacts(contact_id)%s_00, trace)
598 nelectrons_qs_cell0 = nelectrons_qs_cell0 + trace
599 END DO
600
601 ! Use the Fermi level given in the input file or the Fermi level of bulk electrodes as a reference point
602 ! and then refine the Fermi level by using a simple linear interpolation technique
603 CALL get_qs_env(qs_env, energy=energy)
604 IF (negf_control%homo_lumo_gap > 0.0_dp) THEN
605 IF (negf_control%contacts(contact_id)%refine_fermi_level) THEN
606 fermi_level_min = negf_control%contacts(contact_id)%fermi_level
607 ELSE
608 fermi_level_min = energy%efermi
609 END IF
610 fermi_level_max = fermi_level_min + negf_control%homo_lumo_gap
611 ELSE
612 IF (negf_control%contacts(contact_id)%refine_fermi_level) THEN
613 fermi_level_max = negf_control%contacts(contact_id)%fermi_level
614 ELSE
615 fermi_level_max = energy%efermi
616 END IF
617 fermi_level_min = fermi_level_max + negf_control%homo_lumo_gap
618 END IF
619
620 step = 0
621 lbound_cpath = cmplx(negf_control%energy_lbound, negf_control%eta, kind=dp)
622 delta_au = real(negf_control%delta_npoles, kind=dp)*twopi*negf_control%contacts(contact_id)%temperature
623 offset_au = real(negf_control%gamma_kT, kind=dp)*negf_control%contacts(contact_id)%temperature
624 energy_ubound_minus_fermi = -2.0_dp*log(negf_control%conv_density)*negf_control%contacts(contact_id)%temperature
625 t1 = m_walltime()
626
627 DO
628 step = step + 1
629
630 SELECT CASE (step)
631 CASE (1)
632 fermi_level_guess = fermi_level_min
633 CASE (2)
634 fermi_level_guess = fermi_level_max
635 CASE DEFAULT
636 fermi_level_guess = fermi_level_min - (nelectrons_min - nelectrons_qs_cell0)* &
637 (fermi_level_max - fermi_level_min)/(nelectrons_max - nelectrons_min)
638 END SELECT
639
640 negf_control%contacts(contact_id)%fermi_level = fermi_level_guess
641 nelectrons_guess = 0.0_dp
642
643 lbound_lpath = cmplx(fermi_level_guess - offset_au, delta_au, kind=dp)
644 ubound_lpath = cmplx(fermi_level_guess + energy_ubound_minus_fermi, delta_au, kind=dp)
645
646 CALL integration_status_reset(stats)
647
648 DO ispin = 1, nspins
649 CALL negf_init_rho_equiv_residuals(rho_ao_fm=rho_ao_fm, &
650 v_shift=0.0_dp, &
651 ignore_bias=.true., &
652 negf_env=negf_env, &
653 negf_control=negf_control, &
654 sub_env=sub_env, &
655 ispin=ispin, &
656 base_contact=contact_id, &
657 just_contact=contact_id)
658
659 CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_fm, &
660 stats=stats, &
661 v_shift=0.0_dp, &
662 ignore_bias=.true., &
663 negf_env=negf_env, &
664 negf_control=negf_control, &
665 sub_env=sub_env, &
666 ispin=ispin, &
667 base_contact=contact_id, &
668 integr_lbound=lbound_cpath, &
669 integr_ubound=lbound_lpath, &
670 matrix_s_global=matrix_s_fm, &
671 is_circular=.true., &
672 g_surf_cache=g_surf_cache, &
673 just_contact=contact_id)
674 CALL green_functions_cache_release(g_surf_cache)
675
676 CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_fm, &
677 stats=stats, &
678 v_shift=0.0_dp, &
679 ignore_bias=.true., &
680 negf_env=negf_env, &
681 negf_control=negf_control, &
682 sub_env=sub_env, &
683 ispin=ispin, &
684 base_contact=contact_id, &
685 integr_lbound=lbound_lpath, &
686 integr_ubound=ubound_lpath, &
687 matrix_s_global=matrix_s_fm, &
688 is_circular=.false., &
689 g_surf_cache=g_surf_cache, &
690 just_contact=contact_id)
691 CALL green_functions_cache_release(g_surf_cache)
692
693 CALL cp_fm_trace(rho_ao_fm, matrix_s_fm, trace)
694 nelectrons_guess = nelectrons_guess + trace
695 END DO
696 nelectrons_guess = nelectrons_guess*rscale
697
698 t2 = m_walltime()
699
700 IF (log_unit > 0) THEN
701 WRITE (log_unit, '(T2,I5,T12,A,T32,F8.1,T42,F15.8,T60,ES20.5E2)') &
702 step, get_method_description_string(stats, negf_control%integr_method), &
703 t2 - t1, fermi_level_guess, nelectrons_guess - nelectrons_qs_cell0
704 END IF
705
706 IF (abs(nelectrons_qs_cell0 - nelectrons_guess) < negf_control%conv_density) EXIT
707
708 SELECT CASE (step)
709 CASE (1)
710 nelectrons_min = nelectrons_guess
711 CASE (2)
712 nelectrons_max = nelectrons_guess
713 CASE DEFAULT
714 IF (fermi_level_guess < fermi_level_min) THEN
715 fermi_level_max = fermi_level_min
716 nelectrons_max = nelectrons_min
717 fermi_level_min = fermi_level_guess
718 nelectrons_min = nelectrons_guess
719 ELSE IF (fermi_level_guess > fermi_level_max) THEN
720 fermi_level_min = fermi_level_max
721 nelectrons_min = nelectrons_max
722 fermi_level_max = fermi_level_guess
723 nelectrons_max = nelectrons_guess
724 ELSE IF (fermi_level_max - fermi_level_guess < fermi_level_guess - fermi_level_min) THEN
725 fermi_level_max = fermi_level_guess
726 nelectrons_max = nelectrons_guess
727 ELSE
728 fermi_level_min = fermi_level_guess
729 nelectrons_min = nelectrons_guess
730 END IF
731 END SELECT
732
733 t1 = t2
734 END DO
735
736 negf_control%contacts(contact_id)%fermi_level = fermi_level_guess
737
738 IF (sub_env%ngroups > 1) THEN
739 CALL cp_fm_release(matrix_s_fm)
740 DEALLOCATE (matrix_s_fm)
741 END IF
742 CALL cp_fm_release(rho_ao_fm)
743 END IF
744
745 IF (log_unit > 0) THEN
746 WRITE (temperature_str, '(F11.3)') negf_control%contacts(contact_id)%temperature*kelvin
747 WRITE (log_unit, '(/,T2,A,I0)') "NEGF| Contact No. ", contact_id
748 WRITE (log_unit, '(T2,A,T62,F18.8)') "NEGF| Fermi level at "//trim(adjustl(temperature_str))// &
749 " Kelvin (a.u.):", negf_control%contacts(contact_id)%fermi_level
750 WRITE (log_unit, '(T2,A,T62,F18.8)') "NEGF| Electric potential (V):", &
751 negf_control%contacts(contact_id)%v_external*evolt
752 END IF
753
754 CALL timestop(handle)
755 END SUBROUTINE guess_fermi_level
756
757! **************************************************************************************************
758!> \brief Compute shift in Hartree potential
759!> \param negf_env NEGF environment
760!> \param negf_control NEGF control
761!> \param sub_env NEGF parallel (sub)group environment
762!> \param qs_env QuickStep environment
763!> \param base_contact index of the reference contact
764!> \param log_unit output unit
765! **************************************************************************************************
766 SUBROUTINE shift_potential(negf_env, negf_control, sub_env, qs_env, base_contact, log_unit)
767 TYPE(negf_env_type), INTENT(in) :: negf_env
768 TYPE(negf_control_type), POINTER :: negf_control
769 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
770 TYPE(qs_environment_type), POINTER :: qs_env
771 INTEGER, INTENT(in) :: base_contact, log_unit
772
773 CHARACTER(LEN=*), PARAMETER :: routinen = 'shift_potential'
774 TYPE(cp_fm_type), PARAMETER :: fm_dummy = cp_fm_type()
775
776 COMPLEX(kind=dp) :: lbound_cpath, ubound_cpath, ubound_lpath
777 INTEGER :: handle, ispin, iter_count, nao, &
778 ncontacts, nspins
779 LOGICAL :: do_kpoints
780 REAL(kind=dp) :: mu_base, nelectrons_guess, nelectrons_max, nelectrons_min, nelectrons_ref, &
781 t1, t2, temperature, trace, v_shift_guess, v_shift_max, v_shift_min
782 TYPE(cp_blacs_env_type), POINTER :: blacs_env
783 TYPE(cp_fm_struct_type), POINTER :: fm_struct
784 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: rho_ao_fm
785 TYPE(cp_fm_type), POINTER :: matrix_s_fm
786 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_qs_kp
787 TYPE(dft_control_type), POINTER :: dft_control
788 TYPE(green_functions_cache_type), ALLOCATABLE, &
789 DIMENSION(:) :: g_surf_circular, g_surf_linear
790 TYPE(integration_status_type) :: stats
791 TYPE(mp_para_env_type), POINTER :: para_env
792 TYPE(qs_rho_type), POINTER :: rho_struct
793 TYPE(qs_subsys_type), POINTER :: subsys
794
795 ncontacts = SIZE(negf_control%contacts)
796 ! nothing to do
797 IF (.NOT. (ALLOCATED(negf_env%h_s) .AND. ALLOCATED(negf_env%h_sc) .AND. &
798 ASSOCIATED(negf_env%s_s) .AND. ALLOCATED(negf_env%s_sc))) RETURN
799 IF (ncontacts < 2) RETURN
800
801 CALL timeset(routinen, handle)
802
803 CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, dft_control=dft_control, &
804 para_env=para_env, rho=rho_struct, subsys=subsys)
805 cpassert(.NOT. do_kpoints)
806
807 ! apply external NEGF potential
808 t1 = m_walltime()
809
810 ! need a globally distributed overlap matrix in order to compute integration errors
811 IF (sub_env%ngroups > 1) THEN
812 NULLIFY (matrix_s_fm, fm_struct)
813
814 CALL cp_fm_get_info(negf_env%s_s, nrow_global=nao)
815 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
816
817 ALLOCATE (matrix_s_fm)
818 CALL cp_fm_create(matrix_s_fm, fm_struct)
819 CALL cp_fm_struct_release(fm_struct)
820
821 IF (sub_env%group_distribution(sub_env%mepos_global) == 0) THEN
822 CALL cp_fm_copy_general(negf_env%s_s, matrix_s_fm, para_env)
823 ELSE
824 CALL cp_fm_copy_general(fm_dummy, matrix_s_fm, para_env)
825 END IF
826 ELSE
827 matrix_s_fm => negf_env%s_s
828 END IF
829
830 CALL cp_fm_get_info(matrix_s_fm, matrix_struct=fm_struct)
831
832 nspins = SIZE(negf_env%h_s)
833
834 mu_base = negf_control%contacts(base_contact)%fermi_level
835
836 ! keep the initial charge density matrix and Kohn-Sham matrix
837 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_qs_kp)
838
839 ! extract the reference density matrix blocks
840 nelectrons_ref = 0.0_dp
841 ALLOCATE (rho_ao_fm(nspins))
842 DO ispin = 1, nspins
843 CALL cp_fm_create(rho_ao_fm(ispin), fm_struct)
844
845 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_qs_kp(ispin, 1)%matrix, &
846 fm=rho_ao_fm(ispin), &
847 atomlist_row=negf_control%atomlist_S_screening, &
848 atomlist_col=negf_control%atomlist_S_screening, &
849 subsys=subsys, mpi_comm_global=para_env, &
850 do_upper_diag=.true., do_lower=.true.)
851
852 CALL cp_fm_trace(rho_ao_fm(ispin), matrix_s_fm, trace)
853 nelectrons_ref = nelectrons_ref + trace
854 END DO
855
856 IF (log_unit > 0) THEN
857 WRITE (log_unit, '(/,T2,A)') "COMPUTE SHIFT IN HARTREE POTENTIAL"
858 WRITE (log_unit, "( ' ----------------------------------')")
859 WRITE (log_unit, '(/,T2,A,T55,F25.14,/)') "Initial electronic density of the scattering region:", -1.0_dp*nelectrons_ref
860 WRITE (log_unit, '(T3,A)') "Step Integration method Time V shift Convergence (density)"
861 WRITE (log_unit, '(T3,78("-"))')
862 END IF
863
864 temperature = negf_control%contacts(base_contact)%temperature
865
866 ! integration limits: C-path (arch)
867 lbound_cpath = cmplx(negf_control%energy_lbound, negf_control%eta, kind=dp)
868 ubound_cpath = cmplx(mu_base - real(negf_control%gamma_kT, kind=dp)*temperature, &
869 REAL(negf_control%delta_npoles, kind=dp)*twopi*temperature, kind=dp)
870
871 ! integration limits: L-path (linear)
872 ubound_lpath = cmplx(mu_base - log(negf_control%conv_density)*temperature, &
873 REAL(negf_control%delta_npoles, kind=dp)*twopi*temperature, kind=dp)
874
875 v_shift_min = negf_control%v_shift
876 v_shift_max = negf_control%v_shift + negf_control%v_shift_offset
877
878 ALLOCATE (g_surf_circular(nspins), g_surf_linear(nspins))
879
880 DO iter_count = 1, negf_control%v_shift_maxiters
881 SELECT CASE (iter_count)
882 CASE (1)
883 v_shift_guess = v_shift_min
884 CASE (2)
885 v_shift_guess = v_shift_max
886 CASE DEFAULT
887 v_shift_guess = v_shift_min - (nelectrons_min - nelectrons_ref)* &
888 (v_shift_max - v_shift_min)/(nelectrons_max - nelectrons_min)
889 END SELECT
890
891 ! compute an updated density matrix
892 CALL integration_status_reset(stats)
893
894 DO ispin = 1, nspins
895 ! closed contour: residuals
896 CALL negf_init_rho_equiv_residuals(rho_ao_fm=rho_ao_fm(ispin), &
897 v_shift=v_shift_guess, &
898 ignore_bias=.true., &
899 negf_env=negf_env, &
900 negf_control=negf_control, &
901 sub_env=sub_env, &
902 ispin=ispin, &
903 base_contact=base_contact)
904
905 ! closed contour: C-path
906 CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_fm(ispin), &
907 stats=stats, &
908 v_shift=v_shift_guess, &
909 ignore_bias=.true., &
910 negf_env=negf_env, &
911 negf_control=negf_control, &
912 sub_env=sub_env, &
913 ispin=ispin, &
914 base_contact=base_contact, &
915 integr_lbound=lbound_cpath, &
916 integr_ubound=ubound_cpath, &
917 matrix_s_global=matrix_s_fm, &
918 is_circular=.true., &
919 g_surf_cache=g_surf_circular(ispin))
920 IF (negf_control%disable_cache) &
921 CALL green_functions_cache_release(g_surf_circular(ispin))
922
923 ! closed contour: L-path
924 CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_fm(ispin), &
925 stats=stats, &
926 v_shift=v_shift_guess, &
927 ignore_bias=.true., &
928 negf_env=negf_env, &
929 negf_control=negf_control, &
930 sub_env=sub_env, &
931 ispin=ispin, &
932 base_contact=base_contact, &
933 integr_lbound=ubound_cpath, &
934 integr_ubound=ubound_lpath, &
935 matrix_s_global=matrix_s_fm, &
936 is_circular=.false., &
937 g_surf_cache=g_surf_linear(ispin))
938 IF (negf_control%disable_cache) &
939 CALL green_functions_cache_release(g_surf_linear(ispin))
940 END DO
941
942 IF (nspins > 1) THEN
943 DO ispin = 2, nspins
944 CALL cp_fm_scale_and_add(1.0_dp, rho_ao_fm(1), 1.0_dp, rho_ao_fm(ispin))
945 END DO
946 ELSE
947 CALL cp_fm_scale(2.0_dp, rho_ao_fm(1))
948 END IF
949
950 CALL cp_fm_trace(rho_ao_fm(1), matrix_s_fm, nelectrons_guess)
951
952 t2 = m_walltime()
953
954 IF (log_unit > 0) THEN
955 WRITE (log_unit, '(T2,I5,T12,A,T32,F8.1,T42,F15.8,T60,ES20.5E2)') &
956 iter_count, get_method_description_string(stats, negf_control%integr_method), &
957 t2 - t1, v_shift_guess, nelectrons_guess - nelectrons_ref
958 END IF
959
960 IF (abs(nelectrons_guess - nelectrons_ref) < negf_control%conv_scf) EXIT
961
962 ! compute correction
963 SELECT CASE (iter_count)
964 CASE (1)
965 nelectrons_min = nelectrons_guess
966 CASE (2)
967 nelectrons_max = nelectrons_guess
968 CASE DEFAULT
969 IF (v_shift_guess < v_shift_min) THEN
970 v_shift_max = v_shift_min
971 nelectrons_max = nelectrons_min
972 v_shift_min = v_shift_guess
973 nelectrons_min = nelectrons_guess
974 ELSE IF (v_shift_guess > v_shift_max) THEN
975 v_shift_min = v_shift_max
976 nelectrons_min = nelectrons_max
977 v_shift_max = v_shift_guess
978 nelectrons_max = nelectrons_guess
979 ELSE IF (v_shift_max - v_shift_guess < v_shift_guess - v_shift_min) THEN
980 v_shift_max = v_shift_guess
981 nelectrons_max = nelectrons_guess
982 ELSE
983 v_shift_min = v_shift_guess
984 nelectrons_min = nelectrons_guess
985 END IF
986 END SELECT
987
988 t1 = t2
989 END DO
990
991 negf_control%v_shift = v_shift_guess
992
993 IF (log_unit > 0) THEN
994 WRITE (log_unit, '(T2,A,T62,F18.8)') "NEGF| Shift in Hartree potential", negf_control%v_shift
995 END IF
996
997 DO ispin = nspins, 1, -1
998 CALL green_functions_cache_release(g_surf_circular(ispin))
999 CALL green_functions_cache_release(g_surf_linear(ispin))
1000 END DO
1001 DEALLOCATE (g_surf_circular, g_surf_linear)
1002
1003 CALL cp_fm_release(rho_ao_fm)
1004
1005 IF (sub_env%ngroups > 1 .AND. ASSOCIATED(matrix_s_fm)) THEN
1006 CALL cp_fm_release(matrix_s_fm)
1007 DEALLOCATE (matrix_s_fm)
1008 END IF
1009
1010 CALL timestop(handle)
1011 END SUBROUTINE shift_potential
1012
1013! **************************************************************************************************
1014!> \brief Converge electronic density of the scattering region.
1015!> \param negf_env NEGF environment
1016!> \param negf_control NEGF control
1017!> \param sub_env NEGF parallel (sub)group environment
1018!> \param qs_env QuickStep environment
1019!> \param v_shift shift in Hartree potential
1020!> \param base_contact index of the reference contact
1021!> \param log_unit output unit
1022!> \par History
1023!> * 06.2017 created [Sergey Chulkov]
1024! **************************************************************************************************
1025 SUBROUTINE converge_density(negf_env, negf_control, sub_env, qs_env, v_shift, base_contact, log_unit)
1026 TYPE(negf_env_type), INTENT(in) :: negf_env
1027 TYPE(negf_control_type), POINTER :: negf_control
1028 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
1029 TYPE(qs_environment_type), POINTER :: qs_env
1030 REAL(kind=dp), INTENT(in) :: v_shift
1031 INTEGER, INTENT(in) :: base_contact, log_unit
1032
1033 CHARACTER(LEN=*), PARAMETER :: routinen = 'converge_density'
1034 REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*epsilon(0.0_dp)
1035 TYPE(cp_fm_type), PARAMETER :: fm_dummy = cp_fm_type()
1036
1037 COMPLEX(kind=dp) :: lbound_cpath, ubound_cpath, ubound_lpath
1038 INTEGER :: handle, icontact, image, ispin, &
1039 iter_count, nao, ncontacts, nimages, &
1040 nspins
1041 LOGICAL :: do_kpoints
1042 REAL(kind=dp) :: iter_delta, mu_base, nelectrons, &
1043 nelectrons_diff, t1, t2, temperature, &
1044 trace, v_base, v_contact
1045 TYPE(cp_blacs_env_type), POINTER :: blacs_env
1046 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1047 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: rho_ao_delta_fm, rho_ao_new_fm
1048 TYPE(cp_fm_type), POINTER :: matrix_s_fm
1049 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_initial_kp, matrix_ks_qs_kp, &
1050 rho_ao_initial_kp, rho_ao_new_kp, &
1051 rho_ao_qs_kp
1052 TYPE(dft_control_type), POINTER :: dft_control
1053 TYPE(green_functions_cache_type), ALLOCATABLE, &
1054 DIMENSION(:) :: g_surf_circular, g_surf_linear, &
1055 g_surf_nonequiv
1056 TYPE(integration_status_type) :: stats
1057 TYPE(mp_para_env_type), POINTER :: para_env
1058 TYPE(qs_rho_type), POINTER :: rho_struct
1059 TYPE(qs_subsys_type), POINTER :: subsys
1060
1061 ncontacts = SIZE(negf_control%contacts)
1062 ! nothing to do
1063 IF (.NOT. (ALLOCATED(negf_env%h_s) .AND. ALLOCATED(negf_env%h_sc) .AND. &
1064 ASSOCIATED(negf_env%s_s) .AND. ALLOCATED(negf_env%s_sc))) RETURN
1065 IF (ncontacts < 2) RETURN
1066
1067 CALL timeset(routinen, handle)
1068
1069 CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, dft_control=dft_control, &
1070 matrix_ks_kp=matrix_ks_qs_kp, para_env=para_env, rho=rho_struct, subsys=subsys)
1071 cpassert(.NOT. do_kpoints)
1072
1073 ! apply external NEGF potential
1074 t1 = m_walltime()
1075
1076 ! need a globally distributed overlap matrix in order to compute integration errors
1077 IF (sub_env%ngroups > 1) THEN
1078 NULLIFY (matrix_s_fm, fm_struct)
1079
1080 CALL cp_fm_get_info(negf_env%s_s, nrow_global=nao)
1081 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env)
1082
1083 ALLOCATE (matrix_s_fm)
1084 CALL cp_fm_create(matrix_s_fm, fm_struct)
1085 CALL cp_fm_struct_release(fm_struct)
1086
1087 IF (sub_env%group_distribution(sub_env%mepos_global) == 0) THEN
1088 CALL cp_fm_copy_general(negf_env%s_s, matrix_s_fm, para_env)
1089 ELSE
1090 CALL cp_fm_copy_general(fm_dummy, matrix_s_fm, para_env)
1091 END IF
1092 ELSE
1093 matrix_s_fm => negf_env%s_s
1094 END IF
1095
1096 CALL cp_fm_get_info(matrix_s_fm, matrix_struct=fm_struct)
1097
1098 nspins = SIZE(negf_env%h_s)
1099 nimages = dft_control%nimages
1100
1101 v_base = negf_control%contacts(base_contact)%v_external
1102 mu_base = negf_control%contacts(base_contact)%fermi_level + v_base
1103
1104 ! the current subroutine works for the general case as well, but the Poisson solver does not
1105 IF (ncontacts > 2) THEN
1106 cpabort("Poisson solver does not support the general NEGF setup (>2 contacts).")
1107 END IF
1108
1109 ! keep the initial charge density matrix and Kohn-Sham matrix
1110 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_qs_kp)
1111
1112 NULLIFY (matrix_ks_initial_kp, rho_ao_initial_kp, rho_ao_new_kp)
1113 CALL dbcsr_allocate_matrix_set(matrix_ks_initial_kp, nspins, nimages)
1114 CALL dbcsr_allocate_matrix_set(rho_ao_initial_kp, nspins, nimages)
1115 CALL dbcsr_allocate_matrix_set(rho_ao_new_kp, nspins, nimages)
1116
1117 DO image = 1, nimages
1118 DO ispin = 1, nspins
1119 CALL dbcsr_init_p(matrix_ks_initial_kp(ispin, image)%matrix)
1120 CALL dbcsr_copy(matrix_b=matrix_ks_initial_kp(ispin, image)%matrix, matrix_a=matrix_ks_qs_kp(ispin, image)%matrix)
1121
1122 CALL dbcsr_init_p(rho_ao_initial_kp(ispin, image)%matrix)
1123 CALL dbcsr_copy(matrix_b=rho_ao_initial_kp(ispin, image)%matrix, matrix_a=rho_ao_qs_kp(ispin, image)%matrix)
1124
1125 CALL dbcsr_init_p(rho_ao_new_kp(ispin, image)%matrix)
1126 CALL dbcsr_copy(matrix_b=rho_ao_new_kp(ispin, image)%matrix, matrix_a=rho_ao_qs_kp(ispin, image)%matrix)
1127 END DO
1128 END DO
1129
1130 ! extract the reference density matrix blocks
1131 nelectrons = 0.0_dp
1132 ALLOCATE (rho_ao_delta_fm(nspins), rho_ao_new_fm(nspins))
1133 DO ispin = 1, nspins
1134 CALL cp_fm_create(rho_ao_delta_fm(ispin), fm_struct)
1135 CALL cp_fm_create(rho_ao_new_fm(ispin), fm_struct)
1136
1137 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_qs_kp(ispin, 1)%matrix, &
1138 fm=rho_ao_delta_fm(ispin), &
1139 atomlist_row=negf_control%atomlist_S_screening, &
1140 atomlist_col=negf_control%atomlist_S_screening, &
1141 subsys=subsys, mpi_comm_global=para_env, &
1142 do_upper_diag=.true., do_lower=.true.)
1143
1144 CALL cp_fm_trace(rho_ao_delta_fm(ispin), matrix_s_fm, trace)
1145 nelectrons = nelectrons + trace
1146 END DO
1147
1148 ! mixing storage allocation
1149 IF (negf_env%mixing_method >= gspace_mixing_nr) THEN
1150 CALL mixing_allocate(qs_env, negf_env%mixing_method, nspins=nspins, mixing_store=negf_env%mixing_storage)
1151 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
1152 cpabort('TB Code not available')
1153 ELSE IF (dft_control%qs_control%semi_empirical) THEN
1154 cpabort('SE Code not possible')
1155 ELSE
1156 CALL mixing_init(negf_env%mixing_method, rho_struct, negf_env%mixing_storage, para_env)
1157 END IF
1158 END IF
1159
1160 IF (log_unit > 0) THEN
1161 WRITE (log_unit, '(/,T2,A)') "NEGF SELF-CONSISTENT PROCEDURE"
1162 WRITE (log_unit, "( ' ------------------------------')")
1163 WRITE (log_unit, '(/,T2,A,T55,F25.14,/)') "Initial electronic density of the scattering region:", -1.0_dp*nelectrons
1164 WRITE (log_unit, '(T3,A)') "Step Integration method Time Electronic density Convergence"
1165 WRITE (log_unit, '(T3,78("-"))')
1166 END IF
1167
1168 temperature = negf_control%contacts(base_contact)%temperature
1169
1170 DO icontact = 1, ncontacts
1171 IF (icontact /= base_contact) THEN
1172 v_contact = negf_control%contacts(icontact)%v_external
1173
1174 ! integration limits: C-path (arch)
1175 lbound_cpath = cmplx(negf_control%energy_lbound, negf_control%eta, kind=dp)
1176 ubound_cpath = cmplx(mu_base - real(negf_control%gamma_kT, kind=dp)*temperature, &
1177 REAL(negf_control%delta_npoles, kind=dp)*twopi*temperature, kind=dp)
1178
1179 ! integration limits: L-path (linear)
1180 ubound_lpath = cmplx(mu_base - log(negf_control%conv_density)*temperature, &
1181 REAL(negf_control%delta_npoles, kind=dp)*twopi*temperature, kind=dp)
1182
1183 ALLOCATE (g_surf_circular(nspins), g_surf_linear(nspins), g_surf_nonequiv(nspins))
1184
1185 DO iter_count = 1, negf_control%max_scf
1186 ! compute an updated density matrix
1187 CALL integration_status_reset(stats)
1188
1189 DO ispin = 1, nspins
1190 ! closed contour: residuals
1191 CALL negf_init_rho_equiv_residuals(rho_ao_fm=rho_ao_new_fm(ispin), &
1192 v_shift=v_shift, &
1193 ignore_bias=.false., &
1194 negf_env=negf_env, &
1195 negf_control=negf_control, &
1196 sub_env=sub_env, &
1197 ispin=ispin, &
1198 base_contact=base_contact)
1199
1200 ! closed contour: C-path
1201 CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_new_fm(ispin), &
1202 stats=stats, &
1203 v_shift=v_shift, &
1204 ignore_bias=.false., &
1205 negf_env=negf_env, &
1206 negf_control=negf_control, &
1207 sub_env=sub_env, &
1208 ispin=ispin, &
1209 base_contact=base_contact, &
1210 integr_lbound=lbound_cpath, &
1211 integr_ubound=ubound_cpath, &
1212 matrix_s_global=matrix_s_fm, &
1213 is_circular=.true., &
1214 g_surf_cache=g_surf_circular(ispin))
1215 IF (negf_control%disable_cache) &
1216 CALL green_functions_cache_release(g_surf_circular(ispin))
1217
1218 ! closed contour: L-path
1219 CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_new_fm(ispin), &
1220 stats=stats, &
1221 v_shift=v_shift, &
1222 ignore_bias=.false., &
1223 negf_env=negf_env, &
1224 negf_control=negf_control, &
1225 sub_env=sub_env, &
1226 ispin=ispin, &
1227 base_contact=base_contact, &
1228 integr_lbound=ubound_cpath, &
1229 integr_ubound=ubound_lpath, &
1230 matrix_s_global=matrix_s_fm, &
1231 is_circular=.false., &
1232 g_surf_cache=g_surf_linear(ispin))
1233 IF (negf_control%disable_cache) &
1234 CALL green_functions_cache_release(g_surf_linear(ispin))
1235
1236 ! non-equilibrium part
1237 IF (abs(negf_control%contacts(icontact)%v_external - &
1238 negf_control%contacts(base_contact)%v_external) >= threshold) THEN
1239 CALL negf_add_rho_nonequiv(rho_ao_fm=rho_ao_new_fm(ispin), &
1240 stats=stats, &
1241 v_shift=v_shift, &
1242 negf_env=negf_env, &
1243 negf_control=negf_control, &
1244 sub_env=sub_env, &
1245 ispin=ispin, &
1246 base_contact=base_contact, &
1247 matrix_s_global=matrix_s_fm, &
1248 g_surf_cache=g_surf_nonequiv(ispin))
1249 IF (negf_control%disable_cache) &
1250 CALL green_functions_cache_release(g_surf_nonequiv(ispin))
1251 END IF
1252 END DO
1253
1254 IF (nspins == 1) CALL cp_fm_scale(2.0_dp, rho_ao_new_fm(1))
1255
1256 nelectrons = 0.0_dp
1257 nelectrons_diff = 0.0_dp
1258 DO ispin = 1, nspins
1259 CALL cp_fm_trace(rho_ao_new_fm(ispin), matrix_s_fm, trace)
1260 nelectrons = nelectrons + trace
1261
1262 ! rho_ao_delta_fm contains the original (non-mixed) density matrix from the previous iteration
1263 CALL cp_fm_scale_and_add(1.0_dp, rho_ao_delta_fm(ispin), -1.0_dp, rho_ao_new_fm(ispin))
1264 CALL cp_fm_trace(rho_ao_delta_fm(ispin), matrix_s_fm, trace)
1265 nelectrons_diff = nelectrons_diff + trace
1266
1267 ! rho_ao_new_fm -> rho_ao_delta_fm
1268 CALL cp_fm_to_fm(rho_ao_new_fm(ispin), rho_ao_delta_fm(ispin))
1269 END DO
1270
1271 t2 = m_walltime()
1272
1273 IF (log_unit > 0) THEN
1274 WRITE (log_unit, '(T2,I5,T12,A,T32,F8.1,T43,F20.8,T65,ES15.5E2)') &
1275 iter_count, get_method_description_string(stats, negf_control%integr_method), &
1276 t2 - t1, -1.0_dp*nelectrons, nelectrons_diff
1277 END IF
1278
1279 IF (abs(nelectrons_diff) < negf_control%conv_scf) EXIT
1280
1281 t1 = t2
1282
1283 ! mix density matrices
1284 IF (negf_env%mixing_method == direct_mixing_nr) THEN
1285 DO image = 1, nimages
1286 DO ispin = 1, nspins
1287 CALL dbcsr_copy(matrix_b=rho_ao_new_kp(ispin, image)%matrix, &
1288 matrix_a=rho_ao_initial_kp(ispin, image)%matrix)
1289 END DO
1290 END DO
1291
1292 DO ispin = 1, nspins
1293 CALL negf_copy_fm_submat_to_dbcsr(fm=rho_ao_new_fm(ispin), &
1294 matrix=rho_ao_new_kp(ispin, 1)%matrix, &
1295 atomlist_row=negf_control%atomlist_S_screening, &
1296 atomlist_col=negf_control%atomlist_S_screening, &
1297 subsys=subsys)
1298 END DO
1299
1300 CALL scf_env_density_mixing(rho_ao_new_kp, negf_env%mixing_storage, rho_ao_qs_kp, &
1301 para_env, iter_delta, iter_count)
1302
1303 DO image = 1, nimages
1304 DO ispin = 1, nspins
1305 CALL dbcsr_copy(rho_ao_qs_kp(ispin, image)%matrix, rho_ao_new_kp(ispin, image)%matrix)
1306 END DO
1307 END DO
1308 ELSE
1309 ! store the updated density matrix directly into the variable 'rho_ao_qs_kp'
1310 ! (which is qs_env%rho%rho_ao_kp); density mixing will be done on an inverse-space grid
1311 DO image = 1, nimages
1312 DO ispin = 1, nspins
1313 CALL dbcsr_copy(matrix_b=rho_ao_qs_kp(ispin, image)%matrix, &
1314 matrix_a=rho_ao_initial_kp(ispin, image)%matrix)
1315 END DO
1316 END DO
1317
1318 DO ispin = 1, nspins
1319 CALL negf_copy_fm_submat_to_dbcsr(fm=rho_ao_new_fm(ispin), &
1320 matrix=rho_ao_qs_kp(ispin, 1)%matrix, &
1321 atomlist_row=negf_control%atomlist_S_screening, &
1322 atomlist_col=negf_control%atomlist_S_screening, &
1323 subsys=subsys)
1324 END DO
1325 END IF
1326
1327 CALL qs_rho_update_rho(rho_struct, qs_env=qs_env)
1328
1329 IF (negf_env%mixing_method >= gspace_mixing_nr) THEN
1330 CALL gspace_mixing(qs_env, negf_env%mixing_method, negf_env%mixing_storage, &
1331 rho_struct, para_env, iter_count)
1332 END IF
1333
1334 ! update KS-matrix
1335 CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., just_energy=.false.)
1336
1337 ! extract blocks from the updated Kohn-Sham matrix
1338 DO ispin = 1, nspins
1339 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_qs_kp(ispin, 1)%matrix, &
1340 fm=negf_env%h_s(ispin), &
1341 atomlist_row=negf_control%atomlist_S_screening, &
1342 atomlist_col=negf_control%atomlist_S_screening, &
1343 subsys=subsys, mpi_comm_global=para_env, &
1344 do_upper_diag=.true., do_lower=.true.)
1345
1346 END DO
1347 END DO
1348
1349 IF (log_unit > 0) THEN
1350 IF (iter_count <= negf_control%max_scf) THEN
1351 WRITE (log_unit, '(/,T11,1X,A,I0,A)') "*** NEGF run converged in ", iter_count, " iteration(s) ***"
1352 ELSE
1353 WRITE (log_unit, '(/,T11,1X,A,I0,A)') "*** NEGF run did NOT converge after ", iter_count - 1, " iteration(s) ***"
1354 END IF
1355 END IF
1356
1357 DO ispin = nspins, 1, -1
1358 CALL green_functions_cache_release(g_surf_circular(ispin))
1359 CALL green_functions_cache_release(g_surf_linear(ispin))
1360 CALL green_functions_cache_release(g_surf_nonequiv(ispin))
1361 END DO
1362 DEALLOCATE (g_surf_circular, g_surf_linear, g_surf_nonequiv)
1363 END IF
1364 END DO
1365
1366 CALL cp_fm_release(rho_ao_new_fm)
1367 CALL cp_fm_release(rho_ao_delta_fm)
1368
1369 DO image = 1, nimages
1370 DO ispin = 1, nspins
1371 CALL dbcsr_copy(matrix_b=matrix_ks_qs_kp(ispin, image)%matrix, matrix_a=matrix_ks_initial_kp(ispin, image)%matrix)
1372 CALL dbcsr_copy(matrix_b=rho_ao_qs_kp(ispin, image)%matrix, matrix_a=rho_ao_initial_kp(ispin, image)%matrix)
1373
1374 CALL dbcsr_deallocate_matrix(matrix_ks_initial_kp(ispin, image)%matrix)
1375 CALL dbcsr_deallocate_matrix(rho_ao_initial_kp(ispin, image)%matrix)
1376 CALL dbcsr_deallocate_matrix(rho_ao_new_kp(ispin, image)%matrix)
1377 END DO
1378 END DO
1379 DEALLOCATE (matrix_ks_initial_kp, rho_ao_new_kp, rho_ao_initial_kp)
1380
1381 IF (sub_env%ngroups > 1 .AND. ASSOCIATED(matrix_s_fm)) THEN
1382 CALL cp_fm_release(matrix_s_fm)
1383 DEALLOCATE (matrix_s_fm)
1384 END IF
1385
1386 CALL timestop(handle)
1387 END SUBROUTINE converge_density
1388
1389! **************************************************************************************************
1390!> \brief Compute the surface retarded Green's function at a set of points in parallel.
1391!> \param g_surf set of surface Green's functions computed within the given parallel group
1392!> \param omega list of energy points where the surface Green's function need to be computed
1393!> \param h0 diagonal block of the Kohn-Sham matrix (must be Hermitian)
1394!> \param s0 diagonal block of the overlap matrix (must be Hermitian)
1395!> \param h1 off-fiagonal block of the Kohn-Sham matrix
1396!> \param s1 off-fiagonal block of the overlap matrix
1397!> \param sub_env NEGF parallel (sub)group environment
1398!> \param v_external applied electric potential
1399!> \param conv convergence threshold
1400!> \param transp flag which indicates that the matrices h1 and s1 should be transposed
1401!> \par History
1402!> * 07.2017 created [Sergey Chulkov]
1403! **************************************************************************************************
1404 SUBROUTINE negf_surface_green_function_batch(g_surf, omega, h0, s0, h1, s1, sub_env, v_external, conv, transp)
1405 TYPE(cp_cfm_type), DIMENSION(:), INTENT(inout) :: g_surf
1406 COMPLEX(kind=dp), DIMENSION(:), INTENT(in) :: omega
1407 TYPE(cp_fm_type), INTENT(IN) :: h0, s0, h1, s1
1408 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
1409 REAL(kind=dp), INTENT(in) :: v_external, conv
1410 LOGICAL, INTENT(in) :: transp
1411
1412 CHARACTER(len=*), PARAMETER :: routinen = 'negf_surface_green_function_batch'
1413 TYPE(cp_cfm_type), PARAMETER :: cfm_null = cp_cfm_type()
1414
1415 INTEGER :: handle, igroup, ipoint, npoints
1416 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1417 TYPE(sancho_work_matrices_type) :: work
1418
1419 CALL timeset(routinen, handle)
1420 npoints = SIZE(omega)
1421
1422 CALL cp_fm_get_info(s0, matrix_struct=fm_struct)
1423 CALL sancho_work_matrices_create(work, fm_struct)
1424
1425 igroup = sub_env%group_distribution(sub_env%mepos_global)
1426
1427 g_surf(1:npoints) = cfm_null
1428
1429 DO ipoint = igroup + 1, npoints, sub_env%ngroups
1430 IF (debug_this_module) THEN
1431 cpassert(.NOT. ASSOCIATED(g_surf(ipoint)%matrix_struct))
1432 END IF
1433 CALL cp_cfm_create(g_surf(ipoint), fm_struct)
1434
1435 CALL do_sancho(g_surf(ipoint), omega(ipoint) - v_external, &
1436 h0, s0, h1, s1, conv, transp, work)
1437 END DO
1438
1439 CALL sancho_work_matrices_release(work)
1440 CALL timestop(handle)
1441 END SUBROUTINE negf_surface_green_function_batch
1442
1443! **************************************************************************************************
1444!> \brief Compute the retarded Green's function and related properties at a set of points in parallel.
1445!> \param omega list of energy points
1446!> \param v_shift shift in Hartree potential
1447!> \param ignore_bias ignore v_external from negf_control
1448!> \param negf_env NEGF environment
1449!> \param negf_control NEGF control
1450!> \param sub_env (sub)group environment
1451!> \param ispin spin component to compute
1452!> \param g_surf_contacts set of surface Green's functions for every contact that computed
1453!> within the given parallel group
1454!> \param g_ret_s globally distributed matrices to store retarded Green's functions
1455!> \param g_ret_scale scale factor for retarded Green's functions
1456!> \param gamma_contacts 2-D array of globally distributed matrices to store broadening matrices
1457!> for every contact ([n_contacts, npoints])
1458!> \param gret_gamma_gadv 2-D array of globally distributed matrices to store the spectral function:
1459!> g_ret_s * gamma * g_ret_s^C for every contact ([n_contacts, n_points])
1460!> \param dos density of states at 'omega' ([n_points])
1461!> \param transm_coeff transmission coefficients between two contacts 'transm_contact1'
1462!> and 'transm_contact2' computed at points 'omega' ([n_points])
1463!> \param transm_contact1 index of the first contact
1464!> \param transm_contact2 index of the second contact
1465!> \param just_contact if present, compute the retarded Green's function of the system
1466!> lead1 -- device -- lead2. All 3 regions have the same Kohn-Sham
1467!> matrices which are taken from 'negf_env%contacts(just_contact)%h'.
1468!> Useful to apply NEGF procedure a single contact in order to compute
1469!> its Fermi level
1470!> \par History
1471!> * 07.2017 created [Sergey Chulkov]
1472! **************************************************************************************************
1473 SUBROUTINE negf_retarded_green_function_batch(omega, v_shift, ignore_bias, negf_env, negf_control, sub_env, ispin, &
1474 g_surf_contacts, &
1475 g_ret_s, g_ret_scale, gamma_contacts, gret_gamma_gadv, dos, &
1476 transm_coeff, transm_contact1, transm_contact2, just_contact)
1477 COMPLEX(kind=dp), DIMENSION(:), INTENT(in) :: omega
1478 REAL(kind=dp), INTENT(in) :: v_shift
1479 LOGICAL, INTENT(in) :: ignore_bias
1480 TYPE(negf_env_type), INTENT(in) :: negf_env
1481 TYPE(negf_control_type), POINTER :: negf_control
1482 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
1483 INTEGER, INTENT(in) :: ispin
1484 TYPE(cp_cfm_type), DIMENSION(:, :), INTENT(in) :: g_surf_contacts
1485 TYPE(cp_cfm_type), DIMENSION(:), INTENT(in), &
1486 OPTIONAL :: g_ret_s
1487 COMPLEX(kind=dp), DIMENSION(:), INTENT(in), &
1488 OPTIONAL :: g_ret_scale
1489 TYPE(cp_cfm_type), DIMENSION(:, :), INTENT(in), &
1490 OPTIONAL :: gamma_contacts, gret_gamma_gadv
1491 REAL(kind=dp), DIMENSION(:), INTENT(out), OPTIONAL :: dos
1492 COMPLEX(kind=dp), DIMENSION(:), INTENT(out), &
1493 OPTIONAL :: transm_coeff
1494 INTEGER, INTENT(in), OPTIONAL :: transm_contact1, transm_contact2, &
1495 just_contact
1496
1497 CHARACTER(len=*), PARAMETER :: routinen = 'negf_retarded_green_function_batch'
1498
1499 INTEGER :: handle, icontact, igroup, ipoint, &
1500 ncontacts, npoints, nrows
1501 REAL(kind=dp) :: v_external
1502 TYPE(copy_cfm_info_type), ALLOCATABLE, &
1503 DIMENSION(:) :: info1
1504 TYPE(copy_cfm_info_type), ALLOCATABLE, &
1505 DIMENSION(:, :) :: info2
1506 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: g_ret_s_group, self_energy_contacts, &
1507 zwork1_contacts, zwork2_contacts
1508 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:, :) :: gamma_contacts_group, &
1509 gret_gamma_gadv_group
1510 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1511 TYPE(cp_fm_type) :: g_ret_imag
1512 TYPE(cp_fm_type), POINTER :: matrix_s
1513 TYPE(mp_para_env_type), POINTER :: para_env
1514
1515 CALL timeset(routinen, handle)
1516 npoints = SIZE(omega)
1517 ncontacts = SIZE(negf_env%contacts)
1518 cpassert(SIZE(negf_control%contacts) == ncontacts)
1519
1520 IF (PRESENT(just_contact)) THEN
1521 cpassert(just_contact <= ncontacts)
1522 ncontacts = 2
1523 END IF
1524
1525 cpassert(ncontacts >= 2)
1526
1527 IF (ignore_bias) v_external = 0.0_dp
1528
1529 IF (PRESENT(transm_coeff) .OR. PRESENT(transm_contact1) .OR. PRESENT(transm_contact2)) THEN
1530 cpassert(PRESENT(transm_coeff))
1531 cpassert(PRESENT(transm_contact1))
1532 cpassert(PRESENT(transm_contact2))
1533 cpassert(.NOT. PRESENT(just_contact))
1534 END IF
1535
1536 ALLOCATE (self_energy_contacts(ncontacts), zwork1_contacts(ncontacts), zwork2_contacts(ncontacts))
1537
1538 IF (PRESENT(just_contact)) THEN
1539 CALL cp_fm_get_info(negf_env%contacts(just_contact)%s_01, matrix_struct=fm_struct)
1540 DO icontact = 1, ncontacts
1541 CALL cp_cfm_create(zwork1_contacts(icontact), fm_struct)
1542 CALL cp_cfm_create(zwork2_contacts(icontact), fm_struct)
1543 END DO
1544
1545 CALL cp_fm_get_info(negf_env%contacts(just_contact)%s_00, nrow_global=nrows, matrix_struct=fm_struct)
1546 DO icontact = 1, ncontacts
1547 CALL cp_cfm_create(self_energy_contacts(icontact), fm_struct)
1548 END DO
1549 ELSE
1550 DO icontact = 1, ncontacts
1551 CALL cp_fm_get_info(negf_env%s_sc(icontact), matrix_struct=fm_struct)
1552 CALL cp_cfm_create(zwork1_contacts(icontact), fm_struct)
1553 CALL cp_cfm_create(zwork2_contacts(icontact), fm_struct)
1554 END DO
1555
1556 CALL cp_fm_get_info(negf_env%s_s, nrow_global=nrows, matrix_struct=fm_struct)
1557 DO icontact = 1, ncontacts
1558 CALL cp_cfm_create(self_energy_contacts(icontact), fm_struct)
1559 END DO
1560 END IF
1561
1562 IF (PRESENT(g_ret_s) .OR. PRESENT(gret_gamma_gadv) .OR. &
1563 PRESENT(dos) .OR. PRESENT(transm_coeff)) THEN
1564 ALLOCATE (g_ret_s_group(npoints))
1565
1566 IF (sub_env%ngroups <= 1 .AND. PRESENT(g_ret_s)) THEN
1567 g_ret_s_group(1:npoints) = g_ret_s(1:npoints)
1568 END IF
1569 END IF
1570
1571 IF (PRESENT(gamma_contacts) .OR. PRESENT(gret_gamma_gadv) .OR. PRESENT(transm_coeff)) THEN
1572 IF (debug_this_module .AND. PRESENT(gamma_contacts)) THEN
1573 cpassert(SIZE(gamma_contacts, 1) == ncontacts)
1574 END IF
1575
1576 ALLOCATE (gamma_contacts_group(ncontacts, npoints))
1577 IF (sub_env%ngroups <= 1 .AND. PRESENT(gamma_contacts)) THEN
1578 gamma_contacts_group(1:ncontacts, 1:npoints) = gamma_contacts(1:ncontacts, 1:npoints)
1579 END IF
1580 END IF
1581
1582 IF (PRESENT(gret_gamma_gadv)) THEN
1583 IF (debug_this_module .AND. PRESENT(gret_gamma_gadv)) THEN
1584 cpassert(SIZE(gret_gamma_gadv, 1) == ncontacts)
1585 END IF
1586
1587 ALLOCATE (gret_gamma_gadv_group(ncontacts, npoints))
1588 IF (sub_env%ngroups <= 1) THEN
1589 gret_gamma_gadv_group(1:ncontacts, 1:npoints) = gret_gamma_gadv(1:ncontacts, 1:npoints)
1590 END IF
1591 END IF
1592
1593 igroup = sub_env%group_distribution(sub_env%mepos_global)
1594
1595 DO ipoint = 1, npoints
1596 IF (ASSOCIATED(g_surf_contacts(1, ipoint)%matrix_struct)) THEN
1597 IF (sub_env%ngroups > 1 .OR. .NOT. PRESENT(g_ret_s)) THEN
1598 ! create a group-specific matrix to store retarded Green's function if there are
1599 ! at least two parallel groups; otherwise pointers to group-specific matrices have
1600 ! already been initialised and they point to globally distributed matrices
1601 IF (ALLOCATED(g_ret_s_group)) THEN
1602 CALL cp_cfm_create(g_ret_s_group(ipoint), fm_struct)
1603 END IF
1604 END IF
1605
1606 IF (sub_env%ngroups > 1 .OR. .NOT. PRESENT(gamma_contacts)) THEN
1607 IF (ALLOCATED(gamma_contacts_group)) THEN
1608 DO icontact = 1, ncontacts
1609 CALL cp_cfm_create(gamma_contacts_group(icontact, ipoint), fm_struct)
1610 END DO
1611 END IF
1612 END IF
1613
1614 IF (sub_env%ngroups > 1) THEN
1615 IF (ALLOCATED(gret_gamma_gadv_group)) THEN
1616 DO icontact = 1, ncontacts
1617 IF (ASSOCIATED(gret_gamma_gadv(icontact, ipoint)%matrix_struct)) THEN
1618 CALL cp_cfm_create(gret_gamma_gadv_group(icontact, ipoint), fm_struct)
1619 END IF
1620 END DO
1621 END IF
1622 END IF
1623
1624 IF (PRESENT(just_contact)) THEN
1625 ! self energy of the "left" (1) and "right" contacts
1626 DO icontact = 1, ncontacts
1627 CALL negf_contact_self_energy(self_energy_c=self_energy_contacts(icontact), &
1628 omega=omega(ipoint), &
1629 g_surf_c=g_surf_contacts(icontact, ipoint), &
1630 h_sc0=negf_env%contacts(just_contact)%h_01(ispin), &
1631 s_sc0=negf_env%contacts(just_contact)%s_01, &
1632 zwork1=zwork1_contacts(icontact), &
1633 zwork2=zwork2_contacts(icontact), &
1634 transp=(icontact == 1))
1635 END DO
1636 ELSE
1637 ! contact self energies
1638 DO icontact = 1, ncontacts
1639 IF (.NOT. ignore_bias) v_external = negf_control%contacts(icontact)%v_external
1640
1641 CALL negf_contact_self_energy(self_energy_c=self_energy_contacts(icontact), &
1642 omega=omega(ipoint) - v_external, &
1643 g_surf_c=g_surf_contacts(icontact, ipoint), &
1644 h_sc0=negf_env%h_sc(ispin, icontact), &
1645 s_sc0=negf_env%s_sc(icontact), &
1646 zwork1=zwork1_contacts(icontact), &
1647 zwork2=zwork2_contacts(icontact), &
1648 transp=.false.)
1649 END DO
1650 END IF
1651
1652 ! broadening matrices
1653 IF (ALLOCATED(gamma_contacts_group)) THEN
1654 DO icontact = 1, ncontacts
1655 CALL negf_contact_broadening_matrix(gamma_c=gamma_contacts_group(icontact, ipoint), &
1656 self_energy_c=self_energy_contacts(icontact))
1657 END DO
1658 END IF
1659
1660 IF (ALLOCATED(g_ret_s_group)) THEN
1661 ! sum up self energies for all contacts
1662 DO icontact = 2, ncontacts
1663 CALL cp_cfm_scale_and_add(z_one, self_energy_contacts(1), z_one, self_energy_contacts(icontact))
1664 END DO
1665
1666 ! retarded Green's function for the scattering region
1667 IF (PRESENT(just_contact)) THEN
1668 CALL negf_retarded_green_function(g_ret_s=g_ret_s_group(ipoint), &
1669 omega=omega(ipoint) - v_shift, &
1670 self_energy_ret_sum=self_energy_contacts(1), &
1671 h_s=negf_env%contacts(just_contact)%h_00(ispin), &
1672 s_s=negf_env%contacts(just_contact)%s_00)
1673 ELSE IF (ignore_bias) THEN
1674 CALL negf_retarded_green_function(g_ret_s=g_ret_s_group(ipoint), &
1675 omega=omega(ipoint) - v_shift, &
1676 self_energy_ret_sum=self_energy_contacts(1), &
1677 h_s=negf_env%h_s(ispin), &
1678 s_s=negf_env%s_s)
1679 ELSE
1680 CALL negf_retarded_green_function(g_ret_s=g_ret_s_group(ipoint), &
1681 omega=omega(ipoint) - v_shift, &
1682 self_energy_ret_sum=self_energy_contacts(1), &
1683 h_s=negf_env%h_s(ispin), &
1684 s_s=negf_env%s_s, &
1685 v_hartree_s=negf_env%v_hartree_s)
1686 END IF
1687
1688 IF (PRESENT(g_ret_scale)) THEN
1689 IF (g_ret_scale(ipoint) /= z_one) CALL cp_cfm_scale(g_ret_scale(ipoint), g_ret_s_group(ipoint))
1690 END IF
1691 END IF
1692
1693 IF (ALLOCATED(gret_gamma_gadv_group)) THEN
1694 ! we do not need contact self energies any longer, so we can use
1695 ! the array 'self_energy_contacts' as a set of work matrices
1696 DO icontact = 1, ncontacts
1697 IF (ASSOCIATED(gret_gamma_gadv_group(icontact, ipoint)%matrix_struct)) THEN
1698 CALL parallel_gemm('N', 'C', nrows, nrows, nrows, &
1699 z_one, gamma_contacts_group(icontact, ipoint), &
1700 g_ret_s_group(ipoint), &
1701 z_zero, self_energy_contacts(icontact))
1702 CALL parallel_gemm('N', 'N', nrows, nrows, nrows, &
1703 z_one, g_ret_s_group(ipoint), &
1704 self_energy_contacts(icontact), &
1705 z_zero, gret_gamma_gadv_group(icontact, ipoint))
1706 END IF
1707 END DO
1708 END IF
1709 END IF
1710 END DO
1711
1712 ! redistribute locally stored matrices
1713 IF (PRESENT(g_ret_s)) THEN
1714 IF (sub_env%ngroups > 1) THEN
1715 NULLIFY (para_env)
1716 DO ipoint = 1, npoints
1717 IF (ASSOCIATED(g_ret_s(ipoint)%matrix_struct)) THEN
1718 CALL cp_cfm_get_info(g_ret_s(ipoint), para_env=para_env)
1719 EXIT
1720 END IF
1721 END DO
1722
1723 IF (ASSOCIATED(para_env)) THEN
1724 ALLOCATE (info1(npoints))
1725
1726 DO ipoint = 1, npoints
1727 CALL cp_cfm_start_copy_general(g_ret_s_group(ipoint), &
1728 g_ret_s(ipoint), &
1729 para_env, info1(ipoint))
1730 END DO
1731
1732 DO ipoint = 1, npoints
1733 IF (ASSOCIATED(g_ret_s(ipoint)%matrix_struct)) THEN
1734 CALL cp_cfm_finish_copy_general(g_ret_s(ipoint), info1(ipoint))
1735 IF (ASSOCIATED(g_ret_s_group(ipoint)%matrix_struct)) &
1736 CALL cp_cfm_cleanup_copy_general(info1(ipoint))
1737 END IF
1738 END DO
1739
1740 DEALLOCATE (info1)
1741 END IF
1742 END IF
1743 END IF
1744
1745 IF (PRESENT(gamma_contacts)) THEN
1746 IF (sub_env%ngroups > 1) THEN
1747 NULLIFY (para_env)
1748 pnt1: DO ipoint = 1, npoints
1749 DO icontact = 1, ncontacts
1750 IF (ASSOCIATED(gamma_contacts(icontact, ipoint)%matrix_struct)) THEN
1751 CALL cp_cfm_get_info(gamma_contacts(icontact, ipoint), para_env=para_env)
1752 EXIT pnt1
1753 END IF
1754 END DO
1755 END DO pnt1
1756
1757 IF (ASSOCIATED(para_env)) THEN
1758 ALLOCATE (info2(ncontacts, npoints))
1759
1760 DO ipoint = 1, npoints
1761 DO icontact = 1, ncontacts
1762 CALL cp_cfm_start_copy_general(gamma_contacts_group(icontact, ipoint), &
1763 gamma_contacts(icontact, ipoint), &
1764 para_env, info2(icontact, ipoint))
1765 END DO
1766 END DO
1767
1768 DO ipoint = 1, npoints
1769 DO icontact = 1, ncontacts
1770 IF (ASSOCIATED(gamma_contacts(icontact, ipoint)%matrix_struct)) THEN
1771 CALL cp_cfm_finish_copy_general(gamma_contacts(icontact, ipoint), info2(icontact, ipoint))
1772 IF (ASSOCIATED(gamma_contacts_group(icontact, ipoint)%matrix_struct)) THEN
1773 CALL cp_cfm_cleanup_copy_general(info2(icontact, ipoint))
1774 END IF
1775 END IF
1776 END DO
1777 END DO
1778
1779 DEALLOCATE (info2)
1780 END IF
1781 END IF
1782 END IF
1783
1784 IF (PRESENT(gret_gamma_gadv)) THEN
1785 IF (sub_env%ngroups > 1) THEN
1786 NULLIFY (para_env)
1787
1788 pnt2: DO ipoint = 1, npoints
1789 DO icontact = 1, ncontacts
1790 IF (ASSOCIATED(gret_gamma_gadv(icontact, ipoint)%matrix_struct)) THEN
1791 CALL cp_cfm_get_info(gret_gamma_gadv(icontact, ipoint), para_env=para_env)
1792 EXIT pnt2
1793 END IF
1794 END DO
1795 END DO pnt2
1796
1797 IF (ASSOCIATED(para_env)) THEN
1798 ALLOCATE (info2(ncontacts, npoints))
1799
1800 DO ipoint = 1, npoints
1801 DO icontact = 1, ncontacts
1802 CALL cp_cfm_start_copy_general(gret_gamma_gadv_group(icontact, ipoint), &
1803 gret_gamma_gadv(icontact, ipoint), &
1804 para_env, info2(icontact, ipoint))
1805 END DO
1806 END DO
1807
1808 DO ipoint = 1, npoints
1809 DO icontact = 1, ncontacts
1810 IF (ASSOCIATED(gret_gamma_gadv(icontact, ipoint)%matrix_struct)) THEN
1811 CALL cp_cfm_finish_copy_general(gret_gamma_gadv(icontact, ipoint), info2(icontact, ipoint))
1812 IF (ASSOCIATED(gret_gamma_gadv_group(icontact, ipoint)%matrix_struct)) THEN
1813 CALL cp_cfm_cleanup_copy_general(info2(icontact, ipoint))
1814 END IF
1815 END IF
1816 END DO
1817 END DO
1818
1819 DEALLOCATE (info2)
1820 END IF
1821 END IF
1822 END IF
1823
1824 IF (PRESENT(dos)) THEN
1825 dos(:) = 0.0_dp
1826
1827 IF (PRESENT(just_contact)) THEN
1828 matrix_s => negf_env%contacts(just_contact)%s_00
1829 ELSE
1830 matrix_s => negf_env%s_s
1831 END IF
1832
1833 CALL cp_fm_get_info(matrix_s, matrix_struct=fm_struct)
1834 CALL cp_fm_create(g_ret_imag, fm_struct)
1835
1836 DO ipoint = 1, npoints
1837 IF (ASSOCIATED(g_ret_s_group(ipoint)%matrix_struct)) THEN
1838 CALL cp_cfm_to_fm(g_ret_s_group(ipoint), mtargeti=g_ret_imag)
1839 CALL cp_fm_trace(g_ret_imag, matrix_s, dos(ipoint))
1840 IF (sub_env%para_env%mepos /= 0) dos(ipoint) = 0.0_dp
1841 END IF
1842 END DO
1843
1844 CALL cp_fm_release(g_ret_imag)
1845
1846 CALL sub_env%mpi_comm_global%sum(dos)
1847 dos(:) = -1.0_dp/pi*dos(:)
1848 END IF
1849
1850 IF (PRESENT(transm_coeff)) THEN
1851 transm_coeff(:) = z_zero
1852
1853 DO ipoint = 1, npoints
1854 IF (ASSOCIATED(g_ret_s_group(ipoint)%matrix_struct)) THEN
1855 ! gamma_1 * g_adv_s * gamma_2
1856 CALL parallel_gemm('N', 'C', nrows, nrows, nrows, &
1857 z_one, gamma_contacts_group(transm_contact1, ipoint), &
1858 g_ret_s_group(ipoint), &
1859 z_zero, self_energy_contacts(transm_contact1))
1860 CALL parallel_gemm('N', 'N', nrows, nrows, nrows, &
1861 z_one, self_energy_contacts(transm_contact1), &
1862 gamma_contacts_group(transm_contact2, ipoint), &
1863 z_zero, self_energy_contacts(transm_contact2))
1864
1865 ! Trace[ g_ret_s * gamma_1 * g_adv_s * gamma_2 ]
1866 CALL cp_cfm_trace(g_ret_s_group(ipoint), &
1867 self_energy_contacts(transm_contact2), &
1868 transm_coeff(ipoint))
1869 IF (sub_env%para_env%mepos /= 0) transm_coeff(ipoint) = 0.0_dp
1870 END IF
1871 END DO
1872
1873 ! transmission coefficients are scaled by 2/pi
1874 CALL sub_env%mpi_comm_global%sum(transm_coeff)
1875 !transm_coeff(:) = 0.5_dp/pi*transm_coeff(:)
1876 END IF
1877
1878 ! -- deallocate temporary matrices
1879 IF (ALLOCATED(g_ret_s_group)) THEN
1880 DO ipoint = npoints, 1, -1
1881 IF (sub_env%ngroups > 1 .OR. .NOT. PRESENT(g_ret_s)) THEN
1882 CALL cp_cfm_release(g_ret_s_group(ipoint))
1883 END IF
1884 END DO
1885 DEALLOCATE (g_ret_s_group)
1886 END IF
1887
1888 IF (ALLOCATED(gamma_contacts_group)) THEN
1889 DO ipoint = npoints, 1, -1
1890 DO icontact = ncontacts, 1, -1
1891 IF (sub_env%ngroups > 1 .OR. .NOT. PRESENT(gamma_contacts)) THEN
1892 CALL cp_cfm_release(gamma_contacts_group(icontact, ipoint))
1893 END IF
1894 END DO
1895 END DO
1896 DEALLOCATE (gamma_contacts_group)
1897 END IF
1898
1899 IF (ALLOCATED(gret_gamma_gadv_group)) THEN
1900 DO ipoint = npoints, 1, -1
1901 DO icontact = ncontacts, 1, -1
1902 IF (sub_env%ngroups > 1) THEN
1903 CALL cp_cfm_release(gret_gamma_gadv_group(icontact, ipoint))
1904 END IF
1905 END DO
1906 END DO
1907 DEALLOCATE (gret_gamma_gadv_group)
1908 END IF
1909
1910 IF (ALLOCATED(self_energy_contacts)) THEN
1911 DO icontact = ncontacts, 1, -1
1912 CALL cp_cfm_release(self_energy_contacts(icontact))
1913 END DO
1914 DEALLOCATE (self_energy_contacts)
1915 END IF
1916
1917 IF (ALLOCATED(zwork1_contacts)) THEN
1918 DO icontact = ncontacts, 1, -1
1919 CALL cp_cfm_release(zwork1_contacts(icontact))
1920 END DO
1921 DEALLOCATE (zwork1_contacts)
1922 END IF
1923
1924 IF (ALLOCATED(zwork2_contacts)) THEN
1925 DO icontact = ncontacts, 1, -1
1926 CALL cp_cfm_release(zwork2_contacts(icontact))
1927 END DO
1928 DEALLOCATE (zwork2_contacts)
1929 END IF
1930
1931 CALL timestop(handle)
1932 END SUBROUTINE negf_retarded_green_function_batch
1933
1934! **************************************************************************************************
1935!> \brief Fermi function (exp(E/(kT)) + 1) ^ {-1} .
1936!> \param omega 'energy' point on the complex plane
1937!> \param temperature temperature in atomic units
1938!> \return value
1939!> \par History
1940!> * 05.2017 created [Sergey Chulkov]
1941! **************************************************************************************************
1942 PURE FUNCTION fermi_function(omega, temperature) RESULT(val)
1943 COMPLEX(kind=dp), INTENT(in) :: omega
1944 REAL(kind=dp), INTENT(in) :: temperature
1945 COMPLEX(kind=dp) :: val
1946
1947 REAL(kind=dp), PARAMETER :: max_ln_omega_over_t = log(huge(0.0_dp))/16.0_dp
1948
1949 IF (real(omega, kind=dp) <= temperature*max_ln_omega_over_t) THEN
1950 ! exp(omega / T) < huge(0), so EXP() should not return infinity
1951 val = z_one/(exp(omega/temperature) + z_one)
1952 ELSE
1953 val = z_zero
1954 END IF
1955 END FUNCTION fermi_function
1956
1957! **************************************************************************************************
1958!> \brief Compute contribution to the density matrix from the poles of the Fermi function.
1959!> \param rho_ao_fm density matrix (initialised on exit)
1960!> \param v_shift shift in Hartree potential
1961!> \param ignore_bias ignore v_external from negf_control
1962!> \param negf_env NEGF environment
1963!> \param negf_control NEGF control
1964!> \param sub_env NEGF parallel (sub)group environment
1965!> \param ispin spin conponent to proceed
1966!> \param base_contact index of the reference contact
1967!> \param just_contact ...
1968!> \author Sergey Chulkov
1969! **************************************************************************************************
1970 SUBROUTINE negf_init_rho_equiv_residuals(rho_ao_fm, v_shift, ignore_bias, negf_env, &
1971 negf_control, sub_env, ispin, base_contact, just_contact)
1972 TYPE(cp_fm_type), INTENT(IN) :: rho_ao_fm
1973 REAL(kind=dp), INTENT(in) :: v_shift
1974 LOGICAL, INTENT(in) :: ignore_bias
1975 TYPE(negf_env_type), INTENT(in) :: negf_env
1976 TYPE(negf_control_type), POINTER :: negf_control
1977 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
1978 INTEGER, INTENT(in) :: ispin, base_contact
1979 INTEGER, INTENT(in), OPTIONAL :: just_contact
1980
1981 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_init_rho_equiv_residuals'
1982
1983 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: omega
1984 INTEGER :: handle, icontact, ipole, ncontacts, &
1985 npoles
1986 REAL(kind=dp) :: mu_base, pi_temperature, temperature, &
1987 v_external
1988 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: g_ret_s
1989 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1990 TYPE(green_functions_cache_type) :: g_surf_cache
1991 TYPE(mp_para_env_type), POINTER :: para_env
1992
1993 CALL timeset(routinen, handle)
1994
1995 temperature = negf_control%contacts(base_contact)%temperature
1996 IF (ignore_bias) THEN
1997 mu_base = negf_control%contacts(base_contact)%fermi_level
1998 v_external = 0.0_dp
1999 ELSE
2000 mu_base = negf_control%contacts(base_contact)%fermi_level + negf_control%contacts(base_contact)%v_external
2001 END IF
2002
2003 pi_temperature = pi*temperature
2004 npoles = negf_control%delta_npoles
2005
2006 ncontacts = SIZE(negf_env%contacts)
2007 cpassert(base_contact <= ncontacts)
2008 IF (PRESENT(just_contact)) THEN
2009 ncontacts = 2
2010 cpassert(just_contact == base_contact)
2011 END IF
2012
2013 IF (npoles > 0) THEN
2014 CALL cp_fm_get_info(rho_ao_fm, para_env=para_env, matrix_struct=fm_struct)
2015
2016 ALLOCATE (omega(npoles), g_ret_s(npoles))
2017
2018 DO ipole = 1, npoles
2019 CALL cp_cfm_create(g_ret_s(ipole), fm_struct)
2020
2021 omega(ipole) = cmplx(mu_base, real(2*ipole - 1, kind=dp)*pi_temperature, kind=dp)
2022 END DO
2023
2024 CALL green_functions_cache_expand(g_surf_cache, ncontacts, npoles)
2025
2026 IF (PRESENT(just_contact)) THEN
2027 ! do not apply the external potential when computing the Fermi level of a bulk contact.
2028 ! We are using a fictitious electronic device, which identical to the bulk contact in question;
2029 ! icontact == 1 corresponds to the "left" contact, so the matrices h_01 and s_01 needs to be transposed,
2030 ! while icontact == 2 correspond to the "right" contact and we should use the matrices h_01 and s_01 as is.
2031 DO icontact = 1, ncontacts
2032 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, :), &
2033 omega=omega(:), &
2034 h0=negf_env%contacts(just_contact)%h_00(ispin), &
2035 s0=negf_env%contacts(just_contact)%s_00, &
2036 h1=negf_env%contacts(just_contact)%h_01(ispin), &
2037 s1=negf_env%contacts(just_contact)%s_01, &
2038 sub_env=sub_env, v_external=0.0_dp, &
2039 conv=negf_control%conv_green, transp=(icontact == 1))
2040 END DO
2041 ELSE
2042 DO icontact = 1, ncontacts
2043 IF (.NOT. ignore_bias) v_external = negf_control%contacts(icontact)%v_external
2044
2045 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, :), &
2046 omega=omega(:), &
2047 h0=negf_env%contacts(icontact)%h_00(ispin), &
2048 s0=negf_env%contacts(icontact)%s_00, &
2049 h1=negf_env%contacts(icontact)%h_01(ispin), &
2050 s1=negf_env%contacts(icontact)%s_01, &
2051 sub_env=sub_env, &
2052 v_external=v_external, &
2053 conv=negf_control%conv_green, transp=.false.)
2054 END DO
2055 END IF
2056
2057 CALL negf_retarded_green_function_batch(omega=omega(:), &
2058 v_shift=v_shift, &
2059 ignore_bias=ignore_bias, &
2060 negf_env=negf_env, &
2061 negf_control=negf_control, &
2062 sub_env=sub_env, &
2063 ispin=ispin, &
2064 g_surf_contacts=g_surf_cache%g_surf_contacts, &
2065 g_ret_s=g_ret_s, &
2066 just_contact=just_contact)
2067
2068 CALL green_functions_cache_release(g_surf_cache)
2069
2070 DO ipole = 2, npoles
2071 CALL cp_cfm_scale_and_add(z_one, g_ret_s(1), z_one, g_ret_s(ipole))
2072 END DO
2073
2074 !Re(-i * (-2*pi*i*kB*T/(-pi) * [Re(G)+i*Im(G)]) == 2*kB*T * Re(G)
2075 CALL cp_cfm_to_fm(g_ret_s(1), mtargetr=rho_ao_fm)
2076 CALL cp_fm_scale(2.0_dp*temperature, rho_ao_fm)
2077
2078 DO ipole = npoles, 1, -1
2079 CALL cp_cfm_release(g_ret_s(ipole))
2080 END DO
2081 DEALLOCATE (g_ret_s, omega)
2082 END IF
2083
2084 CALL timestop(handle)
2085 END SUBROUTINE negf_init_rho_equiv_residuals
2086
2087! **************************************************************************************************
2088!> \brief Compute equilibrium contribution to the density matrix.
2089!> \param rho_ao_fm density matrix (initialised on exit)
2090!> \param stats integration statistics (updated on exit)
2091!> \param v_shift shift in Hartree potential
2092!> \param ignore_bias ignore v_external from negf_control
2093!> \param negf_env NEGF environment
2094!> \param negf_control NEGF control
2095!> \param sub_env NEGF parallel (sub)group environment
2096!> \param ispin spin conponent to proceed
2097!> \param base_contact index of the reference contact
2098!> \param integr_lbound integration lower bound
2099!> \param integr_ubound integration upper bound
2100!> \param matrix_s_global globally distributed overlap matrix
2101!> \param is_circular compute the integral along the circular path
2102!> \param g_surf_cache set of precomputed surface Green's functions (updated on exit)
2103!> \param just_contact ...
2104!> \author Sergey Chulkov
2105! **************************************************************************************************
2106 SUBROUTINE negf_add_rho_equiv_low(rho_ao_fm, stats, v_shift, ignore_bias, negf_env, negf_control, sub_env, &
2107 ispin, base_contact, integr_lbound, integr_ubound, matrix_s_global, &
2108 is_circular, g_surf_cache, just_contact)
2109 TYPE(cp_fm_type), INTENT(IN) :: rho_ao_fm
2110 TYPE(integration_status_type), INTENT(inout) :: stats
2111 REAL(kind=dp), INTENT(in) :: v_shift
2112 LOGICAL, INTENT(in) :: ignore_bias
2113 TYPE(negf_env_type), INTENT(in) :: negf_env
2114 TYPE(negf_control_type), POINTER :: negf_control
2115 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
2116 INTEGER, INTENT(in) :: ispin, base_contact
2117 COMPLEX(kind=dp), INTENT(in) :: integr_lbound, integr_ubound
2118 TYPE(cp_fm_type), INTENT(IN) :: matrix_s_global
2119 LOGICAL, INTENT(in) :: is_circular
2120 TYPE(green_functions_cache_type), INTENT(inout) :: g_surf_cache
2121 INTEGER, INTENT(in), OPTIONAL :: just_contact
2122
2123 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_add_rho_equiv_low'
2124
2125 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: xnodes, zscale
2126 INTEGER :: handle, icontact, interval_id, ipoint, max_points, min_points, ncontacts, &
2127 npoints, npoints_exist, npoints_tmp, npoints_total, shape_id
2128 LOGICAL :: do_surface_green
2129 REAL(kind=dp) :: conv_integr, mu_base, temperature, &
2130 v_external
2131 TYPE(ccquad_type) :: cc_env
2132 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: zdata, zdata_tmp
2133 TYPE(cp_fm_struct_type), POINTER :: fm_struct
2134 TYPE(cp_fm_type) :: integral_imag
2135 TYPE(mp_para_env_type), POINTER :: para_env
2136 TYPE(simpsonrule_type) :: sr_env
2137
2138 CALL timeset(routinen, handle)
2139
2140 ! convergence criteria for the integral of the retarded Green's function. This integral needs to be
2141 ! computed for both spin-components and needs to be scaled by -1/pi to obtain the electron density.
2142 conv_integr = 0.5_dp*negf_control%conv_density*pi
2143
2144 IF (ignore_bias) THEN
2145 mu_base = negf_control%contacts(base_contact)%fermi_level
2146 v_external = 0.0_dp
2147 ELSE
2148 mu_base = negf_control%contacts(base_contact)%fermi_level + negf_control%contacts(base_contact)%v_external
2149 END IF
2150
2151 min_points = negf_control%integr_min_points
2152 max_points = negf_control%integr_max_points
2153 temperature = negf_control%contacts(base_contact)%temperature
2154
2155 ncontacts = SIZE(negf_env%contacts)
2156 cpassert(base_contact <= ncontacts)
2157 IF (PRESENT(just_contact)) THEN
2158 ncontacts = 2
2159 cpassert(just_contact == base_contact)
2160 END IF
2161
2162 do_surface_green = .NOT. ALLOCATED(g_surf_cache%tnodes)
2163
2164 IF (do_surface_green) THEN
2165 npoints = min_points
2166 ELSE
2167 npoints = SIZE(g_surf_cache%tnodes)
2168 END IF
2169 npoints_total = 0
2170
2171 CALL cp_fm_get_info(rho_ao_fm, para_env=para_env, matrix_struct=fm_struct)
2172 CALL cp_fm_create(integral_imag, fm_struct)
2173
2174 SELECT CASE (negf_control%integr_method)
2175 CASE (negfint_method_cc)
2176 ! Adaptive Clenshaw-Curtis method
2177 ALLOCATE (xnodes(npoints))
2178
2179 IF (is_circular) THEN
2180 shape_id = cc_shape_arc
2181 interval_id = cc_interval_full
2182 ELSE
2183 shape_id = cc_shape_linear
2184 interval_id = cc_interval_half
2185 END IF
2186
2187 IF (do_surface_green) THEN
2188 CALL ccquad_init(cc_env, xnodes, npoints, integr_lbound, integr_ubound, &
2189 interval_id, shape_id, matrix_s_global)
2190 ELSE
2191 CALL ccquad_init(cc_env, xnodes, npoints, integr_lbound, integr_ubound, &
2192 interval_id, shape_id, matrix_s_global, tnodes_restart=g_surf_cache%tnodes)
2193 END IF
2194
2195 ALLOCATE (zdata(npoints))
2196 DO ipoint = 1, npoints
2197 CALL cp_cfm_create(zdata(ipoint), fm_struct)
2198 END DO
2199
2200 DO
2201 IF (do_surface_green) THEN
2202 CALL green_functions_cache_expand(g_surf_cache, ncontacts, npoints)
2203
2204 IF (PRESENT(just_contact)) THEN
2205 ! do not apply the external potential when computing the Fermi level of a bulk contact.
2206 DO icontact = 1, ncontacts
2207 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, npoints_total + 1:), &
2208 omega=xnodes(1:npoints), &
2209 h0=negf_env%contacts(just_contact)%h_00(ispin), &
2210 s0=negf_env%contacts(just_contact)%s_00, &
2211 h1=negf_env%contacts(just_contact)%h_01(ispin), &
2212 s1=negf_env%contacts(just_contact)%s_01, &
2213 sub_env=sub_env, v_external=0.0_dp, &
2214 conv=negf_control%conv_green, transp=(icontact == 1))
2215 END DO
2216 ELSE
2217 DO icontact = 1, ncontacts
2218 IF (.NOT. ignore_bias) v_external = negf_control%contacts(icontact)%v_external
2219
2220 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, npoints_total + 1:), &
2221 omega=xnodes(1:npoints), &
2222 h0=negf_env%contacts(icontact)%h_00(ispin), &
2223 s0=negf_env%contacts(icontact)%s_00, &
2224 h1=negf_env%contacts(icontact)%h_01(ispin), &
2225 s1=negf_env%contacts(icontact)%s_01, &
2226 sub_env=sub_env, &
2227 v_external=v_external, &
2228 conv=negf_control%conv_green, transp=.false.)
2229 END DO
2230 END IF
2231 END IF
2232
2233 ALLOCATE (zscale(npoints))
2234
2235 IF (temperature >= 0.0_dp) THEN
2236 DO ipoint = 1, npoints
2237 zscale(ipoint) = fermi_function(xnodes(ipoint) - mu_base, temperature)
2238 END DO
2239 ELSE
2240 zscale(:) = z_one
2241 END IF
2242
2243 CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints), &
2244 v_shift=v_shift, &
2245 ignore_bias=ignore_bias, &
2246 negf_env=negf_env, &
2247 negf_control=negf_control, &
2248 sub_env=sub_env, &
2249 ispin=ispin, &
2250 g_surf_contacts=g_surf_cache%g_surf_contacts(:, npoints_total + 1:), &
2251 g_ret_s=zdata(1:npoints), &
2252 g_ret_scale=zscale(1:npoints), &
2253 just_contact=just_contact)
2254
2255 DEALLOCATE (xnodes, zscale)
2256 npoints_total = npoints_total + npoints
2257
2258 CALL ccquad_reduce_and_append_zdata(cc_env, zdata)
2259 CALL move_alloc(zdata, zdata_tmp)
2260
2261 CALL ccquad_refine_integral(cc_env)
2262
2263 IF (cc_env%error <= conv_integr) EXIT
2264 IF (2*(npoints_total - 1) + 1 > max_points) EXIT
2265
2266 ! all cached points have been reused at the first iteration;
2267 ! we need to compute surface Green's function at extra points if the integral has not been converged
2268 do_surface_green = .true.
2269
2270 npoints_tmp = npoints
2271 CALL ccquad_double_number_of_points(cc_env, xnodes)
2272 npoints = SIZE(xnodes)
2273
2274 ALLOCATE (zdata(npoints))
2275
2276 npoints_exist = 0
2277 DO ipoint = 1, npoints_tmp
2278 IF (ASSOCIATED(zdata_tmp(ipoint)%matrix_struct)) THEN
2279 npoints_exist = npoints_exist + 1
2280 zdata(npoints_exist) = zdata_tmp(ipoint)
2281 END IF
2282 END DO
2283 DEALLOCATE (zdata_tmp)
2284
2285 DO ipoint = npoints_exist + 1, npoints
2286 CALL cp_cfm_create(zdata(ipoint), fm_struct)
2287 END DO
2288 END DO
2289
2290 ! the obtained integral will be scaled by -1/pi, so scale the error extimate as well
2291 stats%error = stats%error + cc_env%error/pi
2292
2293 DO ipoint = SIZE(zdata_tmp), 1, -1
2294 CALL cp_cfm_release(zdata_tmp(ipoint))
2295 END DO
2296 DEALLOCATE (zdata_tmp)
2297
2298 CALL cp_cfm_to_fm(cc_env%integral, mtargeti=integral_imag)
2299
2300 ! keep the cache
2301 IF (do_surface_green) THEN
2302 CALL green_functions_cache_reorder(g_surf_cache, cc_env%tnodes)
2303 END IF
2304 CALL ccquad_release(cc_env)
2305
2306 CASE (negfint_method_simpson)
2307 ! Adaptive Simpson's rule method
2308 ALLOCATE (xnodes(npoints), zdata(npoints), zscale(npoints))
2309
2310 IF (is_circular) THEN
2311 shape_id = sr_shape_arc
2312 ELSE
2313 shape_id = sr_shape_linear
2314 END IF
2315
2316 IF (do_surface_green) THEN
2317 CALL simpsonrule_init(sr_env, xnodes, npoints, integr_lbound, integr_ubound, &
2318 shape_id, conv_integr, matrix_s_global)
2319 ELSE
2320 CALL simpsonrule_init(sr_env, xnodes, npoints, integr_lbound, integr_ubound, &
2321 shape_id, conv_integr, matrix_s_global, tnodes_restart=g_surf_cache%tnodes)
2322 END IF
2323
2324 DO WHILE (npoints > 0 .AND. npoints_total < max_points)
2325 DO ipoint = 1, npoints
2326 CALL cp_cfm_create(zdata(ipoint), fm_struct)
2327 END DO
2328
2329 IF (do_surface_green) THEN
2330 CALL green_functions_cache_expand(g_surf_cache, ncontacts, npoints)
2331
2332 IF (PRESENT(just_contact)) THEN
2333 ! do not apply the external potential when computing the Fermi level of a bulk contact.
2334 DO icontact = 1, ncontacts
2335 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, npoints_total + 1:), &
2336 omega=xnodes(1:npoints), &
2337 h0=negf_env%contacts(just_contact)%h_00(ispin), &
2338 s0=negf_env%contacts(just_contact)%s_00, &
2339 h1=negf_env%contacts(just_contact)%h_01(ispin), &
2340 s1=negf_env%contacts(just_contact)%s_01, &
2341 sub_env=sub_env, v_external=0.0_dp, &
2342 conv=negf_control%conv_green, transp=(icontact == 1))
2343 END DO
2344 ELSE
2345 DO icontact = 1, ncontacts
2346 IF (.NOT. ignore_bias) v_external = negf_control%contacts(icontact)%v_external
2347
2348 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, npoints_total + 1:), &
2349 omega=xnodes(1:npoints), &
2350 h0=negf_env%contacts(icontact)%h_00(ispin), &
2351 s0=negf_env%contacts(icontact)%s_00, &
2352 h1=negf_env%contacts(icontact)%h_01(ispin), &
2353 s1=negf_env%contacts(icontact)%s_01, &
2354 sub_env=sub_env, &
2355 v_external=v_external, &
2356 conv=negf_control%conv_green, transp=.false.)
2357 END DO
2358 END IF
2359 END IF
2360
2361 IF (temperature >= 0.0_dp) THEN
2362 DO ipoint = 1, npoints
2363 zscale(ipoint) = fermi_function(xnodes(ipoint) - mu_base, temperature)
2364 END DO
2365 ELSE
2366 zscale(:) = z_one
2367 END IF
2368
2369 CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints), &
2370 v_shift=v_shift, &
2371 ignore_bias=ignore_bias, &
2372 negf_env=negf_env, &
2373 negf_control=negf_control, &
2374 sub_env=sub_env, &
2375 ispin=ispin, &
2376 g_surf_contacts=g_surf_cache%g_surf_contacts(:, npoints_total + 1:), &
2377 g_ret_s=zdata(1:npoints), &
2378 g_ret_scale=zscale(1:npoints), &
2379 just_contact=just_contact)
2380
2381 npoints_total = npoints_total + npoints
2382
2383 CALL simpsonrule_refine_integral(sr_env, zdata(1:npoints))
2384
2385 IF (sr_env%error <= conv_integr) EXIT
2386
2387 ! all cached points have been reused at the first iteration;
2388 ! if the integral has not been converged, turn on the 'do_surface_green' flag
2389 ! in order to add more points
2390 do_surface_green = .true.
2391
2392 npoints = max_points - npoints_total
2393 IF (npoints <= 0) EXIT
2394 IF (npoints > SIZE(xnodes)) npoints = SIZE(xnodes)
2395
2396 CALL simpsonrule_get_next_nodes(sr_env, xnodes, npoints)
2397 END DO
2398
2399 ! the obtained integral will be scaled by -1/pi, so scale the error extimate as well
2400 stats%error = stats%error + sr_env%error/pi
2401
2402 CALL cp_cfm_to_fm(sr_env%integral, mtargeti=integral_imag)
2403
2404 ! keep the cache
2405 IF (do_surface_green) THEN
2406 CALL green_functions_cache_reorder(g_surf_cache, sr_env%tnodes)
2407 END IF
2408
2409 CALL simpsonrule_release(sr_env)
2410 DEALLOCATE (xnodes, zdata, zscale)
2411
2412 CASE DEFAULT
2413 cpabort("Unimplemented integration method")
2414 END SELECT
2415
2416 stats%npoints = stats%npoints + npoints_total
2417
2418 CALL cp_fm_scale_and_add(1.0_dp, rho_ao_fm, -1.0_dp/pi, integral_imag)
2419 CALL cp_fm_release(integral_imag)
2420
2421 CALL timestop(handle)
2422 END SUBROUTINE negf_add_rho_equiv_low
2423
2424! **************************************************************************************************
2425!> \brief Compute non-equilibrium contribution to the density matrix.
2426!> \param rho_ao_fm density matrix (initialised on exit)
2427!> \param stats integration statistics (updated on exit)
2428!> \param v_shift shift in Hartree potential
2429!> \param negf_env NEGF environment
2430!> \param negf_control NEGF control
2431!> \param sub_env NEGF parallel (sub)group environment
2432!> \param ispin spin conponent to proceed
2433!> \param base_contact index of the reference contact
2434!> \param matrix_s_global globally distributed overlap matrix
2435!> \param g_surf_cache set of precomputed surface Green's functions (updated on exit)
2436!> \author Sergey Chulkov
2437! **************************************************************************************************
2438 SUBROUTINE negf_add_rho_nonequiv(rho_ao_fm, stats, v_shift, negf_env, negf_control, sub_env, &
2439 ispin, base_contact, matrix_s_global, g_surf_cache)
2440 TYPE(cp_fm_type), INTENT(IN) :: rho_ao_fm
2441 TYPE(integration_status_type), INTENT(inout) :: stats
2442 REAL(kind=dp), INTENT(in) :: v_shift
2443 TYPE(negf_env_type), INTENT(in) :: negf_env
2444 TYPE(negf_control_type), POINTER :: negf_control
2445 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
2446 INTEGER, INTENT(in) :: ispin, base_contact
2447 TYPE(cp_fm_type), INTENT(IN) :: matrix_s_global
2448 TYPE(green_functions_cache_type), INTENT(inout) :: g_surf_cache
2449
2450 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_add_rho_nonequiv'
2451
2452 COMPLEX(kind=dp) :: fermi_base, fermi_contact, &
2453 integr_lbound, integr_ubound
2454 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: xnodes
2455 INTEGER :: handle, icontact, ipoint, jcontact, &
2456 max_points, min_points, ncontacts, &
2457 npoints, npoints_total
2458 LOGICAL :: do_surface_green
2459 REAL(kind=dp) :: conv_density, conv_integr, eta, &
2460 ln_conv_density, mu_base, mu_contact, &
2461 temperature_base, temperature_contact
2462 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:, :) :: zdata
2463 TYPE(cp_fm_struct_type), POINTER :: fm_struct
2464 TYPE(cp_fm_type) :: integral_real
2465 TYPE(simpsonrule_type) :: sr_env
2466
2467 CALL timeset(routinen, handle)
2468
2469 ncontacts = SIZE(negf_env%contacts)
2470 cpassert(base_contact <= ncontacts)
2471
2472 ! the current subroutine works for the general case as well, but the Poisson solver does not
2473 IF (ncontacts > 2) THEN
2474 cpabort("Poisson solver does not support the general NEGF setup (>2 contacts).")
2475 END IF
2476
2477 mu_base = negf_control%contacts(base_contact)%fermi_level + negf_control%contacts(base_contact)%v_external
2478 min_points = negf_control%integr_min_points
2479 max_points = negf_control%integr_max_points
2480 temperature_base = negf_control%contacts(base_contact)%temperature
2481 eta = negf_control%eta
2482 conv_density = negf_control%conv_density
2483 ln_conv_density = log(conv_density)
2484
2485 ! convergence criteria for the integral. This integral needs to be computed for both
2486 ! spin-components and needs to be scaled by -1/pi to obtain the electron density.
2487 conv_integr = 0.5_dp*conv_density*pi
2488
2489 DO icontact = 1, ncontacts
2490 IF (icontact /= base_contact) THEN
2491 mu_contact = negf_control%contacts(icontact)%fermi_level + negf_control%contacts(icontact)%v_external
2492 temperature_contact = negf_control%contacts(icontact)%temperature
2493
2494 integr_lbound = cmplx(min(mu_base + ln_conv_density*temperature_base, &
2495 mu_contact + ln_conv_density*temperature_contact), eta, kind=dp)
2496 integr_ubound = cmplx(max(mu_base - ln_conv_density*temperature_base, &
2497 mu_contact - ln_conv_density*temperature_contact), eta, kind=dp)
2498
2499 do_surface_green = .NOT. ALLOCATED(g_surf_cache%tnodes)
2500
2501 IF (do_surface_green) THEN
2502 npoints = min_points
2503 ELSE
2504 npoints = SIZE(g_surf_cache%tnodes)
2505 END IF
2506 npoints_total = 0
2507
2508 CALL cp_fm_get_info(rho_ao_fm, matrix_struct=fm_struct)
2509
2510 ALLOCATE (xnodes(npoints), zdata(ncontacts, npoints))
2511
2512 IF (do_surface_green) THEN
2513 CALL simpsonrule_init(sr_env, xnodes, npoints, integr_lbound, integr_ubound, &
2514 sr_shape_linear, conv_integr, matrix_s_global)
2515 ELSE
2516 CALL simpsonrule_init(sr_env, xnodes, npoints, integr_lbound, integr_ubound, &
2517 sr_shape_linear, conv_integr, matrix_s_global, tnodes_restart=g_surf_cache%tnodes)
2518 END IF
2519
2520 DO WHILE (npoints > 0 .AND. npoints_total < max_points)
2521 IF (do_surface_green) THEN
2522 CALL green_functions_cache_expand(g_surf_cache, ncontacts, npoints)
2523
2524 DO jcontact = 1, ncontacts
2525 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(jcontact, npoints_total + 1:), &
2526 omega=xnodes(1:npoints), &
2527 h0=negf_env%contacts(jcontact)%h_00(ispin), &
2528 s0=negf_env%contacts(jcontact)%s_00, &
2529 h1=negf_env%contacts(jcontact)%h_01(ispin), &
2530 s1=negf_env%contacts(jcontact)%s_01, &
2531 sub_env=sub_env, &
2532 v_external=negf_control%contacts(jcontact)%v_external, &
2533 conv=negf_control%conv_green, transp=.false.)
2534 END DO
2535 END IF
2536
2537 DO ipoint = 1, npoints
2538 CALL cp_cfm_create(zdata(icontact, ipoint), fm_struct)
2539 END DO
2540
2541 CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints), &
2542 v_shift=v_shift, &
2543 ignore_bias=.false., &
2544 negf_env=negf_env, &
2545 negf_control=negf_control, &
2546 sub_env=sub_env, &
2547 ispin=ispin, &
2548 g_surf_contacts=g_surf_cache%g_surf_contacts(:, npoints_total + 1:), &
2549 gret_gamma_gadv=zdata(:, 1:npoints))
2550
2551 DO ipoint = 1, npoints
2552 fermi_base = fermi_function(cmplx(real(xnodes(ipoint), kind=dp) - mu_base, 0.0_dp, kind=dp), &
2553 temperature_base)
2554 fermi_contact = fermi_function(cmplx(real(xnodes(ipoint), kind=dp) - mu_contact, 0.0_dp, kind=dp), &
2555 temperature_contact)
2556 CALL cp_cfm_scale(fermi_contact - fermi_base, zdata(icontact, ipoint))
2557 END DO
2558
2559 npoints_total = npoints_total + npoints
2560
2561 CALL simpsonrule_refine_integral(sr_env, zdata(icontact, 1:npoints))
2562
2563 IF (sr_env%error <= conv_integr) EXIT
2564
2565 ! not enought cached points to achieve target accuracy
2566 do_surface_green = .true.
2567
2568 npoints = max_points - npoints_total
2569 IF (npoints <= 0) EXIT
2570 IF (npoints > SIZE(xnodes)) npoints = SIZE(xnodes)
2571
2572 CALL simpsonrule_get_next_nodes(sr_env, xnodes, npoints)
2573 END DO
2574
2575 CALL cp_fm_create(integral_real, fm_struct)
2576
2577 CALL cp_cfm_to_fm(sr_env%integral, mtargetr=integral_real)
2578 CALL cp_fm_scale_and_add(1.0_dp, rho_ao_fm, 0.5_dp/pi, integral_real)
2579
2580 CALL cp_fm_release(integral_real)
2581
2582 DEALLOCATE (xnodes, zdata)
2583
2584 stats%error = stats%error + sr_env%error*0.5_dp/pi
2585 stats%npoints = stats%npoints + npoints_total
2586
2587 ! keep the cache
2588 IF (do_surface_green) THEN
2589 CALL green_functions_cache_reorder(g_surf_cache, sr_env%tnodes)
2590 END IF
2591
2592 CALL simpsonrule_release(sr_env)
2593 END IF
2594 END DO
2595
2596 CALL timestop(handle)
2597 END SUBROUTINE negf_add_rho_nonequiv
2598
2599! **************************************************************************************************
2600!> \brief Reset integration statistics.
2601!> \param stats integration statistics
2602!> \author Sergey Chulkov
2603! **************************************************************************************************
2604 ELEMENTAL SUBROUTINE integration_status_reset(stats)
2605 TYPE(integration_status_type), INTENT(out) :: stats
2606
2607 stats%npoints = 0
2608 stats%error = 0.0_dp
2609 END SUBROUTINE integration_status_reset
2610
2611! **************************************************************************************************
2612!> \brief Generate an integration method description string.
2613!> \param stats integration statistics
2614!> \param integration_method integration method used
2615!> \return description string
2616!> \author Sergey Chulkov
2617! **************************************************************************************************
2618 ELEMENTAL FUNCTION get_method_description_string(stats, integration_method) RESULT(method_descr)
2619 TYPE(integration_status_type), INTENT(in) :: stats
2620 INTEGER, INTENT(in) :: integration_method
2621 CHARACTER(len=18) :: method_descr
2622
2623 CHARACTER(len=2) :: method_abbr
2624 CHARACTER(len=6) :: npoints_str
2625
2626 SELECT CASE (integration_method)
2627 CASE (negfint_method_cc)
2628 ! Adaptive Clenshaw-Curtis method
2629 method_abbr = "CC"
2630 CASE (negfint_method_simpson)
2631 ! Adaptive Simpson's rule method
2632 method_abbr = "SR"
2633 CASE DEFAULT
2634 method_abbr = "??"
2635 END SELECT
2636
2637 WRITE (npoints_str, '(I6)') stats%npoints
2638 WRITE (method_descr, '(A2,T4,A,T11,ES8.2E2)') method_abbr, trim(adjustl(npoints_str)), stats%error
2639 END FUNCTION get_method_description_string
2640
2641! **************************************************************************************************
2642!> \brief Compute electric current for one spin-channel through the scattering region.
2643!> \param contact_id1 reference contact
2644!> \param contact_id2 another contact
2645!> \param v_shift shift in Hartree potential
2646!> \param negf_env NEFG environment
2647!> \param negf_control NEGF control
2648!> \param sub_env NEGF parallel (sub)group environment
2649!> \param ispin spin conponent to proceed
2650!> \param blacs_env_global global BLACS environment
2651!> \return electric current in Amper
2652!> \author Sergey Chulkov
2653! **************************************************************************************************
2654 FUNCTION negf_compute_current(contact_id1, contact_id2, v_shift, negf_env, negf_control, sub_env, ispin, &
2655 blacs_env_global) RESULT(current)
2656 INTEGER, INTENT(in) :: contact_id1, contact_id2
2657 REAL(kind=dp), INTENT(in) :: v_shift
2658 TYPE(negf_env_type), INTENT(in) :: negf_env
2659 TYPE(negf_control_type), POINTER :: negf_control
2660 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
2661 INTEGER, INTENT(in) :: ispin
2662 TYPE(cp_blacs_env_type), POINTER :: blacs_env_global
2663 REAL(kind=dp) :: current
2664
2665 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_compute_current'
2666 REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*epsilon(0.0_dp)
2667
2668 COMPLEX(kind=dp) :: fermi_contact1, fermi_contact2, &
2669 integr_lbound, integr_ubound
2670 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: transm_coeff, xnodes
2671 COMPLEX(kind=dp), DIMENSION(1, 1) :: transmission
2672 INTEGER :: handle, icontact, ipoint, max_points, &
2673 min_points, ncontacts, npoints, &
2674 npoints_total
2675 REAL(kind=dp) :: conv_density, energy, eta, ln_conv_density, mu_contact1, mu_contact2, &
2676 temperature_contact1, temperature_contact2, v_contact1, v_contact2
2677 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: zdata
2678 TYPE(cp_fm_struct_type), POINTER :: fm_struct_single
2679 TYPE(cp_fm_type) :: weights
2680 TYPE(green_functions_cache_type) :: g_surf_cache
2681 TYPE(simpsonrule_type) :: sr_env
2682
2683 current = 0.0_dp
2684 ! nothing to do
2685 IF (.NOT. ASSOCIATED(negf_env%s_s)) RETURN
2686
2687 CALL timeset(routinen, handle)
2688
2689 ncontacts = SIZE(negf_env%contacts)
2690 cpassert(contact_id1 <= ncontacts)
2691 cpassert(contact_id2 <= ncontacts)
2692 cpassert(contact_id1 /= contact_id2)
2693
2694 v_contact1 = negf_control%contacts(contact_id1)%v_external
2695 mu_contact1 = negf_control%contacts(contact_id1)%fermi_level + v_contact1
2696 v_contact2 = negf_control%contacts(contact_id2)%v_external
2697 mu_contact2 = negf_control%contacts(contact_id2)%fermi_level + v_contact2
2698
2699 IF (abs(mu_contact1 - mu_contact2) < threshold) THEN
2700 CALL timestop(handle)
2701 RETURN
2702 END IF
2703
2704 min_points = negf_control%integr_min_points
2705 max_points = negf_control%integr_max_points
2706 temperature_contact1 = negf_control%contacts(contact_id1)%temperature
2707 temperature_contact2 = negf_control%contacts(contact_id2)%temperature
2708 eta = negf_control%eta
2709 conv_density = negf_control%conv_density
2710 ln_conv_density = log(conv_density)
2711
2712 integr_lbound = cmplx(min(mu_contact1 + ln_conv_density*temperature_contact1, &
2713 mu_contact2 + ln_conv_density*temperature_contact2), eta, kind=dp)
2714 integr_ubound = cmplx(max(mu_contact1 - ln_conv_density*temperature_contact1, &
2715 mu_contact2 - ln_conv_density*temperature_contact2), eta, kind=dp)
2716
2717 npoints_total = 0
2718 npoints = min_points
2719
2720 NULLIFY (fm_struct_single)
2721 CALL cp_fm_struct_create(fm_struct_single, nrow_global=1, ncol_global=1, context=blacs_env_global)
2722 CALL cp_fm_create(weights, fm_struct_single)
2723 CALL cp_fm_set_all(weights, 1.0_dp)
2724
2725 ALLOCATE (transm_coeff(npoints), xnodes(npoints), zdata(npoints))
2726
2727 CALL simpsonrule_init(sr_env, xnodes, npoints, integr_lbound, integr_ubound, &
2728 sr_shape_linear, negf_control%conv_density, weights)
2729
2730 DO WHILE (npoints > 0 .AND. npoints_total < max_points)
2731 CALL green_functions_cache_expand(g_surf_cache, ncontacts, npoints)
2732
2733 DO icontact = 1, ncontacts
2734 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, 1:npoints), &
2735 omega=xnodes(1:npoints), &
2736 h0=negf_env%contacts(icontact)%h_00(ispin), &
2737 s0=negf_env%contacts(icontact)%s_00, &
2738 h1=negf_env%contacts(icontact)%h_01(ispin), &
2739 s1=negf_env%contacts(icontact)%s_01, &
2740 sub_env=sub_env, &
2741 v_external=negf_control%contacts(icontact)%v_external, &
2742 conv=negf_control%conv_green, transp=.false.)
2743 END DO
2744
2745 CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints), &
2746 v_shift=v_shift, &
2747 ignore_bias=.false., &
2748 negf_env=negf_env, &
2749 negf_control=negf_control, &
2750 sub_env=sub_env, &
2751 ispin=ispin, &
2752 g_surf_contacts=g_surf_cache%g_surf_contacts(:, 1:npoints), &
2753 transm_coeff=transm_coeff(1:npoints), &
2754 transm_contact1=contact_id1, &
2755 transm_contact2=contact_id2)
2756
2757 DO ipoint = 1, npoints
2758 CALL cp_cfm_create(zdata(ipoint), fm_struct_single)
2759
2760 energy = real(xnodes(ipoint), kind=dp)
2761 fermi_contact1 = fermi_function(cmplx(energy - mu_contact1, 0.0_dp, kind=dp), temperature_contact1)
2762 fermi_contact2 = fermi_function(cmplx(energy - mu_contact2, 0.0_dp, kind=dp), temperature_contact2)
2763
2764 transmission(1, 1) = transm_coeff(ipoint)*(fermi_contact1 - fermi_contact2)
2765 CALL cp_cfm_set_submatrix(zdata(ipoint), transmission)
2766 END DO
2767
2768 CALL green_functions_cache_release(g_surf_cache)
2769
2770 npoints_total = npoints_total + npoints
2771
2772 CALL simpsonrule_refine_integral(sr_env, zdata(1:npoints))
2773
2774 IF (sr_env%error <= negf_control%conv_density) EXIT
2775
2776 npoints = max_points - npoints_total
2777 IF (npoints <= 0) EXIT
2778 IF (npoints > SIZE(xnodes)) npoints = SIZE(xnodes)
2779
2780 CALL simpsonrule_get_next_nodes(sr_env, xnodes, npoints)
2781 END DO
2782
2783 CALL cp_cfm_get_submatrix(sr_env%integral, transmission)
2784
2785 current = 0.5_dp/pi*real(transmission(1, 1), kind=dp)*e_charge/seconds
2786
2787 CALL cp_fm_release(weights)
2788 CALL cp_fm_struct_release(fm_struct_single)
2789
2790 CALL simpsonrule_release(sr_env)
2791 DEALLOCATE (transm_coeff, xnodes, zdata)
2792
2793 CALL timestop(handle)
2794 END FUNCTION negf_compute_current
2795
2796! **************************************************************************************************
2797!> \brief Print the Density of States.
2798!> \param log_unit output unit
2799!> \param energy_min energy point to start with
2800!> \param energy_max energy point to end with
2801!> \param npoints number of points to compute
2802!> \param v_shift shift in Hartree potential
2803!> \param negf_env NEFG environment
2804!> \param negf_control NEGF control
2805!> \param sub_env NEGF parallel (sub)group environment
2806!> \param base_contact index of the reference contact
2807!> \param just_contact compute DOS for the given contact rather than for a scattering region
2808!> \param volume unit cell volume
2809!> \author Sergey Chulkov
2810! **************************************************************************************************
2811 SUBROUTINE negf_print_dos(log_unit, energy_min, energy_max, npoints, v_shift, negf_env, &
2812 negf_control, sub_env, base_contact, just_contact, volume)
2813 INTEGER, INTENT(in) :: log_unit
2814 REAL(kind=dp), INTENT(in) :: energy_min, energy_max
2815 INTEGER, INTENT(in) :: npoints
2816 REAL(kind=dp), INTENT(in) :: v_shift
2817 TYPE(negf_env_type), INTENT(in) :: negf_env
2818 TYPE(negf_control_type), POINTER :: negf_control
2819 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
2820 INTEGER, INTENT(in) :: base_contact
2821 INTEGER, INTENT(in), OPTIONAL :: just_contact
2822 REAL(kind=dp), INTENT(in), OPTIONAL :: volume
2823
2824 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_print_dos'
2825
2826 CHARACTER :: uks_str
2827 CHARACTER(len=15) :: units_str
2828 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: xnodes
2829 INTEGER :: handle, icontact, ipoint, ispin, &
2830 ncontacts, npoints_bundle, &
2831 npoints_remain, nspins
2832 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dos
2833 TYPE(green_functions_cache_type) :: g_surf_cache
2834
2835 CALL timeset(routinen, handle)
2836
2837 IF (PRESENT(just_contact)) THEN
2838 nspins = SIZE(negf_env%contacts(just_contact)%h_00)
2839 ELSE
2840 nspins = SIZE(negf_env%h_s)
2841 END IF
2842
2843 IF (log_unit > 0) THEN
2844 IF (PRESENT(volume)) THEN
2845 units_str = ' (angstroms^-3)'
2846 ELSE
2847 units_str = ''
2848 END IF
2849
2850 IF (nspins > 1) THEN
2851 ! [alpha , beta]
2852 uks_str = ','
2853 ELSE
2854 ! [alpha + beta]
2855 uks_str = '+'
2856 END IF
2857
2858 IF (PRESENT(just_contact)) THEN
2859 WRITE (log_unit, '(3A,T70,I11)') "# Density of states", trim(units_str), " for the contact No. ", just_contact
2860 ELSE
2861 WRITE (log_unit, '(3A)') "# Density of states", trim(units_str), " for the scattering region"
2862 END IF
2863
2864 WRITE (log_unit, '(A,T10,A,T43,3A)') "#", "Energy (a.u.)", "Number of states [alpha ", uks_str, " beta]"
2865
2866 WRITE (log_unit, '("#", T3,78("-"))')
2867 END IF
2868
2869 ncontacts = SIZE(negf_env%contacts)
2870 cpassert(base_contact <= ncontacts)
2871 IF (PRESENT(just_contact)) THEN
2872 ncontacts = 2
2873 cpassert(just_contact == base_contact)
2874 END IF
2875 mark_used(base_contact)
2876
2877 npoints_bundle = 4*sub_env%ngroups
2878 IF (npoints_bundle > npoints) npoints_bundle = npoints
2879
2880 ALLOCATE (dos(npoints_bundle, nspins), xnodes(npoints_bundle))
2881
2882 npoints_remain = npoints
2883 DO WHILE (npoints_remain > 0)
2884 IF (npoints_bundle > npoints_remain) npoints_bundle = npoints_remain
2885
2886 IF (npoints > 1) THEN
2887 DO ipoint = 1, npoints_bundle
2888 xnodes(ipoint) = cmplx(energy_min + real(npoints - npoints_remain + ipoint - 1, kind=dp)/ &
2889 REAL(npoints - 1, kind=dp)*(energy_max - energy_min), negf_control%eta, kind=dp)
2890 END DO
2891 ELSE
2892 xnodes(ipoint) = cmplx(energy_min, negf_control%eta, kind=dp)
2893 END IF
2894
2895 DO ispin = 1, nspins
2896 CALL green_functions_cache_expand(g_surf_cache, ncontacts, npoints_bundle)
2897
2898 IF (PRESENT(just_contact)) THEN
2899 DO icontact = 1, ncontacts
2900 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, :), &
2901 omega=xnodes(1:npoints_bundle), &
2902 h0=negf_env%contacts(just_contact)%h_00(ispin), &
2903 s0=negf_env%contacts(just_contact)%s_00, &
2904 h1=negf_env%contacts(just_contact)%h_01(ispin), &
2905 s1=negf_env%contacts(just_contact)%s_01, &
2906 sub_env=sub_env, v_external=0.0_dp, &
2907 conv=negf_control%conv_green, transp=(icontact == 1))
2908 END DO
2909 ELSE
2910 DO icontact = 1, ncontacts
2911 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, :), &
2912 omega=xnodes(1:npoints_bundle), &
2913 h0=negf_env%contacts(icontact)%h_00(ispin), &
2914 s0=negf_env%contacts(icontact)%s_00, &
2915 h1=negf_env%contacts(icontact)%h_01(ispin), &
2916 s1=negf_env%contacts(icontact)%s_01, &
2917 sub_env=sub_env, &
2918 v_external=negf_control%contacts(icontact)%v_external, &
2919 conv=negf_control%conv_green, transp=.false.)
2920 END DO
2921 END IF
2922
2923 CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints_bundle), &
2924 v_shift=v_shift, &
2925 ignore_bias=.false., &
2926 negf_env=negf_env, &
2927 negf_control=negf_control, &
2928 sub_env=sub_env, &
2929 ispin=ispin, &
2930 g_surf_contacts=g_surf_cache%g_surf_contacts, &
2931 dos=dos(1:npoints_bundle, ispin), &
2932 just_contact=just_contact)
2933
2934 CALL green_functions_cache_release(g_surf_cache)
2935 END DO
2936
2937 IF (log_unit > 0) THEN
2938 DO ipoint = 1, npoints_bundle
2939 IF (nspins > 1) THEN
2940 ! spin-polarised calculations: print alpha- and beta-spin components separately
2941 WRITE (log_unit, '(T2,F20.8,T30,2ES25.11E3)') real(xnodes(ipoint), kind=dp), dos(ipoint, 1), dos(ipoint, 2)
2942 ELSE
2943 ! spin-restricted calculations: print alpha- and beta-spin components together
2944 WRITE (log_unit, '(T2,F20.8,T43,ES25.11E3)') real(xnodes(ipoint), kind=dp), 2.0_dp*dos(ipoint, 1)
2945 END IF
2946 END DO
2947 END IF
2948
2949 npoints_remain = npoints_remain - npoints_bundle
2950 END DO
2951
2952 DEALLOCATE (dos, xnodes)
2953 CALL timestop(handle)
2954 END SUBROUTINE negf_print_dos
2955
2956! **************************************************************************************************
2957!> \brief Print the transmission coefficient.
2958!> \param log_unit output unit
2959!> \param energy_min energy point to start with
2960!> \param energy_max energy point to end with
2961!> \param npoints number of points to compute
2962!> \param v_shift shift in Hartree potential
2963!> \param negf_env NEFG environment
2964!> \param negf_control NEGF control
2965!> \param sub_env NEGF parallel (sub)group environment
2966!> \param contact_id1 index of a reference contact
2967!> \param contact_id2 index of another contact
2968!> \author Sergey Chulkov
2969! **************************************************************************************************
2970 SUBROUTINE negf_print_transmission(log_unit, energy_min, energy_max, npoints, v_shift, negf_env, &
2971 negf_control, sub_env, contact_id1, contact_id2)
2972 INTEGER, INTENT(in) :: log_unit
2973 REAL(kind=dp), INTENT(in) :: energy_min, energy_max
2974 INTEGER, INTENT(in) :: npoints
2975 REAL(kind=dp), INTENT(in) :: v_shift
2976 TYPE(negf_env_type), INTENT(in) :: negf_env
2977 TYPE(negf_control_type), POINTER :: negf_control
2978 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
2979 INTEGER, INTENT(in) :: contact_id1, contact_id2
2980
2981 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_print_transmission'
2982
2983 CHARACTER :: uks_str
2984 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: xnodes
2985 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: transm_coeff
2986 INTEGER :: handle, icontact, ipoint, ispin, &
2987 ncontacts, npoints_bundle, &
2988 npoints_remain, nspins
2989 REAL(kind=dp) :: rscale
2990 TYPE(green_functions_cache_type) :: g_surf_cache
2991
2992 CALL timeset(routinen, handle)
2993
2994 nspins = SIZE(negf_env%h_s)
2995
2996 IF (nspins > 1) THEN
2997 ! [alpha , beta]
2998 uks_str = ','
2999 ELSE
3000 ! [alpha + beta]
3001 uks_str = '+'
3002 END IF
3003
3004 IF (log_unit > 0) THEN
3005 WRITE (log_unit, '(A)') "# Transmission coefficient (G0 = 2 e^2/h) for the scattering region"
3006
3007 WRITE (log_unit, '(A,T10,A,T39,3A)') "#", "Energy (a.u.)", "Transmission coefficient [alpha ", uks_str, " beta]"
3008 WRITE (log_unit, '("#", T3,78("-"))')
3009 END IF
3010
3011 ncontacts = SIZE(negf_env%contacts)
3012 cpassert(contact_id1 <= ncontacts)
3013 cpassert(contact_id2 <= ncontacts)
3014
3015 IF (nspins == 1) THEN
3016 rscale = 2.0_dp
3017 ELSE
3018 rscale = 1.0_dp
3019 END IF
3020
3021 ! print transmission coefficients in terms of G0 = 2 * e^2 / h = 1 / pi ;
3022 ! transmission coefficients returned by negf_retarded_green_function_batch() are already multiplied by 2 / pi
3023 rscale = 0.5_dp*rscale
3024
3025 npoints_bundle = 4*sub_env%ngroups
3026 IF (npoints_bundle > npoints) npoints_bundle = npoints
3027
3028 ALLOCATE (transm_coeff(npoints_bundle, nspins), xnodes(npoints_bundle))
3029
3030 npoints_remain = npoints
3031 DO WHILE (npoints_remain > 0)
3032 IF (npoints_bundle > npoints_remain) npoints_bundle = npoints_remain
3033
3034 IF (npoints > 1) THEN
3035 DO ipoint = 1, npoints_bundle
3036 xnodes(ipoint) = cmplx(energy_min + real(npoints - npoints_remain + ipoint - 1, kind=dp)/ &
3037 REAL(npoints - 1, kind=dp)*(energy_max - energy_min), negf_control%eta, kind=dp)
3038 END DO
3039 ELSE
3040 xnodes(ipoint) = cmplx(energy_min, negf_control%eta, kind=dp)
3041 END IF
3042
3043 DO ispin = 1, nspins
3044 CALL green_functions_cache_expand(g_surf_cache, ncontacts, npoints_bundle)
3045
3046 DO icontact = 1, ncontacts
3047 CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, :), &
3048 omega=xnodes(1:npoints_bundle), &
3049 h0=negf_env%contacts(icontact)%h_00(ispin), &
3050 s0=negf_env%contacts(icontact)%s_00, &
3051 h1=negf_env%contacts(icontact)%h_01(ispin), &
3052 s1=negf_env%contacts(icontact)%s_01, &
3053 sub_env=sub_env, &
3054 v_external=negf_control%contacts(icontact)%v_external, &
3055 conv=negf_control%conv_green, transp=.false.)
3056 END DO
3057
3058 CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints_bundle), &
3059 v_shift=v_shift, &
3060 ignore_bias=.false., &
3061 negf_env=negf_env, &
3062 negf_control=negf_control, &
3063 sub_env=sub_env, &
3064 ispin=ispin, &
3065 g_surf_contacts=g_surf_cache%g_surf_contacts, &
3066 transm_coeff=transm_coeff(1:npoints_bundle, ispin), &
3067 transm_contact1=contact_id1, &
3068 transm_contact2=contact_id2)
3069
3070 CALL green_functions_cache_release(g_surf_cache)
3071 END DO
3072
3073 IF (log_unit > 0) THEN
3074 DO ipoint = 1, npoints_bundle
3075 IF (nspins > 1) THEN
3076 ! spin-polarised calculations: print alpha- and beta-spin components separately
3077 WRITE (log_unit, '(T2,F20.8,T30,2ES25.11E3)') &
3078 REAL(xnodes(ipoint), kind=dp), rscale*real(transm_coeff(ipoint, 1:2), kind=dp)
3079 ELSE
3080 ! spin-restricted calculations: print alpha- and beta-spin components together
3081 WRITE (log_unit, '(T2,F20.8,T43,ES25.11E3)') &
3082 REAL(xnodes(ipoint), kind=dp), rscale*real(transm_coeff(ipoint, 1), kind=dp)
3083 END IF
3084 END DO
3085 END IF
3086
3087 npoints_remain = npoints_remain - npoints_bundle
3088 END DO
3089
3090 DEALLOCATE (transm_coeff, xnodes)
3091 CALL timestop(handle)
3092 END SUBROUTINE negf_print_transmission
3093END MODULE negf_methods
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:60
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:74
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::papior2017
integer, save, public papior2017
Definition bibliography.F:36
bibliography::bailey2006
integer, save, public bailey2006
Definition bibliography.F:36
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:240
cp_cfm_basic_linalg::cp_cfm_trace
subroutine, public cp_cfm_trace(matrix_a, matrix_b, trace)
Returns the trace of matrix_a^T matrix_b, i.e sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
Definition cp_cfm_basic_linalg.F:910
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_get_submatrix
subroutine, public cp_cfm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
Extract a sub-matrix from the full matrix: op(target_m)(1:n_rows,1:n_cols) = fm(start_row:start_row+n...
Definition cp_cfm_types.F:329
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:155
cp_cfm_types::cp_cfm_start_copy_general
subroutine, public cp_cfm_start_copy_general(source, destination, para_env, info)
Initiate the copy operation: get distribution data, post MPI isend and irecvs.
Definition cp_cfm_types.F:875
cp_cfm_types::cp_cfm_cleanup_copy_general
subroutine, public cp_cfm_cleanup_copy_general(info)
Complete the copy operation: wait for comms clean up MPI state.
Definition cp_cfm_types.F:1254
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:603
cp_cfm_types::cp_cfm_set_submatrix
subroutine, public cp_cfm_set_submatrix(matrix, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
Set a sub-matrix of the full matrix: matrix(start_row:start_row+n_rows,start_col:start_col+n_cols) = ...
Definition cp_cfm_types.F:466
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:761
cp_cfm_types::cp_cfm_finish_copy_general
subroutine, public cp_cfm_finish_copy_general(destination, info)
Complete the copy operation: wait for comms, unpack, clean up MPI state.
Definition cp_cfm_types.F:1187
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:166
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1661
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition cp_fm_types.F:1537
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1009
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:528
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::debug_print_level
integer, parameter, public debug_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::high_print_level
integer, parameter, public high_print_level
Definition cp_output_handling.F:73
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
global_types
Define type storing the global information of a run. Keep the amount of stored data small....
Definition global_types.F:21
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::negfint_method_simpson
integer, parameter, public negfint_method_simpson
Definition input_constants.F:1273
input_constants::negfint_method_cc
integer, parameter, public negfint_method_cc
Definition input_constants.F:1273
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:365
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:153
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
negf_control_types
Input control types for NEGF based quantum transport calculations.
Definition negf_control_types.F:12
negf_control_types::negf_control_create
subroutine, public negf_control_create(negf_control)
allocate control options for Non-equilibrium Green's Function calculation
Definition negf_control_types.F:132
negf_control_types::read_negf_control
subroutine, public read_negf_control(negf_control, input, subsys)
Read NEGF input parameters.
Definition negf_control_types.F:198
negf_control_types::negf_control_release
subroutine, public negf_control_release(negf_control)
release memory allocated for NEGF control options
Definition negf_control_types.F:153
negf_env_types
Environment for NEGF based quantum transport calculations.
Definition negf_env_types.F:11
negf_env_types::negf_env_release
subroutine, public negf_env_release(negf_env)
Release a NEGF environment variable.
Definition negf_env_types.F:1192
negf_env_types::negf_env_create
subroutine, public negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
Create a new NEGF environment and compute the relevant Kohn-Sham matrices.
Definition negf_env_types.F:155
negf_green_cache
Storage to keep precomputed surface Green's functions.
Definition negf_green_cache.F:11
negf_green_cache::green_functions_cache_reorder
subroutine, public green_functions_cache_reorder(cache, tnodes)
Sort cached items in ascending order.
Definition negf_green_cache.F:87
negf_green_cache::green_functions_cache_release
subroutine, public green_functions_cache_release(cache)
Release storage.
Definition negf_green_cache.F:124
negf_green_cache::green_functions_cache_expand
subroutine, public green_functions_cache_expand(cache, ncontacts, nnodes_extra)
Reallocate storage so it can handle extra 'nnodes_extra' items for each contact.
Definition negf_green_cache.F:53
negf_green_methods
Subroutines to compute Green functions.
Definition negf_green_methods.F:11
negf_green_methods::sancho_work_matrices_create
subroutine, public sancho_work_matrices_create(work, fm_struct)
Create work matrices required for the Lopez-Sancho algorithm.
Definition negf_green_methods.F:73
negf_green_methods::sancho_work_matrices_release
subroutine, public sancho_work_matrices_release(work)
Release work matrices.
Definition negf_green_methods.F:106
negf_green_methods::negf_contact_self_energy
subroutine, public negf_contact_self_energy(self_energy_c, omega, g_surf_c, h_sc0, s_sc0, zwork1, zwork2, transp)
Compute the contact self energy at point 'omega' as self_energy_C = [omega * S_SC0 - KS_SC0] * g_surf...
Definition negf_green_methods.F:270
negf_green_methods::negf_contact_broadening_matrix
subroutine, public negf_contact_broadening_matrix(gamma_c, self_energy_c)
Compute contact broadening matrix as gamma_C = i (self_energy_c^{ret.} - (self_energy_c^{ret....
Definition negf_green_methods.F:315
negf_green_methods::do_sancho
subroutine, public do_sancho(g_surf, omega, h0, s0, h1, s1, conv, transp, work)
Iterative method to compute a retarded surface Green's function at the point omega.
Definition negf_green_methods.F:172
negf_green_methods::negf_retarded_green_function
subroutine, public negf_retarded_green_function(g_ret_s, omega, self_energy_ret_sum, h_s, s_s, v_hartree_s)
Compute the retarded Green's function at point 'omega' as G_S^{ret.} = [ omega * S_S - KS_S - \sum_{c...
Definition negf_green_methods.F:349
negf_integr_cc
Adaptive Clenshaw-Curtis quadrature algorithm to integrate a complex-valued function in a complex pla...
Definition negf_integr_cc.F:14
negf_integr_cc::cc_shape_linear
integer, parameter, public cc_shape_linear
Definition negf_integr_cc.F:52
negf_integr_cc::ccquad_refine_integral
subroutine, public ccquad_refine_integral(cc_env)
Refine approximated integral.
Definition negf_integr_cc.F:429
negf_integr_cc::cc_interval_full
integer, parameter, public cc_interval_full
Definition negf_integr_cc.F:49
negf_integr_cc::ccquad_double_number_of_points
subroutine, public ccquad_double_number_of_points(cc_env, xnodes_next)
Get the next set of points at which the integrand needs to be computed. These points are then can be ...
Definition negf_integr_cc.F:251
negf_integr_cc::ccquad_release
subroutine, public ccquad_release(cc_env)
Release a Clenshaw-Curtis quadrature environment variable.
Definition negf_integr_cc.F:203
negf_integr_cc::cc_interval_half
integer, parameter, public cc_interval_half
Definition negf_integr_cc.F:49
negf_integr_cc::ccquad_init
subroutine, public ccquad_init(cc_env, xnodes, nnodes, a, b, interval_id, shape_id, weights, tnodes_restart)
Initialise a Clenshaw-Curtis quadrature environment variable.
Definition negf_integr_cc.F:114
negf_integr_cc::cc_shape_arc
integer, parameter, public cc_shape_arc
Definition negf_integr_cc.F:52
negf_integr_cc::ccquad_reduce_and_append_zdata
subroutine, public ccquad_reduce_and_append_zdata(cc_env, zdata_next)
Prepare Clenshaw-Curtis environment for the subsequent refinement of the integral.
Definition negf_integr_cc.F:334
negf_integr_simpson
Adaptive Simpson's rule algorithm to integrate a complex-valued function in a complex plane.
Definition negf_integr_simpson.F:11
negf_integr_simpson::sr_shape_arc
integer, parameter, public sr_shape_arc
Definition negf_integr_simpson.F:46
negf_integr_simpson::simpsonrule_refine_integral
subroutine, public simpsonrule_refine_integral(sr_env, zdata_next)
Compute integral using the simpson's rules.
Definition negf_integr_simpson.F:341
negf_integr_simpson::simpsonrule_init
subroutine, public simpsonrule_init(sr_env, xnodes, nnodes, a, b, shape_id, conv, weights, tnodes_restart)
Initialise a Simpson's rule environment variable.
Definition negf_integr_simpson.F:122
negf_integr_simpson::simpsonrule_release
subroutine, public simpsonrule_release(sr_env)
Release a Simpson's rule environment variable.
Definition negf_integr_simpson.F:189
negf_integr_simpson::sr_shape_linear
integer, parameter, public sr_shape_linear
Definition negf_integr_simpson.F:46
negf_integr_simpson::simpsonrule_get_next_nodes
subroutine, public simpsonrule_get_next_nodes(sr_env, xnodes_next, nnodes)
Get the next set of nodes where to compute integrand.
Definition negf_integr_simpson.F:260
negf_matrix_utils
Helper routines to manipulate with matrices.
Definition negf_matrix_utils.F:11
negf_matrix_utils::invert_cell_to_index
subroutine, public invert_cell_to_index(cell_to_index, nimages, index_to_cell)
Invert cell_to_index mapping between unit cells and DBCSR matrix images.
Definition negf_matrix_utils.F:636
negf_matrix_utils::negf_copy_fm_submat_to_dbcsr
subroutine, public negf_copy_fm_submat_to_dbcsr(fm, matrix, atomlist_row, atomlist_col, subsys)
Populate relevant blocks of the DBCSR matrix using data from a ScaLAPACK matrix. Irrelevant blocks of...
Definition negf_matrix_utils.F:105
negf_matrix_utils::negf_copy_sym_dbcsr_to_fm_submat
subroutine, public negf_copy_sym_dbcsr_to_fm_submat(matrix, fm, atomlist_row, atomlist_col, subsys, mpi_comm_global, do_upper_diag, do_lower)
Extract part of the DBCSR matrix based on selected atoms and copy it into a dense matrix.
Definition negf_matrix_utils.F:192
negf_methods
NEGF based quantum transport calculations.
Definition negf_methods.F:11
negf_methods::do_negf
subroutine, public do_negf(force_env)
Perform NEGF calculation.
Definition negf_methods.F:156
negf_subgroup_types
Environment for NEGF based quantum transport calculations.
Definition negf_subgroup_types.F:13
negf_subgroup_types::negf_sub_env_release
subroutine, public negf_sub_env_release(sub_env)
Release a parallel (sub)group environment.
Definition negf_subgroup_types.F:129
negf_subgroup_types::negf_sub_env_create
subroutine, public negf_sub_env_create(sub_env, negf_control, blacs_env_global, blacs_grid_layout, blacs_repeatable)
Split MPI communicator to create a set of parallel (sub)groups.
Definition negf_subgroup_types.F:69
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::e_charge
real(kind=dp), parameter, public e_charge
Definition physcon.F:106
physcon::kelvin
real(kind=dp), parameter, public kelvin
Definition physcon.F:165
physcon::seconds
real(kind=dp), parameter, public seconds
Definition physcon.F:150
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
qs_density_mixing_types
module that contains the definitions of the scf types
Definition qs_density_mixing_types.F:14
qs_density_mixing_types::direct_mixing_nr
integer, parameter, public direct_mixing_nr
Definition qs_density_mixing_types.F:45
qs_density_mixing_types::gspace_mixing_nr
integer, parameter, public gspace_mixing_nr
Definition qs_density_mixing_types.F:45
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:526
qs_gspace_mixing
Definition qs_gspace_mixing.F:9
qs_gspace_mixing::gspace_mixing
subroutine, public gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count)
Driver for the g-space mixing, calls the proper routine given the requested method.
Definition qs_gspace_mixing.F:66
qs_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_methods.F:22
qs_ks_methods::qs_ks_build_kohn_sham_matrix
subroutine, public qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, print_active, ext_ks_matrix)
routine where the real calculations are made: the KS matrix is calculated
Definition qs_ks_methods.F:179
qs_mixing_utils
Definition qs_mixing_utils.F:9
qs_mixing_utils::mixing_init
subroutine, public mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom)
initialiation needed when gspace mixing is used
Definition qs_mixing_utils.F:397
qs_mixing_utils::mixing_allocate
subroutine, public mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
allocation needed when density mixing is used
Definition qs_mixing_utils.F:110
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition qs_rho_methods.F:15
qs_rho_methods::qs_rho_update_rho
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
Definition qs_rho_methods.F:378
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_methods
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
Definition qs_scf_methods.F:19
qs_scf_methods::scf_env_density_mixing
subroutine, public scf_env_density_mixing(p_mix_new, mixing_store, rho_ao, para_env, iter_delta, iter_count, diis, invert)
perform (if requested) a density mixing
Definition qs_scf_methods.F:94
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::integer_to_string
subroutine, public integer_to_string(inumber, string)
Converts an integer number to a string. The WRITE statement will return an error message,...
Definition string_utilities.F:3239
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_cfm_types::copy_cfm_info_type
Stores the state of a copy between cp_cfm_start_copy_general and cp_cfm_finish_copy_general.
Definition cp_cfm_types.F:89
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68
cp_control_types::dft_control_type
Definition cp_control_types.F:622
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
force_env_types::force_env_p_type
allows for the creation of an array of force_env
Definition force_env_types.F:174
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
global_types::global_environment_type
contains the initially parsed file and the initial parallel environment
Definition global_types.F:66
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:150
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
negf_control_types::negf_control_type
Input parameters related to the NEGF run.
Definition negf_control_types.F:72
negf_env_types::negf_env_type
NEGF environment.
Definition negf_env_types.F:113
negf_green_cache::green_functions_cache_type
Storage to keep surface Green's functions.
Definition negf_green_cache.F:34
negf_green_methods::sancho_work_matrices_type
Definition negf_green_methods.F:46
negf_integr_cc::ccquad_type
Adaptive Clenshaw-Curtis environment.
Definition negf_integr_cc.F:66
negf_integr_simpson::simpsonrule_type
A structure to store data needed for adaptive Simpson's rule algorithm.
Definition negf_integr_simpson.F:70
negf_subgroup_types::negf_subgroup_env_type
Parallel (sub)group environment.
Definition negf_subgroup_types.F:36
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56