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