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