(git:e7e05ae)
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 ! **************************************************************************************************
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  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 
402 END 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::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
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::integrate_gf_npbc
subroutine integrate_gf_npbc(grid, xdat, ydat, zdat, bo, zlb, zub, ylb, yub, xlb, xub, force)
...
Definition: mm_collocate_potential.F:91
pw_types
Definition: pw_types.F:24