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