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