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 Calculate the MM potential by collocating the primitive Gaussian
10 : !> functions (pgf)
11 : !> \par History
12 : !> 7.2004 created [tlaino]
13 : !> \author Teodoro Laino
14 : ! **************************************************************************************************
15 : MODULE mm_collocate_potential
16 : USE ao_util, ONLY: exp_radius
17 : USE cell_types, ONLY: cell_type
18 : USE cube_utils, ONLY: cube_info_type,&
19 : return_cube
20 : USE kinds, ONLY: dp
21 : USE pw_types, ONLY: pw_r3d_rs_type
22 : #include "./base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 : PRIVATE
26 :
27 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
28 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mm_collocate_potential'
29 :
30 : PUBLIC :: collocate_gf_rspace_NoPBC, &
31 : integrate_gf_rspace_NoPBC
32 : !***
33 : CONTAINS
34 :
35 : ! **************************************************************************************************
36 : !> \brief ...
37 : !> \param grid ...
38 : !> \param xdat ...
39 : !> \param ydat ...
40 : !> \param zdat ...
41 : !> \param bo1 ...
42 : !> \param bo2 ...
43 : !> \param zlb ...
44 : !> \param zub ...
45 : !> \param ylb ...
46 : !> \param yub ...
47 : !> \param xlb ...
48 : !> \param xub ...
49 : ! **************************************************************************************************
50 77897 : SUBROUTINE collocate_gf_npbc(grid, xdat, ydat, zdat, bo1, bo2, zlb, zub, ylb, yub, xlb, xub)
51 : USE kinds, ONLY: dp
52 : INTEGER, INTENT(IN) :: bo1(2, 3)
53 : REAL(dp), INTENT(INOUT) :: &
54 : grid(bo1(1, 1):bo1(2, 1), bo1(1, 2):bo1(2, 2), bo1(1, 3):bo1(2, 3))
55 : INTEGER, INTENT(IN) :: bo2(2, 3)
56 : REAL(dp), INTENT(IN) :: zdat(bo2(1, 3):bo2(2, 3)), &
57 : ydat(bo2(1, 2):bo2(2, 2)), &
58 : xdat(bo2(1, 1):bo2(2, 1))
59 : INTEGER, INTENT(IN) :: zlb, zub, ylb, yub, xlb, xub
60 :
61 : INTEGER :: ix, iy, iz
62 : REAL(dp) :: tmp1
63 :
64 900071 : DO iz = zlb, zub
65 13679903 : DO iy = ylb, yub
66 12779832 : tmp1 = zdat(iz)*ydat(iy)
67 309548839 : DO ix = xlb, xub
68 308726665 : grid(ix, iy, iz) = grid(ix, iy, iz) + xdat(ix)*tmp1
69 : END DO ! Loop on x
70 : END DO ! Loop on y
71 : END DO ! Loop on z
72 :
73 77897 : END SUBROUTINE collocate_gf_npbc
74 :
75 : ! **************************************************************************************************
76 : !> \brief ...
77 : !> \param grid ...
78 : !> \param xdat ...
79 : !> \param ydat ...
80 : !> \param zdat ...
81 : !> \param bo ...
82 : !> \param zlb ...
83 : !> \param zub ...
84 : !> \param ylb ...
85 : !> \param yub ...
86 : !> \param xlb ...
87 : !> \param xub ...
88 : !> \param force ...
89 : ! **************************************************************************************************
90 41125 : SUBROUTINE integrate_gf_npbc(grid, xdat, ydat, zdat, bo, zlb, zub, ylb, yub, xlb, xub, force)
91 : USE kinds, ONLY: dp
92 : INTEGER, INTENT(IN) :: bo(2, 3)
93 : REAL(dp), INTENT(IN) :: zdat(2, bo(1, 3):bo(2, 3)), &
94 : ydat(2, bo(1, 2):bo(2, 2)), &
95 : xdat(2, bo(1, 1):bo(2, 1))
96 : REAL(dp), INTENT(INOUT) :: grid(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3))
97 : INTEGER, INTENT(IN) :: zlb, zub, ylb, yub, xlb, xub
98 : REAL(dp), INTENT(INOUT) :: force(3)
99 :
100 : INTEGER :: ix, iy, iy2, iz
101 : REAL(dp) :: fx1, fx2, fyz1, fyz2, g1, g2, x1, x2
102 :
103 533632 : DO iz = zlb, zub
104 492507 : iy2 = HUGE(0)
105 : ! unroll by 2
106 492507 : DO iy = ylb, yub - 1, 2
107 3594214 : iy2 = iy + 1
108 3594214 : fx1 = 0.0_dp
109 3594214 : fyz1 = 0.0_dp
110 3594214 : fx2 = 0.0_dp
111 3594214 : fyz2 = 0.0_dp
112 85869874 : DO ix = xlb, xub
113 82275660 : g1 = grid(ix, iy, iz)
114 82275660 : g2 = grid(ix, iy2, iz)
115 82275660 : x1 = xdat(1, ix)
116 82275660 : x2 = xdat(2, ix)
117 82275660 : fyz1 = fyz1 + g1*x1
118 82275660 : fx1 = fx1 + g1*x2
119 82275660 : fyz2 = fyz2 + g2*x1
120 85869874 : fx2 = fx2 + g2*x2
121 : END DO ! Loop on x
122 3594214 : force(1) = force(1) + fx1*zdat(1, iz)*ydat(1, iy)
123 3594214 : force(2) = force(2) + fyz1*zdat(1, iz)*ydat(2, iy)
124 3594214 : force(3) = force(3) + fyz1*zdat(2, iz)*ydat(1, iy)
125 3594214 : force(1) = force(1) + fx2*zdat(1, iz)*ydat(1, iy2)
126 3594214 : force(2) = force(2) + fyz2*zdat(1, iz)*ydat(2, iy2)
127 3604763 : force(3) = force(3) + fyz2*zdat(2, iz)*ydat(1, iy2)
128 : END DO ! Loop on y
129 :
130 : ! cleanup loop: check if the last loop element has done
131 533632 : IF (iy2 .NE. yub) THEN
132 366844 : iy2 = yub
133 366844 : fx2 = 0.0_dp
134 366844 : fyz2 = 0.0_dp
135 5949279 : DO ix = xlb, xub
136 5582435 : g2 = grid(ix, iy2, iz)
137 5582435 : x1 = xdat(1, ix)
138 5582435 : x2 = xdat(2, ix)
139 5582435 : fyz2 = fyz2 + g2*x1
140 5949279 : fx2 = fx2 + g2*x2
141 : END DO ! Loop on x
142 366844 : force(1) = force(1) + fx2*zdat(1, iz)*ydat(1, iy2)
143 366844 : force(2) = force(2) + fyz2*zdat(1, iz)*ydat(2, iy2)
144 366844 : force(3) = force(3) + fyz2*zdat(2, iz)*ydat(1, iy2)
145 : END IF
146 :
147 : END DO ! Loop on z
148 :
149 41125 : END SUBROUTINE integrate_gf_npbc
150 :
151 : ! **************************************************************************************************
152 : !> \brief Main driver to collocate gaussian functions on grid
153 : !> without using periodic boundary conditions (NoPBC)
154 : !> \param zetp ...
155 : !> \param rp ...
156 : !> \param scale ...
157 : !> \param W ...
158 : !> \param pwgrid ...
159 : !> \param cube_info ...
160 : !> \param eps_mm_rspace ...
161 : !> \param xdat ...
162 : !> \param ydat ...
163 : !> \param zdat ...
164 : !> \param bo2 ...
165 : !> \param n_rep_real ...
166 : !> \param mm_cell ...
167 : !> \par History
168 : !> 07.2004 created [tlaino]
169 : !> \author Teodoro Laino
170 : ! **************************************************************************************************
171 24477 : SUBROUTINE collocate_gf_rspace_NoPBC(zetp, rp, scale, W, pwgrid, cube_info, &
172 : eps_mm_rspace, xdat, ydat, zdat, bo2, n_rep_real, mm_cell)
173 : REAL(KIND=dp), INTENT(IN) :: zetp
174 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rp
175 : REAL(KIND=dp), INTENT(IN) :: scale, W
176 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pwgrid
177 : TYPE(cube_info_type), INTENT(IN) :: cube_info
178 : REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace
179 : REAL(KIND=dp), DIMENSION(:), POINTER :: xdat, ydat, zdat
180 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo2
181 : INTEGER, DIMENSION(3), INTENT(IN) :: n_rep_real
182 : TYPE(cell_type), POINTER :: mm_cell
183 :
184 : INTEGER :: ig, ix, iy, iz, xlb, xub, ylb, yub, zlb, &
185 : zub
186 : INTEGER, DIMENSION(2, 3) :: bo, gbo
187 : INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ub_cube
188 24477 : INTEGER, DIMENSION(:), POINTER :: sphere_bounds
189 : REAL(KIND=dp) :: radius, rpg, xap, yap, zap
190 : REAL(KIND=dp), DIMENSION(3) :: dr, my_shift, rpl
191 24477 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: grid
192 :
193 48954 : radius = exp_radius(0, zetp, eps_mm_rspace, scale*W)
194 24477 : IF (radius .EQ. 0.0_dp) THEN
195 53 : RETURN
196 : END IF
197 :
198 : ! *** properties of the grid ***
199 24477 : rpl = rp
200 97908 : dr(:) = pwgrid%pw_grid%dr(:)
201 24477 : grid => pwgrid%array
202 244770 : bo = pwgrid%pw_grid%bounds_local
203 244770 : gbo = pwgrid%pw_grid%bounds
204 :
205 : ! *** get the sub grid properties for the given radius ***
206 24477 : CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
207 :
208 90492 : IF (ALL(n_rep_real == 0)) THEN
209 88020 : cubecenter(:) = FLOOR(rpl(:)/dr(:)) + gbo(1, :)
210 22005 : zub = MIN(bo(2, 3), cubecenter(3) + ub_cube(3))
211 22005 : zlb = MAX(bo(1, 3), cubecenter(3) + lb_cube(3))
212 22005 : yub = MIN(bo(2, 2), cubecenter(2) + ub_cube(2))
213 22005 : ylb = MAX(bo(1, 2), cubecenter(2) + lb_cube(2))
214 22005 : xub = MIN(bo(2, 1), cubecenter(1) + ub_cube(1))
215 22005 : xlb = MAX(bo(1, 1), cubecenter(1) + lb_cube(1))
216 22005 : IF (zlb .GT. zub .OR. ylb .GT. yub .OR. xlb .GT. xub) RETURN
217 465074 : DO ig = zlb, zub
218 443122 : rpg = REAL(ig - gbo(1, 3), dp)*dr(3) - rpl(3)
219 443122 : zap = EXP(-zetp*rpg**2)
220 465074 : zdat(ig) = scale*W*zap
221 : END DO
222 454154 : DO ig = ylb, yub
223 432202 : rpg = REAL(ig - gbo(1, 2), dp)*dr(2) - rpl(2)
224 432202 : yap = EXP(-zetp*rpg**2)
225 454154 : ydat(ig) = yap
226 : END DO
227 500305 : DO ig = xlb, xub
228 478353 : rpg = REAL(ig - gbo(1, 1), dp)*dr(1) - rpl(1)
229 478353 : xap = EXP(-zetp*rpg**2)
230 500305 : xdat(ig) = xap
231 : END DO
232 21952 : CALL collocate_gf_npbc(grid, xdat, ydat, zdat, bo, bo2, zlb, zub, ylb, yub, xlb, xub)
233 : ELSE
234 12432 : DO iz = -n_rep_real(3), n_rep_real(3)
235 9960 : my_shift(3) = mm_cell%hmat(3, 3)*REAL(iz, KIND=dp)
236 57624 : DO iy = -n_rep_real(2), n_rep_real(2)
237 45192 : my_shift(2) = mm_cell%hmat(2, 2)*REAL(iy, KIND=dp)
238 290616 : DO ix = -n_rep_real(1), n_rep_real(1)
239 235464 : my_shift(1) = mm_cell%hmat(1, 1)*REAL(ix, KIND=dp)
240 941856 : rpl = rp + my_shift(:)
241 941856 : cubecenter(:) = FLOOR(rpl(:)/dr(:)) + gbo(1, :)
242 235464 : zub = MIN(bo(2, 3), cubecenter(3) + ub_cube(3))
243 235464 : zlb = MAX(bo(1, 3), cubecenter(3) + lb_cube(3))
244 235464 : yub = MIN(bo(2, 2), cubecenter(2) + ub_cube(2))
245 235464 : ylb = MAX(bo(1, 2), cubecenter(2) + lb_cube(2))
246 235464 : xub = MIN(bo(2, 1), cubecenter(1) + ub_cube(1))
247 235464 : xlb = MAX(bo(1, 1), cubecenter(1) + lb_cube(1))
248 235464 : IF (zlb .GT. zub .OR. ylb .GT. yub .OR. xlb .GT. xub) CYCLE
249 434997 : DO ig = zlb, zub
250 379052 : rpg = REAL(ig - gbo(1, 3), dp)*dr(3) - rpl(3)
251 379052 : zap = EXP(-zetp*rpg**2)
252 434997 : zdat(ig) = scale*W*zap
253 : END DO
254 435469 : DO ig = ylb, yub
255 379524 : rpg = REAL(ig - gbo(1, 2), dp)*dr(2) - rpl(2)
256 379524 : yap = EXP(-zetp*rpg**2)
257 435469 : ydat(ig) = yap
258 : END DO
259 435757 : DO ig = xlb, xub
260 379812 : rpg = REAL(ig - gbo(1, 1), dp)*dr(1) - rpl(1)
261 379812 : xap = EXP(-zetp*rpg**2)
262 435757 : xdat(ig) = xap
263 : END DO
264 101137 : CALL collocate_gf_npbc(grid, xdat, ydat, zdat, bo, bo2, zlb, zub, ylb, yub, xlb, xub)
265 : END DO
266 : END DO
267 : END DO
268 : END IF
269 :
270 24477 : END SUBROUTINE collocate_gf_rspace_NoPBC
271 :
272 : ! **************************************************************************************************
273 : !> \brief Main driver to integrate gaussian functions on a grid function
274 : !> without using periodic boundary conditions (NoPBC)
275 : !> Computes Forces.
276 : !> \param zetp ...
277 : !> \param rp ...
278 : !> \param scale ...
279 : !> \param W ...
280 : !> \param pwgrid ...
281 : !> \param cube_info ...
282 : !> \param eps_mm_rspace ...
283 : !> \param xdat ...
284 : !> \param ydat ...
285 : !> \param zdat ...
286 : !> \param bo ...
287 : !> \param force ...
288 : !> \param n_rep_real ...
289 : !> \param mm_cell ...
290 : !> \par History
291 : !> 07.2004 created [tlaino]
292 : !> \author Teodoro Laino
293 : ! **************************************************************************************************
294 13588 : SUBROUTINE integrate_gf_rspace_NoPBC(zetp, rp, scale, W, pwgrid, cube_info, &
295 13588 : eps_mm_rspace, xdat, ydat, zdat, bo, force, n_rep_real, mm_cell)
296 : REAL(KIND=dp), INTENT(IN) :: zetp
297 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rp
298 : REAL(KIND=dp), INTENT(IN) :: scale, W
299 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pwgrid
300 : TYPE(cube_info_type), INTENT(IN) :: cube_info
301 : REAL(KIND=dp), INTENT(IN) :: eps_mm_rspace
302 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo
303 : REAL(KIND=dp), DIMENSION(2, bo(1, 3):bo(2, 3)) :: zdat
304 : REAL(KIND=dp), DIMENSION(2, bo(1, 2):bo(2, 2)) :: ydat
305 : REAL(KIND=dp), DIMENSION(2, bo(1, 1):bo(2, 1)) :: xdat
306 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: force
307 : INTEGER, DIMENSION(3), INTENT(IN) :: n_rep_real
308 : TYPE(cell_type), POINTER :: mm_cell
309 :
310 : INTEGER :: ig, ix, iy, iz, xlb, xub, ylb, yub, zlb, &
311 : zub
312 : INTEGER, DIMENSION(2, 3) :: gbo
313 : INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ub_cube
314 13588 : INTEGER, DIMENSION(:), POINTER :: sphere_bounds
315 : REAL(KIND=dp) :: radius, rpg, xap, yap, zap
316 : REAL(KIND=dp), DIMENSION(3) :: dr, my_shift, rpl
317 13588 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: grid
318 :
319 13588 : force = 0.0_dp
320 27176 : radius = exp_radius(0, zetp, eps_mm_rspace, scale*W)
321 13641 : IF (radius .EQ. 0.0_dp) RETURN
322 :
323 : ! *** properties of the grid ***
324 13588 : rpl = rp
325 54352 : dr(:) = pwgrid%pw_grid%dr(:)
326 13588 : grid => pwgrid%array
327 135880 : gbo = pwgrid%pw_grid%bounds
328 :
329 : ! *** get the sub grid properties for the given radius ***
330 13588 : CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
331 :
332 50932 : IF (ALL(n_rep_real == 0)) THEN
333 49792 : cubecenter(:) = FLOOR(rpl(:)/dr(:)) + gbo(1, :)
334 12448 : zub = MIN(bo(2, 3), cubecenter(3) + ub_cube(3))
335 12448 : zlb = MAX(bo(1, 3), cubecenter(3) + lb_cube(3))
336 12448 : yub = MIN(bo(2, 2), cubecenter(2) + ub_cube(2))
337 12448 : ylb = MAX(bo(1, 2), cubecenter(2) + lb_cube(2))
338 12448 : xub = MIN(bo(2, 1), cubecenter(1) + ub_cube(1))
339 12448 : xlb = MAX(bo(1, 1), cubecenter(1) + lb_cube(1))
340 12448 : IF (zlb .GT. zub .OR. ylb .GT. yub .OR. xlb .GT. xub) RETURN
341 269039 : DO ig = zlb, zub
342 256644 : rpg = REAL(ig - gbo(1, 3), dp)*dr(3) - rpl(3)
343 256644 : zap = EXP(-zetp*rpg**2)
344 256644 : zdat(1, ig) = scale*W*zap
345 269039 : zdat(2, ig) = rpg*zdat(1, ig)*zetp*2.0_dp
346 : END DO
347 254492 : DO ig = ylb, yub
348 242097 : rpg = REAL(ig - gbo(1, 2), dp)*dr(2) - rpl(2)
349 242097 : yap = EXP(-zetp*rpg**2)
350 242097 : ydat(1, ig) = yap
351 254492 : ydat(2, ig) = rpg*ydat(1, ig)*zetp*2.0_dp
352 : END DO
353 284060 : DO ig = xlb, xub
354 271665 : rpg = REAL(ig - gbo(1, 1), dp)*dr(1) - rpl(1)
355 271665 : xap = EXP(-zetp*rpg**2)
356 271665 : xdat(1, ig) = xap
357 284060 : xdat(2, ig) = rpg*xdat(1, ig)*zetp*2.0_dp
358 : END DO
359 12395 : CALL integrate_gf_npbc(grid, xdat, ydat, zdat, bo, zlb, zub, ylb, yub, xlb, xub, force)
360 : ELSE
361 6456 : DO iz = -n_rep_real(3), n_rep_real(3)
362 5316 : my_shift(3) = mm_cell%hmat(3, 3)*REAL(iz, KIND=dp)
363 31884 : DO iy = -n_rep_real(2), n_rep_real(2)
364 25428 : my_shift(2) = mm_cell%hmat(2, 2)*REAL(iy, KIND=dp)
365 154428 : DO ix = -n_rep_real(1), n_rep_real(1)
366 123684 : my_shift(1) = mm_cell%hmat(1, 1)*REAL(ix, KIND=dp)
367 494736 : rpl = rp + my_shift(:)
368 494736 : cubecenter(:) = FLOOR(rpl(:)/dr(:)) + gbo(1, :)
369 123684 : zub = MIN(bo(2, 3), cubecenter(3) + ub_cube(3))
370 123684 : zlb = MAX(bo(1, 3), cubecenter(3) + lb_cube(3))
371 123684 : yub = MIN(bo(2, 2), cubecenter(2) + ub_cube(2))
372 123684 : ylb = MAX(bo(1, 2), cubecenter(2) + lb_cube(2))
373 123684 : xub = MIN(bo(2, 1), cubecenter(1) + ub_cube(1))
374 123684 : xlb = MAX(bo(1, 1), cubecenter(1) + lb_cube(1))
375 123684 : IF (zlb .GT. zub .OR. ylb .GT. yub .OR. xlb .GT. xub) CYCLE
376 264593 : DO ig = zlb, zub
377 235863 : rpg = REAL(ig - gbo(1, 3), dp)*dr(3) - rpl(3)
378 235863 : zap = EXP(-zetp*rpg**2)
379 235863 : zdat(1, ig) = scale*W*zap
380 264593 : zdat(2, ig) = rpg*zdat(1, ig)*zetp*2.0_dp
381 : END DO
382 264897 : DO ig = ylb, yub
383 236167 : rpg = REAL(ig - gbo(1, 2), dp)*dr(2) - rpl(2)
384 236167 : yap = EXP(-zetp*rpg**2)
385 236167 : ydat(1, ig) = yap
386 264897 : ydat(2, ig) = rpg*ydat(1, ig)*zetp*2.0_dp
387 : END DO
388 264849 : DO ig = xlb, xub
389 236119 : rpg = REAL(ig - gbo(1, 1), dp)*dr(1) - rpl(1)
390 236119 : xap = EXP(-zetp*rpg**2)
391 236119 : xdat(1, ig) = xap
392 264849 : xdat(2, ig) = rpg*xdat(1, ig)*zetp*2.0_dp
393 : END DO
394 : CALL integrate_gf_npbc(grid, xdat, ydat, zdat, bo, &
395 54158 : zlb, zub, ylb, yub, xlb, xub, force)
396 : END DO
397 : END DO
398 : END DO
399 : END IF
400 13588 : END SUBROUTINE integrate_gf_rspace_NoPBC
401 :
402 : END MODULE mm_collocate_potential
|
