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