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