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