LCOV - code coverage report
Current view: top level - src - negf_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 871 1077 80.9 %
Date: 2024-03-28 07:31:50 Functions: 14 16 87.5 %

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

Generated by: LCOV version 1.15