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