(git:374b731)
Loading...
Searching...
No Matches
mm_collocate_potential.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 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! **************************************************************************************************
15MODULE 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!***
33CONTAINS
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 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 DO iz = zlb, zub
65 DO iy = ylb, yub
66 tmp1 = zdat(iz)*ydat(iy)
67 DO ix = xlb, xub
68 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 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 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 DO iz = zlb, zub
104 iy2 = huge(0)
105 ! unroll by 2
106 DO iy = ylb, yub - 1, 2
107 iy2 = iy + 1
108 fx1 = 0.0_dp
109 fyz1 = 0.0_dp
110 fx2 = 0.0_dp
111 fyz2 = 0.0_dp
112 DO ix = xlb, xub
113 g1 = grid(ix, iy, iz)
114 g2 = grid(ix, iy2, iz)
115 x1 = xdat(1, ix)
116 x2 = xdat(2, ix)
117 fyz1 = fyz1 + g1*x1
118 fx1 = fx1 + g1*x2
119 fyz2 = fyz2 + g2*x1
120 fx2 = fx2 + g2*x2
121 END DO ! Loop on x
122 force(1) = force(1) + fx1*zdat(1, iz)*ydat(1, iy)
123 force(2) = force(2) + fyz1*zdat(1, iz)*ydat(2, iy)
124 force(3) = force(3) + fyz1*zdat(2, iz)*ydat(1, iy)
125 force(1) = force(1) + fx2*zdat(1, iz)*ydat(1, iy2)
126 force(2) = force(2) + fyz2*zdat(1, iz)*ydat(2, iy2)
127 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 IF (iy2 .NE. yub) THEN
132 iy2 = yub
133 fx2 = 0.0_dp
134 fyz2 = 0.0_dp
135 DO ix = xlb, xub
136 g2 = grid(ix, iy2, iz)
137 x1 = xdat(1, ix)
138 x2 = xdat(2, ix)
139 fyz2 = fyz2 + g2*x1
140 fx2 = fx2 + g2*x2
141 END DO ! Loop on x
142 force(1) = force(1) + fx2*zdat(1, iz)*ydat(1, iy2)
143 force(2) = force(2) + fyz2*zdat(1, iz)*ydat(2, iy2)
144 force(3) = force(3) + fyz2*zdat(2, iz)*ydat(1, iy2)
145 END IF
146
147 END DO ! Loop on z
148
149 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 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 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 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
192
193 radius = exp_radius(0, zetp, eps_mm_rspace, scale*w)
194 IF (radius .EQ. 0.0_dp) THEN
195 RETURN
196 END IF
197
198! *** properties of the grid ***
199 rpl = rp
200 dr(:) = pwgrid%pw_grid%dr(:)
201 grid => pwgrid%array
202 bo = pwgrid%pw_grid%bounds_local
203 gbo = pwgrid%pw_grid%bounds
204
205! *** get the sub grid properties for the given radius ***
206 CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
207
208 IF (all(n_rep_real == 0)) THEN
209 cubecenter(:) = floor(rpl(:)/dr(:)) + gbo(1, :)
210 zub = min(bo(2, 3), cubecenter(3) + ub_cube(3))
211 zlb = max(bo(1, 3), cubecenter(3) + lb_cube(3))
212 yub = min(bo(2, 2), cubecenter(2) + ub_cube(2))
213 ylb = max(bo(1, 2), cubecenter(2) + lb_cube(2))
214 xub = min(bo(2, 1), cubecenter(1) + ub_cube(1))
215 xlb = max(bo(1, 1), cubecenter(1) + lb_cube(1))
216 IF (zlb .GT. zub .OR. ylb .GT. yub .OR. xlb .GT. xub) RETURN
217 DO ig = zlb, zub
218 rpg = real(ig - gbo(1, 3), dp)*dr(3) - rpl(3)
219 zap = exp(-zetp*rpg**2)
220 zdat(ig) = scale*w*zap
221 END DO
222 DO ig = ylb, yub
223 rpg = real(ig - gbo(1, 2), dp)*dr(2) - rpl(2)
224 yap = exp(-zetp*rpg**2)
225 ydat(ig) = yap
226 END DO
227 DO ig = xlb, xub
228 rpg = real(ig - gbo(1, 1), dp)*dr(1) - rpl(1)
229 xap = exp(-zetp*rpg**2)
230 xdat(ig) = xap
231 END DO
232 CALL collocate_gf_npbc(grid, xdat, ydat, zdat, bo, bo2, zlb, zub, ylb, yub, xlb, xub)
233 ELSE
234 DO iz = -n_rep_real(3), n_rep_real(3)
235 my_shift(3) = mm_cell%hmat(3, 3)*real(iz, kind=dp)
236 DO iy = -n_rep_real(2), n_rep_real(2)
237 my_shift(2) = mm_cell%hmat(2, 2)*real(iy, kind=dp)
238 DO ix = -n_rep_real(1), n_rep_real(1)
239 my_shift(1) = mm_cell%hmat(1, 1)*real(ix, kind=dp)
240 rpl = rp + my_shift(:)
241 cubecenter(:) = floor(rpl(:)/dr(:)) + gbo(1, :)
242 zub = min(bo(2, 3), cubecenter(3) + ub_cube(3))
243 zlb = max(bo(1, 3), cubecenter(3) + lb_cube(3))
244 yub = min(bo(2, 2), cubecenter(2) + ub_cube(2))
245 ylb = max(bo(1, 2), cubecenter(2) + lb_cube(2))
246 xub = min(bo(2, 1), cubecenter(1) + ub_cube(1))
247 xlb = max(bo(1, 1), cubecenter(1) + lb_cube(1))
248 IF (zlb .GT. zub .OR. ylb .GT. yub .OR. xlb .GT. xub) cycle
249 DO ig = zlb, zub
250 rpg = real(ig - gbo(1, 3), dp)*dr(3) - rpl(3)
251 zap = exp(-zetp*rpg**2)
252 zdat(ig) = scale*w*zap
253 END DO
254 DO ig = ylb, yub
255 rpg = real(ig - gbo(1, 2), dp)*dr(2) - rpl(2)
256 yap = exp(-zetp*rpg**2)
257 ydat(ig) = yap
258 END DO
259 DO ig = xlb, xub
260 rpg = real(ig - gbo(1, 1), dp)*dr(1) - rpl(1)
261 xap = exp(-zetp*rpg**2)
262 xdat(ig) = xap
263 END DO
264 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 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 SUBROUTINE integrate_gf_rspace_nopbc(zetp, rp, scale, W, pwgrid, cube_info, &
295 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 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 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
318
319 force = 0.0_dp
320 radius = exp_radius(0, zetp, eps_mm_rspace, scale*w)
321 IF (radius .EQ. 0.0_dp) RETURN
322
323! *** properties of the grid ***
324 rpl = rp
325 dr(:) = pwgrid%pw_grid%dr(:)
326 grid => pwgrid%array
327 gbo = pwgrid%pw_grid%bounds
328
329! *** get the sub grid properties for the given radius ***
330 CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
331
332 IF (all(n_rep_real == 0)) THEN
333 cubecenter(:) = floor(rpl(:)/dr(:)) + gbo(1, :)
334 zub = min(bo(2, 3), cubecenter(3) + ub_cube(3))
335 zlb = max(bo(1, 3), cubecenter(3) + lb_cube(3))
336 yub = min(bo(2, 2), cubecenter(2) + ub_cube(2))
337 ylb = max(bo(1, 2), cubecenter(2) + lb_cube(2))
338 xub = min(bo(2, 1), cubecenter(1) + ub_cube(1))
339 xlb = max(bo(1, 1), cubecenter(1) + lb_cube(1))
340 IF (zlb .GT. zub .OR. ylb .GT. yub .OR. xlb .GT. xub) RETURN
341 DO ig = zlb, zub
342 rpg = real(ig - gbo(1, 3), dp)*dr(3) - rpl(3)
343 zap = exp(-zetp*rpg**2)
344 zdat(1, ig) = scale*w*zap
345 zdat(2, ig) = rpg*zdat(1, ig)*zetp*2.0_dp
346 END DO
347 DO ig = ylb, yub
348 rpg = real(ig - gbo(1, 2), dp)*dr(2) - rpl(2)
349 yap = exp(-zetp*rpg**2)
350 ydat(1, ig) = yap
351 ydat(2, ig) = rpg*ydat(1, ig)*zetp*2.0_dp
352 END DO
353 DO ig = xlb, xub
354 rpg = real(ig - gbo(1, 1), dp)*dr(1) - rpl(1)
355 xap = exp(-zetp*rpg**2)
356 xdat(1, ig) = xap
357 xdat(2, ig) = rpg*xdat(1, ig)*zetp*2.0_dp
358 END DO
359 CALL integrate_gf_npbc(grid, xdat, ydat, zdat, bo, zlb, zub, ylb, yub, xlb, xub, force)
360 ELSE
361 DO iz = -n_rep_real(3), n_rep_real(3)
362 my_shift(3) = mm_cell%hmat(3, 3)*real(iz, kind=dp)
363 DO iy = -n_rep_real(2), n_rep_real(2)
364 my_shift(2) = mm_cell%hmat(2, 2)*real(iy, kind=dp)
365 DO ix = -n_rep_real(1), n_rep_real(1)
366 my_shift(1) = mm_cell%hmat(1, 1)*real(ix, kind=dp)
367 rpl = rp + my_shift(:)
368 cubecenter(:) = floor(rpl(:)/dr(:)) + gbo(1, :)
369 zub = min(bo(2, 3), cubecenter(3) + ub_cube(3))
370 zlb = max(bo(1, 3), cubecenter(3) + lb_cube(3))
371 yub = min(bo(2, 2), cubecenter(2) + ub_cube(2))
372 ylb = max(bo(1, 2), cubecenter(2) + lb_cube(2))
373 xub = min(bo(2, 1), cubecenter(1) + ub_cube(1))
374 xlb = max(bo(1, 1), cubecenter(1) + lb_cube(1))
375 IF (zlb .GT. zub .OR. ylb .GT. yub .OR. xlb .GT. xub) cycle
376 DO ig = zlb, zub
377 rpg = real(ig - gbo(1, 3), dp)*dr(3) - rpl(3)
378 zap = exp(-zetp*rpg**2)
379 zdat(1, ig) = scale*w*zap
380 zdat(2, ig) = rpg*zdat(1, ig)*zetp*2.0_dp
381 END DO
382 DO ig = ylb, yub
383 rpg = real(ig - gbo(1, 2), dp)*dr(2) - rpl(2)
384 yap = exp(-zetp*rpg**2)
385 ydat(1, ig) = yap
386 ydat(2, ig) = rpg*ydat(1, ig)*zetp*2.0_dp
387 END DO
388 DO ig = xlb, xub
389 rpg = real(ig - gbo(1, 1), dp)*dr(1) - rpl(1)
390 xap = exp(-zetp*rpg**2)
391 xdat(1, ig) = xap
392 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 zlb, zub, ylb, yub, xlb, xub, force)
396 END DO
397 END DO
398 END DO
399 END IF
400 END SUBROUTINE integrate_gf_rspace_nopbc
401
402END MODULE mm_collocate_potential
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius
real(kind=dp) function, public exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow)
The radius of a primitive Gaussian function for a given threshold is calculated. g(r) = prefactor*r**...
Definition ao_util.F:96
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cube_utils
for a given dr()/dh(r) this will provide the bounds to be used if one wants to go over a sphere-subre...
Definition cube_utils.F:18
cube_utils::return_cube
subroutine, public return_cube(info, radius, lb_cube, ub_cube, sphere_bounds)
...
Definition cube_utils.F:140
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mm_collocate_potential
Calculate the MM potential by collocating the primitive Gaussian functions (pgf)
Definition mm_collocate_potential.F:15
mm_collocate_potential::integrate_gf_rspace_nopbc
subroutine, public integrate_gf_rspace_nopbc(zetp, rp, scale, w, pwgrid, cube_info, eps_mm_rspace, xdat, ydat, zdat, bo, force, n_rep_real, mm_cell)
Main driver to integrate gaussian functions on a grid function without using periodic boundary condit...
Definition mm_collocate_potential.F:296
mm_collocate_potential::collocate_gf_rspace_nopbc
subroutine, public collocate_gf_rspace_nopbc(zetp, rp, scale, w, pwgrid, cube_info, eps_mm_rspace, xdat, ydat, zdat, bo2, n_rep_real, mm_cell)
Main driver to collocate gaussian functions on grid without using periodic boundary conditions (NoPBC...
Definition mm_collocate_potential.F:173
pw_types
Definition pw_types.F:24
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cube_utils::cube_info_type
Definition cube_utils.F:38
pw_types::pw_r3d_rs_type
Definition pw_types.F:62