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