(git:07c9450)
Loading...
Searching...
No Matches
negf_atom_map.F
Go to the documentation of this file.
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
13MODULE 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
63CONTAINS
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 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set_contact, particle_set_device
96 TYPE(qs_kind_group_type), ALLOCATABLE, &
97 DIMENSION(:) :: kind_groups_contact
98 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set_contact, qs_kind_set_device
99
100 CALL timeset(routinen, handle)
101
102 natoms = SIZE(atom_map, 1)
103 cpassert(SIZE(atom_list) == natoms)
104
105 CALL qs_subsys_get(subsys_device, particle_set=particle_set_device, qs_kind_set=qs_kind_set_device)
106 CALL qs_subsys_get(subsys_contact, cell=cell_contact, particle_set=particle_set_contact, qs_kind_set=qs_kind_set_contact)
107
108 CALL qs_kind_groups_create(kind_groups_contact, particle_set_contact, qs_kind_set_contact)
109 nkinds_contact = SIZE(kind_groups_contact)
110
111 DO iatom = 1, natoms
112 atom_index_device = atom_list(iatom)
113 CALL get_atomic_kind(particle_set_device(atom_index_device)%atomic_kind, kind_number=ikind)
114 CALL get_qs_kind(qs_kind_set_device(ikind), element_symbol=element_device, nsgf=nsgf_device)
115
116 atom_map(iatom)%iatom = 0
117
118 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 IF (kind_groups_contact(ikind_contact)%element_symbol == element_device .AND. &
121 kind_groups_contact(ikind_contact)%nsgf == nsgf_device) THEN
122
123 ! loop over matching atoms
124 DO iatom_kind = 1, kind_groups_contact(ikind_contact)%natoms
125 coords(1:3) = particle_set_device(atom_index_device)%r(1:3) - &
126 kind_groups_contact(ikind_contact)%r(1:3, iatom_kind)
127
128 CALL real_to_scaled(coords_scaled, coords, cell_contact)
129 coords_error = coords_scaled - real(nint(coords_scaled), kind=dp)
130
131 IF (dot_product(coords_error, coords_error) < (eps_geometry*eps_geometry)) THEN
132 atom_map(iatom)%iatom = kind_groups_contact(ikind_contact)%atom_list(iatom_kind)
133 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 IF (atom_map(iatom)%iatom == 0) THEN
141 ! atom has not been found in the corresponding force_env
142 WRITE (atom_str, '(A2,3(1X,F11.6))') element_device, particle_set_device(atom_index_device)%r
143
144 CALL cp_abort(__location__, &
145 "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 CALL qs_kind_groups_release(kind_groups_contact)
150
151 CALL centering_contact_coordinates(subsys=subsys_contact)
152
153 CALL timestop(handle)
154 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 SUBROUTINE centering_contact_coordinates(subsys)
164 TYPE(qs_subsys_type), POINTER :: subsys
165
166 REAL(kind=dp) :: shiftx, shifty, shiftz
167 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
168
169 CALL qs_subsys_get(subsys, particle_set=particle_set)
170
171 shiftx = (maxval(particle_set(:)%r(1)) + minval(particle_set(:)%r(1)))/2.0
172 shifty = (maxval(particle_set(:)%r(2)) + minval(particle_set(:)%r(2)))/2.0
173 shiftz = (maxval(particle_set(:)%r(3)) + minval(particle_set(:)%r(3)))/2.0
174
175 particle_set(:)%r(1) = particle_set(:)%r(1) - shiftx
176 particle_set(:)%r(2) = particle_set(:)%r(2) - shifty
177 particle_set(:)%r(3) = particle_set(:)%r(3) - shiftz
178
179 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 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 CALL timeset(routinen, handle)
202
203 natoms = SIZE(particle_set)
204 nkinds = 0
205
206 DO iatom = 1, natoms
207 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
208 IF (nkinds < ikind) nkinds = ikind
209 END DO
210
211 ALLOCATE (kind_groups(nkinds))
212
213 DO ikind = 1, nkinds
214 kind_groups(ikind)%natoms = 0
215 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 DO iatom = 1, natoms
219 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
220 kind_groups(ikind)%natoms = kind_groups(ikind)%natoms + 1
221 END DO
222
223 DO ikind = 1, nkinds
224 ALLOCATE (kind_groups(ikind)%atom_list(kind_groups(ikind)%natoms))
225 ALLOCATE (kind_groups(ikind)%r(3, kind_groups(ikind)%natoms))
226
227 kind_groups(ikind)%natoms = 0
228 END DO
229
230 DO iatom = 1, natoms
231 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
232 kind_groups(ikind)%natoms = kind_groups(ikind)%natoms + 1
233
234 kind_groups(ikind)%atom_list(kind_groups(ikind)%natoms) = iatom
235 kind_groups(ikind)%r(1:3, kind_groups(ikind)%natoms) = particle_set(iatom)%r(1:3)
236 END DO
237
238 CALL timestop(handle)
239 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 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 CALL timeset(routinen, handle)
256
257 IF (ALLOCATED(kind_groups)) THEN
258 DO ikind = SIZE(kind_groups), 1, -1
259 IF (ALLOCATED(kind_groups(ikind)%atom_list)) DEALLOCATE (kind_groups(ikind)%atom_list)
260 IF (ALLOCATED(kind_groups(ikind)%r)) DEALLOCATE (kind_groups(ikind)%r)
261 END DO
262
263 DEALLOCATE (kind_groups)
264 END IF
265
266 CALL timestop(handle)
267 END SUBROUTINE qs_kind_groups_release
268END MODULE negf_atom_map
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:491
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
negf_atom_map
Map atoms between various force environments.
Definition negf_atom_map.F:13
negf_atom_map::negf_map_atomic_indices
subroutine, public negf_map_atomic_indices(atom_map, atom_list, subsys_device, subsys_contact, eps_geometry)
Map atoms in the cell 'subsys_device' listed in 'atom_list' with the corresponding atoms in the cell ...
Definition negf_atom_map.F:80
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition qs_subsys_types.F:134
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
negf_atom_map::negf_atom_map_type
Structure that maps the given atom in the sourse FORCE_EVAL section with another atom from the target...
Definition negf_atom_map.F:38
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56