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 : !> * 10.2025 Centering of contact coordinates is added in the end. It is necessary to keep the
75 : !> right order of indices in the real space image matrices. It is made only here, because
76 : !> the mapping is based on the comparison of the coordinates. [Dmitry Ryndyk]
77 : !> \note
78 : ! **************************************************************************************************
79 52 : SUBROUTINE negf_map_atomic_indices(atom_map, atom_list, subsys_device, subsys_contact, eps_geometry)
80 : TYPE(negf_atom_map_type), DIMENSION(:), &
81 : INTENT(out) :: atom_map
82 : INTEGER, DIMENSION(:), INTENT(in) :: atom_list
83 : TYPE(qs_subsys_type), POINTER :: subsys_device, subsys_contact
84 : REAL(kind=dp), INTENT(in) :: eps_geometry
85 :
86 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_map_atomic_indices'
87 :
88 : CHARACTER(len=2) :: element_device
89 : CHARACTER(len=default_string_length) :: atom_str
90 : INTEGER :: atom_index_device, handle, iatom, &
91 : iatom_kind, ikind, ikind_contact, &
92 : natoms, nkinds_contact, nsgf_device
93 : REAL(kind=dp), DIMENSION(3) :: coords, coords_error, coords_scaled
94 : TYPE(cell_type), POINTER :: cell_contact
95 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set_contact, particle_set_device
96 : TYPE(qs_kind_group_type), ALLOCATABLE, &
97 4 : DIMENSION(:) :: kind_groups_contact
98 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set_contact, qs_kind_set_device
99 :
100 4 : CALL timeset(routineN, handle)
101 :
102 4 : natoms = SIZE(atom_map, 1)
103 4 : CPASSERT(SIZE(atom_list) == natoms)
104 :
105 4 : CALL qs_subsys_get(subsys_device, particle_set=particle_set_device, qs_kind_set=qs_kind_set_device)
106 4 : CALL qs_subsys_get(subsys_contact, cell=cell_contact, particle_set=particle_set_contact, qs_kind_set=qs_kind_set_contact)
107 :
108 4 : CALL qs_kind_groups_create(kind_groups_contact, particle_set_contact, qs_kind_set_contact)
109 4 : nkinds_contact = SIZE(kind_groups_contact)
110 :
111 36 : DO iatom = 1, natoms
112 32 : atom_index_device = atom_list(iatom)
113 32 : CALL get_atomic_kind(particle_set_device(atom_index_device)%atomic_kind, kind_number=ikind)
114 32 : CALL get_qs_kind(qs_kind_set_device(ikind), element_symbol=element_device, nsgf=nsgf_device)
115 :
116 32 : atom_map(iatom)%iatom = 0
117 :
118 32 : iterate_kind: DO ikind_contact = 1, nkinds_contact
119 : ! looking for an equivalent atomic kind (based on the element symbol and the number of atomic basis functions)
120 32 : IF (kind_groups_contact(ikind_contact)%element_symbol == element_device .AND. &
121 0 : kind_groups_contact(ikind_contact)%nsgf == nsgf_device) THEN
122 :
123 : ! loop over matching atoms
124 80 : DO iatom_kind = 1, kind_groups_contact(ikind_contact)%natoms
125 : coords(1:3) = particle_set_device(atom_index_device)%r(1:3) - &
126 320 : kind_groups_contact(ikind_contact)%r(1:3, iatom_kind)
127 :
128 80 : CALL real_to_scaled(coords_scaled, coords, cell_contact)
129 320 : coords_error = coords_scaled - REAL(NINT(coords_scaled), kind=dp)
130 :
131 320 : IF (DOT_PRODUCT(coords_error, coords_error) < (eps_geometry*eps_geometry)) THEN
132 32 : atom_map(iatom)%iatom = kind_groups_contact(ikind_contact)%atom_list(iatom_kind)
133 128 : atom_map(iatom)%cell = NINT(coords_scaled)
134 : EXIT iterate_kind
135 : END IF
136 : END DO
137 : END IF
138 : END DO iterate_kind
139 :
140 68 : IF (atom_map(iatom)%iatom == 0) THEN
141 : ! atom has not been found in the corresponding force_env
142 0 : WRITE (atom_str, '(A2,3(1X,F11.6))') element_device, particle_set_device(atom_index_device)%r
143 :
144 : CALL cp_abort(__LOCATION__, &
145 0 : "Unable to map the atom ("//TRIM(atom_str)//") onto the atom from the corresponding FORCE_EVAL section")
146 : END IF
147 : END DO
148 :
149 4 : CALL qs_kind_groups_release(kind_groups_contact)
150 :
151 4 : CALL centering_contact_coordinates(subsys=subsys_contact)
152 :
153 4 : CALL timestop(handle)
154 8 : END SUBROUTINE negf_map_atomic_indices
155 :
156 : ! **************************************************************************************************
157 : !> \brief Centering the atom coordinates of the primary unit cell of the bulk electrode.
158 : !> \param subsys ...
159 : !> \par History
160 : !> * 10.2025 created [Dmitry Ryndyk]
161 : !> \note It is necessary to keep the right order of indices in the real space image matrices.
162 : ! **************************************************************************************************
163 4 : SUBROUTINE centering_contact_coordinates(subsys)
164 : TYPE(qs_subsys_type), POINTER :: subsys
165 :
166 : REAL(KIND=dp) :: shiftX, shiftY, shiftZ
167 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
168 :
169 4 : CALL qs_subsys_get(subsys, particle_set=particle_set)
170 :
171 48 : shiftX = (MAXVAL(particle_set(:)%r(1)) + MINVAL(particle_set(:)%r(1)))/2.0
172 48 : shiftY = (MAXVAL(particle_set(:)%r(2)) + MINVAL(particle_set(:)%r(2)))/2.0
173 48 : shiftZ = (MAXVAL(particle_set(:)%r(3)) + MINVAL(particle_set(:)%r(3)))/2.0
174 :
175 20 : particle_set(:)%r(1) = particle_set(:)%r(1) - shiftX
176 20 : particle_set(:)%r(2) = particle_set(:)%r(2) - shiftY
177 20 : particle_set(:)%r(3) = particle_set(:)%r(3) - shiftZ
178 :
179 4 : END SUBROUTINE centering_contact_coordinates
180 :
181 : ! **************************************************************************************************
182 : !> \brief Group particles from 'particle_set' according to their atomic (QS) kind.
183 : !> \param kind_groups kind groups that will be created
184 : !> \param particle_set list of particles
185 : !> \param qs_kind_set list of QS kinds
186 : !> \par History
187 : !> * 08.2017 created [Sergey Chulkov]
188 : !> \note used within the subroutine negf_map_atomic_indices() in order to map atoms in
189 : !> a linear scalling fashion
190 : ! **************************************************************************************************
191 4 : SUBROUTINE qs_kind_groups_create(kind_groups, particle_set, qs_kind_set)
192 : TYPE(qs_kind_group_type), ALLOCATABLE, &
193 : DIMENSION(:), INTENT(inout) :: kind_groups
194 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
195 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
196 :
197 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_kind_groups_create'
198 :
199 : INTEGER :: handle, iatom, ikind, natoms, nkinds
200 :
201 4 : CALL timeset(routineN, handle)
202 :
203 4 : natoms = SIZE(particle_set)
204 4 : nkinds = 0
205 :
206 20 : DO iatom = 1, natoms
207 16 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
208 20 : IF (nkinds < ikind) nkinds = ikind
209 : END DO
210 :
211 16 : ALLOCATE (kind_groups(nkinds))
212 :
213 8 : DO ikind = 1, nkinds
214 4 : kind_groups(ikind)%natoms = 0
215 8 : CALL get_qs_kind(qs_kind_set(ikind), element_symbol=kind_groups(ikind)%element_symbol, nsgf=kind_groups(ikind)%nsgf)
216 : END DO
217 :
218 20 : DO iatom = 1, natoms
219 16 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
220 20 : kind_groups(ikind)%natoms = kind_groups(ikind)%natoms + 1
221 : END DO
222 :
223 8 : DO ikind = 1, nkinds
224 12 : ALLOCATE (kind_groups(ikind)%atom_list(kind_groups(ikind)%natoms))
225 12 : ALLOCATE (kind_groups(ikind)%r(3, kind_groups(ikind)%natoms))
226 :
227 8 : kind_groups(ikind)%natoms = 0
228 : END DO
229 :
230 20 : DO iatom = 1, natoms
231 16 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
232 16 : kind_groups(ikind)%natoms = kind_groups(ikind)%natoms + 1
233 :
234 16 : kind_groups(ikind)%atom_list(kind_groups(ikind)%natoms) = iatom
235 68 : kind_groups(ikind)%r(1:3, kind_groups(ikind)%natoms) = particle_set(iatom)%r(1:3)
236 : END DO
237 :
238 4 : CALL timestop(handle)
239 4 : END SUBROUTINE qs_kind_groups_create
240 :
241 : ! **************************************************************************************************
242 : !> \brief Release groups of particles.
243 : !> \param kind_groups kind groups to release
244 : !> \par History
245 : !> * 08.2017 created [Sergey Chulkov]
246 : ! **************************************************************************************************
247 4 : SUBROUTINE qs_kind_groups_release(kind_groups)
248 : TYPE(qs_kind_group_type), ALLOCATABLE, &
249 : DIMENSION(:), INTENT(inout) :: kind_groups
250 :
251 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_kind_groups_release'
252 :
253 : INTEGER :: handle, ikind
254 :
255 4 : CALL timeset(routineN, handle)
256 :
257 4 : IF (ALLOCATED(kind_groups)) THEN
258 8 : DO ikind = SIZE(kind_groups), 1, -1
259 4 : IF (ALLOCATED(kind_groups(ikind)%atom_list)) DEALLOCATE (kind_groups(ikind)%atom_list)
260 8 : IF (ALLOCATED(kind_groups(ikind)%r)) DEALLOCATE (kind_groups(ikind)%r)
261 : END DO
262 :
263 8 : DEALLOCATE (kind_groups)
264 : END IF
265 :
266 4 : CALL timestop(handle)
267 4 : END SUBROUTINE qs_kind_groups_release
268 0 : END MODULE negf_atom_map
|
