(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
14MODULE 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
56CONTAINS
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
1236END MODULE dielectric_methods
dielectric_methods::dielectric_compute
Definition dielectric_methods.F:49
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_copy
Definition pw_methods.F:145
pw_methods::pw_set
Definition pw_methods.F:126
pw_methods::pw_transfer
Definition pw_methods.F:303
pw_methods::pw_zero
Definition pw_methods.F:82
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
dielectric_types::dielectric_parameters
Definition dielectric_types.F:38
dielectric_types::dielectric_type
Definition dielectric_types.F:62
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:62
pw_pool_types::pw_pool_type
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:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
realspace_grid_types::realspace_grid_type
Definition realspace_grid_types.F:137