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