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