LCOV - code coverage report
Current view: top level - src - negf_env_types.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 92.6 % 421 390
Test Date: 2025-07-25 12:55:17 Functions: 66.7 % 18 12

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Environment for NEGF based quantum transport calculations
      10              : ! **************************************************************************************************
      11              : MODULE negf_env_types
      12              :    USE cell_types,                      ONLY: cell_type,&
      13              :                                               real_to_scaled
      14              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      15              :    USE cp_control_types,                ONLY: dft_control_type
      16              :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      17              :                                               dbcsr_deallocate_matrix,&
      18              :                                               dbcsr_init_p,&
      19              :                                               dbcsr_p_type,&
      20              :                                               dbcsr_set,&
      21              :                                               dbcsr_type
      22              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      23              :                                               cp_fm_struct_release,&
      24              :                                               cp_fm_struct_type
      25              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      26              :                                               cp_fm_release,&
      27              :                                               cp_fm_type
      28              :    USE force_env_types,                 ONLY: force_env_get,&
      29              :                                               force_env_p_type,&
      30              :                                               force_env_type,&
      31              :                                               use_qs_force
      32              :    USE input_section_types,             ONLY: section_vals_type,&
      33              :                                               section_vals_val_get
      34              :    USE kinds,                           ONLY: default_string_length,&
      35              :                                               dp
      36              :    USE kpoint_types,                    ONLY: get_kpoint_env,&
      37              :                                               get_kpoint_info,&
      38              :                                               kpoint_env_p_type,&
      39              :                                               kpoint_type
      40              :    USE message_passing,                 ONLY: mp_para_env_type
      41              :    USE negf_atom_map,                   ONLY: negf_atom_map_type,&
      42              :                                               negf_map_atomic_indices
      43              :    USE negf_control_types,              ONLY: negf_control_contact_type,&
      44              :                                               negf_control_type
      45              :    USE negf_matrix_utils,               ONLY: invert_cell_to_index,&
      46              :                                               negf_copy_contact_matrix,&
      47              :                                               negf_copy_sym_dbcsr_to_fm_submat,&
      48              :                                               negf_reference_contact_matrix,&
      49              :                                               number_of_atomic_orbitals
      50              :    USE negf_subgroup_types,             ONLY: negf_subgroup_env_type
      51              :    USE negf_vectors,                    ONLY: contact_direction_vector,&
      52              :                                               projection_on_direction_vector
      53              :    USE particle_types,                  ONLY: particle_type
      54              :    USE pw_env_types,                    ONLY: pw_env_get,&
      55              :                                               pw_env_type
      56              :    USE pw_pool_types,                   ONLY: pw_pool_type
      57              :    USE pw_types,                        ONLY: pw_r3d_rs_type
      58              :    USE qs_density_mixing_types,         ONLY: mixing_storage_create,&
      59              :                                               mixing_storage_release,&
      60              :                                               mixing_storage_type
      61              :    USE qs_energy,                       ONLY: qs_energies
      62              :    USE qs_energy_init,                  ONLY: qs_energies_init
      63              :    USE qs_environment_types,            ONLY: get_qs_env,&
      64              :                                               qs_environment_type
      65              :    USE qs_integrate_potential,          ONLY: integrate_v_rspace
      66              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      67              :                                               mo_set_type
      68              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      69              :                                               qs_rho_type
      70              :    USE qs_subsys_types,                 ONLY: qs_subsys_get,&
      71              :                                               qs_subsys_type
      72              : #include "./base/base_uses.f90"
      73              : 
      74              :    IMPLICIT NONE
      75              :    PRIVATE
      76              : 
      77              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_env_types'
      78              :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .TRUE.
      79              : 
      80              :    PUBLIC :: negf_env_type, negf_env_contact_type
      81              :    PUBLIC :: negf_env_create, negf_env_release
      82              : 
      83              : ! **************************************************************************************************
      84              : !> \brief  Contact-specific NEGF environment.
      85              : !> \author Sergey Chulkov
      86              : ! **************************************************************************************************
      87              :    TYPE negf_env_contact_type
      88              :       REAL(kind=dp), DIMENSION(3)                        :: direction_vector = -1.0_dp, origin = -1.0_dp
      89              :       REAL(kind=dp), DIMENSION(3)                        :: direction_vector_bias = -1.0_dp, origin_bias = -1.0_dp
      90              :       !> an axis towards the secondary contact unit cell which coincides with the transport direction
      91              :       !> 0 (undefined), 1 (+x), 2 (+y), 3 (+z), -1 (-x), -2 (-y), -3 (-z)
      92              :       INTEGER                                            :: direction_axis = -1
      93              :       !> atoms belonging to a primary contact unit cell
      94              :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atomlist_cell0
      95              :       !> atoms belonging to a secondary contact unit cell (will be removed one day ...)
      96              :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atomlist_cell1
      97              :       !> list of equivalent atoms in an appropriate contact force environment
      98              :       TYPE(negf_atom_map_type), ALLOCATABLE, &
      99              :          DIMENSION(:)                                    :: atom_map_cell0, atom_map_cell1
     100              :       !> energy of the HOMO
     101              :       REAL(kind=dp)                                      :: homo_energy = -1.0_dp
     102              :       !> diagonal (h_00) and off-diagonal (h_01) blocks of the contact Kohn-Sham matrix ([number_of_spins]).
     103              :       !> The matrix h_01 is of the shape [nao_cell0 x nao_cell1]
     104              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)      :: h_00, h_01
     105              :       !> diagonal and off-diagonal blocks of the density matrix
     106              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)      :: rho_00, rho_01
     107              :       !> diagonal and off-diagonal blocks of the overlap matrix
     108              :       TYPE(cp_fm_type), POINTER                          :: s_00 => null(), s_01 => null()
     109              :    END TYPE negf_env_contact_type
     110              : 
     111              : ! **************************************************************************************************
     112              : !> \brief  NEGF environment.
     113              : !> \author Sergey Chulkov
     114              : ! **************************************************************************************************
     115              :    TYPE negf_env_type
     116              :       !> contact-specific NEGF environments
     117              :       TYPE(negf_env_contact_type), ALLOCATABLE, &
     118              :          DIMENSION(:)                                     :: contacts
     119              :       !> Kohn-Sham matrix of the scattering region
     120              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)       :: h_s
     121              :       !> Kohn-Sham matrix of the scattering region -- contact interface ([nspins, ncontacts])
     122              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)    :: h_sc
     123              :       !> overlap matrix of the scattering region
     124              :       TYPE(cp_fm_type), POINTER                           :: s_s => null()
     125              :       !> an external Hartree potential in atomic basis set representation
     126              :       TYPE(cp_fm_type), POINTER                           :: v_hartree_s => null()
     127              :       !> overlap matrix of the scattering region -- contact interface for every contact ([ncontacts])
     128              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)       :: s_sc
     129              :       !> structure needed for density mixing
     130              :       TYPE(mixing_storage_type), POINTER                  :: mixing_storage => NULL()
     131              :       !> density mixing method
     132              :       INTEGER                                             :: mixing_method = -1
     133              :    END TYPE negf_env_type
     134              : 
     135              : ! **************************************************************************************************
     136              : !> \brief  Allocatable list of the type 'negf_atom_map_type'.
     137              : !> \author Sergey Chulkov
     138              : ! **************************************************************************************************
     139              :    TYPE negf_atom_map_contact_type
     140              :       TYPE(negf_atom_map_type), ALLOCATABLE, DIMENSION(:) :: atom_map
     141              :    END TYPE negf_atom_map_contact_type
     142              : 
     143              : CONTAINS
     144              : 
     145              : ! **************************************************************************************************
     146              : !> \brief Create a new NEGF environment and compute the relevant Kohn-Sham matrices.
     147              : !> \param negf_env            NEGF environment to create
     148              : !> \param sub_env             NEGF parallel (sub)group environment
     149              : !> \param negf_control        NEGF control
     150              : !> \param force_env           the primary force environment
     151              : !> \param negf_mixing_section pointer to a mixing section within the NEGF input section
     152              : !> \param log_unit            output unit number
     153              : !> \par History
     154              : !>   * 01.2017 created [Sergey Chulkov]
     155              : ! **************************************************************************************************
     156            4 :    SUBROUTINE negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
     157              :       TYPE(negf_env_type), INTENT(inout)                 :: negf_env
     158              :       TYPE(negf_subgroup_env_type), INTENT(in)           :: sub_env
     159              :       TYPE(negf_control_type), POINTER                   :: negf_control
     160              :       TYPE(force_env_type), POINTER                      :: force_env
     161              :       TYPE(section_vals_type), POINTER                   :: negf_mixing_section
     162              :       INTEGER, INTENT(in)                                :: log_unit
     163              : 
     164              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'negf_env_create'
     165              : 
     166              :       CHARACTER(len=default_string_length)               :: contact_str, force_env_str, &
     167              :                                                             n_force_env_str
     168              :       INTEGER                                            :: handle, icontact, in_use, n_force_env, &
     169              :                                                             ncontacts
     170              :       LOGICAL                                            :: do_kpoints
     171              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     172            4 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_kp, matrix_s_kp
     173              :       TYPE(dft_control_type), POINTER                    :: dft_control
     174            4 :       TYPE(force_env_p_type), DIMENSION(:), POINTER      :: sub_force_env
     175              :       TYPE(negf_atom_map_contact_type), ALLOCATABLE, &
     176            4 :          DIMENSION(:)                                    :: map_contact
     177              :       TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_rspace
     178              :       TYPE(qs_environment_type), POINTER                 :: qs_env, qs_env_contact
     179              :       TYPE(qs_subsys_type), POINTER                      :: subsys, subsys_contact
     180              : 
     181            4 :       CALL timeset(routineN, handle)
     182              : 
     183              :       ! ensure we have Quickstep enabled for all force_env
     184            4 :       NULLIFY (sub_force_env)
     185            4 :       CALL force_env_get(force_env, in_use=in_use, qs_env=qs_env, sub_force_env=sub_force_env)
     186              : 
     187            4 :       IF (ASSOCIATED(sub_force_env)) THEN
     188            2 :          n_force_env = SIZE(sub_force_env)
     189              :       ELSE
     190            2 :          n_force_env = 0
     191              :       END IF
     192              : 
     193            4 :       IF (in_use == use_qs_force) THEN
     194            8 :          DO icontact = 1, n_force_env
     195            4 :             CALL force_env_get(sub_force_env(icontact)%force_env, in_use=in_use)
     196            8 :             IF (in_use /= use_qs_force) EXIT
     197              :          END DO
     198              :       END IF
     199              : 
     200            4 :       IF (in_use /= use_qs_force) THEN
     201            0 :          CPABORT("Quickstep is required for NEGF run.")
     202              :       END IF
     203              : 
     204              :       ! check that all mentioned FORCE_EVAL sections are actually present
     205            4 :       ncontacts = SIZE(negf_control%contacts)
     206              : 
     207           12 :       DO icontact = 1, ncontacts
     208           12 :          IF (negf_control%contacts(icontact)%force_env_index > n_force_env) THEN
     209            0 :             WRITE (contact_str, '(I11)') icontact
     210            0 :             WRITE (force_env_str, '(I11)') negf_control%contacts(icontact)%force_env_index
     211            0 :             WRITE (n_force_env_str, '(I11)') n_force_env
     212              : 
     213              :             CALL cp_abort(__LOCATION__, &
     214              :                           "Contact number "//TRIM(ADJUSTL(contact_str))//" is linked with the FORCE_EVAL section number "// &
     215              :                           TRIM(ADJUSTL(force_env_str))//", however only "//TRIM(ADJUSTL(n_force_env_str))// &
     216              :                           " FORCE_EVAL sections have been found. Note that FORCE_EVAL sections are enumerated from 0"// &
     217            0 :                           " and that the primary (0-th) section must contain all the atoms.")
     218              :          END IF
     219              :       END DO
     220              : 
     221              :       ! create basic matrices and neighbour lists for the primary force_env,
     222              :       ! so we know how matrix elements are actually distributed across CPUs.
     223            4 :       CALL qs_energies_init(qs_env, calc_forces=.FALSE.)
     224              :       CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, &
     225              :                       matrix_s_kp=matrix_s_kp, matrix_ks_kp=matrix_ks_kp, &
     226            4 :                       subsys=subsys, v_hartree_rspace=v_hartree_rspace)
     227              : 
     228            4 :       IF (do_kpoints) THEN
     229            0 :          CPABORT("k-points are currently not supported for device FORCE_EVAL")
     230              :       END IF
     231              : 
     232              :       ! stage 1: map the atoms between the device force_env and all contact force_env-s
     233           80 :       ALLOCATE (negf_env%contacts(ncontacts))
     234           20 :       ALLOCATE (map_contact(ncontacts))
     235              : 
     236           12 :       DO icontact = 1, ncontacts
     237           12 :          IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
     238            4 :             CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
     239            4 :             CALL get_qs_env(qs_env_contact, subsys=subsys_contact)
     240              : 
     241              :             CALL negf_env_contact_init_maps(contact_env=negf_env%contacts(icontact), &
     242              :                                             contact_control=negf_control%contacts(icontact), &
     243              :                                             atom_map=map_contact(icontact)%atom_map, &
     244              :                                             eps_geometry=negf_control%eps_geometry, &
     245              :                                             subsys_device=subsys, &
     246            4 :                                             subsys_contact=subsys_contact)
     247              : 
     248            4 :             IF (negf_env%contacts(icontact)%direction_axis == 0) THEN
     249            0 :                WRITE (contact_str, '(I11)') icontact
     250            0 :                WRITE (force_env_str, '(I11)') negf_control%contacts(icontact)%force_env_index
     251              :                CALL cp_abort(__LOCATION__, &
     252              :                              "One lattice vector of the contact unit cell (FORCE_EVAL section "// &
     253              :                              TRIM(ADJUSTL(force_env_str))//") must be parallel to the direction of the contact "// &
     254            0 :                              TRIM(ADJUSTL(contact_str))//".")
     255              :             END IF
     256              :          END IF
     257              :       END DO
     258              : 
     259              :       ! stage 2: obtain relevant Kohn-Sham matrix blocks for each contact (separate bulk DFT calculation)
     260           12 :       DO icontact = 1, ncontacts
     261           12 :          IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
     262            4 :             IF (log_unit > 0) &
     263            2 :                WRITE (log_unit, '(/,T2,A,T70,I11)') "NEGF| Construct the Kohn-Sham matrix for the contact ", icontact
     264              : 
     265            4 :             CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
     266              : 
     267            4 :             CALL qs_energies(qs_env_contact, consistent_energies=.FALSE., calc_forces=.FALSE.)
     268              : 
     269              :             CALL negf_env_contact_init_matrices(contact_env=negf_env%contacts(icontact), sub_env=sub_env, &
     270            4 :                                                 qs_env_contact=qs_env_contact, matrix_s_device=matrix_s_kp(1, 1)%matrix)
     271              :          END IF
     272              :       END DO
     273              : 
     274              :       ! obtain an initial KS-matrix for the scattering region
     275            4 :       CALL qs_energies(qs_env, consistent_energies=.FALSE., calc_forces=.FALSE.)
     276              : 
     277              :       ! *** obtain relevant Kohn-Sham matrix blocks for each contact with no separate FORCE_ENV ***
     278           12 :       DO icontact = 1, ncontacts
     279           12 :          IF (negf_control%contacts(icontact)%force_env_index <= 0) THEN
     280              :             CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
     281              :                                                       contact_control=negf_control%contacts(icontact), &
     282              :                                                       sub_env=sub_env, qs_env=qs_env, &
     283            4 :                                                       eps_geometry=negf_control%eps_geometry)
     284              :          END IF
     285              :       END DO
     286              : 
     287              :       ! extract device-related matrix blocks
     288            4 :       CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
     289              : 
     290              :       ! electron density mixing;
     291              :       ! the input section below should be consistent with the subroutine create_negf_section()
     292            4 :       NULLIFY (negf_env%mixing_storage)
     293            4 :       CALL section_vals_val_get(negf_mixing_section, "METHOD", i_val=negf_env%mixing_method)
     294              : 
     295            4 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     296           16 :       ALLOCATE (negf_env%mixing_storage)
     297              :       CALL mixing_storage_create(negf_env%mixing_storage, negf_mixing_section, &
     298            4 :                                  negf_env%mixing_method, dft_control%qs_control%cutoff)
     299              : 
     300            4 :       CALL timestop(handle)
     301           16 :    END SUBROUTINE negf_env_create
     302              : 
     303              : ! **************************************************************************************************
     304              : !> \brief Establish mapping between the primary and the contact force environments
     305              : !> \param contact_env         NEGF environment for the given contact (modified on exit)
     306              : !> \param contact_control     NEGF control
     307              : !> \param atom_map            atomic map
     308              : !> \param eps_geometry        accuracy in mapping atoms between different force environments
     309              : !> \param subsys_device       QuickStep subsystem of the device force environment
     310              : !> \param subsys_contact      QuickStep subsystem of the contact force environment
     311              : !> \author Sergey Chulkov
     312              : ! **************************************************************************************************
     313            4 :    SUBROUTINE negf_env_contact_init_maps(contact_env, contact_control, atom_map, &
     314              :                                          eps_geometry, subsys_device, subsys_contact)
     315              :       TYPE(negf_env_contact_type), INTENT(inout)         :: contact_env
     316              :       TYPE(negf_control_contact_type), INTENT(in)        :: contact_control
     317              :       TYPE(negf_atom_map_type), ALLOCATABLE, &
     318              :          DIMENSION(:), INTENT(inout)                     :: atom_map
     319              :       REAL(kind=dp), INTENT(in)                          :: eps_geometry
     320              :       TYPE(qs_subsys_type), POINTER                      :: subsys_device, subsys_contact
     321              : 
     322              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_env_contact_init_maps'
     323              : 
     324              :       INTEGER                                            :: handle, natoms
     325              : 
     326            4 :       CALL timeset(routineN, handle)
     327              : 
     328              :       CALL contact_direction_vector(contact_env%origin, &
     329              :                                     contact_env%direction_vector, &
     330              :                                     contact_env%origin_bias, &
     331              :                                     contact_env%direction_vector_bias, &
     332              :                                     contact_control%atomlist_screening, &
     333              :                                     contact_control%atomlist_bulk, &
     334            4 :                                     subsys_device)
     335              : 
     336            4 :       contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys_contact, eps_geometry)
     337              : 
     338            4 :       IF (contact_env%direction_axis /= 0) THEN
     339            4 :          natoms = SIZE(contact_control%atomlist_bulk)
     340           56 :          ALLOCATE (atom_map(natoms))
     341              : 
     342              :          ! map atom listed in 'contact_control%atomlist_bulk' to the corresponding atom/cell replica from the contact force_env
     343              :          CALL negf_map_atomic_indices(atom_map=atom_map, &
     344              :                                       atom_list=contact_control%atomlist_bulk, &
     345              :                                       subsys_device=subsys_device, &
     346              :                                       subsys_contact=subsys_contact, &
     347            4 :                                       eps_geometry=eps_geometry)
     348              : 
     349              :          ! list atoms from 'contact_control%atomlist_bulk' which belong to
     350              :          ! the primary unit cell of the bulk region for the given contact
     351              :          CALL list_atoms_in_bulk_primary_unit_cell(atomlist_cell0=contact_env%atomlist_cell0, &
     352              :                                                    atom_map_cell0=contact_env%atom_map_cell0, &
     353              :                                                    atomlist_bulk=contact_control%atomlist_bulk, &
     354              :                                                    atom_map=atom_map, &
     355              :                                                    origin=contact_env%origin, &
     356              :                                                    direction_vector=contact_env%direction_vector, &
     357              :                                                    direction_axis=contact_env%direction_axis, &
     358            4 :                                                    subsys_device=subsys_device)
     359              : 
     360              :          ! secondary unit cell
     361              :          CALL list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1=contact_env%atomlist_cell1, &
     362              :                                                      atom_map_cell1=contact_env%atom_map_cell1, &
     363              :                                                      atomlist_bulk=contact_control%atomlist_bulk, &
     364              :                                                      atom_map=atom_map, &
     365              :                                                      origin=contact_env%origin, &
     366              :                                                      direction_vector=contact_env%direction_vector, &
     367              :                                                      direction_axis=contact_env%direction_axis, &
     368            4 :                                                      subsys_device=subsys_device)
     369              :       END IF
     370              : 
     371            4 :       CALL timestop(handle)
     372            4 :    END SUBROUTINE negf_env_contact_init_maps
     373              : 
     374              : ! **************************************************************************************************
     375              : !> \brief Extract relevant matrix blocks for the given contact.
     376              : !> \param contact_env         NEGF environment for the contact (modified on exit)
     377              : !> \param sub_env             NEGF parallel (sub)group environment
     378              : !> \param qs_env_contact      QuickStep environment for the contact force environment
     379              : !> \param matrix_s_device     overlap matrix from device force environment
     380              : !> \author Sergey Chulkov
     381              : ! **************************************************************************************************
     382            4 :    SUBROUTINE negf_env_contact_init_matrices(contact_env, sub_env, qs_env_contact, matrix_s_device)
     383              :       TYPE(negf_env_contact_type), INTENT(inout)         :: contact_env
     384              :       TYPE(negf_subgroup_env_type), INTENT(in)           :: sub_env
     385              :       TYPE(qs_environment_type), POINTER                 :: qs_env_contact
     386              :       TYPE(dbcsr_type), POINTER                          :: matrix_s_device
     387              : 
     388              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_env_contact_init_matrices'
     389              : 
     390              :       INTEGER                                            :: handle, iatom, ispin, nao, natoms, &
     391              :                                                             nimages, nspins
     392            4 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_list0, atom_list1
     393            4 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: index_to_cell, is_same_cell
     394            4 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     395              :       LOGICAL                                            :: do_kpoints
     396              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     397            4 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
     398              :       TYPE(dbcsr_type), POINTER                          :: matrix_s_ref
     399              :       TYPE(dft_control_type), POINTER                    :: dft_control
     400              :       TYPE(kpoint_type), POINTER                         :: kpoints
     401              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     402              :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     403              :       TYPE(qs_subsys_type), POINTER                      :: subsys
     404              : 
     405            4 :       CALL timeset(routineN, handle)
     406              : 
     407              :       CALL get_qs_env(qs_env_contact, &
     408              :                       dft_control=dft_control, &
     409              :                       do_kpoints=do_kpoints, &
     410              :                       kpoints=kpoints, &
     411              :                       matrix_ks_kp=matrix_ks_kp, &
     412              :                       matrix_s_kp=matrix_s_kp, &
     413              :                       para_env=para_env, &
     414              :                       rho=rho_struct, &
     415            4 :                       subsys=subsys)
     416            4 :       CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
     417              : 
     418            4 :       CALL negf_homo_energy_estimate(contact_env%homo_energy, qs_env_contact)
     419              : 
     420            4 :       natoms = SIZE(contact_env%atomlist_cell0)
     421           12 :       ALLOCATE (atom_list0(natoms))
     422           20 :       DO iatom = 1, natoms
     423           16 :          atom_list0(iatom) = contact_env%atom_map_cell0(iatom)%iatom
     424              : 
     425              :          ! with no k-points there is one-to-one correspondence between the primary unit cell
     426              :          ! of the contact force_env and the first contact unit cell of the device force_env
     427           68 :          IF (SUM(ABS(contact_env%atom_map_cell0(iatom)%cell(:))) > 0) THEN
     428            0 :             CPABORT("NEGF K-points are not currently supported")
     429              :          END IF
     430              :       END DO
     431              : 
     432            4 :       CPASSERT(SIZE(contact_env%atomlist_cell1) == natoms)
     433           12 :       ALLOCATE (atom_list1(natoms))
     434           20 :       DO iatom = 1, natoms
     435           20 :          atom_list1(iatom) = contact_env%atom_map_cell1(iatom)%iatom
     436              :       END DO
     437              : 
     438            4 :       nspins = dft_control%nspins
     439            4 :       nimages = dft_control%nimages
     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 :       NULLIFY (fm_struct)
     453            4 :       nao = number_of_atomic_orbitals(subsys, atom_list0)
     454            4 :       CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=sub_env%blacs_env)
     455              : 
     456              :       ! ++ create matrices: s_00, s_01
     457            4 :       ALLOCATE (contact_env%s_00, contact_env%s_01)
     458            4 :       CALL cp_fm_create(contact_env%s_00, fm_struct)
     459            4 :       CALL cp_fm_create(contact_env%s_01, fm_struct)
     460              : 
     461              :       ! ++ create matrices: h_00, h_01, rho_00, rho_01
     462           28 :       ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
     463           28 :       ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
     464            8 :       DO ispin = 1, nspins
     465            4 :          CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
     466            4 :          CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
     467            4 :          CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
     468            8 :          CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
     469              :       END DO
     470              : 
     471            4 :       CALL cp_fm_struct_release(fm_struct)
     472              : 
     473            4 :       NULLIFY (matrix_s_ref)
     474            4 :       CALL dbcsr_init_p(matrix_s_ref)
     475            4 :       CALL dbcsr_copy(matrix_s_ref, matrix_s_kp(1, 1)%matrix)
     476            4 :       CALL dbcsr_set(matrix_s_ref, 0.0_dp)
     477              : 
     478           16 :       ALLOCATE (is_same_cell(natoms, natoms))
     479              : 
     480              :       CALL negf_reference_contact_matrix(matrix_contact=matrix_s_ref, &
     481              :                                          matrix_device=matrix_s_device, &
     482              :                                          atom_list=contact_env%atomlist_cell0, &
     483              :                                          atom_map=contact_env%atom_map_cell0, &
     484            4 :                                          para_env=para_env)
     485              : 
     486              :       ! extract matrices: s_00, s_01
     487              :       CALL negf_copy_contact_matrix(fm_cell0=contact_env%s_00, &
     488              :                                     fm_cell1=contact_env%s_01, &
     489              :                                     direction_axis=contact_env%direction_axis, &
     490              :                                     matrix_kp=matrix_s_kp(1, :), &
     491              :                                     index_to_cell=index_to_cell, &
     492              :                                     atom_list0=atom_list0, atom_list1=atom_list1, &
     493              :                                     subsys=subsys, mpi_comm_global=para_env, &
     494            4 :                                     is_same_cell=is_same_cell, matrix_ref=matrix_s_ref)
     495              : 
     496              :       ! extract matrices: h_00, h_01, rho_00, rho_01
     497            8 :       DO ispin = 1, nspins
     498              :          CALL negf_copy_contact_matrix(fm_cell0=contact_env%h_00(ispin), &
     499              :                                        fm_cell1=contact_env%h_01(ispin), &
     500              :                                        direction_axis=contact_env%direction_axis, &
     501              :                                        matrix_kp=matrix_ks_kp(ispin, :), &
     502              :                                        index_to_cell=index_to_cell, &
     503              :                                        atom_list0=atom_list0, atom_list1=atom_list1, &
     504              :                                        subsys=subsys, mpi_comm_global=para_env, &
     505            4 :                                        is_same_cell=is_same_cell)
     506              : 
     507              :          CALL negf_copy_contact_matrix(fm_cell0=contact_env%rho_00(ispin), &
     508              :                                        fm_cell1=contact_env%rho_01(ispin), &
     509              :                                        direction_axis=contact_env%direction_axis, &
     510              :                                        matrix_kp=rho_ao_kp(ispin, :), &
     511              :                                        index_to_cell=index_to_cell, &
     512              :                                        atom_list0=atom_list0, atom_list1=atom_list1, &
     513              :                                        subsys=subsys, mpi_comm_global=para_env, &
     514            8 :                                        is_same_cell=is_same_cell)
     515              :       END DO
     516              : 
     517            4 :       DEALLOCATE (is_same_cell)
     518            4 :       CALL dbcsr_deallocate_matrix(matrix_s_ref)
     519            4 :       DEALLOCATE (index_to_cell)
     520            4 :       DEALLOCATE (atom_list0, atom_list1)
     521              : 
     522            4 :       CALL timestop(handle)
     523            4 :    END SUBROUTINE negf_env_contact_init_matrices
     524              : 
     525              : ! **************************************************************************************************
     526              : !> \brief Extract relevant matrix blocks for the given contact using the device's force environment.
     527              : !> \param contact_env         NEGF environment for the contact (modified on exit)
     528              : !> \param contact_control     NEGF control for the contact
     529              : !> \param sub_env             NEGF parallel (sub)group environment
     530              : !> \param qs_env              QuickStep environment for the device force environment
     531              : !> \param eps_geometry        accuracy in Cartesian coordinates
     532              : !> \author Sergey Chulkov
     533              : ! **************************************************************************************************
     534            4 :    SUBROUTINE negf_env_contact_init_matrices_gamma(contact_env, contact_control, sub_env, qs_env, eps_geometry)
     535              :       TYPE(negf_env_contact_type), INTENT(inout)         :: contact_env
     536              :       TYPE(negf_control_contact_type), INTENT(in)        :: contact_control
     537              :       TYPE(negf_subgroup_env_type), INTENT(in)           :: sub_env
     538              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     539              :       REAL(kind=dp), INTENT(in)                          :: eps_geometry
     540              : 
     541              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_env_contact_init_matrices_gamma'
     542              : 
     543              :       INTEGER                                            :: handle, iatom, icell, ispin, nao_c, &
     544              :                                                             nspins
     545              :       LOGICAL                                            :: do_kpoints
     546              :       REAL(kind=dp), DIMENSION(2)                        :: r2_origin_cell
     547              :       REAL(kind=dp), DIMENSION(3)                        :: direction_vector, origin
     548              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     549            4 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
     550              :       TYPE(dft_control_type), POINTER                    :: dft_control
     551              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     552            4 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     553              :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     554              :       TYPE(qs_subsys_type), POINTER                      :: subsys
     555              : 
     556            4 :       CALL timeset(routineN, handle)
     557              : 
     558              :       CALL get_qs_env(qs_env, &
     559              :                       dft_control=dft_control, &
     560              :                       do_kpoints=do_kpoints, &
     561              :                       matrix_ks_kp=matrix_ks_kp, &
     562              :                       matrix_s_kp=matrix_s_kp, &
     563              :                       para_env=para_env, &
     564              :                       rho=rho_struct, &
     565            4 :                       subsys=subsys)
     566            4 :       CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
     567              : 
     568            4 :       IF (do_kpoints) THEN
     569              :          CALL cp_abort(__LOCATION__, &
     570            0 :                        "K-points in device region have not been implemented yet.")
     571              :       END IF
     572              : 
     573            4 :       nspins = dft_control%nspins
     574              : 
     575            4 :       nao_c = number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(1)%vector)
     576            4 :       IF (number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(2)%vector) /= nao_c) THEN
     577              :          CALL cp_abort(__LOCATION__, &
     578              :                        "Primary and secondary bulk contact cells should be identical "// &
     579              :                        "in terms of the number of atoms of each kind, and their basis sets. "// &
     580            0 :                        "No single atom, however, can be shared between these two cells.")
     581              :       END IF
     582              : 
     583            4 :       contact_env%homo_energy = 0.0_dp
     584              : 
     585              :       CALL contact_direction_vector(contact_env%origin, &
     586              :                                     contact_env%direction_vector, &
     587              :                                     contact_env%origin_bias, &
     588              :                                     contact_env%direction_vector_bias, &
     589              :                                     contact_control%atomlist_screening, &
     590              :                                     contact_control%atomlist_bulk, &
     591            4 :                                     subsys)
     592              : 
     593            4 :       contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys, eps_geometry)
     594              : 
     595              :       ! choose the primary and secondary contact unit cells
     596            4 :       CALL qs_subsys_get(subsys, particle_set=particle_set)
     597              : 
     598           16 :       origin = particle_set(contact_control%atomlist_screening(1))%r
     599           16 :       DO iatom = 2, SIZE(contact_control%atomlist_screening)
     600           52 :          origin = origin + particle_set(contact_control%atomlist_screening(iatom))%r
     601              :       END DO
     602           16 :       origin = origin/REAL(SIZE(contact_control%atomlist_screening), kind=dp)
     603              : 
     604           12 :       DO icell = 1, 2
     605           32 :          direction_vector = particle_set(contact_control%atomlist_cell(icell)%vector(1))%r
     606           32 :          DO iatom = 2, SIZE(contact_control%atomlist_cell(icell)%vector)
     607          104 :             direction_vector = direction_vector + particle_set(contact_control%atomlist_cell(icell)%vector(iatom))%r
     608              :          END DO
     609           32 :          direction_vector = direction_vector/REAL(SIZE(contact_control%atomlist_cell(icell)%vector), kind=dp)
     610           32 :          direction_vector = direction_vector - origin
     611           36 :          r2_origin_cell(icell) = DOT_PRODUCT(direction_vector, direction_vector)
     612              :       END DO
     613              : 
     614            4 :       IF (ABS(r2_origin_cell(1) - r2_origin_cell(2)) < (eps_geometry*eps_geometry)) THEN
     615              :          ! primary and secondary bulk unit cells should not overlap;
     616              :          ! currently we check that they are different by at least one atom that is, indeed, not sufficient.
     617              :          CALL cp_abort(__LOCATION__, &
     618            0 :                        "Primary and secondary bulk contact cells should not overlap ")
     619            4 :       ELSE IF (r2_origin_cell(1) < r2_origin_cell(2)) THEN
     620           12 :          ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(1)%vector)))
     621           20 :          contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(1)%vector(:)
     622           12 :          ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(2)%vector)))
     623           20 :          contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(2)%vector(:)
     624              :       ELSE
     625            0 :          ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(2)%vector)))
     626            0 :          contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(2)%vector(:)
     627            0 :          ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(1)%vector)))
     628            0 :          contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(1)%vector(:)
     629              :       END IF
     630              : 
     631            4 :       NULLIFY (fm_struct)
     632            4 :       CALL cp_fm_struct_create(fm_struct, nrow_global=nao_c, ncol_global=nao_c, context=sub_env%blacs_env)
     633           28 :       ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
     634           28 :       ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
     635            8 :       DO ispin = 1, nspins
     636            4 :          CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
     637            4 :          CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
     638            4 :          CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
     639            8 :          CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
     640              :       END DO
     641            4 :       ALLOCATE (contact_env%s_00, contact_env%s_01)
     642            4 :       CALL cp_fm_create(contact_env%s_00, fm_struct)
     643            4 :       CALL cp_fm_create(contact_env%s_01, fm_struct)
     644            4 :       CALL cp_fm_struct_release(fm_struct)
     645              : 
     646            8 :       DO ispin = 1, nspins
     647              :          CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
     648              :                                                fm=contact_env%h_00(ispin), &
     649              :                                                atomlist_row=contact_env%atomlist_cell0, &
     650              :                                                atomlist_col=contact_env%atomlist_cell0, &
     651              :                                                subsys=subsys, mpi_comm_global=para_env, &
     652            4 :                                                do_upper_diag=.TRUE., do_lower=.TRUE.)
     653              :          CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
     654              :                                                fm=contact_env%h_01(ispin), &
     655              :                                                atomlist_row=contact_env%atomlist_cell0, &
     656              :                                                atomlist_col=contact_env%atomlist_cell1, &
     657              :                                                subsys=subsys, mpi_comm_global=para_env, &
     658            4 :                                                do_upper_diag=.TRUE., do_lower=.TRUE.)
     659              : 
     660              :          CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
     661              :                                                fm=contact_env%rho_00(ispin), &
     662              :                                                atomlist_row=contact_env%atomlist_cell0, &
     663              :                                                atomlist_col=contact_env%atomlist_cell0, &
     664              :                                                subsys=subsys, mpi_comm_global=para_env, &
     665            4 :                                                do_upper_diag=.TRUE., do_lower=.TRUE.)
     666              :          CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
     667              :                                                fm=contact_env%rho_01(ispin), &
     668              :                                                atomlist_row=contact_env%atomlist_cell0, &
     669              :                                                atomlist_col=contact_env%atomlist_cell1, &
     670              :                                                subsys=subsys, mpi_comm_global=para_env, &
     671            8 :                                                do_upper_diag=.TRUE., do_lower=.TRUE.)
     672              :       END DO
     673              : 
     674              :       CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
     675              :                                             fm=contact_env%s_00, &
     676              :                                             atomlist_row=contact_env%atomlist_cell0, &
     677              :                                             atomlist_col=contact_env%atomlist_cell0, &
     678              :                                             subsys=subsys, mpi_comm_global=para_env, &
     679            4 :                                             do_upper_diag=.TRUE., do_lower=.TRUE.)
     680              :       CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
     681              :                                             fm=contact_env%s_01, &
     682              :                                             atomlist_row=contact_env%atomlist_cell0, &
     683              :                                             atomlist_col=contact_env%atomlist_cell1, &
     684              :                                             subsys=subsys, mpi_comm_global=para_env, &
     685            4 :                                             do_upper_diag=.TRUE., do_lower=.TRUE.)
     686              : 
     687            4 :       CALL timestop(handle)
     688            4 :    END SUBROUTINE negf_env_contact_init_matrices_gamma
     689              : 
     690              : ! **************************************************************************************************
     691              : !> \brief Extract relevant matrix blocks for the scattering region as well as
     692              : !>        all the scattering -- contact interface regions.
     693              : !> \param negf_env            NEGF environment (modified on exit)
     694              : !> \param negf_control        NEGF control
     695              : !> \param sub_env             NEGF parallel (sub)group environment
     696              : !> \param qs_env              Primary QuickStep environment
     697              : !> \author Sergey Chulkov
     698              : ! **************************************************************************************************
     699            4 :    SUBROUTINE negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
     700              :       TYPE(negf_env_type), INTENT(inout)                 :: negf_env
     701              :       TYPE(negf_control_type), POINTER                   :: negf_control
     702              :       TYPE(negf_subgroup_env_type), INTENT(in)           :: sub_env
     703              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     704              : 
     705              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_env_device_init_matrices'
     706              : 
     707              :       INTEGER                                            :: handle, icontact, ispin, nao_c, nao_s, &
     708              :                                                             ncontacts, nspins
     709              :       LOGICAL                                            :: do_kpoints
     710              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     711              :       TYPE(dbcsr_p_type)                                 :: hmat
     712            4 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_kp, matrix_s_kp
     713              :       TYPE(dft_control_type), POINTER                    :: dft_control
     714              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     715              :       TYPE(pw_env_type), POINTER                         :: pw_env
     716              :       TYPE(pw_pool_type), POINTER                        :: pw_pool
     717              :       TYPE(pw_r3d_rs_type)                               :: v_hartree
     718              :       TYPE(qs_subsys_type), POINTER                      :: subsys
     719              : 
     720            4 :       CALL timeset(routineN, handle)
     721              : 
     722            4 :       IF (ALLOCATED(negf_control%atomlist_S_screening)) THEN
     723              :          CALL get_qs_env(qs_env, &
     724              :                          dft_control=dft_control, &
     725              :                          do_kpoints=do_kpoints, &
     726              :                          matrix_ks_kp=matrix_ks_kp, &
     727              :                          matrix_s_kp=matrix_s_kp, &
     728              :                          para_env=para_env, &
     729              :                          pw_env=pw_env, &
     730            4 :                          subsys=subsys)
     731            4 :          CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool)
     732              : 
     733            4 :          IF (do_kpoints) THEN
     734              :             CALL cp_abort(__LOCATION__, &
     735            0 :                           "K-points in device region have not been implemented yet.")
     736              :          END IF
     737              : 
     738            4 :          ncontacts = SIZE(negf_control%contacts)
     739            4 :          nspins = dft_control%nspins
     740              : 
     741            4 :          NULLIFY (fm_struct)
     742            4 :          nao_s = number_of_atomic_orbitals(subsys, negf_control%atomlist_S_screening)
     743              : 
     744              :          ! ++ create matrices: h_s, s_s
     745            4 :          NULLIFY (negf_env%s_s, negf_env%v_hartree_s, fm_struct)
     746           16 :          ALLOCATE (negf_env%h_s(nspins))
     747              : 
     748            4 :          CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_s, context=sub_env%blacs_env)
     749            4 :          ALLOCATE (negf_env%s_s)
     750            4 :          CALL cp_fm_create(negf_env%s_s, fm_struct)
     751            8 :          DO ispin = 1, nspins
     752            8 :             CALL cp_fm_create(negf_env%h_s(ispin), fm_struct)
     753              :          END DO
     754            4 :          ALLOCATE (negf_env%v_hartree_s)
     755            4 :          CALL cp_fm_create(negf_env%v_hartree_s, fm_struct)
     756            4 :          CALL cp_fm_struct_release(fm_struct)
     757              : 
     758              :          ! ++ create matrices: h_sc, s_sc
     759           48 :          ALLOCATE (negf_env%h_sc(nspins, ncontacts), negf_env%s_sc(ncontacts))
     760           12 :          DO icontact = 1, ncontacts
     761            8 :             nao_c = number_of_atomic_orbitals(subsys, negf_env%contacts(icontact)%atomlist_cell0)
     762            8 :             CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_c, context=sub_env%blacs_env)
     763              : 
     764            8 :             CALL cp_fm_create(negf_env%s_sc(icontact), fm_struct)
     765              : 
     766           16 :             DO ispin = 1, nspins
     767           16 :                CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
     768              :             END DO
     769              : 
     770           12 :             CALL cp_fm_struct_release(fm_struct)
     771              :          END DO
     772              : 
     773              :          ! extract matrices: h_s, s_s
     774            8 :          DO ispin = 1, nspins
     775              :             CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
     776              :                                                   fm=negf_env%h_s(ispin), &
     777              :                                                   atomlist_row=negf_control%atomlist_S_screening, &
     778              :                                                   atomlist_col=negf_control%atomlist_S_screening, &
     779              :                                                   subsys=subsys, mpi_comm_global=para_env, &
     780            8 :                                                   do_upper_diag=.TRUE., do_lower=.TRUE.)
     781              :          END DO
     782              : 
     783              :          CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
     784              :                                                fm=negf_env%s_s, &
     785              :                                                atomlist_row=negf_control%atomlist_S_screening, &
     786              :                                                atomlist_col=negf_control%atomlist_S_screening, &
     787              :                                                subsys=subsys, mpi_comm_global=para_env, &
     788            4 :                                                do_upper_diag=.TRUE., do_lower=.TRUE.)
     789              : 
     790              :          ! v_hartree_s
     791            4 :          NULLIFY (hmat%matrix)
     792            4 :          CALL dbcsr_init_p(hmat%matrix)
     793            4 :          CALL dbcsr_copy(matrix_b=hmat%matrix, matrix_a=matrix_s_kp(1, 1)%matrix)
     794            4 :          CALL dbcsr_set(hmat%matrix, 0.0_dp)
     795              : 
     796            4 :          CALL pw_pool%create_pw(v_hartree)
     797            4 :          CALL negf_env_init_v_hartree(v_hartree, negf_env%contacts, negf_control%contacts)
     798              : 
     799              :          CALL integrate_v_rspace(v_rspace=v_hartree, hmat=hmat, qs_env=qs_env, &
     800            4 :                                  calculate_forces=.FALSE., compute_tau=.FALSE., gapw=.FALSE.)
     801              : 
     802            4 :          CALL pw_pool%give_back_pw(v_hartree)
     803              : 
     804              :          CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=hmat%matrix, &
     805              :                                                fm=negf_env%v_hartree_s, &
     806              :                                                atomlist_row=negf_control%atomlist_S_screening, &
     807              :                                                atomlist_col=negf_control%atomlist_S_screening, &
     808              :                                                subsys=subsys, mpi_comm_global=para_env, &
     809            4 :                                                do_upper_diag=.TRUE., do_lower=.TRUE.)
     810              : 
     811            4 :          CALL dbcsr_deallocate_matrix(hmat%matrix)
     812              : 
     813              :          ! extract matrices: h_sc, s_sc
     814           12 :          DO icontact = 1, ncontacts
     815           16 :             DO ispin = 1, nspins
     816              :                CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
     817              :                                                      fm=negf_env%h_sc(ispin, icontact), &
     818              :                                                      atomlist_row=negf_control%atomlist_S_screening, &
     819              :                                                      atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
     820              :                                                      subsys=subsys, mpi_comm_global=para_env, &
     821           16 :                                                      do_upper_diag=.TRUE., do_lower=.TRUE.)
     822              :             END DO
     823              : 
     824              :             CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
     825              :                                                   fm=negf_env%s_sc(icontact), &
     826              :                                                   atomlist_row=negf_control%atomlist_S_screening, &
     827              :                                                   atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
     828              :                                                   subsys=subsys, mpi_comm_global=para_env, &
     829           12 :                                                   do_upper_diag=.TRUE., do_lower=.TRUE.)
     830              :          END DO
     831              :       END IF
     832              : 
     833            4 :       CALL timestop(handle)
     834            4 :    END SUBROUTINE negf_env_device_init_matrices
     835              : 
     836              : ! **************************************************************************************************
     837              : !> \brief Contribution to the Hartree potential related to the external bias voltage.
     838              : !> \param v_hartree        Hartree potential (modified on exit)
     839              : !> \param contact_env      NEGF environment for every contact
     840              : !> \param contact_control  NEGF control for every contact
     841              : !> \author Sergey Chulkov
     842              : ! **************************************************************************************************
     843            4 :    SUBROUTINE negf_env_init_v_hartree(v_hartree, contact_env, contact_control)
     844              :       TYPE(pw_r3d_rs_type), INTENT(IN)                   :: v_hartree
     845              :       TYPE(negf_env_contact_type), DIMENSION(:), &
     846              :          INTENT(in)                                      :: contact_env
     847              :       TYPE(negf_control_contact_type), DIMENSION(:), &
     848              :          INTENT(in)                                      :: contact_control
     849              : 
     850              :       CHARACTER(len=*), PARAMETER :: routineN = 'negf_env_init_v_hartree'
     851              :       REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*EPSILON(0.0_dp)
     852              : 
     853              :       INTEGER                                            :: dx, dy, dz, handle, icontact, ix, iy, &
     854              :                                                             iz, lx, ly, lz, ncontacts, ux, uy, uz
     855              :       REAL(kind=dp)                                      :: dvol, pot, proj, v1, v2
     856              :       REAL(kind=dp), DIMENSION(3)                        :: dirvector_bias, point_coord, &
     857              :                                                             point_indices, vector
     858              : 
     859            4 :       CALL timeset(routineN, handle)
     860              : 
     861            4 :       ncontacts = SIZE(contact_env)
     862            4 :       CPASSERT(SIZE(contact_control) == ncontacts)
     863            4 :       CPASSERT(ncontacts == 2)
     864              : 
     865           16 :       dirvector_bias = contact_env(2)%origin_bias - contact_env(1)%origin_bias
     866            4 :       v1 = contact_control(1)%v_external
     867            4 :       v2 = contact_control(2)%v_external
     868              : 
     869            4 :       lx = v_hartree%pw_grid%bounds_local(1, 1)
     870            4 :       ux = v_hartree%pw_grid%bounds_local(2, 1)
     871            4 :       ly = v_hartree%pw_grid%bounds_local(1, 2)
     872            4 :       uy = v_hartree%pw_grid%bounds_local(2, 2)
     873            4 :       lz = v_hartree%pw_grid%bounds_local(1, 3)
     874            4 :       uz = v_hartree%pw_grid%bounds_local(2, 3)
     875              : 
     876            4 :       dx = v_hartree%pw_grid%npts(1)/2
     877            4 :       dy = v_hartree%pw_grid%npts(2)/2
     878            4 :       dz = v_hartree%pw_grid%npts(3)/2
     879              : 
     880            4 :       dvol = v_hartree%pw_grid%dvol
     881              : 
     882         1936 :       DO iz = lz, uz
     883         1932 :          point_indices(3) = REAL(iz + dz, kind=dp)
     884        59140 :          DO iy = ly, uy
     885        57204 :             point_indices(2) = REAL(iy + dy, kind=dp)
     886              : 
     887       912030 :             DO ix = lx, ux
     888       852894 :                point_indices(1) = REAL(ix + dx, kind=dp)
     889     11087622 :                point_coord(:) = MATMUL(v_hartree%pw_grid%dh, point_indices)
     890              : 
     891      3411576 :                vector = point_coord - contact_env(1)%origin_bias
     892       852894 :                proj = projection_on_direction_vector(vector, dirvector_bias)
     893       852894 :                IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
     894              :                   ! scattering region
     895              :                   ! proj == 0   we are at the first contact boundary
     896              :                   ! proj == 1   we are at the second contact boundary
     897       336454 :                   IF (proj < 0.0_dp) THEN
     898              :                      proj = 0.0_dp
     899       334701 :                   ELSE IF (proj > 1.0_dp) THEN
     900            0 :                      proj = 1.0_dp
     901              :                   END IF
     902       336454 :                   pot = v1 + (v2 - v1)*proj
     903              :                ELSE
     904       818268 :                   pot = 0.0_dp
     905       818268 :                   DO icontact = 1, ncontacts
     906      3156784 :                      vector = point_coord - contact_env(icontact)%origin_bias
     907       789196 :                      proj = projection_on_direction_vector(vector, contact_env(icontact)%direction_vector_bias)
     908              : 
     909       818268 :                      IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
     910       487368 :                         pot = contact_control(icontact)%v_external
     911       487368 :                         EXIT
     912              :                      END IF
     913              :                   END DO
     914              :                END IF
     915              : 
     916       910098 :                v_hartree%array(ix, iy, iz) = pot*dvol
     917              :             END DO
     918              :          END DO
     919              :       END DO
     920              : 
     921            4 :       CALL timestop(handle)
     922            4 :    END SUBROUTINE negf_env_init_v_hartree
     923              : 
     924              : ! **************************************************************************************************
     925              : !> \brief Detect the axis towards secondary unit cell.
     926              : !> \param direction_vector    direction vector
     927              : !> \param subsys_contact      QuickStep subsystem of the contact force environment
     928              : !> \param eps_geometry        accuracy in mapping atoms between different force environments
     929              : !> \return direction axis: 0 (undefined), 1 (x), 2(y), 3 (z)
     930              : !> \par History
     931              : !>   * 08.2017 created [Sergey Chulkov]
     932              : ! **************************************************************************************************
     933            8 :    FUNCTION contact_direction_axis(direction_vector, subsys_contact, eps_geometry) RESULT(direction_axis)
     934              :       REAL(kind=dp), DIMENSION(3), INTENT(in)            :: direction_vector
     935              :       TYPE(qs_subsys_type), POINTER                      :: subsys_contact
     936              :       REAL(kind=dp), INTENT(in)                          :: eps_geometry
     937              :       INTEGER                                            :: direction_axis
     938              : 
     939              :       INTEGER                                            :: i, naxes
     940              :       REAL(kind=dp), DIMENSION(3)                        :: scaled
     941              :       TYPE(cell_type), POINTER                           :: cell
     942              : 
     943            8 :       CALL qs_subsys_get(subsys_contact, cell=cell)
     944            8 :       CALL real_to_scaled(scaled, direction_vector, cell)
     945              : 
     946            8 :       naxes = 0
     947            8 :       direction_axis = 0 ! initialize to make GCC<=6 happy
     948              : 
     949           32 :       DO i = 1, 3
     950           32 :          IF (ABS(scaled(i)) > eps_geometry) THEN
     951            8 :             IF (scaled(i) > 0.0_dp) THEN
     952              :                direction_axis = i
     953              :             ELSE
     954            4 :                direction_axis = -i
     955              :             END IF
     956            8 :             naxes = naxes + 1
     957              :          END IF
     958              :       END DO
     959              : 
     960              :       ! direction_vector is not parallel to one of the unit cell's axis
     961            8 :       IF (naxes /= 1) direction_axis = 0
     962            8 :    END FUNCTION contact_direction_axis
     963              : 
     964              : ! **************************************************************************************************
     965              : !> \brief Estimate energy of the highest spin-alpha occupied molecular orbital.
     966              : !> \param homo_energy  HOMO energy (initialised on exit)
     967              : !> \param qs_env       QuickStep environment
     968              : !> \par History
     969              : !>   * 01.2017 created [Sergey Chulkov]
     970              : ! **************************************************************************************************
     971            4 :    SUBROUTINE negf_homo_energy_estimate(homo_energy, qs_env)
     972              :       REAL(kind=dp), INTENT(out)                         :: homo_energy
     973              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     974              : 
     975              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_homo_energy_estimate'
     976              :       INTEGER, PARAMETER                                 :: gamma_point = 1
     977              : 
     978              :       INTEGER                                            :: handle, homo, ikpgr, ikpoint, imo, &
     979              :                                                             ispin, kplocal, nmo, nspins
     980              :       INTEGER, DIMENSION(2)                              :: kp_range
     981              :       LOGICAL                                            :: do_kpoints
     982              :       REAL(kind=dp)                                      :: my_homo_energy
     983            4 :       REAL(kind=dp), DIMENSION(:), POINTER               :: eigenvalues
     984            4 :       TYPE(kpoint_env_p_type), DIMENSION(:), POINTER     :: kp_env
     985              :       TYPE(kpoint_type), POINTER                         :: kpoints
     986            4 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     987            4 :       TYPE(mo_set_type), DIMENSION(:, :), POINTER        :: mos_kp
     988              :       TYPE(mp_para_env_type), POINTER                    :: para_env, para_env_kp
     989              : 
     990            4 :       CALL timeset(routineN, handle)
     991            4 :       my_homo_energy = 0.0_dp
     992              : 
     993            4 :       CALL get_qs_env(qs_env, para_env=para_env, mos=mos, kpoints=kpoints, do_kpoints=do_kpoints)
     994              : 
     995            4 :       IF (do_kpoints) THEN
     996            4 :          CALL get_kpoint_info(kpoints, kp_env=kp_env, kp_range=kp_range, para_env_kp=para_env_kp)
     997              : 
     998              :          ! looking for a processor that holds the gamma point
     999            4 :          IF (para_env_kp%mepos == 0 .AND. kp_range(1) <= gamma_point .AND. kp_range(2) >= gamma_point) THEN
    1000            2 :             kplocal = kp_range(2) - kp_range(1) + 1
    1001              : 
    1002            2 :             DO ikpgr = 1, kplocal
    1003            2 :                CALL get_kpoint_env(kp_env(ikpgr)%kpoint_env, nkpoint=ikpoint, mos=mos_kp)
    1004              : 
    1005            2 :                IF (ikpoint == gamma_point) THEN
    1006              :                   ! mos_kp(component, spin), where component = 1 (real), or 2 (imaginary)
    1007            2 :                   CALL get_mo_set(mos_kp(1, 1), homo=homo, eigenvalues=eigenvalues) ! mu=fermi_level
    1008              : 
    1009            2 :                   my_homo_energy = eigenvalues(homo)
    1010            2 :                   EXIT
    1011              :                END IF
    1012              :             END DO
    1013              :          END IF
    1014              : 
    1015            4 :          CALL para_env%sum(my_homo_energy)
    1016              :       ELSE
    1017              :          ! Hamiltonian of the bulk contact region has been computed without k-points.
    1018              :          ! Try to obtain the HOMO energy assuming there is no OT. We probably should abort here
    1019              :          ! as we do need a second replica of the bulk contact unit cell along transport
    1020              :          ! direction anyway which is not available without k-points.
    1021              : 
    1022            0 :          CPWARN("It is strongly advised to use k-points within all contact FORCE_EVAL-s")
    1023              : 
    1024            0 :          nspins = SIZE(mos)
    1025              : 
    1026            0 :          spin_loop: DO ispin = 1, nspins
    1027            0 :             CALL get_mo_set(mos(ispin), homo=homo, nmo=nmo, eigenvalues=eigenvalues)
    1028              : 
    1029            0 :             DO imo = nmo, 1, -1
    1030            0 :                IF (eigenvalues(imo) /= 0.0_dp) EXIT spin_loop
    1031              :             END DO
    1032              :          END DO spin_loop
    1033              : 
    1034            0 :          IF (imo == 0) THEN
    1035            0 :             CPABORT("Orbital transformation (OT) for contact FORCE_EVAL-s is not supported")
    1036              :          END IF
    1037              : 
    1038            0 :          my_homo_energy = eigenvalues(homo)
    1039              :       END IF
    1040              : 
    1041            4 :       homo_energy = my_homo_energy
    1042            4 :       CALL timestop(handle)
    1043            4 :    END SUBROUTINE negf_homo_energy_estimate
    1044              : 
    1045              : ! **************************************************************************************************
    1046              : !> \brief List atoms from the contact's primary unit cell.
    1047              : !> \param atomlist_cell0    list of atoms belonging to the contact's primary unit cell
    1048              : !>                          (allocate and initialised on exit)
    1049              : !> \param atom_map_cell0    atomic map of atoms from 'atomlist_cell0' (allocate and initialised on exit)
    1050              : !> \param atomlist_bulk     list of atoms belonging to the bulk contact region
    1051              : !> \param atom_map          atomic map of atoms from 'atomlist_bulk'
    1052              : !> \param origin            origin of the contact
    1053              : !> \param direction_vector  direction vector of the contact
    1054              : !> \param direction_axis    axis towards secondary unit cell
    1055              : !> \param subsys_device     QuickStep subsystem of the device force environment
    1056              : !> \par History
    1057              : !>   * 08.2017 created [Sergey Chulkov]
    1058              : ! **************************************************************************************************
    1059            4 :    SUBROUTINE list_atoms_in_bulk_primary_unit_cell(atomlist_cell0, atom_map_cell0, atomlist_bulk, atom_map, &
    1060              :                                                    origin, direction_vector, direction_axis, subsys_device)
    1061              :       INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout)  :: atomlist_cell0
    1062              :       TYPE(negf_atom_map_type), ALLOCATABLE, &
    1063              :          DIMENSION(:), INTENT(inout)                     :: atom_map_cell0
    1064              :       INTEGER, DIMENSION(:), INTENT(in)                  :: atomlist_bulk
    1065              :       TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
    1066              :       REAL(kind=dp), DIMENSION(3), INTENT(in)            :: origin, direction_vector
    1067              :       INTEGER, INTENT(in)                                :: direction_axis
    1068              :       TYPE(qs_subsys_type), POINTER                      :: subsys_device
    1069              : 
    1070              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'list_atoms_in_bulk_primary_unit_cell'
    1071              : 
    1072              :       INTEGER                                            :: atom_min, dir_axis_min, &
    1073              :                                                             direction_axis_abs, handle, iatom, &
    1074              :                                                             natoms_bulk, natoms_cell0
    1075              :       REAL(kind=dp)                                      :: proj, proj_min
    1076              :       REAL(kind=dp), DIMENSION(3)                        :: vector
    1077            4 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1078              : 
    1079            4 :       CALL timeset(routineN, handle)
    1080            4 :       CALL qs_subsys_get(subsys_device, particle_set=particle_set)
    1081              : 
    1082            4 :       natoms_bulk = SIZE(atomlist_bulk)
    1083            4 :       CPASSERT(SIZE(atom_map, 1) == natoms_bulk)
    1084            4 :       direction_axis_abs = ABS(direction_axis)
    1085              : 
    1086              :       ! looking for the nearest atom from the scattering region
    1087            4 :       proj_min = 1.0_dp
    1088            4 :       atom_min = 1
    1089           36 :       DO iatom = 1, natoms_bulk
    1090          128 :          vector = particle_set(atomlist_bulk(iatom))%r - origin
    1091           32 :          proj = projection_on_direction_vector(vector, direction_vector)
    1092              : 
    1093           36 :          IF (proj < proj_min) THEN
    1094           16 :             proj_min = proj
    1095           16 :             atom_min = iatom
    1096              :          END IF
    1097              :       END DO
    1098              : 
    1099            4 :       dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
    1100              : 
    1101            4 :       natoms_cell0 = 0
    1102           36 :       DO iatom = 1, natoms_bulk
    1103           32 :          IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) &
    1104           20 :             natoms_cell0 = natoms_cell0 + 1
    1105              :       END DO
    1106              : 
    1107           12 :       ALLOCATE (atomlist_cell0(natoms_cell0))
    1108           40 :       ALLOCATE (atom_map_cell0(natoms_cell0))
    1109              : 
    1110            4 :       natoms_cell0 = 0
    1111           36 :       DO iatom = 1, natoms_bulk
    1112           36 :          IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) THEN
    1113           16 :             natoms_cell0 = natoms_cell0 + 1
    1114           16 :             atomlist_cell0(natoms_cell0) = atomlist_bulk(iatom)
    1115           16 :             atom_map_cell0(natoms_cell0) = atom_map(iatom)
    1116              :          END IF
    1117              :       END DO
    1118              : 
    1119            4 :       CALL timestop(handle)
    1120            4 :    END SUBROUTINE list_atoms_in_bulk_primary_unit_cell
    1121              : 
    1122              : ! **************************************************************************************************
    1123              : !> \brief List atoms from the contact's secondary unit cell.
    1124              : !> \param atomlist_cell1    list of atoms belonging to the contact's secondary unit cell
    1125              : !>                          (allocate and initialised on exit)
    1126              : !> \param atom_map_cell1    atomic map of atoms from 'atomlist_cell1'
    1127              : !>                          (allocate and initialised on exit)
    1128              : !> \param atomlist_bulk     list of atoms belonging to the bulk contact region
    1129              : !> \param atom_map          atomic map of atoms from 'atomlist_bulk'
    1130              : !> \param origin            origin of the contact
    1131              : !> \param direction_vector  direction vector of the contact
    1132              : !> \param direction_axis    axis towards the secondary unit cell
    1133              : !> \param subsys_device     QuickStep subsystem of the device force environment
    1134              : !> \par History
    1135              : !>   * 11.2017 created [Sergey Chulkov]
    1136              : !> \note Cloned from list_atoms_in_bulk_primary_unit_cell. Will be removed once we can managed to
    1137              : !>        maintain consistency between real-space matrices from different force_eval sections.
    1138              : ! **************************************************************************************************
    1139            4 :    SUBROUTINE list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1, atom_map_cell1, atomlist_bulk, atom_map, &
    1140              :                                                      origin, direction_vector, direction_axis, subsys_device)
    1141              :       INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout)  :: atomlist_cell1
    1142              :       TYPE(negf_atom_map_type), ALLOCATABLE, &
    1143              :          DIMENSION(:), INTENT(inout)                     :: atom_map_cell1
    1144              :       INTEGER, DIMENSION(:), INTENT(in)                  :: atomlist_bulk
    1145              :       TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
    1146              :       REAL(kind=dp), DIMENSION(3), INTENT(in)            :: origin, direction_vector
    1147              :       INTEGER, INTENT(in)                                :: direction_axis
    1148              :       TYPE(qs_subsys_type), POINTER                      :: subsys_device
    1149              : 
    1150              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'list_atoms_in_bulk_secondary_unit_cell'
    1151              : 
    1152              :       INTEGER                                            :: atom_min, dir_axis_min, &
    1153              :                                                             direction_axis_abs, handle, iatom, &
    1154              :                                                             natoms_bulk, natoms_cell1, offset
    1155              :       REAL(kind=dp)                                      :: proj, proj_min
    1156              :       REAL(kind=dp), DIMENSION(3)                        :: vector
    1157            4 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1158              : 
    1159            4 :       CALL timeset(routineN, handle)
    1160            4 :       CALL qs_subsys_get(subsys_device, particle_set=particle_set)
    1161              : 
    1162            4 :       natoms_bulk = SIZE(atomlist_bulk)
    1163            4 :       CPASSERT(SIZE(atom_map, 1) == natoms_bulk)
    1164            4 :       direction_axis_abs = ABS(direction_axis)
    1165            4 :       offset = SIGN(1, direction_axis)
    1166              : 
    1167              :       ! looking for the nearest atom from the scattering region
    1168            4 :       proj_min = 1.0_dp
    1169            4 :       atom_min = 1
    1170           36 :       DO iatom = 1, natoms_bulk
    1171          128 :          vector = particle_set(atomlist_bulk(iatom))%r - origin
    1172           32 :          proj = projection_on_direction_vector(vector, direction_vector)
    1173              : 
    1174           36 :          IF (proj < proj_min) THEN
    1175           16 :             proj_min = proj
    1176           16 :             atom_min = iatom
    1177              :          END IF
    1178              :       END DO
    1179              : 
    1180            4 :       dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
    1181              : 
    1182            4 :       natoms_cell1 = 0
    1183           36 :       DO iatom = 1, natoms_bulk
    1184           32 :          IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) &
    1185           20 :             natoms_cell1 = natoms_cell1 + 1
    1186              :       END DO
    1187              : 
    1188           12 :       ALLOCATE (atomlist_cell1(natoms_cell1))
    1189           40 :       ALLOCATE (atom_map_cell1(natoms_cell1))
    1190              : 
    1191            4 :       natoms_cell1 = 0
    1192           36 :       DO iatom = 1, natoms_bulk
    1193           36 :          IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) THEN
    1194           16 :             natoms_cell1 = natoms_cell1 + 1
    1195           16 :             atomlist_cell1(natoms_cell1) = atomlist_bulk(iatom)
    1196           16 :             atom_map_cell1(natoms_cell1) = atom_map(iatom)
    1197           16 :             atom_map_cell1(natoms_cell1)%cell(direction_axis_abs) = dir_axis_min
    1198              :          END IF
    1199              :       END DO
    1200              : 
    1201            4 :       CALL timestop(handle)
    1202            4 :    END SUBROUTINE list_atoms_in_bulk_secondary_unit_cell
    1203              : 
    1204              : ! **************************************************************************************************
    1205              : !> \brief Release a NEGF environment variable.
    1206              : !> \param negf_env  NEGF environment to release
    1207              : !> \par History
    1208              : !>   * 01.2017 created [Sergey Chulkov]
    1209              : ! **************************************************************************************************
    1210            4 :    SUBROUTINE negf_env_release(negf_env)
    1211              :       TYPE(negf_env_type), INTENT(inout)                 :: negf_env
    1212              : 
    1213              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'negf_env_release'
    1214              : 
    1215              :       INTEGER                                            :: handle, icontact
    1216              : 
    1217            4 :       CALL timeset(routineN, handle)
    1218              : 
    1219            4 :       IF (ALLOCATED(negf_env%contacts)) THEN
    1220           12 :          DO icontact = SIZE(negf_env%contacts), 1, -1
    1221           12 :             CALL negf_env_contact_release(negf_env%contacts(icontact))
    1222              :          END DO
    1223              : 
    1224           12 :          DEALLOCATE (negf_env%contacts)
    1225              :       END IF
    1226              : 
    1227              :       ! h_s
    1228            4 :       CALL cp_fm_release(negf_env%h_s)
    1229              : 
    1230              :       ! h_sc
    1231            4 :       CALL cp_fm_release(negf_env%h_sc)
    1232              : 
    1233              :       ! s_s
    1234            4 :       IF (ASSOCIATED(negf_env%s_s)) THEN
    1235            4 :          CALL cp_fm_release(negf_env%s_s)
    1236            4 :          DEALLOCATE (negf_env%s_s)
    1237              :          NULLIFY (negf_env%s_s)
    1238              :       END IF
    1239              : 
    1240              :       ! s_sc
    1241            4 :       CALL cp_fm_release(negf_env%s_sc)
    1242              : 
    1243              :       ! v_hartree_s
    1244            4 :       IF (ASSOCIATED(negf_env%v_hartree_s)) THEN
    1245            4 :          CALL cp_fm_release(negf_env%v_hartree_s)
    1246            4 :          DEALLOCATE (negf_env%v_hartree_s)
    1247              :          NULLIFY (negf_env%v_hartree_s)
    1248              :       END IF
    1249              : 
    1250              :       ! density mixing
    1251            4 :       IF (ASSOCIATED(negf_env%mixing_storage)) THEN
    1252            4 :          CALL mixing_storage_release(negf_env%mixing_storage)
    1253            4 :          DEALLOCATE (negf_env%mixing_storage)
    1254              :       END IF
    1255              : 
    1256            4 :       CALL timestop(handle)
    1257            4 :    END SUBROUTINE negf_env_release
    1258              : 
    1259              : ! **************************************************************************************************
    1260              : !> \brief Release a NEGF contact environment variable.
    1261              : !> \param contact_env  NEGF contact environment to release
    1262              : ! **************************************************************************************************
    1263            8 :    SUBROUTINE negf_env_contact_release(contact_env)
    1264              :       TYPE(negf_env_contact_type), INTENT(inout)         :: contact_env
    1265              : 
    1266              :       CHARACTER(len=*), PARAMETER :: routineN = 'negf_env_contact_release'
    1267              : 
    1268              :       INTEGER                                            :: handle
    1269              : 
    1270            8 :       CALL timeset(routineN, handle)
    1271              : 
    1272              :       ! h_00
    1273            8 :       CALL cp_fm_release(contact_env%h_00)
    1274              : 
    1275              :       ! h_01
    1276            8 :       CALL cp_fm_release(contact_env%h_01)
    1277              : 
    1278              :       ! rho_00
    1279            8 :       CALL cp_fm_release(contact_env%rho_00)
    1280              : 
    1281              :       ! rho_01
    1282            8 :       CALL cp_fm_release(contact_env%rho_01)
    1283              : 
    1284              :       ! s_00
    1285            8 :       IF (ASSOCIATED(contact_env%s_00)) THEN
    1286            8 :          CALL cp_fm_release(contact_env%s_00)
    1287            8 :          DEALLOCATE (contact_env%s_00)
    1288              :          NULLIFY (contact_env%s_00)
    1289              :       END IF
    1290              : 
    1291              :       ! s_01
    1292            8 :       IF (ASSOCIATED(contact_env%s_01)) THEN
    1293            8 :          CALL cp_fm_release(contact_env%s_01)
    1294            8 :          DEALLOCATE (contact_env%s_01)
    1295              :          NULLIFY (contact_env%s_01)
    1296              :       END IF
    1297              : 
    1298            8 :       IF (ALLOCATED(contact_env%atomlist_cell0)) DEALLOCATE (contact_env%atomlist_cell0)
    1299            8 :       IF (ALLOCATED(contact_env%atomlist_cell1)) DEALLOCATE (contact_env%atomlist_cell1)
    1300            8 :       IF (ALLOCATED(contact_env%atom_map_cell0)) DEALLOCATE (contact_env%atom_map_cell0)
    1301              : 
    1302            8 :       CALL timestop(handle)
    1303            8 :    END SUBROUTINE negf_env_contact_release
    1304            0 : END MODULE negf_env_types
        

Generated by: LCOV version 2.0-1