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