1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Define the data structure for the particle information.
8!> \par History
9!>      - Atomic kind added in particle_type (MK,08.01.2002)
10!>      - Functionality for particle_type added (MK,14.01.2002)
11!>      - Allow for general coordinate input (MK,13.09.2003)
12!>      - Molecule concept introduced (MK,26.09.2003)
13!>      - Last atom information added (jgh,23.05.2004)
14!>      - particle_type cleaned (MK,03.02.2005)
15!> \author CJM, MK
16! **************************************************************************************************
17MODULE particle_types
18   USE atomic_kind_types,               ONLY: atomic_kind_type
19   USE kinds,                           ONLY: dp
20   USE message_passing,                 ONLY: mp_sum
21#include "../base/base_uses.f90"
22
23   IMPLICIT NONE
24
25   PRIVATE
26
27   ! Global parameters (in this module)
28
29   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'particle_types'
30
31   ! Data types
32! **************************************************************************************************
33   TYPE particle_type
34      TYPE(atomic_kind_type), POINTER       :: atomic_kind => Null() ! atomic kind information
35      REAL(KIND=dp), DIMENSION(3)           :: f = 0.0_dp, & ! force
36                                               r = 0.0_dp, & ! position
37                                               v = 0.0_dp ! velocity
38      ! Particle dependent terms for shell-model
39      INTEGER                               :: atom_index = -1, &
40                                               t_region_index = -1, &
41                                               shell_index = -1
42   END TYPE particle_type
43
44   ! Public data types
45
46   PUBLIC :: particle_type
47
48   ! Public subroutines
49
50   PUBLIC :: allocate_particle_set, &
51             deallocate_particle_set, &
52             update_particle_set, &
53             update_particle_pos_or_vel, &
54             get_particle_pos_or_vel
55
56CONTAINS
57
58! **************************************************************************************************
59!> \brief   Allocate a particle set.
60!> \param particle_set ...
61!> \param nparticle ...
62!> \date    14.01.2002
63!> \author  MK
64!> \version 1.0
65! **************************************************************************************************
66   SUBROUTINE allocate_particle_set(particle_set, nparticle)
67      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
68      INTEGER, INTENT(IN)                                :: nparticle
69
70      CHARACTER(len=*), PARAMETER :: routineN = 'allocate_particle_set', &
71         routineP = moduleN//':'//routineN
72
73      INTEGER                                            :: iparticle
74
75      IF (ASSOCIATED(particle_set)) THEN
76         CALL deallocate_particle_set(particle_set)
77      END IF
78      ALLOCATE (particle_set(nparticle))
79
80      DO iparticle = 1, nparticle
81         NULLIFY (particle_set(iparticle)%atomic_kind)
82         particle_set(iparticle)%f(:) = 0.0_dp
83         particle_set(iparticle)%r(:) = 0.0_dp
84         particle_set(iparticle)%v(:) = 0.0_dp
85         particle_set(iparticle)%shell_index = 0
86         particle_set(iparticle)%atom_index = 0
87         particle_set(iparticle)%t_region_index = 0
88      END DO
89
90   END SUBROUTINE allocate_particle_set
91
92! **************************************************************************************************
93!> \brief   Deallocate a particle set.
94!> \param particle_set ...
95!> \date    14.01.2002
96!> \author  MK
97!> \version 1.0
98! **************************************************************************************************
99   SUBROUTINE deallocate_particle_set(particle_set)
100      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
101
102      CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_particle_set', &
103         routineP = moduleN//':'//routineN
104
105      IF (ASSOCIATED(particle_set)) THEN
106         DEALLOCATE (particle_set)
107      ELSE
108         CALL cp_abort(__LOCATION__, &
109                       "The pointer particle_set is not associated and "// &
110                       "cannot be deallocated")
111      END IF
112
113   END SUBROUTINE deallocate_particle_set
114
115! **************************************************************************************************
116!> \brief ...
117!> \param particle_set ...
118!> \param int_group ...
119!> \param pos ...
120!> \param vel ...
121!> \param for ...
122!> \param add ...
123! **************************************************************************************************
124   SUBROUTINE update_particle_set(particle_set, int_group, pos, vel, for, add)
125
126      TYPE(particle_type), POINTER                       :: particle_set(:)
127      INTEGER, INTENT(IN)                                :: int_group
128      REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: pos(:, :), vel(:, :), for(:, :)
129      LOGICAL, INTENT(IN), OPTIONAL                      :: add
130
131      CHARACTER(len=*), PARAMETER :: routineN = 'update_particle_set', &
132         routineP = moduleN//':'//routineN
133
134      INTEGER                                            :: handle, iparticle, nparticle
135      LOGICAL                                            :: my_add, update_for, update_pos, &
136                                                            update_vel
137
138      CALL timeset(routineN, handle)
139
140      nparticle = SIZE(particle_set)
141      update_pos = PRESENT(pos)
142      update_vel = PRESENT(vel)
143      update_for = PRESENT(for)
144      my_add = .FALSE.
145      IF (PRESENT(add)) my_add = add
146
147      IF (update_pos) THEN
148         CALL mp_sum(pos, int_group)
149         IF (my_add) THEN
150            DO iparticle = 1, nparticle
151               particle_set(iparticle)%r(:) = particle_set(iparticle)%r(:) + pos(:, iparticle)
152            END DO
153         ELSE
154            DO iparticle = 1, nparticle
155               particle_set(iparticle)%r(:) = pos(:, iparticle)
156            END DO
157         END IF
158      END IF
159      IF (update_vel) THEN
160         CALL mp_sum(vel, int_group)
161         IF (my_add) THEN
162            DO iparticle = 1, nparticle
163               particle_set(iparticle)%v(:) = particle_set(iparticle)%v(:) + vel(:, iparticle)
164            END DO
165         ELSE
166            DO iparticle = 1, nparticle
167               particle_set(iparticle)%v(:) = vel(:, iparticle)
168            END DO
169         END IF
170      END IF
171      IF (update_for) THEN
172         CALL mp_sum(for, int_group)
173         IF (my_add) THEN
174            DO iparticle = 1, nparticle
175               particle_set(iparticle)%f(:) = particle_set(iparticle)%f(:) + for(:, iparticle)
176            END DO
177         ELSE
178            DO iparticle = 1, nparticle
179               particle_set(iparticle)%f(:) = for(:, iparticle)
180            END DO
181         END IF
182      END IF
183
184      CALL timestop(handle)
185
186   END SUBROUTINE update_particle_set
187
188! **************************************************************************************************
189!> \brief   Return the atomic position or velocity of atom iatom in x from a
190!>          packed vector even if core-shell particles are present
191!> \param iatom ...
192!> \param particle_set ...
193!> \param vector ...
194!> \return ...
195!> \date    25.11.2010
196!> \author  Matthias Krack
197!> \version 1.0
198! **************************************************************************************************
199   FUNCTION get_particle_pos_or_vel(iatom, particle_set, vector) RESULT(x)
200
201      INTEGER, INTENT(IN)                                :: iatom
202      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
203      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: vector
204      REAL(KIND=dp), DIMENSION(3)                        :: x
205
206      INTEGER                                            :: ic, is
207      REAL(KIND=dp)                                      :: fc, fs, mass
208
209      ic = 3*(iatom - 1)
210      IF (particle_set(iatom)%shell_index == 0) THEN
211         x(1:3) = vector(ic + 1:ic + 3)
212      ELSE
213         is = 3*(SIZE(particle_set) + particle_set(iatom)%shell_index - 1)
214         mass = particle_set(iatom)%atomic_kind%mass
215         fc = particle_set(iatom)%atomic_kind%shell%mass_core/mass
216         fs = particle_set(iatom)%atomic_kind%shell%mass_shell/mass
217         x(1:3) = fc*vector(ic + 1:ic + 3) + fs*vector(is + 1:is + 3)
218      END IF
219
220   END FUNCTION get_particle_pos_or_vel
221
222! **************************************************************************************************
223!> \brief   Update the atomic position or velocity by x and return the updated
224!>          atomic position or velocity in x even if core-shell particles are
225!>          present
226!> \param iatom ...
227!> \param particle_set ...
228!> \param x ...
229!> \param vector ...
230!> \date    26.11.2010
231!> \author  Matthias Krack
232!> \version 1.0
233!> \note    particle-set is not changed, only the positions or velocities in
234!>          the packed vector are updated
235! **************************************************************************************************
236   SUBROUTINE update_particle_pos_or_vel(iatom, particle_set, x, vector)
237
238      INTEGER, INTENT(IN)                                :: iatom
239      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
240      REAL(KIND=dp), DIMENSION(3), INTENT(INOUT)         :: x
241      REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: vector
242
243      INTEGER                                            :: ic, is
244      REAL(KIND=dp)                                      :: fc, fs, mass
245
246      ic = 3*(iatom - 1)
247      IF (particle_set(iatom)%shell_index == 0) THEN
248         vector(ic + 1:ic + 3) = vector(ic + 1:ic + 3) + x(1:3)
249         x(1:3) = vector(ic + 1:ic + 3)
250      ELSE
251         is = 3*(SIZE(particle_set) + particle_set(iatom)%shell_index - 1)
252         mass = particle_set(iatom)%atomic_kind%mass
253         fc = particle_set(iatom)%atomic_kind%shell%mass_core/mass
254         fs = particle_set(iatom)%atomic_kind%shell%mass_shell/mass
255         vector(ic + 1:ic + 3) = vector(ic + 1:ic + 3) + x(1:3)
256         vector(is + 1:is + 3) = vector(is + 1:is + 3) + x(1:3)
257         x(1:3) = fc*vector(ic + 1:ic + 3) + fs*vector(is + 1:is + 3)
258      END IF
259
260   END SUBROUTINE update_particle_pos_or_vel
261
262END MODULE particle_types
263