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