(git:ccc2433)
dielectric_methods.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 methods for evaluating the dielectric constant
10 !> \par History
11 !> 06.2014 created [Hossein Bani-Hashemian]
12 !> \author Mohammad Hossein Bani-Hashemian
13 ! **************************************************************************************************
14 MODULE dielectric_methods
15 
16  USE dct, ONLY: pw_expand
17  USE dielectric_types, ONLY: &
18  derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft, derivative_fft_use_deps, &
19  derivative_fft_use_drho, dielectric_parameters, dielectric_type, rho_dependent, &
20  spatially_dependent, spatially_rho_dependent
21  USE kinds, ONLY: dp
22  USE mathconstants, ONLY: twopi
23  USE pw_grid_types, ONLY: pw_grid_type
24  USE pw_methods, ONLY: pw_axpy, &
25  pw_set, pw_copy, &
26  pw_derive, &
27  pw_transfer, &
28  pw_zero
29  USE pw_pool_types, ONLY: pw_pool_create, &
30  pw_pool_release, &
31  pw_pool_type
32  USE pw_types, ONLY: &
33  pw_c1d_gs_type, &
34  pw_r3d_rs_type
35  USE realspace_grid_types, ONLY: realspace_grid_type
36  USE rs_methods, ONLY: derive_fdm_cd3, &
37  derive_fdm_cd5, &
38  derive_fdm_cd7, &
39  pw_mollifier, &
40  setup_grid_axes
41 #include "../base/base_uses.f90"
42 
43  IMPLICIT NONE
44  PRIVATE
45  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dielectric_methods'
46 
47  PUBLIC :: dielectric_create, dielectric_compute, derive_fft
48 
49  INTERFACE dielectric_compute
50  MODULE PROCEDURE dielectric_compute_periodic_r3d_rs, &
51  dielectric_compute_neumann_r3d_rs
52  MODULE PROCEDURE dielectric_compute_periodic_c1d_gs, &
53  dielectric_compute_neumann_c1d_gs
54  END INTERFACE dielectric_compute
55 
56 CONTAINS
57 
58 ! **************************************************************************************************
59 !> \brief allocates memory for a dielectric data type
60 !> \param dielectric the dielectric data type to be allocated
61 !> \param pw_pool pool of pw grid
62 !> \param dielectric_params dielectric parameters read from input file
63 !> \par History
64 !> 06.2014 created [Hossein Bani-Hashemian]
65 !> \author Mohammad Hossein Bani-Hashemian
66 ! **************************************************************************************************
67  SUBROUTINE dielectric_create(dielectric, pw_pool, dielectric_params)
68  TYPE(dielectric_type), INTENT(INOUT), POINTER :: dielectric
69  TYPE(pw_pool_type), POINTER :: pw_pool
70  TYPE(dielectric_parameters), INTENT(IN) :: dielectric_params
71 
72  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_create'
73 
74  INTEGER :: handle, i
75 
76  CALL timeset(routinen, handle)
77 
78  IF (.NOT. ASSOCIATED(dielectric)) THEN
79  ALLOCATE (dielectric)
80  NULLIFY (dielectric%eps)
81  NULLIFY (dielectric%deps_drho)
82  ALLOCATE (dielectric%eps, dielectric%deps_drho)
83  CALL pw_pool%create_pw(dielectric%eps)
84  CALL pw_pool%create_pw(dielectric%deps_drho)
85  CALL pw_set(dielectric%eps, 1.0_dp)
86  CALL pw_zero(dielectric%deps_drho)
87  DO i = 1, 3
88  CALL pw_pool%create_pw(dielectric%dln_eps(i))
89  CALL pw_zero(dielectric%dln_eps(i))
90  END DO
91  dielectric%params = dielectric_params
92  dielectric%params%times_called = 0
93  END IF
94 
95  CALL timestop(handle)
96 
97  END SUBROUTINE dielectric_create
98 
99 ! **************************************************************************************************
100 !> \brief evaluates the dielectric constant
101 !> \param dielectric the dielectric data type to be initialized
102 !> \param diel_rs_grid real space grid for finite difference derivative
103 !> \param pw_pool pool of plane wave grid
104 !> \param rho electronic density
105 !> \param rho_core core density
106 !> \par History
107 !> 06.2014 created [Hossein Bani-Hashemian]
108 !> 12.2014 added finite difference derivatives [Hossein Bani-Hashemian]
109 !> 07.2015 density-independent dielectric regions [Hossein Bani-Hashemian]
110 !> \author Mohammad Hossein Bani-Hashemian
111 ! **************************************************************************************************
112  SUBROUTINE dielectric_compute_periodic_c1d_gs (dielectric, diel_rs_grid, pw_pool, rho, rho_core)
113 
114  TYPE(dielectric_type), INTENT(INOUT), POINTER :: dielectric
115  TYPE(realspace_grid_type), POINTER :: diel_rs_grid
116  TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
117  TYPE(pw_c1d_gs_type), INTENT(IN) :: rho
118  TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core
119 
120  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_compute_periodic'
121  REAL(dp), PARAMETER :: small_value = epsilon(1.0_dp)
122 
123  INTEGER :: derivative_method, dielec_functiontype, &
124  handle, i, idir, j, k, times_called
125  INTEGER, DIMENSION(3) :: lb, ub
126  REAL(dp) :: eps0, rho_max, rho_min
127  TYPE(pw_r3d_rs_type) :: ln_eps, rho_elec_rs
128  TYPE(pw_r3d_rs_type) :: rho_core_rs
129  TYPE(pw_r3d_rs_type), DIMENSION(3) :: deps, drho
130 
131  CALL timeset(routinen, handle)
132 
133  rho_min = dielectric%params%rho_min
134  rho_max = dielectric%params%rho_max
135  eps0 = dielectric%params%eps0
136  derivative_method = dielectric%params%derivative_method
137  times_called = dielectric%params%times_called
138  dielec_functiontype = dielectric%params%dielec_functiontype
139 
140  IF (.NOT. dielec_functiontype .EQ. rho_dependent .AND. &
141  ((derivative_method .EQ. derivative_fft_use_deps) .OR. &
142  (derivative_method .EQ. derivative_fft_use_deps))) THEN
143  CALL cp_abort(__location__, &
144  "The specified derivative method is not compatible with the type of "// &
145  "the dielectric constant function.")
146  END IF
147 
148  CALL pw_pool%create_pw(rho_elec_rs)
149 
150  ! for evaluating epsilon make sure rho is in the real space
151  CALL pw_transfer(rho, rho_elec_rs)
152 
153  IF (PRESENT(rho_core)) THEN
154  ! make sure rho_core is in the real space
155  CALL pw_pool%create_pw(rho_core_rs)
156  CALL pw_transfer(rho_core, rho_core_rs)
157  IF (dielectric%params%dielec_core_correction) THEN
158  ! use (rho_elec - rho_core) to compute dielectric to avoid obtaining spurious
159  ! epsilon in the core region
160  CALL pw_axpy(rho_core_rs, rho_elec_rs, -2.0_dp)
161  ELSE
162  CALL pw_axpy(rho_core_rs, rho_elec_rs, -1.0_dp)
163  END IF
164  CALL pw_pool%give_back_pw(rho_core_rs)
165  ELSE
166  cpabort("For dielectric constant larger than 1, rho_core has to be present.")
167  END IF
168 
169 ! calculate the dielectric constant
170  SELECT CASE (dielec_functiontype)
171  CASE (rho_dependent)
172  CALL dielectric_constant_sccs(rho_elec_rs, dielectric%eps, dielectric%deps_drho, &
173  eps0, rho_max, rho_min)
174  CASE (spatially_dependent)
175  IF (times_called .EQ. 0) THEN
176  CALL dielectric_constant_spatially_dependent(dielectric%eps, pw_pool, dielectric%params)
177  END IF
178  CASE (spatially_rho_dependent)
179  CALL dielectric_constant_spatially_rho_dependent(rho_elec_rs, dielectric%eps, &
180  dielectric%deps_drho, pw_pool, dielectric%params)
181  END SELECT
182 
183 ! derivatives
184  IF ((dielec_functiontype .EQ. rho_dependent) .OR. &
185  (dielec_functiontype .EQ. spatially_rho_dependent) .OR. &
186  ((dielec_functiontype .EQ. spatially_dependent) .AND. times_called .EQ. 0)) THEN
187  SELECT CASE (derivative_method)
188  CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
189  CALL pw_pool%create_pw(ln_eps)
190  ln_eps%array = log(dielectric%eps%array)
191  CASE (derivative_fft_use_deps)
192  DO i = 1, 3
193  CALL pw_pool%create_pw(deps(i))
194  CALL pw_zero(deps(i))
195  END DO
196  CASE (derivative_fft_use_drho)
197  DO i = 1, 3
198  CALL pw_pool%create_pw(deps(i))
199  CALL pw_pool%create_pw(drho(i))
200  CALL pw_zero(deps(i))
201  CALL pw_zero(drho(i))
202  END DO
203  END SELECT
204 
205  SELECT CASE (derivative_method)
206  CASE (derivative_cd3)
207  CALL derive_fdm_cd3(ln_eps, dielectric%dln_eps, diel_rs_grid)
208  CASE (derivative_cd5)
209  CALL derive_fdm_cd5(ln_eps, dielectric%dln_eps, diel_rs_grid)
210  CASE (derivative_cd7)
211  CALL derive_fdm_cd7(ln_eps, dielectric%dln_eps, diel_rs_grid)
212  CASE (derivative_fft)
213  CALL derive_fft(ln_eps, dielectric%dln_eps, pw_pool)
214  CASE (derivative_fft_use_deps)
215  ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
216  CALL derive_fft(dielectric%eps, deps, pw_pool)
217 
218  lb(1:3) = pw_pool%pw_grid%bounds_local(1, 1:3)
219  ub(1:3) = pw_pool%pw_grid%bounds_local(2, 1:3)
220  ! damp oscillations cuased by fft-based derivative in the region
221  ! where electron density is zero
222  DO idir = 1, 3
223  DO k = lb(3), ub(3)
224  DO j = lb(2), ub(2)
225  DO i = lb(1), ub(1)
226  IF (abs(dielectric%deps_drho%array(i, j, k)) .LE. small_value) THEN
227  deps(idir)%array(i, j, k) = 0.0_dp
228  END IF
229  END DO
230  END DO
231  END DO
232  dielectric%dln_eps(idir)%array = deps(idir)%array/dielectric%eps%array
233  END DO
234  CASE (derivative_fft_use_drho)
235  ! \Nabla \eps = \Nabla \rho \cdot \frac{\partial \eps}{\partial \rho}
236  ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
237  CALL derive_fft(rho_elec_rs, drho, pw_pool)
238  DO i = 1, 3
239  deps(i)%array = drho(i)%array*dielectric%deps_drho%array
240  dielectric%dln_eps(i)%array = deps(i)%array/dielectric%eps%array
241  END DO
242  END SELECT
243 
244  SELECT CASE (derivative_method)
245  CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
246  CALL pw_pool%give_back_pw(ln_eps)
247  CASE (derivative_fft_use_deps)
248  DO i = 1, 3
249  CALL pw_pool%give_back_pw(deps(i))
250  END DO
251  CASE (derivative_fft_use_drho)
252  DO i = 1, 3
253  CALL pw_pool%give_back_pw(drho(i))
254  CALL pw_pool%give_back_pw(deps(i))
255  END DO
256  END SELECT
257  END IF
258 
259  CALL pw_pool%give_back_pw(rho_elec_rs)
260 
261  dielectric%params%times_called = dielectric%params%times_called + 1
262 
263  CALL timestop(handle)
264 
265  END SUBROUTINE dielectric_compute_periodic_c1d_gs
266 
267 ! **************************************************************************************************
268 !> \brief evaluates the dielectric constant for non-periodic (Neumann-type)
269 !> boundaries
270 !> \param dielectric the dielectric data type to be initialized
271 !> \param diel_rs_grid real space grid for finite difference derivative
272 !> \param pw_pool_orig pool of plane wave grid
273 !> \param dct_pw_grid ...
274 !> \param neumann_directions ...
275 !> \param recv_msgs_bnds bounds of the messages to be received (pw_expand)
276 !> \param dests_expand list of the destination processes (pw_expand)
277 !> \param srcs_expand list of the source processes (pw_expand)
278 !> \param flipg_stat flipping status for the received data chunks (pw_expand)
279 !> \param bounds_shftd bounds of the original grid shifted to have g0 in the
280 !> middle of the cell (pw_expand)
281 !> \param rho electronic density
282 !> \param rho_core core density
283 !> \par History
284 !> 11.2015 created [Hossein Bani-Hashemian]
285 !> \author Mohammad Hossein Bani-Hashemian
286 ! **************************************************************************************************
287  SUBROUTINE dielectric_compute_neumann_c1d_gs (dielectric, diel_rs_grid, pw_pool_orig, &
288  dct_pw_grid, neumann_directions, recv_msgs_bnds, &
289  dests_expand, srcs_expand, flipg_stat, bounds_shftd, rho, rho_core)
290 
291  TYPE(dielectric_type), INTENT(INOUT), POINTER :: dielectric
292  TYPE(realspace_grid_type), POINTER :: diel_rs_grid
293  TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool_orig
294  TYPE(pw_grid_type), INTENT(IN), POINTER :: dct_pw_grid
295  INTEGER, INTENT(IN) :: neumann_directions
296  INTEGER, DIMENSION(:, :, :), INTENT(IN), POINTER :: recv_msgs_bnds
297  INTEGER, DIMENSION(:), INTENT(IN), POINTER :: dests_expand, srcs_expand, flipg_stat
298  INTEGER, DIMENSION(2, 3), INTENT(IN) :: bounds_shftd
299  TYPE(pw_c1d_gs_type), INTENT(IN) :: rho
300  TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core
301 
302  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_compute_neumann'
303  REAL(dp), PARAMETER :: small_value = epsilon(1.0_dp)
304 
305  INTEGER :: derivative_method, dielec_functiontype, &
306  handle, i, idir, j, k, times_called
307  INTEGER, DIMENSION(3) :: lb, ub
308  REAL(dp) :: eps0, rho_max, rho_min
309  TYPE(pw_pool_type), POINTER :: pw_pool_xpndd
310  TYPE(pw_r3d_rs_type) :: ln_eps, rho_core_rs, rho_core_rs_xpndd, &
311  rho_elec_rs, rho_elec_rs_xpndd
312  TYPE(pw_r3d_rs_type), DIMENSION(3) :: deps, drho
313 
314  CALL timeset(routinen, handle)
315 
316  rho_min = dielectric%params%rho_min
317  rho_max = dielectric%params%rho_max
318  eps0 = dielectric%params%eps0
319  derivative_method = dielectric%params%derivative_method
320  times_called = dielectric%params%times_called
321  dielec_functiontype = dielectric%params%dielec_functiontype
322 
323  IF (.NOT. dielec_functiontype .EQ. rho_dependent .AND. &
324  ((derivative_method .EQ. derivative_fft_use_deps) .OR. &
325  (derivative_method .EQ. derivative_fft_use_deps))) THEN
326  CALL cp_abort(__location__, &
327  "The specified derivative method is not compatible with the type of "// &
328  "the dielectric constant function.")
329  END IF
330 
331  CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
332 
333  ! make sure rho is in the real space
334  CALL pw_pool_orig%create_pw(rho_elec_rs)
335  CALL pw_transfer(rho, rho_elec_rs)
336  ! expand rho_elec
337  CALL pw_pool_xpndd%create_pw(rho_elec_rs_xpndd)
338  CALL pw_expand(neumann_directions, recv_msgs_bnds, dests_expand, srcs_expand, flipg_stat, bounds_shftd, &
339  rho_elec_rs, rho_elec_rs_xpndd)
340 
341  IF (PRESENT(rho_core)) THEN
342  ! make sure rho_core is in the real space
343  CALL pw_pool_orig%create_pw(rho_core_rs)
344  CALL pw_transfer(rho_core, rho_core_rs)
345  ! expand rho_core
346  CALL pw_pool_xpndd%create_pw(rho_core_rs_xpndd)
347  CALL pw_expand(neumann_directions, recv_msgs_bnds, dests_expand, srcs_expand, flipg_stat, bounds_shftd, &
348  rho_core_rs, rho_core_rs_xpndd)
349 
350  IF (dielectric%params%dielec_core_correction) THEN
351  ! use (rho_elec - rho_core) to compute dielectric
352  CALL pw_axpy(rho_core_rs_xpndd, rho_elec_rs_xpndd, -2.0_dp)
353  ELSE
354  CALL pw_axpy(rho_core_rs_xpndd, rho_elec_rs_xpndd, -1.0_dp)
355  END IF
356  CALL pw_pool_orig%give_back_pw(rho_core_rs)
357  CALL pw_pool_xpndd%give_back_pw(rho_core_rs_xpndd)
358  ELSE
359  cpabort("For dielectric constant larger than 1, rho_core has to be present.")
360  END IF
361 
362 ! calculate the dielectric constant
363  SELECT CASE (dielec_functiontype)
364  CASE (rho_dependent)
365  CALL dielectric_constant_sccs(rho_elec_rs_xpndd, dielectric%eps, dielectric%deps_drho, &
366  eps0, rho_max, rho_min)
367  CASE (spatially_dependent)
368  IF (times_called .EQ. 0) THEN
369  CALL dielectric_constant_spatially_dependent(dielectric%eps, pw_pool_xpndd, dielectric%params)
370  END IF
371  CASE (spatially_rho_dependent)
372  CALL dielectric_constant_spatially_rho_dependent(rho_elec_rs_xpndd, dielectric%eps, &
373  dielectric%deps_drho, pw_pool_xpndd, dielectric%params)
374  END SELECT
375 
376 ! derivatives
377  IF ((dielec_functiontype .EQ. rho_dependent) .OR. &
378  (dielec_functiontype .EQ. spatially_rho_dependent) .OR. &
379  ((dielec_functiontype .EQ. spatially_dependent) .AND. times_called .EQ. 0)) THEN
380  SELECT CASE (derivative_method)
381  CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
382  CALL pw_pool_xpndd%create_pw(ln_eps)
383  ln_eps%array = log(dielectric%eps%array)
384  CASE (derivative_fft_use_deps)
385  DO i = 1, 3
386  CALL pw_pool_xpndd%create_pw(deps(i))
387  CALL pw_zero(deps(i))
388  END DO
389  CASE (derivative_fft_use_drho)
390  DO i = 1, 3
391  CALL pw_pool_xpndd%create_pw(deps(i))
392  CALL pw_pool_xpndd%create_pw(drho(i))
393  CALL pw_zero(deps(i))
394  CALL pw_zero(drho(i))
395  END DO
396  END SELECT
397 
398  SELECT CASE (derivative_method)
399  CASE (derivative_cd3)
400  CALL derive_fdm_cd3(ln_eps, dielectric%dln_eps, diel_rs_grid)
401  CASE (derivative_cd5)
402  CALL derive_fdm_cd5(ln_eps, dielectric%dln_eps, diel_rs_grid)
403  CASE (derivative_cd7)
404  CALL derive_fdm_cd7(ln_eps, dielectric%dln_eps, diel_rs_grid)
405  CASE (derivative_fft)
406  CALL derive_fft(ln_eps, dielectric%dln_eps, pw_pool_xpndd)
407  CASE (derivative_fft_use_deps)
408  ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
409  CALL derive_fft(dielectric%eps, deps, pw_pool_xpndd)
410 
411  lb(1:3) = pw_pool_xpndd%pw_grid%bounds_local(1, 1:3)
412  ub(1:3) = pw_pool_xpndd%pw_grid%bounds_local(2, 1:3)
413  ! damp oscillations cuased by fft-based derivative in the region
414  ! where electron density is zero
415  DO idir = 1, 3
416  DO k = lb(3), ub(3)
417  DO j = lb(2), ub(2)
418  DO i = lb(1), ub(1)
419  IF (abs(dielectric%deps_drho%array(i, j, k)) .LE. small_value) THEN
420  deps(idir)%array(i, j, k) = 0.0_dp
421  END IF
422  END DO
423  END DO
424  END DO
425  dielectric%dln_eps(idir)%array = deps(idir)%array/dielectric%eps%array
426  END DO
427  CASE (derivative_fft_use_drho)
428  ! \Nabla \eps = \Nabla \rho \cdot \frac{\partial \eps}{\partial \rho}
429  ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
430  CALL derive_fft(rho_elec_rs_xpndd, drho, pw_pool_xpndd)
431  DO i = 1, 3
432  deps(i)%array = drho(i)%array*dielectric%deps_drho%array
433  dielectric%dln_eps(i)%array = deps(i)%array/dielectric%eps%array
434  END DO
435  END SELECT
436 
437  SELECT CASE (derivative_method)
438  CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
439  CALL pw_pool_xpndd%give_back_pw(ln_eps)
440  CASE (derivative_fft_use_deps)
441  DO i = 1, 3
442  CALL pw_pool_xpndd%give_back_pw(deps(i))
443  END DO
444  CASE (derivative_fft_use_drho)
445  DO i = 1, 3
446  CALL pw_pool_xpndd%give_back_pw(drho(i))
447  CALL pw_pool_xpndd%give_back_pw(deps(i))
448  END DO
449  END SELECT
450  END IF
451 
452  CALL pw_pool_orig%give_back_pw(rho_elec_rs)
453  CALL pw_pool_xpndd%give_back_pw(rho_elec_rs_xpndd)
454  CALL pw_pool_release(pw_pool_xpndd)
455 
456  dielectric%params%times_called = dielectric%params%times_called + 1
457 
458  CALL timestop(handle)
459 
460  END SUBROUTINE dielectric_compute_neumann_c1d_gs
461 ! **************************************************************************************************
462 !> \brief evaluates the dielectric constant
463 !> \param dielectric the dielectric data type to be initialized
464 !> \param diel_rs_grid real space grid for finite difference derivative
465 !> \param pw_pool pool of plane wave grid
466 !> \param rho electronic density
467 !> \param rho_core core density
468 !> \par History
469 !> 06.2014 created [Hossein Bani-Hashemian]
470 !> 12.2014 added finite difference derivatives [Hossein Bani-Hashemian]
471 !> 07.2015 density-independent dielectric regions [Hossein Bani-Hashemian]
472 !> \author Mohammad Hossein Bani-Hashemian
473 ! **************************************************************************************************
474  SUBROUTINE dielectric_compute_periodic_r3d_rs (dielectric, diel_rs_grid, pw_pool, rho, rho_core)
475 
476  TYPE(dielectric_type), INTENT(INOUT), POINTER :: dielectric
477  TYPE(realspace_grid_type), POINTER :: diel_rs_grid
478  TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
479  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
480  TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core
481 
482  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_compute_periodic'
483  REAL(dp), PARAMETER :: small_value = epsilon(1.0_dp)
484 
485  INTEGER :: derivative_method, dielec_functiontype, &
486  handle, i, idir, j, k, times_called
487  INTEGER, DIMENSION(3) :: lb, ub
488  REAL(dp) :: eps0, rho_max, rho_min
489  TYPE(pw_r3d_rs_type) :: ln_eps, rho_elec_rs
490  TYPE(pw_r3d_rs_type) :: rho_core_rs
491  TYPE(pw_r3d_rs_type), DIMENSION(3) :: deps, drho
492 
493  CALL timeset(routinen, handle)
494 
495  rho_min = dielectric%params%rho_min
496  rho_max = dielectric%params%rho_max
497  eps0 = dielectric%params%eps0
498  derivative_method = dielectric%params%derivative_method
499  times_called = dielectric%params%times_called
500  dielec_functiontype = dielectric%params%dielec_functiontype
501 
502  IF (.NOT. dielec_functiontype .EQ. rho_dependent .AND. &
503  ((derivative_method .EQ. derivative_fft_use_deps) .OR. &
504  (derivative_method .EQ. derivative_fft_use_deps))) THEN
505  CALL cp_abort(__location__, &
506  "The specified derivative method is not compatible with the type of "// &
507  "the dielectric constant function.")
508  END IF
509 
510  CALL pw_pool%create_pw(rho_elec_rs)
511 
512  ! for evaluating epsilon make sure rho is in the real space
513  CALL pw_transfer(rho, rho_elec_rs)
514 
515  IF (PRESENT(rho_core)) THEN
516  ! make sure rho_core is in the real space
517  CALL pw_pool%create_pw(rho_core_rs)
518  CALL pw_transfer(rho_core, rho_core_rs)
519  IF (dielectric%params%dielec_core_correction) THEN
520  ! use (rho_elec - rho_core) to compute dielectric to avoid obtaining spurious
521  ! epsilon in the core region
522  CALL pw_axpy(rho_core_rs, rho_elec_rs, -2.0_dp)
523  ELSE
524  CALL pw_axpy(rho_core_rs, rho_elec_rs, -1.0_dp)
525  END IF
526  CALL pw_pool%give_back_pw(rho_core_rs)
527  ELSE
528  cpabort("For dielectric constant larger than 1, rho_core has to be present.")
529  END IF
530 
531 ! calculate the dielectric constant
532  SELECT CASE (dielec_functiontype)
533  CASE (rho_dependent)
534  CALL dielectric_constant_sccs(rho_elec_rs, dielectric%eps, dielectric%deps_drho, &
535  eps0, rho_max, rho_min)
536  CASE (spatially_dependent)
537  IF (times_called .EQ. 0) THEN
538  CALL dielectric_constant_spatially_dependent(dielectric%eps, pw_pool, dielectric%params)
539  END IF
540  CASE (spatially_rho_dependent)
541  CALL dielectric_constant_spatially_rho_dependent(rho_elec_rs, dielectric%eps, &
542  dielectric%deps_drho, pw_pool, dielectric%params)
543  END SELECT
544 
545 ! derivatives
546  IF ((dielec_functiontype .EQ. rho_dependent) .OR. &
547  (dielec_functiontype .EQ. spatially_rho_dependent) .OR. &
548  ((dielec_functiontype .EQ. spatially_dependent) .AND. times_called .EQ. 0)) THEN
549  SELECT CASE (derivative_method)
550  CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
551  CALL pw_pool%create_pw(ln_eps)
552  ln_eps%array = log(dielectric%eps%array)
553  CASE (derivative_fft_use_deps)
554  DO i = 1, 3
555  CALL pw_pool%create_pw(deps(i))
556  CALL pw_zero(deps(i))
557  END DO
558  CASE (derivative_fft_use_drho)
559  DO i = 1, 3
560  CALL pw_pool%create_pw(deps(i))
561  CALL pw_pool%create_pw(drho(i))
562  CALL pw_zero(deps(i))
563  CALL pw_zero(drho(i))
564  END DO
565  END SELECT
566 
567  SELECT CASE (derivative_method)
568  CASE (derivative_cd3)
569  CALL derive_fdm_cd3(ln_eps, dielectric%dln_eps, diel_rs_grid)
570  CASE (derivative_cd5)
571  CALL derive_fdm_cd5(ln_eps, dielectric%dln_eps, diel_rs_grid)
572  CASE (derivative_cd7)
573  CALL derive_fdm_cd7(ln_eps, dielectric%dln_eps, diel_rs_grid)
574  CASE (derivative_fft)
575  CALL derive_fft(ln_eps, dielectric%dln_eps, pw_pool)
576  CASE (derivative_fft_use_deps)
577  ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
578  CALL derive_fft(dielectric%eps, deps, pw_pool)
579 
580  lb(1:3) = pw_pool%pw_grid%bounds_local(1, 1:3)
581  ub(1:3) = pw_pool%pw_grid%bounds_local(2, 1:3)
582  ! damp oscillations cuased by fft-based derivative in the region
583  ! where electron density is zero
584  DO idir = 1, 3
585  DO k = lb(3), ub(3)
586  DO j = lb(2), ub(2)
587  DO i = lb(1), ub(1)
588  IF (abs(dielectric%deps_drho%array(i, j, k)) .LE. small_value) THEN
589  deps(idir)%array(i, j, k) = 0.0_dp
590  END IF
591  END DO
592  END DO
593  END DO
594  dielectric%dln_eps(idir)%array = deps(idir)%array/dielectric%eps%array
595  END DO
596  CASE (derivative_fft_use_drho)
597  ! \Nabla \eps = \Nabla \rho \cdot \frac{\partial \eps}{\partial \rho}
598  ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
599  CALL derive_fft(rho_elec_rs, drho, pw_pool)
600  DO i = 1, 3
601  deps(i)%array = drho(i)%array*dielectric%deps_drho%array
602  dielectric%dln_eps(i)%array = deps(i)%array/dielectric%eps%array
603  END DO
604  END SELECT
605 
606  SELECT CASE (derivative_method)
607  CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
608  CALL pw_pool%give_back_pw(ln_eps)
609  CASE (derivative_fft_use_deps)
610  DO i = 1, 3
611  CALL pw_pool%give_back_pw(deps(i))
612  END DO
613  CASE (derivative_fft_use_drho)
614  DO i = 1, 3
615  CALL pw_pool%give_back_pw(drho(i))
616  CALL pw_pool%give_back_pw(deps(i))
617  END DO
618  END SELECT
619  END IF
620 
621  CALL pw_pool%give_back_pw(rho_elec_rs)
622 
623  dielectric%params%times_called = dielectric%params%times_called + 1
624 
625  CALL timestop(handle)
626 
627  END SUBROUTINE dielectric_compute_periodic_r3d_rs
628 
629 ! **************************************************************************************************
630 !> \brief evaluates the dielectric constant for non-periodic (Neumann-type)
631 !> boundaries
632 !> \param dielectric the dielectric data type to be initialized
633 !> \param diel_rs_grid real space grid for finite difference derivative
634 !> \param pw_pool_orig pool of plane wave grid
635 !> \param dct_pw_grid ...
636 !> \param neumann_directions ...
637 !> \param recv_msgs_bnds bounds of the messages to be received (pw_expand)
638 !> \param dests_expand list of the destination processes (pw_expand)
639 !> \param srcs_expand list of the source processes (pw_expand)
640 !> \param flipg_stat flipping status for the received data chunks (pw_expand)
641 !> \param bounds_shftd bounds of the original grid shifted to have g0 in the
642 !> middle of the cell (pw_expand)
643 !> \param rho electronic density
644 !> \param rho_core core density
645 !> \par History
646 !> 11.2015 created [Hossein Bani-Hashemian]
647 !> \author Mohammad Hossein Bani-Hashemian
648 ! **************************************************************************************************
649  SUBROUTINE dielectric_compute_neumann_r3d_rs (dielectric, diel_rs_grid, pw_pool_orig, &
650  dct_pw_grid, neumann_directions, recv_msgs_bnds, &
651  dests_expand, srcs_expand, flipg_stat, bounds_shftd, rho, rho_core)
652 
653  TYPE(dielectric_type), INTENT(INOUT), POINTER :: dielectric
654  TYPE(realspace_grid_type), POINTER :: diel_rs_grid
655  TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool_orig
656  TYPE(pw_grid_type), INTENT(IN), POINTER :: dct_pw_grid
657  INTEGER, INTENT(IN) :: neumann_directions
658  INTEGER, DIMENSION(:, :, :), INTENT(IN), POINTER :: recv_msgs_bnds
659  INTEGER, DIMENSION(:), INTENT(IN), POINTER :: dests_expand, srcs_expand, flipg_stat
660  INTEGER, DIMENSION(2, 3), INTENT(IN) :: bounds_shftd
661  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho
662  TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core
663 
664  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_compute_neumann'
665  REAL(dp), PARAMETER :: small_value = epsilon(1.0_dp)
666 
667  INTEGER :: derivative_method, dielec_functiontype, &
668  handle, i, idir, j, k, times_called
669  INTEGER, DIMENSION(3) :: lb, ub
670  REAL(dp) :: eps0, rho_max, rho_min
671  TYPE(pw_pool_type), POINTER :: pw_pool_xpndd
672  TYPE(pw_r3d_rs_type) :: ln_eps, rho_core_rs, rho_core_rs_xpndd, &
673  rho_elec_rs, rho_elec_rs_xpndd
674  TYPE(pw_r3d_rs_type), DIMENSION(3) :: deps, drho
675 
676  CALL timeset(routinen, handle)
677 
678  rho_min = dielectric%params%rho_min
679  rho_max = dielectric%params%rho_max
680  eps0 = dielectric%params%eps0
681  derivative_method = dielectric%params%derivative_method
682  times_called = dielectric%params%times_called
683  dielec_functiontype = dielectric%params%dielec_functiontype
684 
685  IF (.NOT. dielec_functiontype .EQ. rho_dependent .AND. &
686  ((derivative_method .EQ. derivative_fft_use_deps) .OR. &
687  (derivative_method .EQ. derivative_fft_use_deps))) THEN
688  CALL cp_abort(__location__, &
689  "The specified derivative method is not compatible with the type of "// &
690  "the dielectric constant function.")
691  END IF
692 
693  CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
694 
695  ! make sure rho is in the real space
696  CALL pw_pool_orig%create_pw(rho_elec_rs)
697  CALL pw_transfer(rho, rho_elec_rs)
698  ! expand rho_elec
699  CALL pw_pool_xpndd%create_pw(rho_elec_rs_xpndd)
700  CALL pw_expand(neumann_directions, recv_msgs_bnds, dests_expand, srcs_expand, flipg_stat, bounds_shftd, &
701  rho_elec_rs, rho_elec_rs_xpndd)
702 
703  IF (PRESENT(rho_core)) THEN
704  ! make sure rho_core is in the real space
705  CALL pw_pool_orig%create_pw(rho_core_rs)
706  CALL pw_transfer(rho_core, rho_core_rs)
707  ! expand rho_core
708  CALL pw_pool_xpndd%create_pw(rho_core_rs_xpndd)
709  CALL pw_expand(neumann_directions, recv_msgs_bnds, dests_expand, srcs_expand, flipg_stat, bounds_shftd, &
710  rho_core_rs, rho_core_rs_xpndd)
711 
712  IF (dielectric%params%dielec_core_correction) THEN
713  ! use (rho_elec - rho_core) to compute dielectric
714  CALL pw_axpy(rho_core_rs_xpndd, rho_elec_rs_xpndd, -2.0_dp)
715  ELSE
716  CALL pw_axpy(rho_core_rs_xpndd, rho_elec_rs_xpndd, -1.0_dp)
717  END IF
718  CALL pw_pool_orig%give_back_pw(rho_core_rs)
719  CALL pw_pool_xpndd%give_back_pw(rho_core_rs_xpndd)
720  ELSE
721  cpabort("For dielectric constant larger than 1, rho_core has to be present.")
722  END IF
723 
724 ! calculate the dielectric constant
725  SELECT CASE (dielec_functiontype)
726  CASE (rho_dependent)
727  CALL dielectric_constant_sccs(rho_elec_rs_xpndd, dielectric%eps, dielectric%deps_drho, &
728  eps0, rho_max, rho_min)
729  CASE (spatially_dependent)
730  IF (times_called .EQ. 0) THEN
731  CALL dielectric_constant_spatially_dependent(dielectric%eps, pw_pool_xpndd, dielectric%params)
732  END IF
733  CASE (spatially_rho_dependent)
734  CALL dielectric_constant_spatially_rho_dependent(rho_elec_rs_xpndd, dielectric%eps, &
735  dielectric%deps_drho, pw_pool_xpndd, dielectric%params)
736  END SELECT
737 
738 ! derivatives
739  IF ((dielec_functiontype .EQ. rho_dependent) .OR. &
740  (dielec_functiontype .EQ. spatially_rho_dependent) .OR. &
741  ((dielec_functiontype .EQ. spatially_dependent) .AND. times_called .EQ. 0)) THEN
742  SELECT CASE (derivative_method)
743  CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
744  CALL pw_pool_xpndd%create_pw(ln_eps)
745  ln_eps%array = log(dielectric%eps%array)
746  CASE (derivative_fft_use_deps)
747  DO i = 1, 3
748  CALL pw_pool_xpndd%create_pw(deps(i))
749  CALL pw_zero(deps(i))
750  END DO
751  CASE (derivative_fft_use_drho)
752  DO i = 1, 3
753  CALL pw_pool_xpndd%create_pw(deps(i))
754  CALL pw_pool_xpndd%create_pw(drho(i))
755  CALL pw_zero(deps(i))
756  CALL pw_zero(drho(i))
757  END DO
758  END SELECT
759 
760  SELECT CASE (derivative_method)
761  CASE (derivative_cd3)
762  CALL derive_fdm_cd3(ln_eps, dielectric%dln_eps, diel_rs_grid)
763  CASE (derivative_cd5)
764  CALL derive_fdm_cd5(ln_eps, dielectric%dln_eps, diel_rs_grid)
765  CASE (derivative_cd7)
766  CALL derive_fdm_cd7(ln_eps, dielectric%dln_eps, diel_rs_grid)
767  CASE (derivative_fft)
768  CALL derive_fft(ln_eps, dielectric%dln_eps, pw_pool_xpndd)
769  CASE (derivative_fft_use_deps)
770  ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
771  CALL derive_fft(dielectric%eps, deps, pw_pool_xpndd)
772 
773  lb(1:3) = pw_pool_xpndd%pw_grid%bounds_local(1, 1:3)
774  ub(1:3) = pw_pool_xpndd%pw_grid%bounds_local(2, 1:3)
775  ! damp oscillations cuased by fft-based derivative in the region
776  ! where electron density is zero
777  DO idir = 1, 3
778  DO k = lb(3), ub(3)
779  DO j = lb(2), ub(2)
780  DO i = lb(1), ub(1)
781  IF (abs(dielectric%deps_drho%array(i, j, k)) .LE. small_value) THEN
782  deps(idir)%array(i, j, k) = 0.0_dp
783  END IF
784  END DO
785  END DO
786  END DO
787  dielectric%dln_eps(idir)%array = deps(idir)%array/dielectric%eps%array
788  END DO
789  CASE (derivative_fft_use_drho)
790  ! \Nabla \eps = \Nabla \rho \cdot \frac{\partial \eps}{\partial \rho}
791  ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
792  CALL derive_fft(rho_elec_rs_xpndd, drho, pw_pool_xpndd)
793  DO i = 1, 3
794  deps(i)%array = drho(i)%array*dielectric%deps_drho%array
795  dielectric%dln_eps(i)%array = deps(i)%array/dielectric%eps%array
796  END DO
797  END SELECT
798 
799  SELECT CASE (derivative_method)
800  CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
801  CALL pw_pool_xpndd%give_back_pw(ln_eps)
802  CASE (derivative_fft_use_deps)
803  DO i = 1, 3
804  CALL pw_pool_xpndd%give_back_pw(deps(i))
805  END DO
806  CASE (derivative_fft_use_drho)
807  DO i = 1, 3
808  CALL pw_pool_xpndd%give_back_pw(drho(i))
809  CALL pw_pool_xpndd%give_back_pw(deps(i))
810  END DO
811  END SELECT
812  END IF
813 
814  CALL pw_pool_orig%give_back_pw(rho_elec_rs)
815  CALL pw_pool_xpndd%give_back_pw(rho_elec_rs_xpndd)
816  CALL pw_pool_release(pw_pool_xpndd)
817 
818  dielectric%params%times_called = dielectric%params%times_called + 1
819 
820  CALL timestop(handle)
821 
822  END SUBROUTINE dielectric_compute_neumann_r3d_rs
823 
824 ! **************************************************************************************************
825 !> \brief calculates the dielectric constant as a function of the electronic density
826 !> [see O. Andreussi, I. Dabo, and N. Marzari, J. Chem. Phys., 136, 064102 (2012)]
827 !> \param rho electron density
828 !> \param eps dielectric constant
829 !> \param deps_drho derivative of the dielectric constant wrt the density
830 !> \param eps0 dielectric constant in the bulk of the solvent
831 !> \param rho_max upper density threshold
832 !> \param rho_min lower density threshold
833 !> \par History
834 !> 06.2014 created [Hossein Bani-Hashemian]
835 !> \author Mohammad Hossein Bani-Hashemian
836 ! **************************************************************************************************
837  SUBROUTINE dielectric_constant_sccs(rho, eps, deps_drho, eps0, rho_max, rho_min)
838 
839  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho, eps, deps_drho
840  REAL(kind=dp), INTENT(IN) :: eps0, rho_max, rho_min
841 
842  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_constant_sccs'
843 
844  INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
845  ub2, ub3
846  INTEGER, DIMENSION(2, 3) :: bounds_local
847  REAL(kind=dp) :: denom, t
848 
849  CALL timeset(routinen, handle)
850 
851  IF (eps0 .LT. 1.0_dp) THEN
852  cpabort("The dielectric constant has to be greater than or equal to 1.")
853  END IF
854 
855  bounds_local = rho%pw_grid%bounds_local
856  lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
857  lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
858  lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
859 
860  denom = log(rho_max) - log(rho_min)
861  DO k = lb3, ub3
862  DO j = lb2, ub2
863  DO i = lb1, ub1
864  IF (rho%array(i, j, k) .LT. rho_min) THEN
865  eps%array(i, j, k) = eps0
866  deps_drho%array(i, j, k) = 0.0_dp
867  ELSE IF (rho%array(i, j, k) .GT. rho_max) THEN
868  eps%array(i, j, k) = 1.0_dp
869  deps_drho%array(i, j, k) = 0.0_dp
870  ELSE
871  t = twopi*(log(rho_max) - log(rho%array(i, j, k)))/denom
872  eps%array(i, j, k) = exp(log(eps0)*(t - sin(t))/twopi)
873  deps_drho%array(i, j, k) = -eps%array(i, j, k)*log(eps0)*(1.0_dp - cos(t))/(denom*rho%array(i, j, k))
874  END IF
875  END DO
876  END DO
877  END DO
878 
879  CALL timestop(handle)
880 
881  END SUBROUTINE dielectric_constant_sccs
882 
883 ! **************************************************************************************************
884 !> \brief creates an axis-aligned cuboidal dielectric region
885 !> \param eps dielectric constant function
886 !> \param dielec_const dielectric constant inside the region
887 !> \param pw_pool pool of planewave grid
888 !> \param zeta the mollifier's width
889 !> \param x_xtnt the x extent of the cuboidal region
890 !> \param y_xtnt the y extent of the cuboidal region
891 !> \param z_xtnt the z extent of the cuboidal region
892 !> \param x_glbl x grid vector of the simulation box
893 !> \param y_glbl y grid vector of the simulation box
894 !> \param z_glbl z grid vector of the simulation box
895 !> \param x_locl x grid vector of the simulation box local to this process
896 !> \param y_locl y grid vector of the simulation box local to this process
897 !> \param z_locl z grid vector of the simulation box local to this process
898 !> \par History
899 !> 07.2015 created [Hossein Bani-Hashemian]
900 !> \author Mohammad Hossein Bani-Hashemian
901 ! **************************************************************************************************
902  SUBROUTINE dielectric_constant_aa_cuboidal(eps, dielec_const, pw_pool, zeta, &
903  x_xtnt, y_xtnt, z_xtnt, &
904  x_glbl, y_glbl, z_glbl, &
905  x_locl, y_locl, z_locl)
906 
907  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: eps
908  REAL(kind=dp), INTENT(IN) :: dielec_const
909  TYPE(pw_pool_type), POINTER :: pw_pool
910  REAL(kind=dp), INTENT(IN) :: zeta
911  REAL(dp), DIMENSION(2), INTENT(IN) :: x_xtnt, y_xtnt, z_xtnt
912  REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_glbl, y_glbl, z_glbl, x_locl, y_locl, &
913  z_locl
914 
915  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_constant_aa_cuboidal'
916  LOGICAL, DIMENSION(6), PARAMETER :: test_forb_xtnts = .true.
917 
918  INTEGER :: handle, i, j, k, lb1, lb2, lb3, &
919  n_forb_xtnts, ub1, ub2, ub3
920  INTEGER, DIMENSION(2, 3) :: bounds_local
921  LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
922  forb_xtnt4, forb_xtnt5, forb_xtnt6
923  REAL(kind=dp) :: dx, dy, dz
924  TYPE(pw_grid_type), POINTER :: pw_grid
925  TYPE(pw_r3d_rs_type) :: eps_tmp
926 
927  CALL timeset(routinen, handle)
928 
929  IF (dielec_const .LT. 1.0_dp) THEN
930  cpabort("The dielectric constant has to be greater than or equal to 1.")
931  END IF
932 
933  pw_grid => eps%pw_grid
934 
935  dx = pw_grid%dr(1); dy = pw_grid%dr(2); dz = pw_grid%dr(3)
936 
937  forb_xtnt1 = x_xtnt(1) .LT. x_glbl(lbound(x_glbl, 1))
938  forb_xtnt2 = x_xtnt(2) .GT. x_glbl(ubound(x_glbl, 1)) + dx
939  forb_xtnt3 = y_xtnt(1) .LT. y_glbl(lbound(y_glbl, 1))
940  forb_xtnt4 = y_xtnt(2) .GT. y_glbl(ubound(y_glbl, 1)) + dy
941  forb_xtnt5 = z_xtnt(1) .LT. z_glbl(lbound(z_glbl, 1))
942  forb_xtnt6 = z_xtnt(2) .GT. z_glbl(ubound(z_glbl, 1)) + dz
943  n_forb_xtnts = count((/forb_xtnt1, forb_xtnt2, forb_xtnt3, &
944  forb_xtnt4, forb_xtnt5, forb_xtnt6/) .EQV. test_forb_xtnts)
945  IF (n_forb_xtnts .GT. 0) THEN
946  CALL cp_abort(__location__, &
947  "The given extents for the dielectric region are outside the range of "// &
948  "the simulation cell.")
949  END IF
950 
951  CALL pw_pool%create_pw(eps_tmp)
952  CALL pw_copy(eps, eps_tmp)
953 
954  bounds_local = pw_grid%bounds_local
955  lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
956  lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
957  lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
958 
959  DO k = lb3, ub3
960  DO j = lb2, ub2
961  DO i = lb1, ub1
962  IF ((x_locl(i) .GE. x_xtnt(1)) .AND. (x_locl(i) .LE. x_xtnt(2)) .AND. &
963  (y_locl(j) .GE. y_xtnt(1)) .AND. (y_locl(j) .LE. y_xtnt(2)) .AND. &
964  (z_locl(k) .GE. z_xtnt(1)) .AND. (z_locl(k) .LE. z_xtnt(2))) THEN
965  eps_tmp%array(i, j, k) = dielec_const
966  END IF
967  END DO
968  END DO
969  END DO
970 
971  CALL pw_mollifier(pw_pool, zeta, x_glbl, y_glbl, z_glbl, eps_tmp, eps)
972  CALL pw_pool%give_back_pw(eps_tmp)
973 
974  CALL timestop(handle)
975 
976  END SUBROUTINE dielectric_constant_aa_cuboidal
977 
978 ! **************************************************************************************************
979 !> \brief creates an x-axis aligned annular dielectric region
980 !> \param eps dielectric constant function
981 !> \param dielec_const dielectric constant inside the region
982 !> \param pw_pool pool of planewave grid
983 !> \param zeta the mollifier's width
984 !> \param x_xtnt the x extent of the annular region
985 !> \param base_center the y and z coordinates of the cylinder's base
986 !> \param base_radii the radii of the annulus' base
987 !> \param x_glbl x grid vector of the simulation box
988 !> \param y_glbl y grid vector of the simulation box
989 !> \param z_glbl z grid vector of the simulation box
990 !> \param x_locl x grid vector of the simulation box local to this process
991 !> \param y_locl y grid vector of the simulation box local to this process
992 !> \param z_locl z grid vector of the simulation box local to this process
993 !> \par History
994 !> 07.2015 created [Hossein Bani-Hashemian]
995 !> \author Mohammad Hossein Bani-Hashemian
996 ! **************************************************************************************************
997  SUBROUTINE dielectric_constant_xaa_annular(eps, dielec_const, pw_pool, zeta, &
998  x_xtnt, base_center, base_radii, &
999  x_glbl, y_glbl, z_glbl, &
1000  x_locl, y_locl, z_locl)
1001 
1002  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: eps
1003  REAL(kind=dp), INTENT(IN) :: dielec_const
1004  TYPE(pw_pool_type), POINTER :: pw_pool
1005  REAL(dp), INTENT(IN) :: zeta
1006  REAL(dp), DIMENSION(2), INTENT(IN) :: x_xtnt, base_center, base_radii
1007  REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_glbl, y_glbl, z_glbl, x_locl, y_locl, &
1008  z_locl
1009 
1010  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_constant_xaa_annular'
1011  LOGICAL, DIMENSION(6), PARAMETER :: test_forb_xtnts = .true.
1012 
1013  INTEGER :: handle, i, j, k, lb1, lb2, lb3, &
1014  n_forb_xtnts, ub1, ub2, ub3
1015  INTEGER, DIMENSION(2, 3) :: bounds_local
1016  LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
1017  forb_xtnt4, forb_xtnt5, forb_xtnt6
1018  REAL(kind=dp) :: bctry, bctrz, distsq, dx, dy, dz
1019  TYPE(pw_grid_type), POINTER :: pw_grid
1020  TYPE(pw_r3d_rs_type) :: eps_tmp
1021 
1022  CALL timeset(routinen, handle)
1023 
1024  IF (dielec_const .LT. 1.0_dp) THEN
1025  cpabort("The dielectric constant has to be greater than or equal to 1.")
1026  END IF
1027 
1028  pw_grid => eps%pw_grid
1029 
1030  bctry = base_center(1); bctrz = base_center(2)
1031  dx = pw_grid%dr(1); dy = pw_grid%dr(2); dz = pw_grid%dr(3)
1032 
1033  forb_xtnt1 = x_xtnt(1) .LT. x_glbl(lbound(x_glbl, 1))
1034  forb_xtnt2 = x_xtnt(2) .GT. x_glbl(ubound(x_glbl, 1)) + dx
1035  forb_xtnt3 = bctry - maxval(base_radii) .LT. y_glbl(lbound(y_glbl, 1))
1036  forb_xtnt4 = bctry + maxval(base_radii) .GT. y_glbl(ubound(y_glbl, 1)) + dy
1037  forb_xtnt5 = bctrz - maxval(base_radii) .LT. z_glbl(lbound(z_glbl, 1))
1038  forb_xtnt6 = bctrz + maxval(base_radii) .GT. z_glbl(ubound(z_glbl, 1)) + dz
1039  n_forb_xtnts = count((/forb_xtnt1, forb_xtnt2, forb_xtnt3, &
1040  forb_xtnt4, forb_xtnt5, forb_xtnt6/) .EQV. test_forb_xtnts)
1041  IF (n_forb_xtnts .GT. 0) THEN
1042  CALL cp_abort(__location__, &
1043  "The given extents for the dielectric region are outside the range of "// &
1044  "the simulation cell.")
1045  END IF
1046 
1047  CALL pw_pool%create_pw(eps_tmp)
1048  CALL pw_copy(eps, eps_tmp)
1049 
1050  bounds_local = pw_grid%bounds_local
1051  lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
1052  lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
1053  lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
1054 
1055  DO k = lb3, ub3
1056  DO j = lb2, ub2
1057  DO i = lb1, ub1
1058  distsq = (y_locl(j) - bctry)**2 + (z_locl(k) - bctrz)**2
1059  IF ((x_locl(i) .GE. x_xtnt(1)) .AND. (x_locl(i) .LE. x_xtnt(2)) .AND. &
1060  (distsq .GE. minval(base_radii)**2) .AND. (distsq .LE. maxval(base_radii)**2)) THEN
1061  eps_tmp%array(i, j, k) = dielec_const
1062  END IF
1063  END DO
1064  END DO
1065  END DO
1066 
1067  CALL pw_mollifier(pw_pool, zeta, x_glbl, y_glbl, z_glbl, eps_tmp, eps)
1068  CALL pw_pool%give_back_pw(eps_tmp)
1069 
1070  CALL timestop(handle)
1071 
1072  END SUBROUTINE dielectric_constant_xaa_annular
1073 
1074 ! **************************************************************************************************
1075 !> \brief Sets up density independent dielectric regions
1076 !> \param eps dielectric constant function
1077 !> \param pw_pool pool of planewave grid
1078 !> \param dielectric_params dielectric parameters read from input file
1079 !> \par History
1080 !> 07.2015 created [Hossein Bani-Hashemian]
1081 !> \author Mohammad Hossein Bani-Hashemian
1082 ! **************************************************************************************************
1083  SUBROUTINE dielectric_constant_spatially_dependent(eps, pw_pool, dielectric_params)
1084 
1085  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: eps
1086  TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1087  TYPE(dielectric_parameters), INTENT(IN) :: dielectric_params
1088 
1089  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_constant_spatially_dependent'
1090 
1091  INTEGER :: handle, j, n_aa_cuboidal, &
1092  n_dielectric_region, n_xaa_annular
1093  REAL(dp) :: dielec_const, zeta
1094  REAL(dp), ALLOCATABLE, DIMENSION(:) :: x_glbl, x_locl, y_glbl, y_locl, z_glbl, &
1095  z_locl
1096  REAL(dp), DIMENSION(2) :: base_center, base_radii
1097  TYPE(pw_grid_type), POINTER :: pw_grid
1098 
1099  CALL timeset(routinen, handle)
1100 
1101  eps%array = dielectric_params%eps0
1102 
1103  n_aa_cuboidal = dielectric_params%n_aa_cuboidal
1104  n_xaa_annular = dielectric_params%n_xaa_annular
1105  pw_grid => pw_pool%pw_grid
1106  CALL setup_grid_axes(pw_grid, x_glbl, y_glbl, z_glbl, x_locl, y_locl, z_locl)
1107  n_dielectric_region = n_aa_cuboidal + n_xaa_annular
1108 
1109  IF (n_dielectric_region .EQ. 0) THEN
1110  cpabort("No density independent dielectric region is defined.")
1111  END IF
1112 
1113  DO j = 1, n_aa_cuboidal
1114  dielec_const = dielectric_params%aa_cuboidal_eps(j)
1115  zeta = dielectric_params%aa_cuboidal_zeta(j)
1116 
1117  CALL dielectric_constant_aa_cuboidal(eps, dielec_const, pw_pool, zeta, &
1118  dielectric_params%aa_cuboidal_xxtnt(:, j), &
1119  dielectric_params%aa_cuboidal_yxtnt(:, j), &
1120  dielectric_params%aa_cuboidal_zxtnt(:, j), &
1121  x_glbl, y_glbl, z_glbl, &
1122  x_locl, y_locl, z_locl)
1123  END DO
1124 
1125  DO j = 1, n_xaa_annular
1126  base_center = dielectric_params%xaa_annular_bctr(:, j)
1127  base_radii = dielectric_params%xaa_annular_brad(:, j)
1128  dielec_const = dielectric_params%xaa_annular_eps(j)
1129  zeta = dielectric_params%xaa_annular_zeta(j)
1130 
1131  CALL dielectric_constant_xaa_annular(eps, dielec_const, pw_pool, zeta, &
1132  dielectric_params%xaa_annular_xxtnt(:, j), &
1133  base_center, base_radii, &
1134  x_glbl, y_glbl, z_glbl, &
1135  x_locl, y_locl, z_locl)
1136  END DO
1137 
1138  CALL timestop(handle)
1139 
1140  END SUBROUTINE dielectric_constant_spatially_dependent
1141 
1142 ! **************************************************************************************************
1143 !> \brief Sets up various dielectric constant regions. Using the sccs switching
1144 !> function the dielectric constant smoothly varies between 1, where the density
1145 !> is present and the values inside the dielectric regions, otherwise.
1146 !> \param rho electron density
1147 !> \param eps dielectric constant function
1148 !> \param deps_drho derivative of the dielectric constant wrt the density
1149 !> \param pw_pool pool of planewave grid
1150 !> \param dielectric_params dielectric parameters read from input file
1151 !> \par History
1152 !> 07.2015 created [Hossein Bani-Hashemian]
1153 !> \author Mohammad Hossein Bani-Hashemian
1154 ! **************************************************************************************************
1155  SUBROUTINE dielectric_constant_spatially_rho_dependent(rho, eps, deps_drho, &
1156  pw_pool, dielectric_params)
1157 
1158  TYPE(pw_r3d_rs_type), INTENT(IN) :: rho, eps, deps_drho
1159  TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
1160  TYPE(dielectric_parameters), INTENT(IN) :: dielectric_params
1161 
1162  CHARACTER(LEN=*), PARAMETER :: routinen = 'dielectric_constant_spatially_rho_dependent'
1163 
1164  INTEGER :: handle
1165  TYPE(pw_r3d_rs_type) :: dswch_func_drho, eps_sptldep, swch_func
1166 
1167  CALL timeset(routinen, handle)
1168 
1169  IF (dielectric_params%eps0 .LT. 1.0_dp) THEN
1170  cpabort("The dielectric constant has to be greater than or equal to 1.")
1171  END IF
1172 
1173  CALL pw_pool%create_pw(eps_sptldep)
1174  CALL pw_pool%create_pw(swch_func)
1175  CALL pw_pool%create_pw(dswch_func_drho)
1176  CALL pw_zero(eps_sptldep)
1177  CALL pw_zero(swch_func)
1178  CALL pw_zero(dswch_func_drho)
1179 
1180  CALL dielectric_constant_spatially_dependent(eps_sptldep, pw_pool, dielectric_params)
1181  CALL dielectric_constant_sccs(rho, swch_func, dswch_func_drho, 2.0_dp, &
1182  dielectric_params%rho_max, dielectric_params%rho_min)
1183 
1184  eps%array = ((swch_func%array - 1.0_dp)*(eps_sptldep%array - 1.0_dp)) + 1.0_dp
1185  deps_drho%array = dswch_func_drho%array*(eps_sptldep%array - 1.0_dp)
1186 
1187  CALL pw_pool%give_back_pw(dswch_func_drho)
1188  CALL pw_pool%give_back_pw(swch_func)
1189  CALL pw_pool%give_back_pw(eps_sptldep)
1190 
1191  CALL timestop(handle)
1192 
1193  END SUBROUTINE dielectric_constant_spatially_rho_dependent
1194 
1195 ! **************************************************************************************************
1196 !> \brief computes the derivative of a function using FFT
1197 !> \param f input funcition
1198 !> \param df derivative of f
1199 !> \param pw_pool pool of plane-wave grid
1200 ! **************************************************************************************************
1201  SUBROUTINE derive_fft(f, df, pw_pool)
1202 
1203  TYPE(pw_r3d_rs_type), INTENT(IN) :: f
1204  TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: df
1205  TYPE(pw_pool_type), POINTER :: pw_pool
1206 
1207  CHARACTER(LEN=*), PARAMETER :: routinen = 'derive_fft_r3d'
1208 
1209  INTEGER :: handle, i
1210  INTEGER, DIMENSION(3) :: nd
1211  TYPE(pw_c1d_gs_type), DIMENSION(2) :: work_gs
1212 
1213  CALL timeset(routinen, handle)
1214 
1215  DO i = 1, 2
1216  CALL pw_pool%create_pw(work_gs(i))
1217  END DO
1218 
1219  CALL pw_transfer(f, work_gs(1))
1220  DO i = 1, 3
1221  nd(:) = 0
1222  nd(i) = 1
1223  CALL pw_copy(work_gs(1), work_gs(2))
1224  CALL pw_derive(work_gs(2), nd(:))
1225  CALL pw_transfer(work_gs(2), df(i))
1226  END DO
1227 
1228  DO i = 1, 2
1229  CALL pw_pool%give_back_pw(work_gs(i))
1230  END DO
1231 
1232  CALL timestop(handle)
1233 
1234  END SUBROUTINE derive_fft
1235 
1236 END MODULE dielectric_methods
dct
the type I Discrete Cosine Transform (DCT-I)
Definition: dct.F:16
dct::pw_expand
subroutine, public pw_expand(neumann_directions, recv_msgs_bnds, dests_expand, srcs_expand, flipg_stat, bounds_shftd, pw_in, pw_expanded)
expands a pw_r3d_rs_type data to an evenly symmetric pw_r3d_rs_type data that is 8 times larger than ...
Definition: dct.F:518
dielectric_methods
methods for evaluating the dielectric constant
Definition: dielectric_methods.F:14
dielectric_methods::dielectric_create
subroutine, public dielectric_create(dielectric, pw_pool, dielectric_params)
allocates memory for a dielectric data type
Definition: dielectric_methods.F:68
dielectric_methods::derive_fft
subroutine, public derive_fft(f, df, pw_pool)
computes the derivative of a function using FFT
Definition: dielectric_methods.F:1202
dielectric_types
dielectric constant data type
Definition: dielectric_types.F:14
dielectric_types::derivative_fft_use_drho
integer, parameter, public derivative_fft_use_drho
Definition: dielectric_types.F:27
dielectric_types::derivative_fft_use_deps
integer, parameter, public derivative_fft_use_deps
Definition: dielectric_types.F:27
dielectric_types::derivative_fft
integer, parameter, public derivative_fft
Definition: dielectric_types.F:27
dielectric_types::derivative_cd5
integer, parameter, public derivative_cd5
Definition: dielectric_types.F:27
dielectric_types::spatially_rho_dependent
integer, parameter, public spatially_rho_dependent
Definition: dielectric_types.F:34
dielectric_types::derivative_cd3
integer, parameter, public derivative_cd3
Definition: dielectric_types.F:27
dielectric_types::spatially_dependent
integer, parameter, public spatially_dependent
Definition: dielectric_types.F:34
dielectric_types::rho_dependent
integer, parameter, public rho_dependent
Definition: dielectric_types.F:34
dielectric_types::derivative_cd7
integer, parameter, public derivative_cd7
Definition: dielectric_types.F:27
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
pw_grid_types
Definition: pw_grid_types.F:13
pw_methods
Definition: pw_methods.F:25
pw_methods::pw_derive
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition: pw_methods.F:10106
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_pool_types::pw_pool_release
subroutine, public pw_pool_release(pool)
releases the given pool (see cp2k/doc/ReferenceCounting.html)
Definition: pw_pool_types.F:259
pw_pool_types::pw_pool_create
subroutine, public pw_pool_create(pool, pw_grid, max_cache)
creates a pool for pw
Definition: pw_pool_types.F:173
pw_types
Definition: pw_types.F:24
realspace_grid_types
Definition: realspace_grid_types.F:18
rs_methods
numerical operations on real-space grid
Definition: rs_methods.F:14
rs_methods::setup_grid_axes
subroutine, public setup_grid_axes(pw_grid, x_glbl, y_glbl, z_glbl, x_locl, y_locl, z_locl)
returns the global axes and the portion of the axes that are local to the current mpi rank
Definition: rs_methods.F:274
rs_methods::derive_fdm_cd7
subroutine, public derive_fdm_cd7(f, df, rs_grid)
6th order finite difference derivative of a function on realspace grid
Definition: rs_methods.F:198
rs_methods::pw_mollifier
subroutine, public pw_mollifier(pw_pool, zeta, x_glbl, y_glbl, z_glbl, pw_in, pw_out)
convolutes a function with a smoothing kernel K_\zeta v * K_\zeta K_\zeta is the standard mollifier d...
Definition: rs_methods.F:366
rs_methods::derive_fdm_cd5
subroutine, public derive_fdm_cd5(f, df, rs_grid)
4th order finite difference derivative of a function on realspace grid
Definition: rs_methods.F:129
rs_methods::derive_fdm_cd3
subroutine, public derive_fdm_cd3(f, df, rs_grid)
2nd order finite difference derivative of a function on realspace grid
Definition: rs_methods.F:60