1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Input control types for NEGF based quantum transport calculations 8! ************************************************************************************************** 9 10MODULE negf_control_types 11 USE cp_subsys_types, ONLY: cp_subsys_get,& 12 cp_subsys_type 13 USE input_constants, ONLY: negf_run 14 USE input_section_types, ONLY: section_vals_get,& 15 section_vals_get_subs_vals,& 16 section_vals_type,& 17 section_vals_val_get 18 USE kinds, ONLY: default_string_length,& 19 dp 20 USE mathconstants, ONLY: pi 21 USE molecule_kind_types, ONLY: get_molecule_kind,& 22 molecule_kind_type 23 USE molecule_types, ONLY: get_molecule,& 24 molecule_type 25 USE negf_alloc_types, ONLY: negf_allocatable_ivector 26 USE particle_types, ONLY: particle_type 27 USE physcon, ONLY: kelvin 28 USE string_utilities, ONLY: integer_to_string 29 USE util, ONLY: sort 30#include "./base/base_uses.f90" 31 32 IMPLICIT NONE 33 PRIVATE 34 35 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_control_types' 36 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE. 37 38 PUBLIC :: negf_control_type, negf_control_contact_type 39 PUBLIC :: negf_control_create, negf_control_release, read_negf_control 40 41! ************************************************************************************************** 42!> \brief Input parameters related to a single contact. 43!> \author Sergey Chulkov 44! ************************************************************************************************** 45 TYPE negf_control_contact_type 46 !> atoms belonging to bulk and screening regions 47 INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_bulk, atomlist_screening 48 !> atom belonging to the primary and secondary bulk unit cells 49 TYPE(negf_allocatable_ivector), ALLOCATABLE, & 50 DIMENSION(:) :: atomlist_cell 51 !> index of the sub_force_env which should be used for bulk calculation 52 INTEGER :: force_env_index 53 !> contact Fermi level needs to be computed prior NEGF run 54 LOGICAL :: compute_fermi_level 55 !> when computing contact Fermi level, use the energy given in 'fermi_level' (insted of HOMO) 56 !> (insted of the HOMO energy ) as a starting point 57 LOGICAL :: refine_fermi_level 58 !> Fermi level 59 REAL(kind=dp) :: fermi_level 60 !> temperature [in a.u.] 61 REAL(kind=dp) :: temperature 62 !> applied electric potential 63 REAL(kind=dp) :: v_external 64 END TYPE negf_control_contact_type 65 66! ************************************************************************************************** 67!> \brief Input parameters related to the NEGF run. 68!> \author Sergey Chulkov 69! ************************************************************************************************** 70 TYPE negf_control_type 71 !> input options for every contact 72 TYPE(negf_control_contact_type), ALLOCATABLE, & 73 DIMENSION(:) :: contacts 74 !> atoms belonging to the scattering region 75 INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_S 76 !> atoms belonging to the scattering region as well as atoms belonging to 77 !> screening regions of all the contacts 78 INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_S_screening 79 !> do not keep contact self-energy matrices 80 LOGICAL :: disable_cache 81 !> convergence criteria for adaptive integration methods 82 REAL(kind=dp) :: conv_density 83 !> convergence criteria for iterative Lopez-Sancho algorithm 84 REAL(kind=dp) :: conv_green 85 !> convergence criteria for self-consistent iterations 86 REAL(kind=dp) :: conv_scf 87 !> accuracy in mapping atoms between different force environments 88 REAL(kind=dp) :: eps_geometry 89 !> applied bias [in a.u.] 90 REAL(kind=dp) :: v_bias 91 !> integration lower bound [in a.u.] 92 REAL(kind=dp) :: energy_lbound 93 !> infinitesimal offset along the imaginary axis [in a.u.] 94 REAL(kind=dp) :: eta 95 !> initial guess to determine the actual Fermi level of bulk contacts [in a.u.] 96 REAL(kind=dp) :: homo_lumo_gap 97 !> number of residuals (poles of the Fermi function) 98 INTEGER :: delta_npoles 99 !> offset along the x-axis away from the poles of the Fermi function [in units of kT] 100 INTEGER :: gamma_kT 101 !> integration method 102 INTEGER :: integr_method 103 !> minimal number of grid points along the closed contour 104 INTEGER :: integr_min_points 105 !> maximal number of grid points along the closed contour 106 INTEGER :: integr_max_points 107 !> maximal number of SCF iterations 108 INTEGER :: max_scf 109 !> minimal number of MPI processes to be used to compute Green's function per energy point 110 INTEGER :: nprocs 111 !> shift in Hartree potential [in a.u.] 112 REAL(kind=dp) :: v_shift 113 !> initial offset to determine the correct shift in Hartree potential [in a.u.] 114 REAL(kind=dp) :: v_shift_offset 115 !> maximal number of iteration to determine the shift in Hartree potential 116 INTEGER :: v_shift_maxiters 117 END TYPE negf_control_type 118 119 PRIVATE :: read_negf_atomlist 120 121CONTAINS 122 123! ************************************************************************************************** 124!> \brief allocate control options for Non-equilibrium Green's Function calculation 125!> \param negf_control an object to create 126!> \par History 127!> * 02.2017 created [Sergey Chulkov] 128! ************************************************************************************************** 129 SUBROUTINE negf_control_create(negf_control) 130 TYPE(negf_control_type), POINTER :: negf_control 131 132 CHARACTER(len=*), PARAMETER :: routineN = 'negf_control_create', & 133 routineP = moduleN//':'//routineN 134 135 INTEGER :: handle 136 137 CPASSERT(.NOT. ASSOCIATED(negf_control)) 138 CALL timeset(routineN, handle) 139 140 ALLOCATE (negf_control) 141 142 CALL timestop(handle) 143 END SUBROUTINE negf_control_create 144 145! ************************************************************************************************** 146!> \brief release memory allocated for NEGF control options 147!> \param negf_control an object to release 148!> \par History 149!> * 02.2017 created [Sergey Chulkov] 150! ************************************************************************************************** 151 SUBROUTINE negf_control_release(negf_control) 152 TYPE(negf_control_type), POINTER :: negf_control 153 154 CHARACTER(len=*), PARAMETER :: routineN = 'negf_control_release', & 155 routineP = moduleN//':'//routineN 156 157 INTEGER :: handle, i, j 158 159 CALL timeset(routineN, handle) 160 161 IF (ASSOCIATED(negf_control)) THEN 162 IF (ALLOCATED(negf_control%atomlist_S)) DEALLOCATE (negf_control%atomlist_S) 163 IF (ALLOCATED(negf_control%atomlist_S_screening)) DEALLOCATE (negf_control%atomlist_S_screening) 164 165 IF (ALLOCATED(negf_control%contacts)) THEN 166 DO i = SIZE(negf_control%contacts), 1, -1 167 IF (ALLOCATED(negf_control%contacts(i)%atomlist_bulk)) & 168 DEALLOCATE (negf_control%contacts(i)%atomlist_bulk) 169 170 IF (ALLOCATED(negf_control%contacts(i)%atomlist_screening)) & 171 DEALLOCATE (negf_control%contacts(i)%atomlist_screening) 172 173 IF (ALLOCATED(negf_control%contacts(i)%atomlist_cell)) THEN 174 DO j = SIZE(negf_control%contacts(i)%atomlist_cell), 1, -1 175 IF (ALLOCATED(negf_control%contacts(i)%atomlist_cell(j)%vector)) & 176 DEALLOCATE (negf_control%contacts(i)%atomlist_cell(j)%vector) 177 END DO 178 DEALLOCATE (negf_control%contacts(i)%atomlist_cell) 179 END IF 180 END DO 181 182 DEALLOCATE (negf_control%contacts) 183 END IF 184 185 DEALLOCATE (negf_control) 186 END IF 187 188 CALL timestop(handle) 189 END SUBROUTINE negf_control_release 190 191! ************************************************************************************************** 192!> \brief Read NEGF input parameters. 193!> \param negf_control NEGF control parameters 194!> \param input root input section 195!> \param subsys subsystem environment 196! ************************************************************************************************** 197 SUBROUTINE read_negf_control(negf_control, input, subsys) 198 TYPE(negf_control_type), POINTER :: negf_control 199 TYPE(section_vals_type), POINTER :: input 200 TYPE(cp_subsys_type), POINTER :: subsys 201 202 CHARACTER(len=*), PARAMETER :: routineN = 'read_negf_control', & 203 routineP = moduleN//':'//routineN 204 205 CHARACTER(len=default_string_length) :: contact_id_str, eta_current_str, eta_max_str, & 206 npoles_current_str, npoles_min_str, temp_current_str, temp_min_str 207 INTEGER :: delta_npoles_min, handle, i2_rep, i_rep, & 208 n2_rep, n_rep, natoms_current, & 209 natoms_total, run_type 210 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds 211 LOGICAL :: do_negf, is_explicit 212 REAL(kind=dp) :: eta_max, temp_current, temp_min 213 TYPE(section_vals_type), POINTER :: cell_section, contact_section, & 214 negf_section, region_section 215 216 CALL timeset(routineN, handle) 217 218 CALL section_vals_val_get(input, "GLOBAL%RUN_TYPE", i_val=run_type) 219 do_negf = run_type == negf_run 220 221 negf_section => section_vals_get_subs_vals(input, "NEGF") 222 223 contact_section => section_vals_get_subs_vals(negf_section, "CONTACT") 224 CALL section_vals_get(contact_section, n_repetition=n_rep, explicit=is_explicit) 225 IF ((.NOT. is_explicit) .AND. do_negf) THEN 226 CALL cp_abort(__LOCATION__, & 227 "At least one contact is needed for NEGF calculation.") 228 END IF 229 230 ALLOCATE (negf_control%contacts(n_rep)) 231 DO i_rep = 1, n_rep 232 region_section => section_vals_get_subs_vals(contact_section, "SCREENING_REGION", i_rep_section=i_rep) 233 CALL section_vals_get(region_section, explicit=is_explicit) 234 235 IF ((.NOT. is_explicit) .AND. do_negf) THEN 236 WRITE (contact_id_str, '(I11)') i_rep 237 CALL cp_abort(__LOCATION__, & 238 "The screening region must be defined for the contact "//TRIM(ADJUSTL(contact_id_str))//".") 239 END IF 240 241 IF (is_explicit) THEN 242 CALL read_negf_atomlist(negf_control%contacts(i_rep)%atomlist_screening, region_section, 1, subsys) 243 END IF 244 245 region_section => section_vals_get_subs_vals(contact_section, "BULK_REGION", i_rep_section=i_rep) 246 247 CALL section_vals_get(region_section, explicit=is_explicit) 248 249 IF ((.NOT. is_explicit) .AND. do_negf) THEN 250 WRITE (contact_id_str, '(I11)') i_rep 251 CALL cp_abort(__LOCATION__, & 252 "The bulk region must be defined for the contact "//TRIM(ADJUSTL(contact_id_str))//".") 253 END IF 254 255 IF (is_explicit) THEN 256 CALL read_negf_atomlist(negf_control%contacts(i_rep)%atomlist_bulk, region_section, 1, subsys) 257 END IF 258 259 CALL section_vals_val_get(contact_section, "FORCE_EVAL_SECTION", & 260 i_val=negf_control%contacts(i_rep)%force_env_index, & 261 i_rep_section=i_rep) 262 263 cell_section => section_vals_get_subs_vals(region_section, "CELL") 264 CALL section_vals_get(cell_section, n_repetition=n2_rep, explicit=is_explicit) 265 266 IF (((.NOT. is_explicit) .OR. n2_rep /= 2) .AND. negf_control%contacts(i_rep)%force_env_index <= 0 .AND. do_negf) THEN 267 WRITE (contact_id_str, '(I11)') i_rep 268 CALL cp_abort(__LOCATION__, & 269 "You must either provide indices of atoms belonging to two adjacent bulk unit cells "// & 270 "(BULK_REGION/CELL) for the contact, or the index of the FORCE_EVAL section (FORCE_EVAL_SECTION) "// & 271 "which will be used to construct Kohn-Sham matrix for the bulk contact "// & 272 TRIM(ADJUSTL(contact_id_str))//".") 273 END IF 274 275 IF (is_explicit .AND. n2_rep > 0) THEN 276 ALLOCATE (negf_control%contacts(i_rep)%atomlist_cell(n2_rep)) 277 278 DO i2_rep = 1, n2_rep 279 CALL read_negf_atomlist(negf_control%contacts(i_rep)%atomlist_cell(i2_rep)%vector, cell_section, i2_rep, subsys) 280 END DO 281 END IF 282 283 CALL section_vals_val_get(contact_section, "REFINE_FERMI_LEVEL", & 284 l_val=negf_control%contacts(i_rep)%refine_fermi_level, & 285 i_rep_section=i_rep) 286 287 CALL section_vals_val_get(contact_section, "FERMI_LEVEL", & 288 r_val=negf_control%contacts(i_rep)%fermi_level, & 289 i_rep_section=i_rep, explicit=is_explicit) 290 negf_control%contacts(i_rep)%compute_fermi_level = (.NOT. is_explicit) .OR. & 291 negf_control%contacts(i_rep)%refine_fermi_level 292 293 IF (do_negf .AND. negf_control%contacts(i_rep)%force_env_index <= 0 .AND. & 294 (.NOT. (is_explicit .OR. negf_control%contacts(i_rep)%refine_fermi_level))) THEN 295 WRITE (contact_id_str, '(I11)') i_rep 296 CALL cp_warn(__LOCATION__, & 297 "There is no way to reasonably guess the Fermi level for the bulk contact "// & 298 TRIM(ADJUSTL(contact_id_str))//" without running a separate bulk DFT calculation first. "// & 299 "Therefore, 0.0 Hartree will be used as an initial guess. It is strongly advised to enable "// & 300 "the REFINE_FERMI_LEVEL switch and to provide an initial guess using the FERMI_LEVEL keyword. "// & 301 "Alternatively, a bulk FORCE_EVAL_SECTION can be set up.") 302 END IF 303 304 CALL section_vals_val_get(contact_section, "TEMPERATURE", & 305 r_val=negf_control%contacts(i_rep)%temperature, & 306 i_rep_section=i_rep) 307 IF (negf_control%contacts(i_rep)%temperature <= 0.0_dp) THEN 308 CALL cp_abort(__LOCATION__, "Electronic temperature must be > 0") 309 END IF 310 311 CALL section_vals_val_get(contact_section, "ELECTRIC_POTENTIAL", & 312 r_val=negf_control%contacts(i_rep)%v_external, & 313 i_rep_section=i_rep) 314 END DO 315 316 region_section => section_vals_get_subs_vals(negf_section, "SCATTERING_REGION") 317 CALL section_vals_get(region_section, explicit=is_explicit) 318 IF (is_explicit) THEN 319 CALL read_negf_atomlist(negf_control%atomlist_S, region_section, 1, subsys) 320 END IF 321 322 CALL section_vals_val_get(negf_section, "DISABLE_CACHE", l_val=negf_control%disable_cache) 323 324 CALL section_vals_val_get(negf_section, "EPS_DENSITY", r_val=negf_control%conv_density) 325 CALL section_vals_val_get(negf_section, "EPS_GREEN", r_val=negf_control%conv_green) 326 CALL section_vals_val_get(negf_section, "EPS_SCF", r_val=negf_control%conv_scf) 327 328 CALL section_vals_val_get(negf_section, "EPS_GEO", r_val=negf_control%eps_geometry) 329 330 CALL section_vals_val_get(negf_section, "ENERGY_LBOUND", r_val=negf_control%energy_lbound) 331 CALL section_vals_val_get(negf_section, "ETA", r_val=negf_control%eta) 332 CALL section_vals_val_get(negf_section, "HOMO_LUMO_GAP", r_val=negf_control%homo_lumo_gap) 333 CALL section_vals_val_get(negf_section, "DELTA_NPOLES", i_val=negf_control%delta_npoles) 334 CALL section_vals_val_get(negf_section, "GAMMA_KT", i_val=negf_control%gamma_kT) 335 336 CALL section_vals_val_get(negf_section, "INTEGRATION_METHOD", i_val=negf_control%integr_method) 337 CALL section_vals_val_get(negf_section, "INTEGRATION_MIN_POINTS", i_val=negf_control%integr_min_points) 338 CALL section_vals_val_get(negf_section, "INTEGRATION_MAX_POINTS", i_val=negf_control%integr_max_points) 339 340 IF (negf_control%integr_max_points < negf_control%integr_min_points) & 341 negf_control%integr_max_points = negf_control%integr_min_points 342 343 CALL section_vals_val_get(negf_section, "MAX_SCF", i_val=negf_control%max_scf) 344 345 CALL section_vals_val_get(negf_section, "NPROC_POINT", i_val=negf_control%nprocs) 346 347 CALL section_vals_val_get(negf_section, "V_SHIFT", r_val=negf_control%v_shift) 348 CALL section_vals_val_get(negf_section, "V_SHIFT_OFFSET", r_val=negf_control%v_shift_offset) 349 CALL section_vals_val_get(negf_section, "V_SHIFT_MAX_ITERS", i_val=negf_control%v_shift_maxiters) 350 351 ! check consistency 352 IF (negf_control%eta < 0.0_dp) THEN 353 CALL cp_abort(__LOCATION__, "ETA must be >= 0") 354 END IF 355 356 IF (n_rep > 0) THEN 357 delta_npoles_min = NINT(0.5_dp*(negf_control%eta/(pi*MAXVAL(negf_control%contacts(:)%temperature)) + 1.0_dp)) 358 ELSE 359 delta_npoles_min = 1 360 END IF 361 362 IF (negf_control%delta_npoles < delta_npoles_min) THEN 363 IF (n_rep > 0) THEN 364 eta_max = REAL(2*negf_control%delta_npoles - 1, kind=dp)*pi*MAXVAL(negf_control%contacts(:)%temperature) 365 temp_current = MAXVAL(negf_control%contacts(:)%temperature)*kelvin 366 temp_min = negf_control%eta/(pi*REAL(2*negf_control%delta_npoles - 1, kind=dp))*kelvin 367 368 WRITE (eta_current_str, '(ES11.4E2)') negf_control%eta 369 WRITE (eta_max_str, '(ES11.4E2)') eta_max 370 WRITE (npoles_current_str, '(I11)') negf_control%delta_npoles 371 WRITE (npoles_min_str, '(I11)') delta_npoles_min 372 WRITE (temp_current_str, '(F11.3)') temp_current 373 WRITE (temp_min_str, '(F11.3)') temp_min 374 375 CALL cp_abort(__LOCATION__, & 376 "Parameter DELTA_NPOLES must be at least "//TRIM(ADJUSTL(npoles_min_str))// & 377 " (instead of "//TRIM(ADJUSTL(npoles_current_str))// & 378 ") for given TEMPERATURE ("//TRIM(ADJUSTL(temp_current_str))// & 379 " K) and ETA ("//TRIM(ADJUSTL(eta_current_str))// & 380 "). Alternatively you can increase TEMPERATURE above "//TRIM(ADJUSTL(temp_min_str))// & 381 " K, or decrease ETA below "//TRIM(ADJUSTL(eta_max_str))// & 382 ". Please keep in mind that very tight ETA may result in dramatical precision loss"// & 383 " due to inversion of ill-conditioned matrices.") 384 ELSE 385 ! no leads have been defined, so calculation will abort anyway 386 negf_control%delta_npoles = delta_npoles_min 387 END IF 388 END IF 389 390 ! expand scattering region by adding atoms from contact screening regions 391 n_rep = SIZE(negf_control%contacts) 392 IF (ALLOCATED(negf_control%atomlist_S)) THEN 393 natoms_total = SIZE(negf_control%atomlist_S) 394 ELSE 395 natoms_total = 0 396 END IF 397 398 DO i_rep = 1, n_rep 399 IF (ALLOCATED(negf_control%contacts(i_rep)%atomlist_screening)) THEN 400 IF (ALLOCATED(negf_control%contacts(i_rep)%atomlist_screening)) & 401 natoms_total = natoms_total + SIZE(negf_control%contacts(i_rep)%atomlist_screening) 402 END IF 403 END DO 404 405 IF (natoms_total > 0) THEN 406 ALLOCATE (negf_control%atomlist_S_screening(natoms_total)) 407 IF (ALLOCATED(negf_control%atomlist_S)) THEN 408 natoms_total = SIZE(negf_control%atomlist_S) 409 negf_control%atomlist_S_screening(1:natoms_total) = negf_control%atomlist_S(1:natoms_total) 410 ELSE 411 natoms_total = 0 412 END IF 413 414 DO i_rep = 1, n_rep 415 IF (ALLOCATED(negf_control%contacts(i_rep)%atomlist_screening)) THEN 416 natoms_current = SIZE(negf_control%contacts(i_rep)%atomlist_screening) 417 418 negf_control%atomlist_S_screening(natoms_total + 1:natoms_total + natoms_current) = & 419 negf_control%contacts(i_rep)%atomlist_screening(1:natoms_current) 420 421 natoms_total = natoms_total + natoms_current 422 END IF 423 END DO 424 425 ! sort and remove duplicated atoms 426 ALLOCATE (inds(natoms_total)) 427 CALL sort(negf_control%atomlist_S_screening, natoms_total, inds) 428 DEALLOCATE (inds) 429 430 natoms_current = 1 431 DO i_rep = natoms_current + 1, natoms_total 432 IF (negf_control%atomlist_S_screening(i_rep) /= negf_control%atomlist_S_screening(natoms_current)) THEN 433 natoms_current = natoms_current + 1 434 negf_control%atomlist_S_screening(natoms_current) = negf_control%atomlist_S_screening(i_rep) 435 END IF 436 END DO 437 438 IF (natoms_current < natoms_total) THEN 439 CALL MOVE_ALLOC(negf_control%atomlist_S_screening, inds) 440 441 ALLOCATE (negf_control%atomlist_S_screening(natoms_current)) 442 negf_control%atomlist_S_screening(1:natoms_current) = inds(1:natoms_current) 443 DEALLOCATE (inds) 444 END IF 445 END IF 446 447 IF (do_negf .AND. SIZE(negf_control%contacts) > 2) THEN 448 CALL cp_abort(__LOCATION__, & 449 "General case (> 2 contacts) has not been implemented yet") 450 END IF 451 452 CALL timestop(handle) 453 END SUBROUTINE read_negf_control 454 455! ************************************************************************************************** 456!> \brief Read region-specific list of atoms. 457!> \param atomlist list of atoms 458!> \param input_section input section which contains 'LIST' and 'MOLNAME' keywords 459!> \param i_rep_section repetition index of the input_section 460!> \param subsys subsystem environment 461! ************************************************************************************************** 462 SUBROUTINE read_negf_atomlist(atomlist, input_section, i_rep_section, subsys) 463 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(out) :: atomlist 464 TYPE(section_vals_type), POINTER :: input_section 465 INTEGER, INTENT(in) :: i_rep_section 466 TYPE(cp_subsys_type), POINTER :: subsys 467 468 CHARACTER(len=*), PARAMETER :: routineN = 'read_negf_atomlist', & 469 routineP = moduleN//':'//routineN 470 471 CHARACTER(len=default_string_length) :: index_str, natoms_str 472 CHARACTER(len=default_string_length), & 473 DIMENSION(:), POINTER :: cptr 474 INTEGER :: first_atom, handle, iatom, ikind, imol, iname, irep, last_atom, natoms_current, & 475 natoms_max, natoms_total, nkinds, nmols, nnames, nrep_list, nrep_molname 476 INTEGER, ALLOCATABLE, DIMENSION(:) :: inds 477 INTEGER, DIMENSION(:), POINTER :: iptr 478 LOGICAL :: is_list, is_molname 479 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 480 TYPE(molecule_kind_type), POINTER :: molecule_kind 481 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 482 TYPE(molecule_type), POINTER :: molecule 483 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 484 485 CALL timeset(routineN, handle) 486 487 CALL cp_subsys_get(subsys, particle_set=particle_set, & 488 molecule_set=molecule_set, & 489 molecule_kind_set=molecule_kind_set) 490 natoms_max = SIZE(particle_set) 491 nkinds = SIZE(molecule_kind_set) 492 493 CALL section_vals_val_get(input_section, "LIST", i_rep_section=i_rep_section, & 494 n_rep_val=nrep_list, explicit=is_list) 495 CALL section_vals_val_get(input_section, "MOLNAME", i_rep_section=i_rep_section, & 496 n_rep_val=nrep_molname, explicit=is_molname) 497 498 ! compute the number of atoms in the NEGF region, and check the validity of giben atomic indices 499 natoms_total = 0 500 IF (is_list .AND. nrep_list > 0) THEN 501 DO irep = 1, nrep_list 502 CALL section_vals_val_get(input_section, "LIST", i_rep_section=i_rep_section, i_rep_val=irep, i_vals=iptr) 503 504 natoms_current = SIZE(iptr) 505 DO iatom = 1, natoms_current 506 IF (iptr(iatom) > natoms_max) THEN 507 CALL integer_to_string(iptr(iatom), index_str) 508 CALL integer_to_string(natoms_max, natoms_str) 509 CALL cp_abort(__LOCATION__, & 510 "NEGF: Atomic index "//TRIM(index_str)//" given in section "// & 511 TRIM(input_section%section%name)//" exceeds the maximum number of atoms ("// & 512 TRIM(natoms_str)//").") 513 END IF 514 END DO 515 516 natoms_total = natoms_total + natoms_current 517 END DO 518 END IF 519 520 IF (is_molname .AND. nrep_molname > 0) THEN 521 DO irep = 1, nrep_molname 522 CALL section_vals_val_get(input_section, "MOLNAME", i_rep_section=i_rep_section, i_rep_val=irep, c_vals=cptr) 523 nnames = SIZE(cptr) 524 525 DO iname = 1, nnames 526 DO ikind = 1, nkinds 527 IF (molecule_kind_set(ikind)%name .EQ. cptr(iname)) EXIT 528 END DO 529 530 IF (ikind <= nkinds) THEN 531 molecule_kind => molecule_kind_set(ikind) 532 CALL get_molecule_kind(molecule_kind, nmolecule=nmols, molecule_list=iptr) 533 534 DO imol = 1, nmols 535 molecule => molecule_set(iptr(imol)) 536 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom) 537 natoms_current = last_atom - first_atom + 1 538 natoms_total = natoms_total + natoms_current 539 END DO 540 ELSE 541 CALL cp_abort(__LOCATION__, & 542 "NEGF: A molecule with the name '"//TRIM(cptr(iname))//"' mentioned in section "// & 543 TRIM(input_section%section%name)//" has not been defined. Note that names are case sensitive.") 544 END IF 545 END DO 546 END DO 547 END IF 548 549 ! create a list of atomic indices 550 IF (natoms_total > 0) THEN 551 ALLOCATE (atomlist(natoms_total)) 552 553 natoms_total = 0 554 555 IF (is_list .AND. nrep_list > 0) THEN 556 DO irep = 1, nrep_list 557 CALL section_vals_val_get(input_section, "LIST", i_rep_section=i_rep_section, i_rep_val=irep, i_vals=iptr) 558 559 natoms_current = SIZE(iptr) 560 atomlist(natoms_total + 1:natoms_total + natoms_current) = iptr(1:natoms_current) 561 natoms_total = natoms_total + natoms_current 562 END DO 563 END IF 564 565 IF (is_molname .AND. nrep_molname > 0) THEN 566 DO irep = 1, nrep_molname 567 CALL section_vals_val_get(input_section, "MOLNAME", i_rep_section=i_rep_section, i_rep_val=irep, c_vals=cptr) 568 nnames = SIZE(cptr) 569 570 DO iname = 1, nnames 571 DO ikind = 1, nkinds 572 IF (molecule_kind_set(ikind)%name .EQ. cptr(iname)) EXIT 573 END DO 574 575 IF (ikind <= nkinds) THEN 576 molecule_kind => molecule_kind_set(ikind) 577 CALL get_molecule_kind(molecule_kind, nmolecule=nmols, molecule_list=iptr) 578 579 DO imol = 1, nmols 580 molecule => molecule_set(iptr(imol)) 581 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom) 582 583 DO natoms_current = first_atom, last_atom 584 natoms_total = natoms_total + 1 585 atomlist(natoms_total) = natoms_current 586 END DO 587 END DO 588 END IF 589 END DO 590 END DO 591 END IF 592 593 ! remove duplicated atoms 594 ALLOCATE (inds(natoms_total)) 595 CALL sort(atomlist, natoms_total, inds) 596 DEALLOCATE (inds) 597 598 natoms_current = 1 599 DO iatom = natoms_current + 1, natoms_total 600 IF (atomlist(iatom) /= atomlist(natoms_current)) THEN 601 natoms_current = natoms_current + 1 602 atomlist(natoms_current) = atomlist(iatom) 603 END IF 604 END DO 605 606 IF (natoms_current < natoms_total) THEN 607 CALL MOVE_ALLOC(atomlist, inds) 608 609 ALLOCATE (atomlist(natoms_current)) 610 atomlist(1:natoms_current) = inds(1:natoms_current) 611 DEALLOCATE (inds) 612 END IF 613 END IF 614 615 CALL timestop(handle) 616 END SUBROUTINE read_negf_atomlist 617END MODULE negf_control_types 618