Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 (DOT_PRODUCT(coords_error, coords_error) < (eps_geometry*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
|
