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