Line data Source code
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 54 : 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 54 : CALL timeset(routineN, handle)
77 :
78 54 : IF (.NOT. ASSOCIATED(dielectric)) THEN
79 270 : ALLOCATE (dielectric)
80 : NULLIFY (dielectric%eps)
81 : NULLIFY (dielectric%deps_drho)
82 54 : ALLOCATE (dielectric%eps, dielectric%deps_drho)
83 54 : CALL pw_pool%create_pw(dielectric%eps)
84 54 : CALL pw_pool%create_pw(dielectric%deps_drho)
85 54 : CALL pw_set(dielectric%eps, 1.0_dp)
86 54 : CALL pw_zero(dielectric%deps_drho)
87 216 : DO i = 1, 3
88 162 : CALL pw_pool%create_pw(dielectric%dln_eps(i))
89 216 : CALL pw_zero(dielectric%dln_eps(i))
90 : END DO
91 54 : dielectric%params = dielectric_params
92 54 : dielectric%params%times_called = 0
93 : END IF
94 :
95 54 : CALL timestop(handle)
96 :
97 54 : END SUBROUTINE dielectric_create
98 :
99 : #:for kind in ["c1d_gs", "r3d_rs"]
100 : ! **************************************************************************************************
101 : !> \brief evaluates the dielectric constant
102 : !> \param dielectric the dielectric data type to be initialized
103 : !> \param diel_rs_grid real space grid for finite difference derivative
104 : !> \param pw_pool pool of plane wave grid
105 : !> \param rho electronic density
106 : !> \param rho_core core density
107 : !> \par History
108 : !> 06.2014 created [Hossein Bani-Hashemian]
109 : !> 12.2014 added finite difference derivatives [Hossein Bani-Hashemian]
110 : !> 07.2015 density-independent dielectric regions [Hossein Bani-Hashemian]
111 : !> \author Mohammad Hossein Bani-Hashemian
112 : ! **************************************************************************************************
113 288 : SUBROUTINE dielectric_compute_periodic_${kind}$ (dielectric, diel_rs_grid, pw_pool, rho, rho_core)
114 :
115 : TYPE(dielectric_type), INTENT(INOUT), POINTER :: dielectric
116 : TYPE(realspace_grid_type), POINTER :: diel_rs_grid
117 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
118 : TYPE(pw_${kind}$_type), INTENT(IN) :: rho
119 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core
120 :
121 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dielectric_compute_periodic'
122 : REAL(dp), PARAMETER :: small_value = EPSILON(1.0_dp)
123 :
124 : INTEGER :: derivative_method, dielec_functiontype, &
125 : handle, i, idir, j, k, times_called
126 : INTEGER, DIMENSION(3) :: lb, ub
127 : REAL(dp) :: eps0, rho_max, rho_min
128 : TYPE(pw_r3d_rs_type) :: ln_eps, rho_elec_rs
129 : TYPE(pw_r3d_rs_type) :: rho_core_rs
130 2016 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: deps, drho
131 :
132 288 : CALL timeset(routineN, handle)
133 :
134 288 : rho_min = dielectric%params%rho_min
135 288 : rho_max = dielectric%params%rho_max
136 288 : eps0 = dielectric%params%eps0
137 288 : derivative_method = dielectric%params%derivative_method
138 288 : times_called = dielectric%params%times_called
139 288 : dielec_functiontype = dielectric%params%dielec_functiontype
140 :
141 288 : IF (.NOT. dielec_functiontype .EQ. rho_dependent .AND. &
142 : ((derivative_method .EQ. derivative_fft_use_deps) .OR. &
143 : (derivative_method .EQ. derivative_fft_use_deps))) THEN
144 : CALL cp_abort(__LOCATION__, &
145 : "The specified derivative method is not compatible with the type of "// &
146 0 : "the dielectric constant function.")
147 : END IF
148 :
149 288 : CALL pw_pool%create_pw(rho_elec_rs)
150 :
151 : ! for evaluating epsilon make sure rho is in the real space
152 288 : CALL pw_transfer(rho, rho_elec_rs)
153 :
154 288 : IF (PRESENT(rho_core)) THEN
155 : ! make sure rho_core is in the real space
156 288 : CALL pw_pool%create_pw(rho_core_rs)
157 288 : CALL pw_transfer(rho_core, rho_core_rs)
158 288 : IF (dielectric%params%dielec_core_correction) THEN
159 : ! use (rho_elec - rho_core) to compute dielectric to avoid obtaining spurious
160 : ! epsilon in the core region
161 288 : CALL pw_axpy(rho_core_rs, rho_elec_rs, -2.0_dp)
162 : ELSE
163 0 : CALL pw_axpy(rho_core_rs, rho_elec_rs, -1.0_dp)
164 : END IF
165 288 : CALL pw_pool%give_back_pw(rho_core_rs)
166 : ELSE
167 0 : CPABORT("For dielectric constant larger than 1, rho_core has to be present.")
168 : END IF
169 :
170 : ! calculate the dielectric constant
171 242 : SELECT CASE (dielec_functiontype)
172 : CASE (rho_dependent)
173 : CALL dielectric_constant_sccs(rho_elec_rs, dielectric%eps, dielectric%deps_drho, &
174 242 : eps0, rho_max, rho_min)
175 : CASE (spatially_dependent)
176 22 : IF (times_called .EQ. 0) THEN
177 2 : CALL dielectric_constant_spatially_dependent(dielectric%eps, pw_pool, dielectric%params)
178 : END IF
179 : CASE (spatially_rho_dependent)
180 : CALL dielectric_constant_spatially_rho_dependent(rho_elec_rs, dielectric%eps, &
181 288 : dielectric%deps_drho, pw_pool, dielectric%params)
182 : END SELECT
183 :
184 : ! derivatives
185 : IF ((dielec_functiontype .EQ. rho_dependent) .OR. &
186 288 : (dielec_functiontype .EQ. spatially_rho_dependent) .OR. &
187 : ((dielec_functiontype .EQ. spatially_dependent) .AND. times_called .EQ. 0)) THEN
188 258 : SELECT CASE (derivative_method)
189 : CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
190 258 : CALL pw_pool%create_pw(ln_eps)
191 41625606 : ln_eps%array = LOG(dielectric%eps%array)
192 : CASE (derivative_fft_use_deps)
193 40 : DO i = 1, 3
194 30 : CALL pw_pool%create_pw(deps(i))
195 40 : CALL pw_zero(deps(i))
196 : END DO
197 : CASE (derivative_fft_use_drho)
198 268 : DO i = 1, 3
199 0 : CALL pw_pool%create_pw(deps(i))
200 0 : CALL pw_pool%create_pw(drho(i))
201 0 : CALL pw_zero(deps(i))
202 0 : CALL pw_zero(drho(i))
203 : END DO
204 : END SELECT
205 :
206 72 : SELECT CASE (derivative_method)
207 : CASE (derivative_cd3)
208 72 : CALL derive_fdm_cd3(ln_eps, dielectric%dln_eps, diel_rs_grid)
209 : CASE (derivative_cd5)
210 114 : CALL derive_fdm_cd5(ln_eps, dielectric%dln_eps, diel_rs_grid)
211 : CASE (derivative_cd7)
212 72 : CALL derive_fdm_cd7(ln_eps, dielectric%dln_eps, diel_rs_grid)
213 : CASE (derivative_fft)
214 0 : CALL derive_fft(ln_eps, dielectric%dln_eps, pw_pool)
215 : CASE (derivative_fft_use_deps)
216 : ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
217 10 : CALL derive_fft(dielectric%eps, deps, pw_pool)
218 :
219 40 : lb(1:3) = pw_pool%pw_grid%bounds_local(1, 1:3)
220 40 : ub(1:3) = pw_pool%pw_grid%bounds_local(2, 1:3)
221 : ! damp oscillations cuased by fft-based derivative in the region
222 : ! where electron density is zero
223 40 : DO idir = 1, 3
224 2190 : DO k = lb(3), ub(3)
225 157710 : DO j = lb(2), ub(2)
226 5756400 : DO i = lb(1), ub(1)
227 5754240 : IF (ABS(dielectric%deps_drho%array(i, j, k)) .LE. small_value) THEN
228 5193252 : deps(idir)%array(i, j, k) = 0.0_dp
229 : END IF
230 : END DO
231 : END DO
232 : END DO
233 5756440 : dielectric%dln_eps(idir)%array = deps(idir)%array/dielectric%eps%array
234 : END DO
235 : CASE (derivative_fft_use_drho)
236 : ! \Nabla \eps = \Nabla \rho \cdot \frac{\partial \eps}{\partial \rho}
237 : ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
238 0 : CALL derive_fft(rho_elec_rs, drho, pw_pool)
239 268 : DO i = 1, 3
240 0 : deps(i)%array = drho(i)%array*dielectric%deps_drho%array
241 0 : dielectric%dln_eps(i)%array = deps(i)%array/dielectric%eps%array
242 : END DO
243 : END SELECT
244 :
245 258 : SELECT CASE (derivative_method)
246 : CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
247 258 : CALL pw_pool%give_back_pw(ln_eps)
248 : CASE (derivative_fft_use_deps)
249 40 : DO i = 1, 3
250 40 : CALL pw_pool%give_back_pw(deps(i))
251 : END DO
252 : CASE (derivative_fft_use_drho)
253 268 : DO i = 1, 3
254 0 : CALL pw_pool%give_back_pw(drho(i))
255 0 : CALL pw_pool%give_back_pw(deps(i))
256 : END DO
257 : END SELECT
258 : END IF
259 :
260 288 : CALL pw_pool%give_back_pw(rho_elec_rs)
261 :
262 288 : dielectric%params%times_called = dielectric%params%times_called + 1
263 :
264 288 : CALL timestop(handle)
265 :
266 288 : END SUBROUTINE dielectric_compute_periodic_${kind}$
267 :
268 : ! **************************************************************************************************
269 : !> \brief evaluates the dielectric constant for non-periodic (Neumann-type)
270 : !> boundaries
271 : !> \param dielectric the dielectric data type to be initialized
272 : !> \param diel_rs_grid real space grid for finite difference derivative
273 : !> \param pw_pool_orig pool of plane wave grid
274 : !> \param dct_pw_grid ...
275 : !> \param neumann_directions ...
276 : !> \param recv_msgs_bnds bounds of the messages to be received (pw_expand)
277 : !> \param dests_expand list of the destination processes (pw_expand)
278 : !> \param srcs_expand list of the source processes (pw_expand)
279 : !> \param flipg_stat flipping status for the received data chunks (pw_expand)
280 : !> \param bounds_shftd bounds of the original grid shifted to have g0 in the
281 : !> middle of the cell (pw_expand)
282 : !> \param rho electronic density
283 : !> \param rho_core core density
284 : !> \par History
285 : !> 11.2015 created [Hossein Bani-Hashemian]
286 : !> \author Mohammad Hossein Bani-Hashemian
287 : ! **************************************************************************************************
288 160 : SUBROUTINE dielectric_compute_neumann_${kind}$ (dielectric, diel_rs_grid, pw_pool_orig, &
289 : dct_pw_grid, neumann_directions, recv_msgs_bnds, &
290 : dests_expand, srcs_expand, flipg_stat, bounds_shftd, rho, rho_core)
291 :
292 : TYPE(dielectric_type), INTENT(INOUT), POINTER :: dielectric
293 : TYPE(realspace_grid_type), POINTER :: diel_rs_grid
294 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool_orig
295 : TYPE(pw_grid_type), INTENT(IN), POINTER :: dct_pw_grid
296 : INTEGER, INTENT(IN) :: neumann_directions
297 : INTEGER, DIMENSION(:, :, :), INTENT(IN), POINTER :: recv_msgs_bnds
298 : INTEGER, DIMENSION(:), INTENT(IN), POINTER :: dests_expand, srcs_expand, flipg_stat
299 : INTEGER, DIMENSION(2, 3), INTENT(IN) :: bounds_shftd
300 : TYPE(pw_${kind}$_type), INTENT(IN) :: rho
301 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: rho_core
302 :
303 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dielectric_compute_neumann'
304 : REAL(dp), PARAMETER :: small_value = EPSILON(1.0_dp)
305 :
306 : INTEGER :: derivative_method, dielec_functiontype, &
307 : handle, i, idir, j, k, times_called
308 : INTEGER, DIMENSION(3) :: lb, ub
309 : REAL(dp) :: eps0, rho_max, rho_min
310 : TYPE(pw_pool_type), POINTER :: pw_pool_xpndd
311 : TYPE(pw_r3d_rs_type) :: ln_eps, rho_core_rs, rho_core_rs_xpndd, &
312 : rho_elec_rs, rho_elec_rs_xpndd
313 1120 : TYPE(pw_r3d_rs_type), DIMENSION(3) :: deps, drho
314 :
315 160 : CALL timeset(routineN, handle)
316 :
317 160 : rho_min = dielectric%params%rho_min
318 160 : rho_max = dielectric%params%rho_max
319 160 : eps0 = dielectric%params%eps0
320 160 : derivative_method = dielectric%params%derivative_method
321 160 : times_called = dielectric%params%times_called
322 160 : dielec_functiontype = dielectric%params%dielec_functiontype
323 :
324 160 : IF (.NOT. dielec_functiontype .EQ. rho_dependent .AND. &
325 : ((derivative_method .EQ. derivative_fft_use_deps) .OR. &
326 : (derivative_method .EQ. derivative_fft_use_deps))) THEN
327 : CALL cp_abort(__LOCATION__, &
328 : "The specified derivative method is not compatible with the type of "// &
329 0 : "the dielectric constant function.")
330 : END IF
331 :
332 160 : CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid)
333 :
334 : ! make sure rho is in the real space
335 160 : CALL pw_pool_orig%create_pw(rho_elec_rs)
336 160 : CALL pw_transfer(rho, rho_elec_rs)
337 : ! expand rho_elec
338 160 : CALL pw_pool_xpndd%create_pw(rho_elec_rs_xpndd)
339 : CALL pw_expand(neumann_directions, recv_msgs_bnds, dests_expand, srcs_expand, flipg_stat, bounds_shftd, &
340 160 : rho_elec_rs, rho_elec_rs_xpndd)
341 :
342 160 : IF (PRESENT(rho_core)) THEN
343 : ! make sure rho_core is in the real space
344 160 : CALL pw_pool_orig%create_pw(rho_core_rs)
345 160 : CALL pw_transfer(rho_core, rho_core_rs)
346 : ! expand rho_core
347 160 : CALL pw_pool_xpndd%create_pw(rho_core_rs_xpndd)
348 : CALL pw_expand(neumann_directions, recv_msgs_bnds, dests_expand, srcs_expand, flipg_stat, bounds_shftd, &
349 160 : rho_core_rs, rho_core_rs_xpndd)
350 :
351 160 : IF (dielectric%params%dielec_core_correction) THEN
352 : ! use (rho_elec - rho_core) to compute dielectric
353 160 : CALL pw_axpy(rho_core_rs_xpndd, rho_elec_rs_xpndd, -2.0_dp)
354 : ELSE
355 0 : CALL pw_axpy(rho_core_rs_xpndd, rho_elec_rs_xpndd, -1.0_dp)
356 : END IF
357 160 : CALL pw_pool_orig%give_back_pw(rho_core_rs)
358 160 : CALL pw_pool_xpndd%give_back_pw(rho_core_rs_xpndd)
359 : ELSE
360 0 : CPABORT("For dielectric constant larger than 1, rho_core has to be present.")
361 : END IF
362 :
363 : ! calculate the dielectric constant
364 160 : SELECT CASE (dielec_functiontype)
365 : CASE (rho_dependent)
366 : CALL dielectric_constant_sccs(rho_elec_rs_xpndd, dielectric%eps, dielectric%deps_drho, &
367 160 : eps0, rho_max, rho_min)
368 : CASE (spatially_dependent)
369 0 : IF (times_called .EQ. 0) THEN
370 0 : CALL dielectric_constant_spatially_dependent(dielectric%eps, pw_pool_xpndd, dielectric%params)
371 : END IF
372 : CASE (spatially_rho_dependent)
373 : CALL dielectric_constant_spatially_rho_dependent(rho_elec_rs_xpndd, dielectric%eps, &
374 160 : dielectric%deps_drho, pw_pool_xpndd, dielectric%params)
375 : END SELECT
376 :
377 : ! derivatives
378 : IF ((dielec_functiontype .EQ. rho_dependent) .OR. &
379 160 : (dielec_functiontype .EQ. spatially_rho_dependent) .OR. &
380 : ((dielec_functiontype .EQ. spatially_dependent) .AND. times_called .EQ. 0)) THEN
381 152 : SELECT CASE (derivative_method)
382 : CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
383 152 : CALL pw_pool_xpndd%create_pw(ln_eps)
384 62009272 : ln_eps%array = LOG(dielectric%eps%array)
385 : CASE (derivative_fft_use_deps)
386 32 : DO i = 1, 3
387 24 : CALL pw_pool_xpndd%create_pw(deps(i))
388 32 : CALL pw_zero(deps(i))
389 : END DO
390 : CASE (derivative_fft_use_drho)
391 160 : DO i = 1, 3
392 0 : CALL pw_pool_xpndd%create_pw(deps(i))
393 0 : CALL pw_pool_xpndd%create_pw(drho(i))
394 0 : CALL pw_zero(deps(i))
395 0 : CALL pw_zero(drho(i))
396 : END DO
397 : END SELECT
398 :
399 128 : SELECT CASE (derivative_method)
400 : CASE (derivative_cd3)
401 128 : CALL derive_fdm_cd3(ln_eps, dielectric%dln_eps, diel_rs_grid)
402 : CASE (derivative_cd5)
403 0 : CALL derive_fdm_cd5(ln_eps, dielectric%dln_eps, diel_rs_grid)
404 : CASE (derivative_cd7)
405 24 : CALL derive_fdm_cd7(ln_eps, dielectric%dln_eps, diel_rs_grid)
406 : CASE (derivative_fft)
407 0 : CALL derive_fft(ln_eps, dielectric%dln_eps, pw_pool_xpndd)
408 : CASE (derivative_fft_use_deps)
409 : ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
410 8 : CALL derive_fft(dielectric%eps, deps, pw_pool_xpndd)
411 :
412 32 : lb(1:3) = pw_pool_xpndd%pw_grid%bounds_local(1, 1:3)
413 32 : ub(1:3) = pw_pool_xpndd%pw_grid%bounds_local(2, 1:3)
414 : ! damp oscillations cuased by fft-based derivative in the region
415 : ! where electron density is zero
416 32 : DO idir = 1, 3
417 2904 : DO k = lb(3), ub(3)
418 348504 : DO j = lb(2), ub(2)
419 21084480 : DO i = lb(1), ub(1)
420 21081600 : IF (ABS(dielectric%deps_drho%array(i, j, k)) .LE. small_value) THEN
421 18486864 : deps(idir)%array(i, j, k) = 0.0_dp
422 : END IF
423 : END DO
424 : END DO
425 : END DO
426 21084512 : dielectric%dln_eps(idir)%array = deps(idir)%array/dielectric%eps%array
427 : END DO
428 : CASE (derivative_fft_use_drho)
429 : ! \Nabla \eps = \Nabla \rho \cdot \frac{\partial \eps}{\partial \rho}
430 : ! \Nabla ln(\eps) = \frac{\Nabla \eps}{\eps}
431 0 : CALL derive_fft(rho_elec_rs_xpndd, drho, pw_pool_xpndd)
432 160 : DO i = 1, 3
433 0 : deps(i)%array = drho(i)%array*dielectric%deps_drho%array
434 0 : dielectric%dln_eps(i)%array = deps(i)%array/dielectric%eps%array
435 : END DO
436 : END SELECT
437 :
438 152 : SELECT CASE (derivative_method)
439 : CASE (derivative_cd3, derivative_cd5, derivative_cd7, derivative_fft)
440 152 : CALL pw_pool_xpndd%give_back_pw(ln_eps)
441 : CASE (derivative_fft_use_deps)
442 32 : DO i = 1, 3
443 32 : CALL pw_pool_xpndd%give_back_pw(deps(i))
444 : END DO
445 : CASE (derivative_fft_use_drho)
446 160 : DO i = 1, 3
447 0 : CALL pw_pool_xpndd%give_back_pw(drho(i))
448 0 : CALL pw_pool_xpndd%give_back_pw(deps(i))
449 : END DO
450 : END SELECT
451 : END IF
452 :
453 160 : CALL pw_pool_orig%give_back_pw(rho_elec_rs)
454 160 : CALL pw_pool_xpndd%give_back_pw(rho_elec_rs_xpndd)
455 160 : CALL pw_pool_release(pw_pool_xpndd)
456 :
457 160 : dielectric%params%times_called = dielectric%params%times_called + 1
458 :
459 160 : CALL timestop(handle)
460 :
461 160 : END SUBROUTINE dielectric_compute_neumann_${kind}$
462 : #:endfor
463 :
464 : ! **************************************************************************************************
465 : !> \brief calculates the dielectric constant as a function of the electronic density
466 : !> [see O. Andreussi, I. Dabo, and N. Marzari, J. Chem. Phys., 136, 064102 (2012)]
467 : !> \param rho electron density
468 : !> \param eps dielectric constant
469 : !> \param deps_drho derivative of the dielectric constant wrt the density
470 : !> \param eps0 dielectric constant in the bulk of the solvent
471 : !> \param rho_max upper density threshold
472 : !> \param rho_min lower density threshold
473 : !> \par History
474 : !> 06.2014 created [Hossein Bani-Hashemian]
475 : !> \author Mohammad Hossein Bani-Hashemian
476 : ! **************************************************************************************************
477 426 : SUBROUTINE dielectric_constant_sccs(rho, eps, deps_drho, eps0, rho_max, rho_min)
478 :
479 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho, eps, deps_drho
480 : REAL(KIND=dp), INTENT(IN) :: eps0, rho_max, rho_min
481 :
482 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dielectric_constant_sccs'
483 :
484 : INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
485 : ub2, ub3
486 : INTEGER, DIMENSION(2, 3) :: bounds_local
487 : REAL(KIND=dp) :: denom, t
488 :
489 426 : CALL timeset(routineN, handle)
490 :
491 426 : IF (eps0 .LT. 1.0_dp) THEN
492 0 : CPABORT("The dielectric constant has to be greater than or equal to 1.")
493 : END IF
494 :
495 4260 : bounds_local = rho%pw_grid%bounds_local
496 426 : lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
497 426 : lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
498 426 : lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
499 :
500 426 : denom = LOG(rho_max) - LOG(rho_min)
501 32538 : DO k = lb3, ub3
502 2564714 : DO j = lb2, ub2
503 112197668 : DO i = lb1, ub1
504 112165556 : IF (rho%array(i, j, k) .LT. rho_min) THEN
505 93348789 : eps%array(i, j, k) = eps0
506 93348789 : deps_drho%array(i, j, k) = 0.0_dp
507 16284591 : ELSE IF (rho%array(i, j, k) .GT. rho_max) THEN
508 7611877 : eps%array(i, j, k) = 1.0_dp
509 7611877 : deps_drho%array(i, j, k) = 0.0_dp
510 : ELSE
511 8672714 : t = twopi*(LOG(rho_max) - LOG(rho%array(i, j, k)))/denom
512 8672714 : eps%array(i, j, k) = EXP(LOG(eps0)*(t - SIN(t))/twopi)
513 8672714 : deps_drho%array(i, j, k) = -eps%array(i, j, k)*LOG(eps0)*(1.0_dp - COS(t))/(denom*rho%array(i, j, k))
514 : END IF
515 : END DO
516 : END DO
517 : END DO
518 :
519 426 : CALL timestop(handle)
520 :
521 426 : END SUBROUTINE dielectric_constant_sccs
522 :
523 : ! **************************************************************************************************
524 : !> \brief creates an axis-aligned cuboidal dielectric region
525 : !> \param eps dielectric constant function
526 : !> \param dielec_const dielectric constant inside the region
527 : !> \param pw_pool pool of planewave grid
528 : !> \param zeta the mollifier's width
529 : !> \param x_xtnt the x extent of the cuboidal region
530 : !> \param y_xtnt the y extent of the cuboidal region
531 : !> \param z_xtnt the z extent of the cuboidal region
532 : !> \param x_glbl x grid vector of the simulation box
533 : !> \param y_glbl y grid vector of the simulation box
534 : !> \param z_glbl z grid vector of the simulation box
535 : !> \param x_locl x grid vector of the simulation box local to this process
536 : !> \param y_locl y grid vector of the simulation box local to this process
537 : !> \param z_locl z grid vector of the simulation box local to this process
538 : !> \par History
539 : !> 07.2015 created [Hossein Bani-Hashemian]
540 : !> \author Mohammad Hossein Bani-Hashemian
541 : ! **************************************************************************************************
542 2 : SUBROUTINE dielectric_constant_aa_cuboidal(eps, dielec_const, pw_pool, zeta, &
543 : x_xtnt, y_xtnt, z_xtnt, &
544 : x_glbl, y_glbl, z_glbl, &
545 : x_locl, y_locl, z_locl)
546 :
547 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: eps
548 : REAL(KIND=dp), INTENT(IN) :: dielec_const
549 : TYPE(pw_pool_type), POINTER :: pw_pool
550 : REAL(KIND=dp), INTENT(IN) :: zeta
551 : REAL(dp), DIMENSION(2), INTENT(IN) :: x_xtnt, y_xtnt, z_xtnt
552 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_glbl, y_glbl, z_glbl, x_locl, y_locl, &
553 : z_locl
554 :
555 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dielectric_constant_aa_cuboidal'
556 : LOGICAL, DIMENSION(6), PARAMETER :: test_forb_xtnts = .TRUE.
557 :
558 : INTEGER :: handle, i, j, k, lb1, lb2, lb3, &
559 : n_forb_xtnts, ub1, ub2, ub3
560 : INTEGER, DIMENSION(2, 3) :: bounds_local
561 : LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
562 : forb_xtnt4, forb_xtnt5, forb_xtnt6
563 : REAL(KIND=dp) :: dx, dy, dz
564 : TYPE(pw_grid_type), POINTER :: pw_grid
565 : TYPE(pw_r3d_rs_type) :: eps_tmp
566 :
567 2 : CALL timeset(routineN, handle)
568 :
569 2 : IF (dielec_const .LT. 1.0_dp) THEN
570 0 : CPABORT("The dielectric constant has to be greater than or equal to 1.")
571 : END IF
572 :
573 2 : pw_grid => eps%pw_grid
574 :
575 2 : dx = pw_grid%dr(1); dy = pw_grid%dr(2); dz = pw_grid%dr(3)
576 :
577 2 : forb_xtnt1 = x_xtnt(1) .LT. x_glbl(LBOUND(x_glbl, 1))
578 2 : forb_xtnt2 = x_xtnt(2) .GT. x_glbl(UBOUND(x_glbl, 1)) + dx
579 2 : forb_xtnt3 = y_xtnt(1) .LT. y_glbl(LBOUND(y_glbl, 1))
580 2 : forb_xtnt4 = y_xtnt(2) .GT. y_glbl(UBOUND(y_glbl, 1)) + dy
581 2 : forb_xtnt5 = z_xtnt(1) .LT. z_glbl(LBOUND(z_glbl, 1))
582 2 : forb_xtnt6 = z_xtnt(2) .GT. z_glbl(UBOUND(z_glbl, 1)) + dz
583 : n_forb_xtnts = COUNT((/forb_xtnt1, forb_xtnt2, forb_xtnt3, &
584 14 : forb_xtnt4, forb_xtnt5, forb_xtnt6/) .EQV. test_forb_xtnts)
585 2 : IF (n_forb_xtnts .GT. 0) THEN
586 : CALL cp_abort(__LOCATION__, &
587 : "The given extents for the dielectric region are outside the range of "// &
588 0 : "the simulation cell.")
589 : END IF
590 :
591 2 : CALL pw_pool%create_pw(eps_tmp)
592 2 : CALL pw_copy(eps, eps_tmp)
593 :
594 20 : bounds_local = pw_grid%bounds_local
595 2 : lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
596 2 : lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
597 2 : lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
598 :
599 146 : DO k = lb3, ub3
600 10514 : DO j = lb2, ub2
601 383760 : DO i = lb1, ub1
602 : IF ((x_locl(i) .GE. x_xtnt(1)) .AND. (x_locl(i) .LE. x_xtnt(2)) .AND. &
603 : (y_locl(j) .GE. y_xtnt(1)) .AND. (y_locl(j) .LE. y_xtnt(2)) .AND. &
604 383616 : (z_locl(k) .GE. z_xtnt(1)) .AND. (z_locl(k) .LE. z_xtnt(2))) THEN
605 35721 : eps_tmp%array(i, j, k) = dielec_const
606 : END IF
607 : END DO
608 : END DO
609 : END DO
610 :
611 2 : CALL pw_mollifier(pw_pool, zeta, x_glbl, y_glbl, z_glbl, eps_tmp, eps)
612 2 : CALL pw_pool%give_back_pw(eps_tmp)
613 :
614 2 : CALL timestop(handle)
615 :
616 2 : END SUBROUTINE dielectric_constant_aa_cuboidal
617 :
618 : ! **************************************************************************************************
619 : !> \brief creates an x-axis aligned annular dielectric region
620 : !> \param eps dielectric constant function
621 : !> \param dielec_const dielectric constant inside the region
622 : !> \param pw_pool pool of planewave grid
623 : !> \param zeta the mollifier's width
624 : !> \param x_xtnt the x extent of the annular region
625 : !> \param base_center the y and z coordinates of the cylinder's base
626 : !> \param base_radii the radii of the annulus' base
627 : !> \param x_glbl x grid vector of the simulation box
628 : !> \param y_glbl y grid vector of the simulation box
629 : !> \param z_glbl z grid vector of the simulation box
630 : !> \param x_locl x grid vector of the simulation box local to this process
631 : !> \param y_locl y grid vector of the simulation box local to this process
632 : !> \param z_locl z grid vector of the simulation box local to this process
633 : !> \par History
634 : !> 07.2015 created [Hossein Bani-Hashemian]
635 : !> \author Mohammad Hossein Bani-Hashemian
636 : ! **************************************************************************************************
637 26 : SUBROUTINE dielectric_constant_xaa_annular(eps, dielec_const, pw_pool, zeta, &
638 : x_xtnt, base_center, base_radii, &
639 : x_glbl, y_glbl, z_glbl, &
640 : x_locl, y_locl, z_locl)
641 :
642 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: eps
643 : REAL(KIND=dp), INTENT(IN) :: dielec_const
644 : TYPE(pw_pool_type), POINTER :: pw_pool
645 : REAL(dp), INTENT(IN) :: zeta
646 : REAL(dp), DIMENSION(2), INTENT(IN) :: x_xtnt, base_center, base_radii
647 : REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_glbl, y_glbl, z_glbl, x_locl, y_locl, &
648 : z_locl
649 :
650 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dielectric_constant_xaa_annular'
651 : LOGICAL, DIMENSION(6), PARAMETER :: test_forb_xtnts = .TRUE.
652 :
653 : INTEGER :: handle, i, j, k, lb1, lb2, lb3, &
654 : n_forb_xtnts, ub1, ub2, ub3
655 : INTEGER, DIMENSION(2, 3) :: bounds_local
656 : LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
657 : forb_xtnt4, forb_xtnt5, forb_xtnt6
658 : REAL(KIND=dp) :: bctry, bctrz, distsq, dx, dy, dz
659 : TYPE(pw_grid_type), POINTER :: pw_grid
660 : TYPE(pw_r3d_rs_type) :: eps_tmp
661 :
662 26 : CALL timeset(routineN, handle)
663 :
664 26 : IF (dielec_const .LT. 1.0_dp) THEN
665 0 : CPABORT("The dielectric constant has to be greater than or equal to 1.")
666 : END IF
667 :
668 26 : pw_grid => eps%pw_grid
669 :
670 26 : bctry = base_center(1); bctrz = base_center(2)
671 26 : dx = pw_grid%dr(1); dy = pw_grid%dr(2); dz = pw_grid%dr(3)
672 :
673 26 : forb_xtnt1 = x_xtnt(1) .LT. x_glbl(LBOUND(x_glbl, 1))
674 26 : forb_xtnt2 = x_xtnt(2) .GT. x_glbl(UBOUND(x_glbl, 1)) + dx
675 104 : forb_xtnt3 = bctry - MAXVAL(base_radii) .LT. y_glbl(LBOUND(y_glbl, 1))
676 104 : forb_xtnt4 = bctry + MAXVAL(base_radii) .GT. y_glbl(UBOUND(y_glbl, 1)) + dy
677 104 : forb_xtnt5 = bctrz - MAXVAL(base_radii) .LT. z_glbl(LBOUND(z_glbl, 1))
678 104 : forb_xtnt6 = bctrz + MAXVAL(base_radii) .GT. z_glbl(UBOUND(z_glbl, 1)) + dz
679 : n_forb_xtnts = COUNT((/forb_xtnt1, forb_xtnt2, forb_xtnt3, &
680 182 : forb_xtnt4, forb_xtnt5, forb_xtnt6/) .EQV. test_forb_xtnts)
681 26 : IF (n_forb_xtnts .GT. 0) THEN
682 : CALL cp_abort(__LOCATION__, &
683 : "The given extents for the dielectric region are outside the range of "// &
684 0 : "the simulation cell.")
685 : END IF
686 :
687 26 : CALL pw_pool%create_pw(eps_tmp)
688 26 : CALL pw_copy(eps, eps_tmp)
689 :
690 260 : bounds_local = pw_grid%bounds_local
691 26 : lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
692 26 : lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
693 26 : lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
694 :
695 1898 : DO k = lb3, ub3
696 136682 : DO j = lb2, ub2
697 4988880 : DO i = lb1, ub1
698 4852224 : distsq = (y_locl(j) - bctry)**2 + (z_locl(k) - bctrz)**2
699 : IF ((x_locl(i) .GE. x_xtnt(1)) .AND. (x_locl(i) .LE. x_xtnt(2)) .AND. &
700 34100352 : (distsq .GE. MINVAL(base_radii)**2) .AND. (distsq .LE. MAXVAL(base_radii)**2)) THEN
701 390159 : eps_tmp%array(i, j, k) = dielec_const
702 : END IF
703 : END DO
704 : END DO
705 : END DO
706 :
707 26 : CALL pw_mollifier(pw_pool, zeta, x_glbl, y_glbl, z_glbl, eps_tmp, eps)
708 26 : CALL pw_pool%give_back_pw(eps_tmp)
709 :
710 26 : CALL timestop(handle)
711 :
712 26 : END SUBROUTINE dielectric_constant_xaa_annular
713 :
714 : ! **************************************************************************************************
715 : !> \brief Sets up density independent dielectric regions
716 : !> \param eps dielectric constant function
717 : !> \param pw_pool pool of planewave grid
718 : !> \param dielectric_params dielectric parameters read from input file
719 : !> \par History
720 : !> 07.2015 created [Hossein Bani-Hashemian]
721 : !> \author Mohammad Hossein Bani-Hashemian
722 : ! **************************************************************************************************
723 26 : SUBROUTINE dielectric_constant_spatially_dependent(eps, pw_pool, dielectric_params)
724 :
725 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: eps
726 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
727 : TYPE(dielectric_parameters), INTENT(IN) :: dielectric_params
728 :
729 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dielectric_constant_spatially_dependent'
730 :
731 : INTEGER :: handle, j, n_aa_cuboidal, &
732 : n_dielectric_region, n_xaa_annular
733 : REAL(dp) :: dielec_const, zeta
734 26 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: x_glbl, x_locl, y_glbl, y_locl, z_glbl, &
735 26 : z_locl
736 : REAL(dp), DIMENSION(2) :: base_center, base_radii
737 : TYPE(pw_grid_type), POINTER :: pw_grid
738 :
739 26 : CALL timeset(routineN, handle)
740 :
741 4988906 : eps%array = dielectric_params%eps0
742 :
743 26 : n_aa_cuboidal = dielectric_params%n_aa_cuboidal
744 26 : n_xaa_annular = dielectric_params%n_xaa_annular
745 26 : pw_grid => pw_pool%pw_grid
746 : CALL setup_grid_axes(pw_grid, x_glbl, y_glbl, z_glbl, x_locl, y_locl, z_locl)
747 26 : n_dielectric_region = n_aa_cuboidal + n_xaa_annular
748 :
749 26 : IF (n_dielectric_region .EQ. 0) THEN
750 0 : CPABORT("No density independent dielectric region is defined.")
751 : END IF
752 :
753 28 : DO j = 1, n_aa_cuboidal
754 2 : dielec_const = dielectric_params%aa_cuboidal_eps(j)
755 2 : zeta = dielectric_params%aa_cuboidal_zeta(j)
756 :
757 : CALL dielectric_constant_aa_cuboidal(eps, dielec_const, pw_pool, zeta, &
758 : dielectric_params%aa_cuboidal_xxtnt(:, j), &
759 : dielectric_params%aa_cuboidal_yxtnt(:, j), &
760 : dielectric_params%aa_cuboidal_zxtnt(:, j), &
761 : x_glbl, y_glbl, z_glbl, &
762 28 : x_locl, y_locl, z_locl)
763 : END DO
764 :
765 52 : DO j = 1, n_xaa_annular
766 78 : base_center = dielectric_params%xaa_annular_bctr(:, j)
767 78 : base_radii = dielectric_params%xaa_annular_brad(:, j)
768 26 : dielec_const = dielectric_params%xaa_annular_eps(j)
769 26 : zeta = dielectric_params%xaa_annular_zeta(j)
770 :
771 : CALL dielectric_constant_xaa_annular(eps, dielec_const, pw_pool, zeta, &
772 : dielectric_params%xaa_annular_xxtnt(:, j), &
773 : base_center, base_radii, &
774 : x_glbl, y_glbl, z_glbl, &
775 52 : x_locl, y_locl, z_locl)
776 : END DO
777 :
778 26 : CALL timestop(handle)
779 :
780 26 : END SUBROUTINE dielectric_constant_spatially_dependent
781 :
782 : ! **************************************************************************************************
783 : !> \brief Sets up various dielectric constant regions. Using the sccs switching
784 : !> function the dielectric constant smoothly varies between 1, where the density
785 : !> is present and the values inside the dielectric regions, otherwise.
786 : !> \param rho electron density
787 : !> \param eps dielectric constant function
788 : !> \param deps_drho derivative of the dielectric constant wrt the density
789 : !> \param pw_pool pool of planewave grid
790 : !> \param dielectric_params dielectric parameters read from input file
791 : !> \par History
792 : !> 07.2015 created [Hossein Bani-Hashemian]
793 : !> \author Mohammad Hossein Bani-Hashemian
794 : ! **************************************************************************************************
795 24 : SUBROUTINE dielectric_constant_spatially_rho_dependent(rho, eps, deps_drho, &
796 : pw_pool, dielectric_params)
797 :
798 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho, eps, deps_drho
799 : TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
800 : TYPE(dielectric_parameters), INTENT(IN) :: dielectric_params
801 :
802 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dielectric_constant_spatially_rho_dependent'
803 :
804 : INTEGER :: handle
805 : TYPE(pw_r3d_rs_type) :: dswch_func_drho, eps_sptldep, swch_func
806 :
807 24 : CALL timeset(routineN, handle)
808 :
809 24 : IF (dielectric_params%eps0 .LT. 1.0_dp) THEN
810 0 : CPABORT("The dielectric constant has to be greater than or equal to 1.")
811 : END IF
812 :
813 24 : CALL pw_pool%create_pw(eps_sptldep)
814 24 : CALL pw_pool%create_pw(swch_func)
815 24 : CALL pw_pool%create_pw(dswch_func_drho)
816 24 : CALL pw_zero(eps_sptldep)
817 24 : CALL pw_zero(swch_func)
818 24 : CALL pw_zero(dswch_func_drho)
819 :
820 24 : CALL dielectric_constant_spatially_dependent(eps_sptldep, pw_pool, dielectric_params)
821 : CALL dielectric_constant_sccs(rho, swch_func, dswch_func_drho, 2.0_dp, &
822 24 : dielectric_params%rho_max, dielectric_params%rho_min)
823 :
824 4605144 : eps%array = ((swch_func%array - 1.0_dp)*(eps_sptldep%array - 1.0_dp)) + 1.0_dp
825 4605144 : deps_drho%array = dswch_func_drho%array*(eps_sptldep%array - 1.0_dp)
826 :
827 24 : CALL pw_pool%give_back_pw(dswch_func_drho)
828 24 : CALL pw_pool%give_back_pw(swch_func)
829 24 : CALL pw_pool%give_back_pw(eps_sptldep)
830 :
831 24 : CALL timestop(handle)
832 :
833 24 : END SUBROUTINE dielectric_constant_spatially_rho_dependent
834 :
835 : ! **************************************************************************************************
836 : !> \brief computes the derivative of a function using FFT
837 : !> \param f input funcition
838 : !> \param df derivative of f
839 : !> \param pw_pool pool of plane-wave grid
840 : ! **************************************************************************************************
841 5108 : SUBROUTINE derive_fft(f, df, pw_pool)
842 :
843 : TYPE(pw_r3d_rs_type), INTENT(IN) :: f
844 : TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(INOUT) :: df
845 : TYPE(pw_pool_type), POINTER :: pw_pool
846 :
847 : CHARACTER(LEN=*), PARAMETER :: routineN = 'derive_fft_r3d'
848 :
849 : INTEGER :: handle, i
850 : INTEGER, DIMENSION(3) :: nd
851 15324 : TYPE(pw_c1d_gs_type), DIMENSION(2) :: work_gs
852 :
853 5108 : CALL timeset(routineN, handle)
854 :
855 15324 : DO i = 1, 2
856 15324 : CALL pw_pool%create_pw(work_gs(i))
857 : END DO
858 :
859 5108 : CALL pw_transfer(f, work_gs(1))
860 20432 : DO i = 1, 3
861 15324 : nd(:) = 0
862 15324 : nd(i) = 1
863 15324 : CALL pw_copy(work_gs(1), work_gs(2))
864 15324 : CALL pw_derive(work_gs(2), nd(:))
865 20432 : CALL pw_transfer(work_gs(2), df(i))
866 : END DO
867 :
868 15324 : DO i = 1, 2
869 15324 : CALL pw_pool%give_back_pw(work_gs(i))
870 : END DO
871 :
872 5108 : CALL timestop(handle)
873 :
874 5108 : END SUBROUTINE derive_fft
875 :
876 : END MODULE dielectric_methods
|
