Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Map atoms between various force environments.
      10             : !> \author Sergey Chulkov
      11             : ! **************************************************************************************************
      12             : 
      13             : MODULE negf_atom_map
      14             :    USE atomic_kind_types,               ONLY: get_atomic_kind
      15             :    USE cell_types,                      ONLY: cell_type,&
      16             :                                               real_to_scaled
      17             :    USE kinds,                           ONLY: default_string_length,&
      18             :                                               dp
      19             :    USE particle_types,                  ONLY: particle_type
      20             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      21             :                                               qs_kind_type
      22             :    USE qs_subsys_types,                 ONLY: qs_subsys_get,&
      23             :                                               qs_subsys_type
      24             : #include "./base/base_uses.f90"
      25             : 
      26             :    IMPLICIT NONE
      27             :    PRIVATE
      28             : 
      29             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_atom_map'
      30             :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .TRUE.
      31             : 
      32             :    PUBLIC :: negf_atom_map_type, negf_map_atomic_indices
      33             : 
      34             : ! **************************************************************************************************
      35             : !> \brief Structure that maps the given atom in the sourse FORCE_EVAL section with another atom
      36             : !>        from the target FORCE_EVAL section.
      37             : ! **************************************************************************************************
      38             :    TYPE negf_atom_map_type
      39             :       !> atomic index within the target FORCE_EVAL
      40             :       INTEGER                                            :: iatom = -1
      41             :       !> cell replica
      42             :       INTEGER, DIMENSION(3)                              :: cell = -1
      43             :    END TYPE negf_atom_map_type
      44             : 
      45             :    PRIVATE :: qs_kind_group_type, qs_kind_groups_create, qs_kind_groups_release
      46             : 
      47             : ! **************************************************************************************************
      48             : !> \brief List of equivalent atoms.
      49             : ! **************************************************************************************************
      50             :    TYPE qs_kind_group_type
      51             :       !> atomic symbol
      52             :       CHARACTER(len=2)                                   :: element_symbol = ""
      53             :       !> number of atoms of this kind in 'atom_list'
      54             :       INTEGER                                            :: natoms = -1
      55             :       !> number of spherical Gaussian functions per atom
      56             :       INTEGER                                            :: nsgf = -1
      57             :       !> list of atomic indices
      58             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_list
      59             :       !> atomic coordinates [3 x natoms]
      60             :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: r
      61             :    END TYPE qs_kind_group_type
      62             : 
      63             : CONTAINS
      64             : ! **************************************************************************************************
      65             : !> \brief Map atoms in the cell 'subsys_device' listed in 'atom_list' with the corresponding
      66             : !>        atoms in the cell 'subsys_contact'.
      67             : !> \param atom_map        list of atoms in the cell 'subsys_contact' (initialised on exit)
      68             : !> \param atom_list       atomic indices of selected atoms in the cell 'subsys_device'
      69             : !> \param subsys_device   QuickStep subsystem of the device force environment
      70             : !> \param subsys_contact  QuickStep subsystem of the contact force environment
      71             : !> \param eps_geometry    accuracy in mapping atoms based on their Cartesian coordinates
      72             : !> \par History
      73             : !>   * 08.2017 created [Sergey Chulkov]
      74             : ! **************************************************************************************************
      75          52 :    SUBROUTINE negf_map_atomic_indices(atom_map, atom_list, subsys_device, subsys_contact, eps_geometry)
      76             :       TYPE(negf_atom_map_type), DIMENSION(:), &
      77             :          INTENT(out)                                     :: atom_map
      78             :       INTEGER, DIMENSION(:), INTENT(in)                  :: atom_list
      79             :       TYPE(qs_subsys_type), POINTER                      :: subsys_device, subsys_contact
      80             :       REAL(kind=dp), INTENT(in)                          :: eps_geometry
      81             : 
      82             :       CHARACTER(len=*), PARAMETER :: routineN = 'negf_map_atomic_indices'
      83             : 
      84             :       CHARACTER(len=2)                                   :: element_device
      85             :       CHARACTER(len=default_string_length)               :: atom_str
      86             :       INTEGER                                            :: atom_index_device, handle, iatom, &
      87             :                                                             iatom_kind, ikind, ikind_contact, &
      88             :                                                             natoms, nkinds_contact, nsgf_device
      89             :       REAL(kind=dp), DIMENSION(3)                        :: coords, coords_error, coords_scaled
      90             :       TYPE(cell_type), POINTER                           :: cell_contact
      91           4 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set_contact, particle_set_device
      92             :       TYPE(qs_kind_group_type), ALLOCATABLE, &
      93           4 :          DIMENSION(:)                                    :: kind_groups_contact
      94           4 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set_contact, qs_kind_set_device
      95             : 
      96           4 :       CALL timeset(routineN, handle)
      97             : 
      98           4 :       natoms = SIZE(atom_map, 1)
      99           4 :       CPASSERT(SIZE(atom_list) == natoms)
     100             : 
     101           4 :       CALL qs_subsys_get(subsys_device, particle_set=particle_set_device, qs_kind_set=qs_kind_set_device)
     102           4 :       CALL qs_subsys_get(subsys_contact, cell=cell_contact, particle_set=particle_set_contact, qs_kind_set=qs_kind_set_contact)
     103             : 
     104           4 :       CALL qs_kind_groups_create(kind_groups_contact, particle_set_contact, qs_kind_set_contact)
     105           4 :       nkinds_contact = SIZE(kind_groups_contact)
     106             : 
     107          36 :       DO iatom = 1, natoms
     108          32 :          atom_index_device = atom_list(iatom)
     109          32 :          CALL get_atomic_kind(particle_set_device(atom_index_device)%atomic_kind, kind_number=ikind)
     110          32 :          CALL get_qs_kind(qs_kind_set_device(ikind), element_symbol=element_device, nsgf=nsgf_device)
     111             : 
     112          32 :          atom_map(iatom)%iatom = 0
     113             : 
     114          32 :          iterate_kind: DO ikind_contact = 1, nkinds_contact
     115             :             ! looking for an equivalent atomic kind (based on the element symbol and the number of atomic basis functions)
     116          32 :             IF (kind_groups_contact(ikind_contact)%element_symbol == element_device .AND. &
     117           0 :                 kind_groups_contact(ikind_contact)%nsgf == nsgf_device) THEN
     118             : 
     119             :                ! loop over matching atoms
     120          80 :                DO iatom_kind = 1, kind_groups_contact(ikind_contact)%natoms
     121             :                   coords(1:3) = particle_set_device(atom_index_device)%r(1:3) - &
     122         320 :                                 kind_groups_contact(ikind_contact)%r(1:3, iatom_kind)
     123             : 
     124          80 :                   CALL real_to_scaled(coords_scaled, coords, cell_contact)
     125         320 :                   coords_error = coords_scaled - REAL(NINT(coords_scaled), kind=dp)
     126             : 
     127         320 :                   IF (SQRT(DOT_PRODUCT(coords_error, coords_error)) < eps_geometry) THEN
     128          32 :                      atom_map(iatom)%iatom = kind_groups_contact(ikind_contact)%atom_list(iatom_kind)
     129         128 :                      atom_map(iatom)%cell = NINT(coords_scaled)
     130             :                      EXIT iterate_kind
     131             :                   END IF
     132             :                END DO
     133             :             END IF
     134             :          END DO iterate_kind
     135             : 
     136          68 :          IF (atom_map(iatom)%iatom == 0) THEN
     137             :             ! atom has not been found in the corresponding force_env
     138           0 :             WRITE (atom_str, '(A2,3(1X,F11.6))') element_device, particle_set_device(atom_index_device)%r
     139             : 
     140             :             CALL cp_abort(__LOCATION__, &
     141           0 :                           "Unable to map the atom ("//TRIM(atom_str)//") onto the atom from the corresponding FORCE_EVAL section")
     142             :          END IF
     143             :       END DO
     144             : 
     145           4 :       CALL qs_kind_groups_release(kind_groups_contact)
     146             : 
     147           4 :       CALL timestop(handle)
     148           8 :    END SUBROUTINE negf_map_atomic_indices
     149             : 
     150             : ! **************************************************************************************************
     151             : !> \brief Group particles from 'particle_set' according to their atomic (QS) kind.
     152             : !> \param kind_groups   kind groups that will be created
     153             : !> \param particle_set  list of particles
     154             : !> \param qs_kind_set   list of QS kinds
     155             : !> \par History
     156             : !>   * 08.2017 created [Sergey Chulkov]
     157             : !> \note used within the subroutine negf_map_atomic_indices() in order to map atoms in
     158             : !>       a linear scalling fashion
     159             : ! **************************************************************************************************
     160           4 :    SUBROUTINE qs_kind_groups_create(kind_groups, particle_set, qs_kind_set)
     161             :       TYPE(qs_kind_group_type), ALLOCATABLE, &
     162             :          DIMENSION(:), INTENT(inout)                     :: kind_groups
     163             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     164             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     165             : 
     166             :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_kind_groups_create'
     167             : 
     168             :       INTEGER                                            :: handle, iatom, ikind, natoms, nkinds
     169             : 
     170           4 :       CALL timeset(routineN, handle)
     171             : 
     172           4 :       natoms = SIZE(particle_set)
     173           4 :       nkinds = 0
     174             : 
     175          20 :       DO iatom = 1, natoms
     176          16 :          CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     177          20 :          IF (nkinds < ikind) nkinds = ikind
     178             :       END DO
     179             : 
     180          16 :       ALLOCATE (kind_groups(nkinds))
     181             : 
     182           8 :       DO ikind = 1, nkinds
     183           4 :          kind_groups(ikind)%natoms = 0
     184           8 :          CALL get_qs_kind(qs_kind_set(ikind), element_symbol=kind_groups(ikind)%element_symbol, nsgf=kind_groups(ikind)%nsgf)
     185             :       END DO
     186             : 
     187          20 :       DO iatom = 1, natoms
     188          16 :          CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     189          20 :          kind_groups(ikind)%natoms = kind_groups(ikind)%natoms + 1
     190             :       END DO
     191             : 
     192           8 :       DO ikind = 1, nkinds
     193          12 :          ALLOCATE (kind_groups(ikind)%atom_list(kind_groups(ikind)%natoms))
     194          12 :          ALLOCATE (kind_groups(ikind)%r(3, kind_groups(ikind)%natoms))
     195             : 
     196           8 :          kind_groups(ikind)%natoms = 0
     197             :       END DO
     198             : 
     199          20 :       DO iatom = 1, natoms
     200          16 :          CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
     201          16 :          kind_groups(ikind)%natoms = kind_groups(ikind)%natoms + 1
     202             : 
     203          16 :          kind_groups(ikind)%atom_list(kind_groups(ikind)%natoms) = iatom
     204          68 :          kind_groups(ikind)%r(1:3, kind_groups(ikind)%natoms) = particle_set(iatom)%r(1:3)
     205             :       END DO
     206             : 
     207           4 :       CALL timestop(handle)
     208           4 :    END SUBROUTINE qs_kind_groups_create
     209             : 
     210             : ! **************************************************************************************************
     211             : !> \brief Release groups of particles.
     212             : !> \param kind_groups  kind groups to release
     213             : !> \par History
     214             : !>   * 08.2017 created [Sergey Chulkov]
     215             : ! **************************************************************************************************
     216           4 :    SUBROUTINE qs_kind_groups_release(kind_groups)
     217             :       TYPE(qs_kind_group_type), ALLOCATABLE, &
     218             :          DIMENSION(:), INTENT(inout)                     :: kind_groups
     219             : 
     220             :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_kind_groups_release'
     221             : 
     222             :       INTEGER                                            :: handle, ikind
     223             : 
     224           4 :       CALL timeset(routineN, handle)
     225             : 
     226           4 :       IF (ALLOCATED(kind_groups)) THEN
     227           8 :          DO ikind = SIZE(kind_groups), 1, -1
     228           4 :             IF (ALLOCATED(kind_groups(ikind)%atom_list)) DEALLOCATE (kind_groups(ikind)%atom_list)
     229           8 :             IF (ALLOCATED(kind_groups(ikind)%r)) DEALLOCATE (kind_groups(ikind)%r)
     230             :          END DO
     231             : 
     232           8 :          DEALLOCATE (kind_groups)
     233             :       END IF
     234             : 
     235           4 :       CALL timestop(handle)
     236           4 :    END SUBROUTINE qs_kind_groups_release
     237           0 : END MODULE negf_atom_map