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