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