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