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