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