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