(git:34ef472)
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-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
41  !> cell replica
42  INTEGER, DIMENSION(3) :: cell
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
55  !> number of spherical Gaussian functions per atom
56  INTEGER :: nsgf
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  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  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set_contact, particle_set_device
92  TYPE(qs_kind_group_type), ALLOCATABLE, &
93  DIMENSION(:) :: kind_groups_contact
94  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set_contact, qs_kind_set_device
95 
96  CALL timeset(routinen, handle)
97 
98  natoms = SIZE(atom_map, 1)
99  cpassert(SIZE(atom_list) == natoms)
100 
101  CALL qs_subsys_get(subsys_device, particle_set=particle_set_device, qs_kind_set=qs_kind_set_device)
102  CALL qs_subsys_get(subsys_contact, cell=cell_contact, particle_set=particle_set_contact, qs_kind_set=qs_kind_set_contact)
103 
104  CALL qs_kind_groups_create(kind_groups_contact, particle_set_contact, qs_kind_set_contact)
105  nkinds_contact = SIZE(kind_groups_contact)
106 
107  DO iatom = 1, natoms
108  atom_index_device = atom_list(iatom)
109  CALL get_atomic_kind(particle_set_device(atom_index_device)%atomic_kind, kind_number=ikind)
110  CALL get_qs_kind(qs_kind_set_device(ikind), element_symbol=element_device, nsgf=nsgf_device)
111 
112  atom_map(iatom)%iatom = 0
113 
114  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  IF (kind_groups_contact(ikind_contact)%element_symbol == element_device .AND. &
117  kind_groups_contact(ikind_contact)%nsgf == nsgf_device) THEN
118 
119  ! loop over matching atoms
120  DO iatom_kind = 1, kind_groups_contact(ikind_contact)%natoms
121  coords(1:3) = particle_set_device(atom_index_device)%r(1:3) - &
122  kind_groups_contact(ikind_contact)%r(1:3, iatom_kind)
123 
124  CALL real_to_scaled(coords_scaled, coords, cell_contact)
125  coords_error = coords_scaled - real(nint(coords_scaled), kind=dp)
126 
127  IF (sqrt(dot_product(coords_error, coords_error)) < eps_geometry) THEN
128  atom_map(iatom)%iatom = kind_groups_contact(ikind_contact)%atom_list(iatom_kind)
129  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  IF (atom_map(iatom)%iatom == 0) THEN
137  ! atom has not been found in the corresponding force_env
138  WRITE (atom_str, '(A2,3(1X,F11.6))') element_device, particle_set_device(atom_index_device)%r
139 
140  CALL cp_abort(__location__, &
141  "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  CALL qs_kind_groups_release(kind_groups_contact)
146 
147  CALL timestop(handle)
148  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  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  CALL timeset(routinen, handle)
171 
172  natoms = SIZE(particle_set)
173  nkinds = 0
174 
175  DO iatom = 1, natoms
176  CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
177  IF (nkinds < ikind) nkinds = ikind
178  END DO
179 
180  ALLOCATE (kind_groups(nkinds))
181 
182  DO ikind = 1, nkinds
183  kind_groups(ikind)%natoms = 0
184  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  DO iatom = 1, natoms
188  CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
189  kind_groups(ikind)%natoms = kind_groups(ikind)%natoms + 1
190  END DO
191 
192  DO ikind = 1, nkinds
193  ALLOCATE (kind_groups(ikind)%atom_list(kind_groups(ikind)%natoms))
194  ALLOCATE (kind_groups(ikind)%r(3, kind_groups(ikind)%natoms))
195 
196  kind_groups(ikind)%natoms = 0
197  END DO
198 
199  DO iatom = 1, natoms
200  CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
201  kind_groups(ikind)%natoms = kind_groups(ikind)%natoms + 1
202 
203  kind_groups(ikind)%atom_list(kind_groups(ikind)%natoms) = iatom
204  kind_groups(ikind)%r(1:3, kind_groups(ikind)%natoms) = particle_set(iatom)%r(1:3)
205  END DO
206 
207  CALL timestop(handle)
208  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  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  CALL timeset(routinen, handle)
225 
226  IF (ALLOCATED(kind_groups)) THEN
227  DO ikind = SIZE(kind_groups), 1, -1
228  IF (ALLOCATED(kind_groups(ikind)%atom_list)) DEALLOCATE (kind_groups(ikind)%atom_list)
229  IF (ALLOCATED(kind_groups(ikind)%r)) DEALLOCATE (kind_groups(ikind)%r)
230  END DO
231 
232  DEALLOCATE (kind_groups)
233  END IF
234 
235  CALL timestop(handle)
236  END SUBROUTINE qs_kind_groups_release
237 END 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:486
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:76
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, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, 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_r3d_rs_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, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
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