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 Tools common both to PME and SPME
10 : !> \par History
11 : !> JGH (03-May-2001) : first correctly working version
12 : !> teo (Feb-2007) : Merging common routines to spme and pme
13 : !> \author JGH (21-Mar-2001)
14 : ! **************************************************************************************************
15 : MODULE pme_tools
16 :
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind
19 : USE cell_types, ONLY: cell_type,&
20 : real_to_scaled
21 : USE kinds, ONLY: dp
22 : USE particle_types, ONLY: particle_type
23 : USE realspace_grid_types, ONLY: realspace_grid_type
24 : #include "./base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 : PUBLIC :: get_center, set_list
30 :
31 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pme_tools'
32 :
33 : CONTAINS
34 :
35 : ! **************************************************************************************************
36 : !> \brief ...
37 : !> \param part ...
38 : !> \param npart ...
39 : !> \param center ...
40 : !> \param p1 ...
41 : !> \param rs ...
42 : !> \param ipart ...
43 : !> \param core_center ...
44 : ! **************************************************************************************************
45 9415363 : SUBROUTINE set_list(part, npart, center, p1, rs, ipart, core_center)
46 :
47 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: part
48 : INTEGER, INTENT(IN) :: npart
49 : INTEGER, DIMENSION(:, :), INTENT(IN) :: center
50 : INTEGER, INTENT(OUT) :: p1
51 : TYPE(realspace_grid_type), INTENT(IN) :: rs
52 : INTEGER, INTENT(INOUT) :: ipart
53 : INTEGER, DIMENSION(:, :), OPTIONAL, POINTER :: core_center
54 :
55 : INTEGER :: ndim, npos
56 : INTEGER, DIMENSION(3) :: lb, ub
57 : REAL(KIND=dp) :: charge
58 : TYPE(atomic_kind_type), POINTER :: atomic_kind
59 :
60 9415363 : p1 = 0
61 37661452 : lb = rs%lb_real
62 37661452 : ub = rs%ub_real
63 :
64 16753852 : DO
65 16931213 : ipart = ipart + 1
66 16931213 : IF (ipart > npart) EXIT
67 16753852 : atomic_kind => part(ipart)%atomic_kind
68 16753852 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=charge)
69 16753852 : IF (charge == 0.0_dp .AND. part(ipart)%shell_index == 0) CYCLE
70 16843973 : IF (rs%desc%parallel) THEN
71 : ! check if the rs grid is distributed or not
72 57357536 : IF (ALL(rs%desc%group_dim == 1)) THEN
73 13821548 : ndim = rs%desc%group_size
74 13821548 : npos = rs%desc%my_pos
75 : ! All processors work on the same grid
76 13821548 : IF (MOD(ipart, ndim) == npos) THEN
77 6910774 : p1 = ipart
78 6910774 : EXIT
79 : END IF
80 : ELSE
81 : ! First check if this atom is on my grid
82 1035672 : IF (part(ipart)%shell_index /= 0 .AND. PRESENT(core_center)) THEN
83 0 : IF (in_slice(core_center(:, part(ipart)%shell_index), lb, ub)) THEN
84 0 : p1 = ipart
85 : END IF
86 : ELSE
87 1035672 : IF (in_slice(center(:, ipart), lb, ub)) THEN
88 517836 : p1 = ipart
89 : EXIT
90 : END IF
91 : END IF
92 : END IF
93 : ELSE
94 1809392 : p1 = ipart
95 1809392 : EXIT
96 : END IF
97 : END DO
98 :
99 9415363 : END SUBROUTINE set_list
100 :
101 : ! **************************************************************************************************
102 : !> \brief ...
103 : !> \param pos ...
104 : !> \param lb ...
105 : !> \param ub ...
106 : !> \return ...
107 : ! **************************************************************************************************
108 1035672 : FUNCTION in_slice(pos, lb, ub) RESULT(internal)
109 :
110 : INTEGER, DIMENSION(3), INTENT(IN) :: pos, lb, ub
111 : LOGICAL :: internal
112 :
113 7249704 : IF (ALL(pos >= lb) .AND. ALL(pos <= ub)) THEN
114 : internal = .TRUE.
115 : ELSE
116 517836 : internal = .FALSE.
117 : END IF
118 :
119 1035672 : END FUNCTION in_slice
120 :
121 : ! **************************************************************************************************
122 : !> \brief ...
123 : !> \param part ...
124 : !> \param box ...
125 : !> \param center ...
126 : !> \param delta ...
127 : !> \param npts ...
128 : !> \param n ...
129 : ! **************************************************************************************************
130 97518 : SUBROUTINE get_center(part, box, center, delta, npts, n)
131 :
132 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: part
133 : TYPE(cell_type), POINTER :: box
134 : INTEGER, DIMENSION(:, :), INTENT(OUT) :: center
135 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: delta
136 : INTEGER, DIMENSION(:), INTENT(IN) :: npts
137 : INTEGER, INTENT(IN) :: n
138 :
139 : INTEGER :: ipart, mp
140 : REAL(KIND=dp) :: rmp
141 : REAL(KIND=dp), DIMENSION(3) :: ca, gp, s
142 :
143 : ! The pbc algorithm is sensitive to numeric noise and compiler optimization because of ANINT.
144 : ! Therefore center and delta have to be computed simultaneously to ensure they are consistent.
145 390072 : mp = MAXVAL(npts(:))
146 97518 : rmp = REAL(mp, KIND=dp)
147 8507114 : DO ipart = 1, SIZE(part)
148 : ! compute the scaled coordinate of atom ipart
149 8409596 : CALL real_to_scaled(s, part(ipart)%r, box)
150 33638384 : s = s - ANINT(s)
151 : ! find the continuous ``grid'' point
152 33638384 : gp = REAL(npts, KIND=dp)*s
153 : ! find the closest grid point (on big grid)
154 8409596 : IF (MOD(n, 2) == 0) THEN
155 33326368 : center(:, ipart) = INT(gp + rmp) - mp
156 33326368 : ca(:) = REAL(center(:, ipart), KIND=dp) + 0.5_dp
157 : ELSE
158 312016 : center(:, ipart) = NINT(gp)
159 312016 : ca(:) = REAL(center(:, ipart), KIND=dp)
160 : END IF
161 33638384 : center(:, ipart) = center(:, ipart) + npts(:)/2
162 33638384 : center(:, ipart) = MODULO(center(:, ipart), npts(:))
163 33638384 : center(:, ipart) = center(:, ipart) - npts(:)/2
164 : ! find the distance vector
165 33735902 : delta(:, ipart) = gp - ca(:)
166 : END DO
167 :
168 97518 : END SUBROUTINE get_center
169 :
170 : END MODULE pme_tools
171 :
|
