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