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