(git:374b731)
Loading...
Searching...
No Matches
dirichlet_bc_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 subroutines for defining and creating Dirichlet type subdomains
10!> \par History
11!> 08.2014 created [Hossein Bani-Hashemian]
12!> 10.2015 completely revised [Hossein Bani-Hashemian]
13!> \author Mohammad Hossein Bani-Hashemian
14! **************************************************************************************************
15MODULE dirichlet_bc_methods
16
17 USE cp_log_handling, ONLY: cp_get_default_logger,&
18 cp_logger_get_default_unit_nr,&
19 cp_logger_type
20 USE dirichlet_bc_types, ONLY: &
21 aa_cuboidal, aa_planar, circumscribed, cylindrical, inscribed, planar, &
22 dbc_parameters_dealloc, dirichlet_bc_p_type, dirichlet_bc_type, tile_p_type, x_axis, &
23 xy_plane, xz_plane, y_axis, yz_plane, z_axis
24 USE kinds, ONLY: dp
25 USE mathconstants, ONLY: pi,&
26 twopi
27 USE mathlib, ONLY: inv_3x3,&
28 vector_product
29 USE physcon, ONLY: angstrom
30 USE ps_implicit_types, ONLY: mixed_bc,&
31 mixed_periodic_bc,&
32 neumann_bc,&
33 periodic_bc
34 USE pw_grid_types, ONLY: pw_grid_type
35 USE pw_methods, ONLY: pw_axpy,&
36 pw_zero
37 USE pw_poisson_types, ONLY: pw_poisson_parameter_type
38 USE pw_pool_types, ONLY: pw_pool_type
39 USE pw_types, ONLY: pw_r3d_rs_type
40 USE rs_methods, ONLY: setup_grid_axes
41#include "../base/base_uses.f90"
42
43 IMPLICIT NONE
44 PRIVATE
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dirichlet_bc_methods'
46
47 PUBLIC dirichlet_boundary_region_setup
48
49 REAL(dp), PARAMETER, PRIVATE :: small_value = 1.0e-8_dp
50 INTEGER, PARAMETER, PRIVATE :: FWROT = 1, &
51 bwrot = -1
52
53CONTAINS
54
55! **************************************************************************************************
56!> \brief Sets up the Dirichlet boundary condition
57!> \param pw_pool pool of plane wave grid
58!> \param poisson_params poisson_env parameters
59!> \param dbcs the DBC region to be created
60!> \par History
61!> 10.2014 created [Hossein Bani-Hashemian]
62!> \author Mohammad Hossein Bani-Hashemian
63! **************************************************************************************************
64 SUBROUTINE dirichlet_boundary_region_setup(pw_pool, poisson_params, dbcs)
65
66 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
67 TYPE(pw_poisson_parameter_type), INTENT(INOUT) :: poisson_params
68 TYPE(dirichlet_bc_p_type), ALLOCATABLE, &
69 DIMENSION(:), INTENT(INOUT) :: dbcs
70
71 CHARACTER(LEN=*), PARAMETER :: routinen = 'dirichlet_boundary_region_setup'
72
73 INTEGER :: apx_type, dbc_id, handle, ind_end, ind_start, j, l, n_aa_cuboidal, &
74 n_aa_cylindrical, n_aa_planar, n_dbcs, n_planar, parallel_axis, parallel_plane, unit_nr
75 INTEGER, DIMENSION(3) :: n_prtn, npts
76 INTEGER, DIMENSION(:), POINTER :: aa_cylindrical_nsides
77 LOGICAL :: is_periodic, verbose
78 REAL(dp) :: base_radius, delta_alpha, frequency, &
79 oscillating_fraction, phase, sigma, &
80 thickness, v_d
81 REAL(dp), ALLOCATABLE, DIMENSION(:) :: x_glbl, x_locl, y_glbl, y_locl, z_glbl, &
82 z_locl
83 REAL(dp), DIMENSION(2) :: base_center
84 REAL(dp), DIMENSION(2, 3) :: cell_xtnts
85 REAL(dp), DIMENSION(3) :: dr
86 TYPE(cp_logger_type), POINTER :: logger
87 TYPE(pw_grid_type), POINTER :: pw_grid
88
89 CALL timeset(routinen, handle)
90
91 logger => cp_get_default_logger()
92 IF (logger%para_env%is_source()) THEN
93 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
94 ELSE
95 unit_nr = -1
96 END IF
97
98 dr = pw_pool%pw_grid%dr
99 npts = pw_pool%pw_grid%npts
100 cell_xtnts = 0.0_dp
101 cell_xtnts(2, 1) = npts(1)*dr(1)
102 cell_xtnts(2, 2) = npts(2)*dr(2)
103 cell_xtnts(2, 3) = npts(3)*dr(3)
104
105 verbose = poisson_params%dbc_params%verbose_output
106
107 n_aa_planar = poisson_params%dbc_params%n_aa_planar
108 n_aa_cuboidal = poisson_params%dbc_params%n_aa_cuboidal
109 n_planar = poisson_params%dbc_params%n_planar
110 n_aa_cylindrical = poisson_params%dbc_params%n_aa_cylindrical
111 aa_cylindrical_nsides => poisson_params%dbc_params%aa_cylindrical_nsides
112 pw_grid => pw_pool%pw_grid
113 CALL setup_grid_axes(pw_grid, x_glbl, y_glbl, z_glbl, x_locl, y_locl, z_locl)
114 n_dbcs = n_aa_planar + n_aa_cuboidal + n_planar + sum(aa_cylindrical_nsides)
115
116 SELECT CASE (poisson_params%ps_implicit_params%boundary_condition)
117 CASE (mixed_bc, mixed_periodic_bc)
118 IF (n_dbcs .EQ. 0) THEN
119 cpabort("No Dirichlet region is defined.")
120 END IF
121
122 ALLOCATE (dbcs(n_dbcs))
123 IF (unit_nr .GT. 0) THEN
124 WRITE (unit_nr, '(/,T3,A,A,/,T3,A)') "POISSON| IMPLICIT (GENERALIZED) SOLVER ", repeat('-', 39), &
125 "POISSON| Preparing Dirichlet regions ..."
126 END IF
127
128 DO j = 1, n_aa_planar
129 ALLOCATE (dbcs(j)%dirichlet_bc)
130 n_prtn = poisson_params%dbc_params%aa_planar_nprtn(:, j)
131 dbc_id = aa_planar + j
132 v_d = poisson_params%dbc_params%aa_planar_vD(j)
133 frequency = poisson_params%dbc_params%aa_planar_frequency(j)
134 oscillating_fraction = poisson_params%dbc_params%aa_planar_osc_frac(j)
135 phase = poisson_params%dbc_params%aa_planar_phase(j)
136 sigma = poisson_params%dbc_params%aa_planar_sigma(j)
137 thickness = poisson_params%dbc_params%aa_planar_thickness(j)
138 parallel_plane = poisson_params%dbc_params%aa_planar_pplane(j)
139 is_periodic = poisson_params%dbc_params%aa_planar_is_periodic(j)
140
141 IF (unit_nr .GT. 0) THEN
142 WRITE (unit_nr, '(T3,A,I5)') "POISSON| Dirichlet region", j
143 WRITE (unit_nr, '(T3,A)') "POISSON| type : axis-aligned planar"
144 WRITE (unit_nr, '(T3,A,E13.4,2X,A)') "POISSON| applied potential :", v_d, "[Eh/e]"
145 END IF
146 CALL aa_planar_dbc_setup(cell_xtnts, parallel_plane, &
147 poisson_params%dbc_params%aa_planar_xxtnt(:, j), &
148 poisson_params%dbc_params%aa_planar_yxtnt(:, j), &
149 poisson_params%dbc_params%aa_planar_zxtnt(:, j), &
150 thickness, sigma, v_d, oscillating_fraction, frequency, &
151 phase, dbc_id, verbose, dbcs(j)%dirichlet_bc)
152
153 CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
154 x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
155 END DO
156
157 l = n_aa_planar
158 DO j = l + 1, l + n_aa_cuboidal
159 ALLOCATE (dbcs(j)%dirichlet_bc)
160 n_prtn = poisson_params%dbc_params%aa_cuboidal_nprtn(:, j - l)
161 dbc_id = aa_cuboidal + j - l
162 v_d = poisson_params%dbc_params%aa_cuboidal_vD(j - l)
163 frequency = poisson_params%dbc_params%aa_cuboidal_frequency(j - l)
164 oscillating_fraction = poisson_params%dbc_params%aa_cuboidal_osc_frac(j - l)
165 phase = poisson_params%dbc_params%aa_cuboidal_phase(j - l)
166 sigma = poisson_params%dbc_params%aa_cuboidal_sigma(j - l)
167 is_periodic = poisson_params%dbc_params%aa_cuboidal_is_periodic(j - l)
168
169 IF (unit_nr .GT. 0) THEN
170 WRITE (unit_nr, '(T3,A,I5)') "POISSON| Dirichlet region", j
171 WRITE (unit_nr, '(T3,A)') "POISSON| type : axis-aligned cuboidal"
172 WRITE (unit_nr, '(T3,A,E13.4,2X,A)') "POISSON| applied potential :", v_d, "[Eh/e]"
173 END IF
174 CALL aa_cuboidal_dbc_setup(cell_xtnts, &
175 poisson_params%dbc_params%aa_cuboidal_xxtnt(:, j - l), &
176 poisson_params%dbc_params%aa_cuboidal_yxtnt(:, j - l), &
177 poisson_params%dbc_params%aa_cuboidal_zxtnt(:, j - l), &
178 sigma, v_d, oscillating_fraction, frequency, &
179 phase, dbc_id, verbose, dbcs(j)%dirichlet_bc)
180
181 CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
182 x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
183 END DO
184
185 l = n_aa_planar + n_aa_cuboidal
186 DO j = l + 1, l + n_planar
187 ALLOCATE (dbcs(j)%dirichlet_bc)
188 n_prtn = 1
189 n_prtn(1:2) = poisson_params%dbc_params%planar_nprtn(:, j - l)
190 dbc_id = planar + j - l
191 v_d = poisson_params%dbc_params%planar_vD(j - l)
192 frequency = poisson_params%dbc_params%planar_frequency(j - l)
193 oscillating_fraction = poisson_params%dbc_params%planar_osc_frac(j - l)
194 phase = poisson_params%dbc_params%planar_phase(j - l)
195 sigma = poisson_params%dbc_params%planar_sigma(j - l)
196 thickness = poisson_params%dbc_params%planar_thickness(j - l)
197 is_periodic = poisson_params%dbc_params%planar_is_periodic(j - l)
198
199 IF (unit_nr .GT. 0) THEN
200 WRITE (unit_nr, '(T3,A,I5)') "POISSON| Dirichlet region", j
201 WRITE (unit_nr, '(T3,A)') "POISSON| type : planar"
202 WRITE (unit_nr, '(T3,A,E13.4,2X,A)') "POISSON| applied potential :", v_d, "[Eh/e]"
203 END IF
204 CALL arbitrary_planar_dbc_setup(cell_xtnts, &
205 poisson_params%dbc_params%planar_Avtx(:, j - l), &
206 poisson_params%dbc_params%planar_Bvtx(:, j - l), &
207 poisson_params%dbc_params%planar_Cvtx(:, j - l), &
208 thickness, sigma, v_d, oscillating_fraction, frequency, &
209 phase, dbc_id, verbose, dbcs(j)%dirichlet_bc)
210
211 CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
212 x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
213 END DO
214
215 l = n_aa_planar + n_aa_cuboidal + n_planar
216 DO j = 1, n_aa_cylindrical
217 ind_start = l + 1
218 ind_end = l + aa_cylindrical_nsides(j)
219
220 n_prtn = 1
221 n_prtn(1:2) = poisson_params%dbc_params%aa_cylindrical_nprtn(:, j)
222 parallel_axis = poisson_params%dbc_params%aa_cylindrical_paxis(j)
223 apx_type = poisson_params%dbc_params%aa_cylindrical_apxtyp(j)
224 base_center = poisson_params%dbc_params%aa_cylindrical_bctr(:, j)
225 base_radius = poisson_params%dbc_params%aa_cylindrical_brad(j)
226 thickness = poisson_params%dbc_params%aa_cylindrical_thickness(j)
227 delta_alpha = poisson_params%dbc_params%aa_cylindrical_sgap(j)
228 sigma = poisson_params%dbc_params%aa_cylindrical_sigma(j)
229 v_d = poisson_params%dbc_params%aa_cylindrical_vD(j)
230 frequency = poisson_params%dbc_params%aa_cylindrical_frequency(j)
231 oscillating_fraction = poisson_params%dbc_params%aa_cylindrical_osc_frac(j)
232 phase = poisson_params%dbc_params%aa_cylindrical_phase(j)
233 is_periodic = poisson_params%dbc_params%aa_cylindrical_is_periodic(j)
234
235 IF (unit_nr .GT. 0) THEN
236 WRITE (unit_nr, '(T3,A,I5)') "POISSON| Dirichlet region", l + j
237 WRITE (unit_nr, '(T3,A)') "POISSON| type : axis-aligned cylindrical"
238 WRITE (unit_nr, '(T3,A,E13.4,2X,A)') "POISSON| applied potential :", v_d, "[Eh/e]"
239 END IF
240 CALL aa_cylindrical_dbc_setup(pw_pool, cell_xtnts, x_locl, y_locl, z_locl, apx_type, parallel_axis, &
241 poisson_params%dbc_params%aa_cylindrical_xtnt(:, j), &
242 base_center, base_radius, thickness, delta_alpha, sigma, v_d, &
243 oscillating_fraction, frequency, phase, &
244 n_prtn, is_periodic, verbose, dbcs(ind_start:ind_end))
245
246 l = l + aa_cylindrical_nsides(j)
247 END DO
248
249 CASE (periodic_bc, neumann_bc)
250 ! do nothing
251 END SELECT
252
253! we won't need parameters anymore so deallocate them
254 CALL dbc_parameters_dealloc(poisson_params%dbc_params)
255
256 CALL timestop(handle)
257
258 END SUBROUTINE dirichlet_boundary_region_setup
259
260! **************************************************************************************************
261!> \brief Partitions a dirichlet_bc_type into tiles
262!> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
263!> \param n_prtn vector of size 3 specifying the number of times that the x, y and
264!> z interval (defining the region) should be partitioned into
265!> \param pw_pool pool of plane-wave grid
266!> \param cell_xtnts extents of the simulation cell
267!> \param x_locl x grid vector of the simulation box local to this process
268!> \param y_locl y grid vector of the simulation box local to this process
269!> \param z_locl z grid vector of the simulation box local to this process
270!> \param is_periodic whether or not the tile is periodic with a periodicity of one unit cell
271!> \param verbose whether or not to print out the coordinates of the vertices
272!> \param dirichlet_bc the dirichlet_bc object to be partitioned
273!> \par History
274!> 10.2014 created [Hossein Bani-Hashemian]
275!> \author Mohammad Hossein Bani-Hashemian
276! **************************************************************************************************
277 SUBROUTINE dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
278 x_locl, y_locl, z_locl, is_periodic, verbose, dirichlet_bc)
279
280 REAL(dp), INTENT(IN) :: sigma
281 INTEGER, DIMENSION(3), INTENT(IN) :: n_prtn
282 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool
283 REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
284 REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_locl, y_locl, z_locl
285 LOGICAL, INTENT(IN) :: is_periodic, verbose
286 TYPE(dirichlet_bc_type), INTENT(INOUT), POINTER :: dirichlet_bc
287
288 CHARACTER(LEN=*), PARAMETER :: routinen = 'dirichlet_bc_partition'
289
290 INTEGER :: handle, i, k, n_tiles, unit_nr
291 REAL(dp) :: cyl_maxval, tile_volume
292 REAL(dp), ALLOCATABLE, DIMENSION(:) :: tile_maxvals
293 TYPE(cp_logger_type), POINTER :: logger
294 TYPE(pw_r3d_rs_type), POINTER :: cylinder_pw
295 TYPE(tile_p_type), DIMENSION(:), POINTER :: tiles
296
297 CALL timeset(routinen, handle)
298
299 logger => cp_get_default_logger()
300 IF (logger%para_env%is_source()) THEN
301 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
302 ELSE
303 unit_nr = -1
304 END IF
305
306 SELECT CASE (dirichlet_bc%dbc_geom)
307 CASE (planar)
308
309 CALL arbitrary_dbc_partition(dirichlet_bc%vertices, n_prtn(1:2), tiles)
310 n_tiles = SIZE(tiles)
311 dirichlet_bc%n_tiles = n_tiles
312 ALLOCATE (dirichlet_bc%tiles(n_tiles))
313 dirichlet_bc%tiles = tiles
314 IF ((unit_nr .GT. 0) .AND. verbose) WRITE (unit_nr, '(T7,A,I5)') "number of partitions : ", n_tiles
315
316 ALLOCATE (tile_maxvals(n_tiles))
317 tile_maxvals = 0.0_dp
318 DO k = 1, n_tiles
319 ALLOCATE (dirichlet_bc%tiles(k)%tile%tile_pw)
320 CALL pw_pool%create_pw(dirichlet_bc%tiles(k)%tile%tile_pw)
321 CALL pw_zero(dirichlet_bc%tiles(k)%tile%tile_pw)
322
323 IF ((unit_nr .GT. 0) .AND. verbose) THEN
324 WRITE (unit_nr, '(T7,A,I5)') "tile", k
325 DO i = 1, 8
326 WRITE (unit_nr, '(T10,A,I1,3F10.3)') &
327 " vertex ", i, angstrom*dirichlet_bc%tiles(k)%tile%vertices(:, i)
328 END DO
329 END IF
330
331 CALL arbitrary_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, dirichlet_bc%tiles(k)%tile%vertices, &
332 sigma, dirichlet_bc%tiles(k)%tile%tile_pw)
333 tile_volume = sum(dirichlet_bc%tiles(k)%tile%tile_pw%array)
334 CALL pw_pool%pw_grid%para%group%sum(tile_volume)
335 dirichlet_bc%tiles(k)%tile%volume = tile_volume
336 END DO
337
338 IF ((unit_nr .GT. 0) .AND. verbose) WRITE (unit_nr, '(T3,A)') repeat('=', 78)
339
340 DEALLOCATE (tiles)
341
342 CASE (cylindrical)
343
344 CALL arbitrary_dbc_partition(dirichlet_bc%vertices, n_prtn(1:2), tiles)
345 n_tiles = SIZE(tiles)
346 dirichlet_bc%n_tiles = n_tiles
347 ALLOCATE (dirichlet_bc%tiles(n_tiles))
348 dirichlet_bc%tiles = tiles
349 IF ((unit_nr .GT. 0) .AND. verbose) WRITE (unit_nr, '(T7,A,I5)') "number of partitions : ", n_tiles
350
351 NULLIFY (cylinder_pw)
352 ALLOCATE (cylinder_pw)
353 CALL pw_pool%create_pw(cylinder_pw)
354 CALL pw_zero(cylinder_pw)
355
356 DO k = 1, n_tiles
357 ALLOCATE (dirichlet_bc%tiles(k)%tile%tile_pw)
358 CALL pw_pool%create_pw(dirichlet_bc%tiles(k)%tile%tile_pw)
359 CALL pw_zero(dirichlet_bc%tiles(k)%tile%tile_pw)
360
361 IF ((unit_nr .GT. 0) .AND. verbose) THEN
362 WRITE (unit_nr, '(T7,A,I5)') "tile", k
363 DO i = 1, 8
364 WRITE (unit_nr, '(T10,A,I1,3F10.3)') &
365 " vertex ", i, angstrom*dirichlet_bc%tiles(k)%tile%vertices(:, i)
366 END DO
367 END IF
368
369 CALL arbitrary_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, dirichlet_bc%tiles(k)%tile%vertices, &
370 sigma, dirichlet_bc%tiles(k)%tile%tile_pw)
371 CALL pw_axpy(dirichlet_bc%tiles(k)%tile%tile_pw, cylinder_pw)
372 tile_volume = sum(dirichlet_bc%tiles(k)%tile%tile_pw%array)
373 CALL pw_pool%pw_grid%para%group%sum(tile_volume)
374 dirichlet_bc%tiles(k)%tile%volume = tile_volume
375 END DO
376
377 cyl_maxval = maxval(cylinder_pw%array)
378 CALL pw_pool%pw_grid%para%group%max(cyl_maxval)
379 IF (cyl_maxval .GT. 1.01_dp) THEN
380 CALL cp_abort(__location__, &
381 "Inaccurate approximation of unit step function at cylindrical Dirichlet region. "// &
382 "Increase the number of sides (N_SIDES) and/or the base radius. Alternatively, "// &
383 "one can increase DELTA_ALPHA (see the reference manual).")
384 END IF
385
386 CALL pw_pool%give_back_pw(cylinder_pw)
387 DEALLOCATE (cylinder_pw)
388
389 IF ((unit_nr .GT. 0) .AND. verbose) WRITE (unit_nr, '(T3,A)') repeat('=', 78)
390
391 DEALLOCATE (tiles)
392
393 CASE (aa_planar, aa_cuboidal)
394
395 CALL aa_dbc_partition(dirichlet_bc%vertices, n_prtn, tiles)
396 n_tiles = SIZE(tiles)
397 dirichlet_bc%n_tiles = n_tiles
398 ALLOCATE (dirichlet_bc%tiles(n_tiles))
399 dirichlet_bc%tiles = tiles
400 IF ((unit_nr .GT. 0) .AND. verbose) WRITE (unit_nr, '(T7,A,I5)') "number of partitions : ", n_tiles
401
402 DO k = 1, n_tiles
403 ALLOCATE (dirichlet_bc%tiles(k)%tile%tile_pw)
404 CALL pw_pool%create_pw(dirichlet_bc%tiles(k)%tile%tile_pw)
405 CALL pw_zero(dirichlet_bc%tiles(k)%tile%tile_pw)
406
407 IF ((unit_nr .GT. 0) .AND. verbose) THEN
408 WRITE (unit_nr, '(T7,A,I5)') "tile", k
409 DO i = 1, 8
410 WRITE (unit_nr, '(T10,A,I1,3F10.3)') &
411 " vertex ", i, angstrom*dirichlet_bc%tiles(k)%tile%vertices(:, i)
412 END DO
413 END IF
414
415 CALL aa_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, dirichlet_bc%tiles(k)%tile%vertices, &
416 sigma, is_periodic, dirichlet_bc%tiles(k)%tile%tile_pw)
417 tile_volume = sum(dirichlet_bc%tiles(k)%tile%tile_pw%array)
418 CALL pw_pool%pw_grid%para%group%sum(tile_volume)
419 dirichlet_bc%tiles(k)%tile%volume = tile_volume
420 END DO
421
422 IF ((unit_nr .GT. 0) .AND. verbose) WRITE (unit_nr, '(T3,A)') repeat('=', 78)
423
424 DEALLOCATE (tiles)
425
426 END SELECT
427
428 CALL timestop(handle)
429
430 END SUBROUTINE dirichlet_bc_partition
431
432! **************************************************************************************************
433!> \brief Creates an axis-aligned planar Dirichlet region.
434!> \param cell_xtnts extents of the simulation cell
435!> \param parallel_plane the coordinate plane that the region is parallel to
436!> \param x_xtnt the x extent of the planar region
437!> \param y_xtnt the y extent of the planar region
438!> \param z_xtnt the z extent of the planar region
439!> \param thickness the thickness of the region
440!> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
441!> \param v_D potential imposed at the Dirichlet part of the boundary (gate voltage)
442!> \param osc_frac ...
443!> \param frequency ...
444!> \param phase ...
445!> \param dbc_id unique ID for the Dirichlet region
446!> \param verbose whether or not to print out the coordinates of the vertices
447!> \param dirichlet_bc the dirichlet_bc object to be created
448!> \par History
449!> 08.2014 created [Hossein Bani-Hashemian]
450!> \author Mohammad Hossein Bani-Hashemian
451! **************************************************************************************************
452 SUBROUTINE aa_planar_dbc_setup(cell_xtnts, parallel_plane, x_xtnt, y_xtnt, z_xtnt, &
453 thickness, sigma, v_D, osc_frac, frequency, phase, dbc_id, verbose, dirichlet_bc)
454
455 REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
456 INTEGER, INTENT(IN) :: parallel_plane
457 REAL(dp), DIMENSION(2), INTENT(IN) :: x_xtnt, y_xtnt, z_xtnt
458 REAL(dp), INTENT(IN) :: thickness, sigma, v_d, osc_frac, &
459 frequency, phase
460 INTEGER, INTENT(IN) :: dbc_id
461 LOGICAL, INTENT(IN) :: verbose
462 TYPE(dirichlet_bc_type), INTENT(INOUT), POINTER :: dirichlet_bc
463
464 CHARACTER(LEN=*), PARAMETER :: routinen = 'aa_planar_dbc_setup'
465 LOGICAL, DIMENSION(6), PARAMETER :: test_forb_xtnts = .true.
466
467 INTEGER :: handle, i, n_forb_xtnts, unit_nr
468 LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
469 forb_xtnt4, forb_xtnt5, forb_xtnt6
470 REAL(dp) :: xlb, xub, ylb, yub, zlb, zub
471 TYPE(cp_logger_type), POINTER :: logger
472
473 CALL timeset(routinen, handle)
474
475 logger => cp_get_default_logger()
476 IF (logger%para_env%is_source()) THEN
477 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
478 ELSE
479 unit_nr = -1
480 END IF
481
482 xlb = x_xtnt(1); xub = x_xtnt(2)
483 ylb = y_xtnt(1); yub = y_xtnt(2)
484 zlb = z_xtnt(1); zub = z_xtnt(2)
485
486 SELECT CASE (parallel_plane)
487 CASE (xy_plane)
488 zlb = zlb - thickness*0.5_dp
489 zub = zub + thickness*0.5_dp
490 CASE (xz_plane)
491 ylb = ylb - thickness*0.5_dp
492 yub = yub + thickness*0.5_dp
493 CASE (yz_plane)
494 xlb = xlb - thickness*0.5_dp
495 xub = xub + thickness*0.5_dp
496 END SELECT
497
498 forb_xtnt1 = xlb .LT. cell_xtnts(1, 1)
499 forb_xtnt2 = xub .GT. cell_xtnts(2, 1)
500 forb_xtnt3 = ylb .LT. cell_xtnts(1, 2)
501 forb_xtnt4 = yub .GT. cell_xtnts(2, 2)
502 forb_xtnt5 = zlb .LT. cell_xtnts(1, 3)
503 forb_xtnt6 = zub .GT. cell_xtnts(2, 3)
504 n_forb_xtnts = count((/forb_xtnt1, forb_xtnt2, forb_xtnt3, &
505 forb_xtnt4, forb_xtnt5, forb_xtnt6/) .EQV. test_forb_xtnts)
506 IF (n_forb_xtnts .GT. 0) THEN
507 CALL cp_abort(__location__, &
508 "The given extents for the Dirichlet region are outside the range of "// &
509 "the simulation cell.")
510 END IF
511
512 dirichlet_bc%v_D = v_d
513 dirichlet_bc%osc_frac = osc_frac
514 dirichlet_bc%frequency = frequency
515 dirichlet_bc%phase = phase
516 dirichlet_bc%dbc_id = dbc_id
517 dirichlet_bc%dbc_geom = aa_planar
518 dirichlet_bc%smoothing_width = sigma
519 dirichlet_bc%n_tiles = 1
520
521 dirichlet_bc%vertices(1:3, 1) = (/xlb, ylb, zlb/)
522 dirichlet_bc%vertices(1:3, 2) = (/xub, ylb, zlb/)
523 dirichlet_bc%vertices(1:3, 3) = (/xub, yub, zlb/)
524 dirichlet_bc%vertices(1:3, 4) = (/xlb, yub, zlb/)
525 dirichlet_bc%vertices(1:3, 5) = (/xlb, yub, zub/)
526 dirichlet_bc%vertices(1:3, 6) = (/xlb, ylb, zub/)
527 dirichlet_bc%vertices(1:3, 7) = (/xub, ylb, zub/)
528 dirichlet_bc%vertices(1:3, 8) = (/xub, yub, zub/)
529
530 IF ((unit_nr .GT. 0) .AND. verbose) THEN
531 WRITE (unit_nr, '(T3,A,A)') "======== verbose ", repeat('=', 61)
532 DO i = 1, 8
533 WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, angstrom*dirichlet_bc%vertices(:, i)
534 END DO
535 END IF
536
537 CALL timestop(handle)
538
539 END SUBROUTINE aa_planar_dbc_setup
540
541! **************************************************************************************************
542!> \brief Creates an axis-aligned cuboidal Dirichlet region.
543!> \param cell_xtnts extents of the simulation cell
544!> \param x_xtnt the x extent of the cuboidal region
545!> \param y_xtnt the y extent of the cuboidal region
546!> \param z_xtnt the z extent of the cuboidal region
547!> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
548!> \param v_D potential imposed at the Dirichlet part of the boundary (gate voltage)
549!> \param osc_frac ...
550!> \param frequency ...
551!> \param phase ...
552!> \param dbc_id unique ID for the Dirichlet region
553!> \param verbose whether or not to print out the coordinates of the vertices
554!> \param dirichlet_bc the dirichlet_bc object to be created
555!> \par History
556!> 12.2014 created [Hossein Bani-Hashemian]
557!> \author Mohammad Hossein Bani-Hashemian
558! **************************************************************************************************
559 SUBROUTINE aa_cuboidal_dbc_setup(cell_xtnts, x_xtnt, y_xtnt, z_xtnt, &
560 sigma, v_D, osc_frac, frequency, phase, dbc_id, verbose, dirichlet_bc)
561
562 REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
563 REAL(dp), DIMENSION(2), INTENT(IN) :: x_xtnt, y_xtnt, z_xtnt
564 REAL(dp), INTENT(IN) :: sigma, v_d, osc_frac, frequency, phase
565 INTEGER, INTENT(IN) :: dbc_id
566 LOGICAL, INTENT(IN) :: verbose
567 TYPE(dirichlet_bc_type), INTENT(INOUT), POINTER :: dirichlet_bc
568
569 CHARACTER(LEN=*), PARAMETER :: routinen = 'aa_cuboidal_dbc_setup'
570 LOGICAL, DIMENSION(6), PARAMETER :: test_forb_xtnts = .true.
571
572 INTEGER :: handle, i, n_forb_xtnts, unit_nr
573 LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
574 forb_xtnt4, forb_xtnt5, forb_xtnt6
575 TYPE(cp_logger_type), POINTER :: logger
576
577 CALL timeset(routinen, handle)
578
579 logger => cp_get_default_logger()
580 IF (logger%para_env%is_source()) THEN
581 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
582 ELSE
583 unit_nr = -1
584 END IF
585
586 forb_xtnt1 = x_xtnt(1) .LT. cell_xtnts(1, 1)
587 forb_xtnt2 = x_xtnt(2) .GT. cell_xtnts(2, 1)
588 forb_xtnt3 = y_xtnt(1) .LT. cell_xtnts(1, 2)
589 forb_xtnt4 = y_xtnt(2) .GT. cell_xtnts(2, 2)
590 forb_xtnt5 = z_xtnt(1) .LT. cell_xtnts(1, 3)
591 forb_xtnt6 = z_xtnt(2) .GT. cell_xtnts(2, 3)
592 n_forb_xtnts = count((/forb_xtnt1, forb_xtnt2, forb_xtnt3, &
593 forb_xtnt4, forb_xtnt5, forb_xtnt6/) .EQV. test_forb_xtnts)
594 IF (n_forb_xtnts .GT. 0) THEN
595 CALL cp_abort(__location__, &
596 "The given extents for the Dirichlet region are outside the range of "// &
597 "the simulation cell.")
598 END IF
599
600 dirichlet_bc%v_D = v_d
601 dirichlet_bc%osc_frac = osc_frac
602 dirichlet_bc%frequency = frequency
603 dirichlet_bc%phase = phase
604 dirichlet_bc%dbc_id = dbc_id
605 dirichlet_bc%dbc_geom = aa_cuboidal
606 dirichlet_bc%smoothing_width = sigma
607 dirichlet_bc%n_tiles = 1
608
609 dirichlet_bc%vertices(1:3, 1) = (/x_xtnt(1), y_xtnt(1), z_xtnt(1)/)
610 dirichlet_bc%vertices(1:3, 2) = (/x_xtnt(2), y_xtnt(1), z_xtnt(1)/)
611 dirichlet_bc%vertices(1:3, 3) = (/x_xtnt(2), y_xtnt(2), z_xtnt(1)/)
612 dirichlet_bc%vertices(1:3, 4) = (/x_xtnt(1), y_xtnt(2), z_xtnt(1)/)
613 dirichlet_bc%vertices(1:3, 5) = (/x_xtnt(1), y_xtnt(2), z_xtnt(2)/)
614 dirichlet_bc%vertices(1:3, 6) = (/x_xtnt(1), y_xtnt(1), z_xtnt(2)/)
615 dirichlet_bc%vertices(1:3, 7) = (/x_xtnt(2), y_xtnt(1), z_xtnt(2)/)
616 dirichlet_bc%vertices(1:3, 8) = (/x_xtnt(2), y_xtnt(2), z_xtnt(2)/)
617
618 IF ((unit_nr .GT. 0) .AND. verbose) THEN
619 WRITE (unit_nr, '(T3,A,A)') "======== verbose ", repeat('=', 61)
620 DO i = 1, 8
621 WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, angstrom*dirichlet_bc%vertices(:, i)
622 END DO
623 END IF
624
625 CALL timestop(handle)
626
627 END SUBROUTINE aa_cuboidal_dbc_setup
628
629! **************************************************************************************************
630!> \brief Creates an arbitrary (rectangular-shaped) planar Dirichlet region given
631!> the coordinates of its vertices.
632!> \param cell_xtnts extents of the simulation cell
633!> \param A coordinates of the vertex A
634!> \param B coordinates of the vertex B
635!> \param C coordinates of the vertex C
636!> \param thickness the thickness of the region
637!> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
638!> \param v_D potential imposed at the Dirichlet part of the boundary (gate voltage)
639!> \param osc_frac ...
640!> \param frequency ...
641!> \param phase ...
642!> \param dbc_id unique ID for the Dirichlet region
643!> \param verbose whether or not to print out the coordinates of the vertices
644!> \param dirichlet_bc the dirichlet_bc object to be created
645!> \par History
646!> 08.2014 created [Hossein Bani-Hashemian]
647!> \author Mohammad Hossein Bani-Hashemian
648! **************************************************************************************************
649 SUBROUTINE arbitrary_planar_dbc_setup(cell_xtnts, A, B, C, thickness, &
650 sigma, v_D, osc_frac, frequency, phase, dbc_id, verbose, dirichlet_bc)
651
652 REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
653 REAL(dp), DIMENSION(3), INTENT(IN) :: a, b, c
654 REAL(dp), INTENT(IN) :: thickness, sigma, v_d, osc_frac, &
655 frequency, phase
656 INTEGER, INTENT(IN) :: dbc_id
657 LOGICAL, INTENT(IN) :: verbose
658 TYPE(dirichlet_bc_type), INTENT(INOUT), POINTER :: dirichlet_bc
659
660 CHARACTER(LEN=*), PARAMETER :: routinen = 'arbitrary_planar_dbc_setup'
661
662 INTEGER :: handle, i, unit_nr
663 LOGICAL :: a_is_inside, are_coplanar, &
664 are_orthogonal, b_is_inside, &
665 c_is_inside, d_is_inside, is_rectangle
666 REAL(dp) :: dist1, dist2, dist3, dist4
667 REAL(dp), DIMENSION(3) :: ab, ac, ad, bc, cm, d, normal_vector, &
668 unit_normal
669 TYPE(cp_logger_type), POINTER :: logger
670
671 CALL timeset(routinen, handle)
672
673 logger => cp_get_default_logger()
674 IF (logger%para_env%is_source()) THEN
675 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
676 ELSE
677 unit_nr = -1
678 END IF
679
680 d = a + (c - b)
681 ab = b - a; ac = c - a; ad = d - a; bc = c - b
682 are_orthogonal = abs(dot_product(ab, bc)/(sqrt(sum(ab**2))*sqrt(sum(bc**2)))) .LE. small_value
683 IF (.NOT. are_orthogonal) THEN
684 CALL cp_abort(__location__, "The given vertices for defining a "// &
685 "planar Dirichlet region do not form orthogonal edges.")
686 END IF
687
688 normal_vector = vector_product(ab, ac)
689 unit_normal = normal_vector/sqrt(sum(normal_vector**2))
690
691 a_is_inside = (a(1) .GT. cell_xtnts(1, 1) .AND. a(1) .LT. cell_xtnts(2, 1)) .AND. &
692 (a(2) .GT. cell_xtnts(1, 2) .AND. a(2) .LT. cell_xtnts(2, 2)) .AND. &
693 (a(3) .GT. cell_xtnts(1, 3) .AND. a(3) .LT. cell_xtnts(2, 3))
694 b_is_inside = (b(1) .GT. cell_xtnts(1, 1) .AND. b(1) .LT. cell_xtnts(2, 1)) .AND. &
695 (b(2) .GT. cell_xtnts(1, 2) .AND. b(2) .LT. cell_xtnts(2, 2)) .AND. &
696 (b(3) .GT. cell_xtnts(1, 3) .AND. b(3) .LT. cell_xtnts(2, 3))
697 c_is_inside = (c(1) .GT. cell_xtnts(1, 1) .AND. c(1) .LT. cell_xtnts(2, 1)) .AND. &
698 (c(2) .GT. cell_xtnts(1, 2) .AND. c(2) .LT. cell_xtnts(2, 2)) .AND. &
699 (c(3) .GT. cell_xtnts(1, 3) .AND. c(3) .LT. cell_xtnts(2, 3))
700 d_is_inside = (d(1) .GT. cell_xtnts(1, 1) .AND. d(1) .LT. cell_xtnts(2, 1)) .AND. &
701 (d(2) .GT. cell_xtnts(1, 2) .AND. d(2) .LT. cell_xtnts(2, 2)) .AND. &
702 (d(3) .GT. cell_xtnts(1, 3) .AND. d(3) .LT. cell_xtnts(2, 3))
703 IF (.NOT. (a_is_inside .AND. b_is_inside .AND. c_is_inside .AND. d_is_inside)) THEN
704 CALL cp_abort(__location__, &
705 "At least one of the given vertices for defining a planar Dirichlet "// &
706 "region is outside the simulation box.")
707 END IF
708
709 are_coplanar = abs(dot_product(vector_product(ab, ac), ad)) .LE. small_value
710 IF (.NOT. are_coplanar) THEN
711 cpabort("The given vertices for defining a planar Dirichlet region are not coplanar.")
712 END IF
713
714 cm = (a + b + c + d)/4.0_dp
715 dist1 = sum((a - cm)**2)
716 dist2 = sum((b - cm)**2)
717 dist3 = sum((c - cm)**2)
718 dist4 = sum((d - cm)**2)
719 is_rectangle = (abs(dist1 - dist2) .LE. small_value) .AND. &
720 (abs(dist1 - dist3) .LE. small_value) .AND. &
721 (abs(dist1 - dist4) .LE. small_value)
722 IF (.NOT. is_rectangle) THEN
723 CALL cp_abort(__location__, "The given vertices for defining "// &
724 "a planar Dirichlet region do not form a rectangle.")
725 END IF
726
727 dirichlet_bc%v_D = v_d
728 dirichlet_bc%osc_frac = osc_frac
729 dirichlet_bc%frequency = frequency
730 dirichlet_bc%phase = phase
731 dirichlet_bc%dbc_id = dbc_id
732 dirichlet_bc%dbc_geom = planar
733 dirichlet_bc%smoothing_width = sigma
734 dirichlet_bc%n_tiles = 1
735
736 dirichlet_bc%vertices(1:3, 1) = a - 0.5_dp*thickness*unit_normal
737 dirichlet_bc%vertices(1:3, 2) = b - 0.5_dp*thickness*unit_normal
738 dirichlet_bc%vertices(1:3, 3) = c - 0.5_dp*thickness*unit_normal
739 dirichlet_bc%vertices(1:3, 4) = d - 0.5_dp*thickness*unit_normal
740 dirichlet_bc%vertices(1:3, 5) = d + 0.5_dp*thickness*unit_normal
741 dirichlet_bc%vertices(1:3, 6) = a + 0.5_dp*thickness*unit_normal
742 dirichlet_bc%vertices(1:3, 7) = b + 0.5_dp*thickness*unit_normal
743 dirichlet_bc%vertices(1:3, 8) = c + 0.5_dp*thickness*unit_normal
744
745 IF ((unit_nr .GT. 0) .AND. verbose) THEN
746 WRITE (unit_nr, '(T3,A,A)') "======== verbose ", repeat('=', 61)
747 DO i = 1, 8
748 WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, angstrom*dirichlet_bc%vertices(:, i)
749 END DO
750 END IF
751
752 CALL timestop(handle)
753
754 END SUBROUTINE arbitrary_planar_dbc_setup
755
756! **************************************************************************************************
757!> \brief Creates an x-axis-aligned cylindrical Dirichlet region. The cylindrical
758!> shape is approximated by an n-gonal (circumscribed or inscribed) uniform prism-
759!> like object. The circular base is approximated by a polygon. A second polygon is
760!> then created by rotating the first one by delta_alpha radians. Let's name the
761!> vertices of the first and the second polygon as v1, v2, v3, ... vn-1, vn and
762!> w1, w2, w3, ... wn-1, wn, respectively. Then w1v2, w2v3, ... wn-1vn, wnv1 form
763!> the edges of the cross-section of the approximating object. This way the dbcs
764!> do not share any edge.
765!> \param pw_pool pool of plane wave grid
766!> \param cell_xtnts extents of the simulation cell
767!> \param x_locl x grid vector of the simulation box local to this process
768!> \param y_locl y grid vector of the simulation box local to this process
769!> \param z_locl z grid vector of the simulation box local to this process
770!> \param apx_type the type of the n-gonal prism approximating the cylinder
771!> \param parallel_axis the coordinate axis that the cylindrical region extends along
772!> \param xtnt the x extent of the cylindrical region along its central axis
773!> \param base_center the y and z coordinates of the cylinder's base
774!> \param base_radius the radius of the cylinder's base
775!> \param thickness the thickness of the region
776!> \param delta_alpha ...
777!> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
778!> \param v_D potential imposed at the Dirichlet part of the boundary (gate voltage)
779!> \param osc_frac ...
780!> \param frequency ...
781!> \param phase ...
782!> \param n_prtn number of partitions along the edges of the faces of the prism
783!> \param is_periodic whether or not the tile is periodic with a periodicity of one unit cell
784!> \param verbose whether or not to print out the coordinates of the vertices
785!> \param dbcs the x-axis-aligned cylindrical gate region to be created
786!> \par History
787!> 08.2014 created [Hossein Bani-Hashemian]
788!> \author Mohammad Hossein Bani-Hashemian
789! **************************************************************************************************
790 SUBROUTINE aa_cylindrical_dbc_setup(pw_pool, cell_xtnts, x_locl, y_locl, z_locl, apx_type, parallel_axis, &
791 xtnt, base_center, base_radius, thickness, delta_alpha, &
792 sigma, v_D, osc_frac, frequency, phase, n_prtn, is_periodic, verbose, dbcs)
793
794 TYPE(pw_pool_type), POINTER :: pw_pool
795 REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
796 REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_locl, y_locl, z_locl
797 INTEGER, INTENT(IN), OPTIONAL :: apx_type
798 INTEGER, INTENT(IN) :: parallel_axis
799 REAL(dp), DIMENSION(2), INTENT(IN) :: xtnt, base_center
800 REAL(dp), INTENT(IN) :: base_radius, thickness, delta_alpha, &
801 sigma, v_d, osc_frac, frequency, phase
802 INTEGER, DIMENSION(3), INTENT(IN) :: n_prtn
803 LOGICAL, INTENT(IN) :: is_periodic, verbose
804 TYPE(dirichlet_bc_p_type), DIMENSION(:), &
805 INTENT(INOUT) :: dbcs
806
807 CHARACTER(LEN=*), PARAMETER :: routinen = 'aa_cylindrical_dbc_setup'
808
809 INTEGER :: handle, i, intern_apx_type, j, n_dbcs, &
810 unit_nr
811 LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
812 forb_xtnt4
813 REAL(dp) :: cntrlaxis_lb, cntrlaxis_ub, h, lx, ly, &
814 lz, theta
815 REAL(dp), DIMENSION(3) :: a, b, c, d, normal_vec, unit_normal
816 TYPE(cp_logger_type), POINTER :: logger
817 REAL(dp), DIMENSION(SIZE(dbcs)+1) :: alpha, alpha_rotd, xlb, xub, ylb, yub, &
818 zlb, zub
819
820 CALL timeset(routinen, handle)
821
822 logger => cp_get_default_logger()
823 IF (logger%para_env%is_source()) THEN
824 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
825 ELSE
826 unit_nr = -1
827 END IF
828
829 lx = cell_xtnts(2, 1) - cell_xtnts(1, 1)
830 ly = cell_xtnts(2, 2) - cell_xtnts(1, 2)
831 lz = cell_xtnts(2, 3) - cell_xtnts(1, 3)
832
833 SELECT CASE (parallel_axis)
834 CASE (x_axis)
835 IF (xtnt(1) .LT. cell_xtnts(1, 1) .OR. xtnt(2) .GT. cell_xtnts(2, 1)) THEN
836 CALL cp_abort(__location__, &
837 "The length of the cylindrical Dirichlet region is larger than the "// &
838 "x range of the simulation cell.")
839 END IF
840 forb_xtnt1 = base_center(1) - base_radius .LT. cell_xtnts(1, 2)
841 forb_xtnt2 = base_center(1) + base_radius .GT. cell_xtnts(2, 2)
842 forb_xtnt3 = base_center(2) - base_radius .LT. cell_xtnts(1, 3)
843 forb_xtnt4 = base_center(2) + base_radius .GT. cell_xtnts(2, 3)
844 IF (forb_xtnt1 .OR. forb_xtnt2 .OR. forb_xtnt3 .OR. forb_xtnt4) THEN
845 cpabort("The cylinder does not fit entirely inside the simulation cell.")
846 END IF
847 CASE (y_axis)
848 IF (xtnt(1) .LT. cell_xtnts(1, 2) .OR. xtnt(2) .GT. cell_xtnts(2, 2)) THEN
849 CALL cp_abort(__location__, &
850 "The length of the cylindrical Dirichlet region is larger than the "// &
851 "y range of the simulation cell.")
852 END IF
853 forb_xtnt1 = base_center(1) - base_radius .LT. cell_xtnts(1, 1)
854 forb_xtnt2 = base_center(1) + base_radius .GT. cell_xtnts(2, 1)
855 forb_xtnt3 = base_center(2) - base_radius .LT. cell_xtnts(1, 3)
856 forb_xtnt4 = base_center(2) + base_radius .GT. cell_xtnts(2, 3)
857 IF (forb_xtnt1 .OR. forb_xtnt2 .OR. forb_xtnt3 .OR. forb_xtnt4) THEN
858 cpabort("The cylinder does not fit entirely inside the simulation cell.")
859 END IF
860 CASE (z_axis)
861 IF (xtnt(1) .LT. cell_xtnts(1, 3) .OR. xtnt(2) .GT. cell_xtnts(2, 3)) THEN
862 CALL cp_abort(__location__, &
863 "The length of the cylindrical Dirichlet region is larger than the "// &
864 "z range of the simulation cell.")
865 END IF
866 forb_xtnt1 = base_center(1) - base_radius .LT. cell_xtnts(1, 1)
867 forb_xtnt2 = base_center(1) + base_radius .GT. cell_xtnts(2, 1)
868 forb_xtnt3 = base_center(2) - base_radius .LT. cell_xtnts(1, 2)
869 forb_xtnt4 = base_center(2) + base_radius .GT. cell_xtnts(2, 2)
870 IF (forb_xtnt1 .OR. forb_xtnt2 .OR. forb_xtnt3 .OR. forb_xtnt4) THEN
871 cpabort("The cylinder does not fit entirely inside the simulation cell.")
872 END IF
873 END SELECT
874
875 intern_apx_type = circumscribed
876 IF (PRESENT(apx_type)) intern_apx_type = apx_type
877
878 n_dbcs = SIZE(dbcs)
879
880 cntrlaxis_lb = xtnt(1)
881 cntrlaxis_ub = xtnt(2)
882
883 theta = twopi/n_dbcs
884
885 IF (intern_apx_type .EQ. inscribed) THEN
886 h = base_radius ! inscribed uniform prism
887 ELSE IF (intern_apx_type .EQ. circumscribed) THEN
888 h = base_radius/cos(0.5*theta) ! circumscribed uniform prism
889 ELSE
890 cpabort("Unknown approximation type for cylinder.")
891 END IF
892 SELECT CASE (parallel_axis)
893 CASE (x_axis)
894 IF (h .GT. minval((/ly, lz/))) THEN
895 cpabort("Reduce the base radius!")
896 END IF
897 CASE (y_axis)
898 IF (h .GT. minval((/lx, lz/))) THEN
899 cpabort("Reduce the base radius!")
900 END IF
901 CASE (z_axis)
902 IF (h .GT. minval((/lx, ly/))) THEN
903 cpabort("Reduce the base radius!")
904 END IF
905 END SELECT
906
907 alpha = linspace(0.0_dp, 2*pi, n_dbcs + 1)
908 alpha_rotd = alpha + delta_alpha;
909 SELECT CASE (parallel_axis)
910 CASE (x_axis)
911 DO j = 1, n_dbcs
912 ylb(j) = base_center(1) + h*sin(alpha(j))
913 zlb(j) = base_center(2) + h*cos(alpha(j))
914 yub(j) = base_center(1) + h*sin(alpha_rotd(j))
915 zub(j) = base_center(2) + h*cos(alpha_rotd(j))
916 END DO
917 ylb(n_dbcs + 1) = ylb(1)
918 yub(n_dbcs + 1) = yub(1)
919 zlb(n_dbcs + 1) = zlb(1)
920 zub(n_dbcs + 1) = zub(1)
921 DO j = 1, n_dbcs
922 ALLOCATE (dbcs(j)%dirichlet_bc)
923 dbcs(j)%dirichlet_bc%dbc_geom = cylindrical
924 dbcs(j)%dirichlet_bc%v_D = v_d
925 dbcs(j)%dirichlet_bc%osc_frac = osc_frac
926 dbcs(j)%dirichlet_bc%frequency = frequency
927 dbcs(j)%dirichlet_bc%phase = phase
928 dbcs(j)%dirichlet_bc%dbc_id = cylindrical + j
929 dbcs(j)%dirichlet_bc%smoothing_width = sigma
930
931 a = (/cntrlaxis_lb, yub(j), zub(j)/)
932 b = (/cntrlaxis_lb, ylb(j + 1), zlb(j + 1)/)
933 c = (/cntrlaxis_ub, ylb(j + 1), zlb(j + 1)/)
934 d = (/cntrlaxis_ub, yub(j), zub(j)/)
935 normal_vec = vector_product((a - c), (d - b))
936 unit_normal = normal_vec/sqrt(sum(normal_vec**2))
937
938 dbcs(j)%dirichlet_bc%vertices(1:3, 1) = a
939 dbcs(j)%dirichlet_bc%vertices(1:3, 2) = b
940 dbcs(j)%dirichlet_bc%vertices(1:3, 3) = c
941 dbcs(j)%dirichlet_bc%vertices(1:3, 4) = d
942 dbcs(j)%dirichlet_bc%vertices(1:3, 5) = d + thickness*unit_normal
943 dbcs(j)%dirichlet_bc%vertices(1:3, 6) = a + thickness*unit_normal
944 dbcs(j)%dirichlet_bc%vertices(1:3, 7) = b + thickness*unit_normal
945 dbcs(j)%dirichlet_bc%vertices(1:3, 8) = c + thickness*unit_normal
946
947 dbcs(j)%dirichlet_bc%n_tiles = 1
948
949 IF ((unit_nr .GT. 0) .AND. verbose) THEN
950 WRITE (unit_nr, '(T3,A,A)') "======== verbose ", repeat('=', 61)
951 WRITE (unit_nr, '(T7,A,I5,T20,A,I5,A)') "edge", j, "of the", n_dbcs, &
952 "-gonal prism approximating the cylinder"
953 DO i = 1, 8
954 WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, &
955 angstrom*dbcs(j)%dirichlet_bc%vertices(:, i)
956 END DO
957 END IF
958
959 CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
960 x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
961 END DO
962 CASE (y_axis)
963 DO j = 1, n_dbcs
964 xlb(j) = base_center(1) + h*sin(alpha(j))
965 zlb(j) = base_center(2) + h*cos(alpha(j))
966 xub(j) = base_center(1) + h*sin(alpha_rotd(j))
967 zub(j) = base_center(2) + h*cos(alpha_rotd(j))
968 END DO
969 xlb(n_dbcs + 1) = xlb(1)
970 xub(n_dbcs + 1) = xub(1)
971 zlb(n_dbcs + 1) = zlb(1)
972 zub(n_dbcs + 1) = zub(1)
973 DO j = 1, n_dbcs
974 ALLOCATE (dbcs(j)%dirichlet_bc)
975 dbcs(j)%dirichlet_bc%dbc_geom = cylindrical
976 dbcs(j)%dirichlet_bc%v_D = v_d
977 dbcs(j)%dirichlet_bc%osc_frac = osc_frac
978 dbcs(j)%dirichlet_bc%frequency = frequency
979 dbcs(j)%dirichlet_bc%phase = phase
980 dbcs(j)%dirichlet_bc%dbc_id = cylindrical + j
981 dbcs(j)%dirichlet_bc%smoothing_width = sigma
982
983 a = (/xub(j), cntrlaxis_lb, zub(j)/)
984 b = (/xlb(j + 1), cntrlaxis_lb, zlb(j + 1)/)
985 c = (/xlb(j + 1), cntrlaxis_ub, zlb(j + 1)/)
986 d = (/xub(j), cntrlaxis_ub, zub(j)/)
987 normal_vec = vector_product((a - c), (d - b))
988 unit_normal = normal_vec/sqrt(sum(normal_vec**2))
989
990 dbcs(j)%dirichlet_bc%vertices(1:3, 1) = a
991 dbcs(j)%dirichlet_bc%vertices(1:3, 2) = b
992 dbcs(j)%dirichlet_bc%vertices(1:3, 3) = c
993 dbcs(j)%dirichlet_bc%vertices(1:3, 4) = d
994 dbcs(j)%dirichlet_bc%vertices(1:3, 5) = d + thickness*unit_normal
995 dbcs(j)%dirichlet_bc%vertices(1:3, 6) = a + thickness*unit_normal
996 dbcs(j)%dirichlet_bc%vertices(1:3, 7) = b + thickness*unit_normal
997 dbcs(j)%dirichlet_bc%vertices(1:3, 8) = c + thickness*unit_normal
998
999 dbcs(j)%dirichlet_bc%n_tiles = 1
1000
1001 IF ((unit_nr .GT. 0) .AND. verbose) THEN
1002 WRITE (unit_nr, '(T3,A,A)') "======== verbose ", repeat('=', 61)
1003 WRITE (unit_nr, '(T7,A,I5,T20,A,I5,A)') "edge", j, "of the", n_dbcs, &
1004 "-gonal prism approximating the cylinder"
1005 DO i = 1, 8
1006 WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, &
1007 angstrom*dbcs(j)%dirichlet_bc%vertices(:, i)
1008 END DO
1009 END IF
1010
1011 CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
1012 x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
1013 END DO
1014 CASE (z_axis)
1015 DO j = 1, n_dbcs
1016 xlb(j) = base_center(1) + h*sin(alpha(j))
1017 ylb(j) = base_center(2) + h*cos(alpha(j))
1018 xub(j) = base_center(1) + h*sin(alpha_rotd(j))
1019 yub(j) = base_center(2) + h*cos(alpha_rotd(j))
1020 END DO
1021 xlb(n_dbcs + 1) = xlb(1)
1022 xub(n_dbcs + 1) = xub(1)
1023 ylb(n_dbcs + 1) = ylb(1)
1024 yub(n_dbcs + 1) = yub(1)
1025 DO j = 1, n_dbcs
1026 ALLOCATE (dbcs(j)%dirichlet_bc)
1027 dbcs(j)%dirichlet_bc%dbc_geom = cylindrical
1028 dbcs(j)%dirichlet_bc%v_D = v_d
1029 dbcs(j)%dirichlet_bc%osc_frac = osc_frac
1030 dbcs(j)%dirichlet_bc%frequency = frequency
1031 dbcs(j)%dirichlet_bc%phase = phase
1032 dbcs(j)%dirichlet_bc%dbc_id = cylindrical + j
1033 dbcs(j)%dirichlet_bc%smoothing_width = sigma
1034
1035 a = (/xub(j), yub(j), cntrlaxis_lb/)
1036 b = (/xlb(j + 1), ylb(j + 1), cntrlaxis_lb/)
1037 c = (/xlb(j + 1), ylb(j + 1), cntrlaxis_ub/)
1038 d = (/xub(j), yub(j), cntrlaxis_ub/)
1039 normal_vec = vector_product((a - c), (d - b))
1040 unit_normal = normal_vec/sqrt(sum(normal_vec**2))
1041
1042 dbcs(j)%dirichlet_bc%vertices(1:3, 1) = a
1043 dbcs(j)%dirichlet_bc%vertices(1:3, 2) = b
1044 dbcs(j)%dirichlet_bc%vertices(1:3, 3) = c
1045 dbcs(j)%dirichlet_bc%vertices(1:3, 4) = d
1046 dbcs(j)%dirichlet_bc%vertices(1:3, 5) = d + thickness*unit_normal
1047 dbcs(j)%dirichlet_bc%vertices(1:3, 6) = a + thickness*unit_normal
1048 dbcs(j)%dirichlet_bc%vertices(1:3, 7) = b + thickness*unit_normal
1049 dbcs(j)%dirichlet_bc%vertices(1:3, 8) = c + thickness*unit_normal
1050
1051 dbcs(j)%dirichlet_bc%n_tiles = 1
1052
1053 IF ((unit_nr .GT. 0) .AND. verbose) THEN
1054 WRITE (unit_nr, '(T3,A,A)') "======== verbose ", repeat('=', 61)
1055 WRITE (unit_nr, '(T7,A,I5,T20,A,I5,A)') "edge", j, "of the", n_dbcs, &
1056 "-gonal prism approximating the cylinder"
1057 DO i = 1, 8
1058 WRITE (unit_nr, '(T10,A,I1,3F10.3)') " vertex ", i, &
1059 angstrom*dbcs(j)%dirichlet_bc%vertices(:, i)
1060 END DO
1061 END IF
1062
1063 CALL dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
1064 x_locl, y_locl, z_locl, is_periodic, verbose, dbcs(j)%dirichlet_bc)
1065 END DO
1066 END SELECT
1067
1068 CALL timestop(handle)
1069
1070 END SUBROUTINE aa_cylindrical_dbc_setup
1071
1072! **************************************************************************************************
1073!> \brief Creteas an evenly-spaced 1D array between two values
1074!> \param a starting value
1075!> \param b end value
1076!> \param N number of evenly-spaced points between a and b
1077!> \return array containg N evenly-spaced points between a and b
1078!> \par History
1079!> 07.2014 created [Hossein Bani-Hashemian]
1080!> \author Mohammad Hossein Bani-Hashemian
1081! **************************************************************************************************
1082 FUNCTION linspace(a, b, N) RESULT(arr)
1083
1084 REAL(dp), INTENT(IN) :: a, b
1085 INTEGER, INTENT(IN) :: n
1086 REAL(dp), DIMENSION(N) :: arr
1087
1088 INTEGER :: i
1089 REAL(dp) :: dx
1090
1091 dx = (b - a)/(n - 1)
1092 arr = (/(a + (i - 1)*dx, i=1, n)/)
1093
1094 END FUNCTION linspace
1095
1096! **************************************************************************************************
1097!> \brief rotates and translates a cuboid
1098!> \param cuboid_vtx vertices of the cuboid
1099!> \param Rmat rotation matrix
1100!> \param Tpnt translation vector (translates the origin to the center of the cuboid)
1101!> \param cuboid_transfd_vtx vertices of the transformed cuboid
1102!> \par History
1103!> 11.2015 created [Hossein Bani-Hashemian]
1104!> \author Mohammad Hossein Bani-Hashemian
1105! **************************************************************************************************
1106 SUBROUTINE rotate_translate_cuboid(cuboid_vtx, Rmat, Tpnt, cuboid_transfd_vtx)
1107
1108 REAL(dp), DIMENSION(3, 8), INTENT(IN) :: cuboid_vtx
1109 REAL(dp), DIMENSION(3, 3), INTENT(OUT) :: rmat
1110 REAL(dp), DIMENSION(3), INTENT(OUT) :: tpnt
1111 REAL(dp), DIMENSION(3, 8), INTENT(OUT) :: cuboid_transfd_vtx
1112
1113 CHARACTER(LEN=*), PARAMETER :: routinen = 'rotate_translate_cuboid'
1114
1115 INTEGER :: handle
1116 REAL(dp), DIMENSION(3) :: a, a2, aa2, ab, ad, b, b2, c, c2, d, d2, &
1117 e1_locl, e2_locl, e3_locl
1118
1119 CALL timeset(routinen, handle)
1120
1121 rmat = 0.0_dp
1122
1123 a = cuboid_vtx(1:3, 1); a2 = cuboid_vtx(1:3, 6)
1124 b = cuboid_vtx(1:3, 2); b2 = cuboid_vtx(1:3, 7)
1125 c = cuboid_vtx(1:3, 3); c2 = cuboid_vtx(1:3, 8)
1126 d = cuboid_vtx(1:3, 4); d2 = cuboid_vtx(1:3, 5)
1127
1128 tpnt = (a + c2)/2.0_dp
1129
1130 ab = b - a
1131 ad = d - a
1132 aa2 = a2 - a
1133
1134! unit vectors generating the local coordinate system
1135 e1_locl = ab/sqrt(sum(ab**2))
1136 e2_locl = ad/sqrt(sum(ad**2))
1137 e3_locl = aa2/sqrt(sum(aa2**2))
1138! rotation matrix
1139 rmat(1:3, 1) = e1_locl
1140 rmat(1:3, 2) = e2_locl
1141 rmat(1:3, 3) = e3_locl
1142
1143 cuboid_transfd_vtx(1:3, 1) = rotate_translate_vector(rmat, tpnt, bwrot, a)
1144 cuboid_transfd_vtx(1:3, 2) = rotate_translate_vector(rmat, tpnt, bwrot, b)
1145 cuboid_transfd_vtx(1:3, 3) = rotate_translate_vector(rmat, tpnt, bwrot, c)
1146 cuboid_transfd_vtx(1:3, 4) = rotate_translate_vector(rmat, tpnt, bwrot, d)
1147 cuboid_transfd_vtx(1:3, 5) = rotate_translate_vector(rmat, tpnt, bwrot, d2)
1148 cuboid_transfd_vtx(1:3, 6) = rotate_translate_vector(rmat, tpnt, bwrot, a2)
1149 cuboid_transfd_vtx(1:3, 7) = rotate_translate_vector(rmat, tpnt, bwrot, b2)
1150 cuboid_transfd_vtx(1:3, 8) = rotate_translate_vector(rmat, tpnt, bwrot, c2)
1151
1152 CALL timestop(handle)
1153
1154 END SUBROUTINE rotate_translate_cuboid
1155
1156! **************************************************************************************************
1157!> \brief transforms a position vector according to a given rotation matrix and a
1158!> translation vector
1159!> \param Rmat rotation matrix
1160!> \param Tp image of the origin under the translation
1161!> \param direction forward or backward transformation
1162!> \param vec input vector
1163!> \return transformed vector
1164!> \par History
1165!> 11.2015 created [Hossein Bani-Hashemian]
1166!> \author Mohammad Hossein Bani-Hashemian
1167! **************************************************************************************************
1168 FUNCTION rotate_translate_vector(Rmat, Tp, direction, vec) RESULT(vec_transfd)
1169 REAL(dp), DIMENSION(3, 3), INTENT(IN) :: rmat
1170 REAL(dp), DIMENSION(3), INTENT(IN) :: tp
1171 INTEGER, INTENT(IN) :: direction
1172 REAL(dp), DIMENSION(3), INTENT(IN) :: vec
1173 REAL(dp), DIMENSION(3) :: vec_transfd
1174
1175 REAL(dp), DIMENSION(3) :: tpoint, vec_tmp
1176 REAL(dp), DIMENSION(3, 3) :: rmat_inv
1177
1178 tpoint = 0.0_dp
1179
1180 IF (direction .EQ. fwrot) THEN
1181 vec_tmp = matmul(rmat, vec)
1182 vec_transfd = vec_tmp + tp
1183 ELSEIF (direction .EQ. bwrot) THEN
1184 rmat_inv = inv_3x3(rmat)
1185 tpoint(1) = rmat_inv(1, 1)*tp(1) + rmat_inv(1, 2)*tp(2) + rmat_inv(1, 3)*tp(3)
1186 tpoint(2) = rmat_inv(2, 1)*tp(1) + rmat_inv(2, 2)*tp(2) + rmat_inv(2, 3)*tp(3)
1187 tpoint(3) = rmat_inv(3, 1)*tp(1) + rmat_inv(3, 2)*tp(2) + rmat_inv(3, 3)*tp(3)
1188 vec_tmp = matmul(rmat_inv, vec)
1189 vec_transfd = vec_tmp - tpoint
1190 END IF
1191
1192 END FUNCTION rotate_translate_vector
1193
1194! **************************************************************************************************
1195!> \brief Partitions an axis-aligned cuboid into tiles.
1196!> \param dbc_vertices vertices of the axis-aligned Dirichlet region to be partitioned
1197!> \param n_prtn vector of size 3 specifying the number of times that the x, y and
1198!> z interval (defining the region) should be partitioned into
1199!> \param tiles the obtained tiles
1200!> \par History
1201!> 10.2015 created [Hossein Bani-Hashemian]
1202!> \author Mohammad Hossein Bani-Hashemian
1203! **************************************************************************************************
1204 SUBROUTINE aa_dbc_partition(dbc_vertices, n_prtn, tiles)
1205
1206 REAL(dp), DIMENSION(3, 8), INTENT(IN) :: dbc_vertices
1207 INTEGER, DIMENSION(3), INTENT(IN) :: n_prtn
1208 TYPE(tile_p_type), DIMENSION(:), INTENT(OUT), &
1209 POINTER :: tiles
1210
1211 CHARACTER(LEN=*), PARAMETER :: routinen = 'aa_dbc_partition'
1212
1213 INTEGER :: handle, i, ii, jj, k, kk
1214 REAL(dp) :: step
1215 REAL(dp), ALLOCATABLE, DIMENSION(:) :: xprtn_lb, xprtn_ub, yprtn_lb, yprtn_ub, &
1216 zprtn_lb, zprtn_ub
1217 REAL(dp), DIMENSION(2) :: x_xtnt, y_xtnt, z_xtnt
1218
1219 CALL timeset(routinen, handle)
1220
1221! vertices are labeled as:
1222! 6-------5 A2------D2
1223! /| /| /| /|
1224! 7-|-----8 | B2-|----C2 |
1225! | 1-----|-4 or | A-----|-D
1226! |/ |/ |/ |/
1227! 2-------3 B-------C
1228
1229 ALLOCATE (xprtn_lb(n_prtn(1)), xprtn_ub(n_prtn(1)))
1230 ALLOCATE (yprtn_lb(n_prtn(2)), yprtn_ub(n_prtn(2)))
1231 ALLOCATE (zprtn_lb(n_prtn(3)), zprtn_ub(n_prtn(3)))
1232
1233! find the extents of the cuboid in all three directions
1234 x_xtnt(1) = minval(dbc_vertices(1, :)); x_xtnt(2) = maxval(dbc_vertices(1, :))
1235 y_xtnt(1) = minval(dbc_vertices(2, :)); y_xtnt(2) = maxval(dbc_vertices(2, :))
1236 z_xtnt(1) = minval(dbc_vertices(3, :)); z_xtnt(2) = maxval(dbc_vertices(3, :))
1237
1238! divide the x, y and z extents into n_prtn partitions
1239 step = (x_xtnt(2) - x_xtnt(1))/real(n_prtn(1), kind=dp)
1240 xprtn_lb(:) = x_xtnt(1) + (/(i, i=0, n_prtn(1) - 1)/)*step ! lower bounds
1241 xprtn_ub(:) = x_xtnt(1) + (/(i, i=1, n_prtn(1))/)*step ! upper bounds
1242
1243 step = (y_xtnt(2) - y_xtnt(1))/real(n_prtn(2), kind=dp)
1244 yprtn_lb(:) = y_xtnt(1) + (/(i, i=0, n_prtn(2) - 1)/)*step
1245 yprtn_ub(:) = y_xtnt(1) + (/(i, i=1, n_prtn(2))/)*step
1246
1247 step = (z_xtnt(2) - z_xtnt(1))/real(n_prtn(3), kind=dp)
1248 zprtn_lb(:) = z_xtnt(1) + (/(i, i=0, n_prtn(3) - 1)/)*step
1249 zprtn_ub(:) = z_xtnt(1) + (/(i, i=1, n_prtn(3))/)*step
1250
1251 ALLOCATE (tiles(n_prtn(1)*n_prtn(2)*n_prtn(3)))
1252 k = 1
1253 DO kk = 1, n_prtn(3)
1254 DO jj = 1, n_prtn(2)
1255 DO ii = 1, n_prtn(1)
1256 ALLOCATE (tiles(k)%tile)
1257 tiles(k)%tile%tile_id = k
1258
1259 tiles(k)%tile%vertices(1:3, 1) = (/xprtn_lb(ii), yprtn_lb(jj), zprtn_lb(kk)/)
1260 tiles(k)%tile%vertices(1:3, 2) = (/xprtn_ub(ii), yprtn_lb(jj), zprtn_lb(kk)/)
1261 tiles(k)%tile%vertices(1:3, 3) = (/xprtn_ub(ii), yprtn_ub(jj), zprtn_lb(kk)/)
1262 tiles(k)%tile%vertices(1:3, 4) = (/xprtn_lb(ii), yprtn_ub(jj), zprtn_lb(kk)/)
1263 tiles(k)%tile%vertices(1:3, 5) = (/xprtn_lb(ii), yprtn_ub(jj), zprtn_ub(kk)/)
1264 tiles(k)%tile%vertices(1:3, 6) = (/xprtn_lb(ii), yprtn_lb(jj), zprtn_ub(kk)/)
1265 tiles(k)%tile%vertices(1:3, 7) = (/xprtn_ub(ii), yprtn_lb(jj), zprtn_ub(kk)/)
1266 tiles(k)%tile%vertices(1:3, 8) = (/xprtn_ub(ii), yprtn_ub(jj), zprtn_ub(kk)/)
1267
1268 tiles(k)%tile%volume = 0.0_dp
1269 k = k + 1
1270 END DO
1271 END DO
1272 END DO
1273
1274 CALL timestop(handle)
1275
1276 END SUBROUTINE aa_dbc_partition
1277
1278! **************************************************************************************************
1279!> \brief Partitions an arbitrary cuboid into tiles.
1280!> \param dbc_vertices vertices of the axis-aligned Dirichlet region to be partitioned
1281!> \param n_prtn number of partitions along the edges
1282!> \param tiles the obtained tiles
1283!> \par History
1284!> 10.2015 created [Hossein Bani-Hashemian]
1285!> \author Mohammad Hossein Bani-Hashemian
1286! **************************************************************************************************
1287 SUBROUTINE arbitrary_dbc_partition(dbc_vertices, n_prtn, tiles)
1288
1289 REAL(dp), DIMENSION(3, 8), INTENT(IN) :: dbc_vertices
1290 INTEGER, DIMENSION(2), INTENT(IN) :: n_prtn
1291 TYPE(tile_p_type), DIMENSION(:), INTENT(OUT), &
1292 POINTER :: tiles
1293
1294 CHARACTER(LEN=*), PARAMETER :: routinen = 'arbitrary_dbc_partition'
1295
1296 INTEGER :: handle, i, k, l, np1, np2
1297 REAL(dp) :: ablength, adlength, step, thickness
1298 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: end_points, x
1299 REAL(dp), DIMENSION(3) :: a, ab, ac, ad, b, c, d, d2, &
1300 normal_vector, point1, point2, &
1301 unit_normal
1302
1303 CALL timeset(routinen, handle)
1304
1305 a = dbc_vertices(:, 1)
1306 b = dbc_vertices(:, 2)
1307 c = dbc_vertices(:, 3)
1308 d = dbc_vertices(:, 4)
1309 d2 = dbc_vertices(:, 5)
1310
1311 ab = b - a
1312 ac = c - a
1313 ad = d - a
1314 normal_vector = vector_product(ab, ac)
1315 unit_normal = normal_vector/sqrt(sum(normal_vector**2))
1316 thickness = sqrt(sum((d - d2)**2))
1317
1318! the larger n_prtn is assigned to the longer edge
1319 ablength = sqrt(sum(ab**2))
1320 adlength = sqrt(sum(ad**2))
1321 IF (adlength .GE. ablength) THEN
1322 np1 = max(n_prtn(1) + 1, n_prtn(2) + 1)
1323 np2 = min(n_prtn(1) + 1, n_prtn(2) + 1)
1324 ELSE
1325 np1 = min(n_prtn(1) + 1, n_prtn(2) + 1)
1326 np2 = max(n_prtn(1) + 1, n_prtn(2) + 1)
1327 END IF
1328
1329 ALLOCATE (x(np1, np2, 3))
1330 ALLOCATE (end_points(2, np1, 3))
1331
1332! partition AD and BC
1333 DO l = 1, np1
1334 step = real((l - 1), kind=dp)/real((np1 - 1), kind=dp)
1335 end_points(1, l, :) = a*(1.0_dp - step) + d*step
1336 end_points(2, l, :) = b*(1.0_dp - step) + c*step
1337 END DO
1338
1339! partition in the second direction along the line segments with endpoints from
1340! the previous partitioning step
1341 DO l = 1, np1
1342 point1(:) = end_points(1, l, :)
1343 point2(:) = end_points(2, l, :)
1344 DO i = 1, np2
1345 step = real((i - 1), kind=dp)/real((np2 - 1), kind=dp)
1346 x(l, i, :) = point1*(1.0_dp - step) + point2*step
1347 END DO
1348 END DO
1349
1350 ALLOCATE (tiles((np1 - 1)*(np2 - 1)))
1351 k = 1
1352 DO l = 1, np1 - 1
1353 DO i = 1, np2 - 1
1354 ALLOCATE (tiles(k)%tile)
1355 tiles(k)%tile%tile_id = k
1356
1357 tiles(k)%tile%vertices(1:3, 1) = x(l, i, :)
1358 tiles(k)%tile%vertices(1:3, 2) = x(l, i + 1, :)
1359 tiles(k)%tile%vertices(1:3, 3) = x(l + 1, i + 1, :)
1360 tiles(k)%tile%vertices(1:3, 4) = x(l + 1, i, :)
1361 tiles(k)%tile%vertices(1:3, 5) = x(l + 1, i, :) + thickness*unit_normal
1362 tiles(k)%tile%vertices(1:3, 6) = x(l, i, :) + thickness*unit_normal
1363 tiles(k)%tile%vertices(1:3, 7) = x(l, i + 1, :) + thickness*unit_normal
1364 tiles(k)%tile%vertices(1:3, 8) = x(l + 1, i + 1, :) + thickness*unit_normal
1365
1366 tiles(k)%tile%volume = 0.0_dp
1367 k = k + 1
1368 END DO
1369 END DO
1370
1371 CALL timestop(handle)
1372
1373 END SUBROUTINE arbitrary_dbc_partition
1374
1375! **************************************************************************************************
1376!> \brief computes the pw data corresponding to an axis-aligned tile
1377!> \param cell_xtnts extents of the simulation cell
1378!> \param x_locl x grid vector of the simulation box local to this process
1379!> \param y_locl y grid vector of the simulation box local to this process
1380!> \param z_locl z grid vector of the simulation box local to this process
1381!> \param tile_vertices coordinates of the vertices of the input tile
1382!> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
1383!> \param is_periodic whether or not the tile is periodic with a periodicity of one unit cell
1384!> \param tile_pw the output pw data
1385!> \par History
1386!> 10.2015 created [Hossein Bani-Hashemian]
1387!> \author Mohammad Hossein Bani-Hashemian
1388! **************************************************************************************************
1389 SUBROUTINE aa_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, tile_vertices, &
1390 sigma, is_periodic, tile_pw)
1391
1392 REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
1393 REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_locl, y_locl, z_locl
1394 REAL(dp), DIMENSION(3, 8), INTENT(IN) :: tile_vertices
1395 REAL(dp), INTENT(IN) :: sigma
1396 LOGICAL, INTENT(IN) :: is_periodic
1397 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: tile_pw
1398
1399 CHARACTER(LEN=*), PARAMETER :: routinen = 'aa_tile_pw_compute'
1400
1401 INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
1402 ub2, ub3
1403 INTEGER, DIMENSION(2, 3) :: bounds_local
1404 REAL(dp) :: fx, fy, fz, gx, gy, gz, hx, hy, hz, xi0, &
1405 xil, xir, yj0, yjl, yjr, zk0, zkl, zkr
1406 REAL(dp), ALLOCATABLE, DIMENSION(:) :: xlocll, xloclr, ylocll, yloclr, zlocll, &
1407 zloclr
1408 REAL(dp), DIMENSION(2) :: x_xtnt, y_xtnt, z_xtnt
1409 REAL(dp), DIMENSION(3) :: per
1410 REAL(dp), DIMENSION(3, 3) :: prefactor
1411
1412 CALL timeset(routinen, handle)
1413
1414 bounds_local = tile_pw%pw_grid%bounds_local
1415 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
1416 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
1417 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
1418
1419 ALLOCATE (xlocll(lb1:ub1), ylocll(lb2:ub2), zlocll(lb3:ub3))
1420 ALLOCATE (xloclr(lb1:ub1), yloclr(lb2:ub2), zloclr(lb3:ub3))
1421
1422! periodicities
1423 per = cell_xtnts(2, :) - cell_xtnts(1, :)
1424
1425 x_xtnt(1) = minval(tile_vertices(1, :)); x_xtnt(2) = maxval(tile_vertices(1, :))
1426 y_xtnt(1) = minval(tile_vertices(2, :)); y_xtnt(2) = maxval(tile_vertices(2, :))
1427 z_xtnt(1) = minval(tile_vertices(3, :)); z_xtnt(2) = maxval(tile_vertices(3, :))
1428
1429! prefactor: columns: left, original, right - rows: x , y, z directions
1430 prefactor = 0.5_dp
1431
1432 IF (is_periodic) THEN
1433 DO k = lb3, ub3
1434 DO j = lb2, ub2
1435 DO i = lb1, ub1
1436 xi0 = x_locl(i); yj0 = y_locl(j); zk0 = z_locl(k)
1437 xil = x_locl(i) - per(1); yjl = y_locl(j) - per(2); zkl = z_locl(k) - per(3)
1438 xir = x_locl(i) + per(1); yjr = y_locl(j) + per(2); zkr = z_locl(k) + per(3)
1439
1440 ! points from the original cell local to the current processor
1441 fx = prefactor(1, 2)*(erf((xi0 - x_xtnt(1))/sigma) - erf((xi0 - x_xtnt(2))/sigma))
1442 fy = prefactor(2, 2)*(erf((yj0 - y_xtnt(1))/sigma) - erf((yj0 - y_xtnt(2))/sigma))
1443 fz = prefactor(3, 2)*(erf((zk0 - z_xtnt(1))/sigma) - erf((zk0 - z_xtnt(2))/sigma))
1444
1445 ! periodically replicated cell on the left, points local to the current processor
1446 gx = prefactor(1, 1)*(erf((xil - x_xtnt(1))/sigma) - erf((xil - x_xtnt(2))/sigma))
1447 gy = prefactor(2, 1)*(erf((yjl - y_xtnt(1))/sigma) - erf((yjl - y_xtnt(2))/sigma))
1448 gz = prefactor(3, 1)*(erf((zkl - z_xtnt(1))/sigma) - erf((zkl - z_xtnt(2))/sigma))
1449
1450 ! periodically replicated cell on the right, points local to the current processor
1451 hx = prefactor(1, 3)*(erf((xir - x_xtnt(1))/sigma) - erf((xir - x_xtnt(2))/sigma))
1452 hy = prefactor(2, 3)*(erf((yjr - y_xtnt(1))/sigma) - erf((yjr - y_xtnt(2))/sigma))
1453 hz = prefactor(3, 3)*(erf((zkr - z_xtnt(1))/sigma) - erf((zkr - z_xtnt(2))/sigma))
1454
1455 tile_pw%array(i, j, k) = (fx + gx + hx)*(fy + gy + hy)*(fz + gz + hz)
1456 END DO
1457 END DO
1458 END DO
1459 ELSE
1460 DO k = lb3, ub3
1461 DO j = lb2, ub2
1462 DO i = lb1, ub1
1463 xi0 = x_locl(i); yj0 = y_locl(j); zk0 = z_locl(k)
1464
1465 fx = 0.5_dp*(erf((xi0 - x_xtnt(1))/sigma) - erf((xi0 - x_xtnt(2))/sigma))
1466 fy = 0.5_dp*(erf((yj0 - y_xtnt(1))/sigma) - erf((yj0 - y_xtnt(2))/sigma))
1467 fz = 0.5_dp*(erf((zk0 - z_xtnt(1))/sigma) - erf((zk0 - z_xtnt(2))/sigma))
1468
1469 tile_pw%array(i, j, k) = fx*fy*fz
1470 END DO
1471 END DO
1472 END DO
1473 END IF
1474
1475 CALL timestop(handle)
1476
1477 END SUBROUTINE aa_tile_pw_compute
1478
1479! **************************************************************************************************
1480!> \brief computes the pw data corresponding to an arbitrary (rectangular-shaped) tile
1481!> \param cell_xtnts extents of the simulation cell
1482!> \param x_locl x grid vector of the simulation box local to this process
1483!> \param y_locl y grid vector of the simulation box local to this process
1484!> \param z_locl z grid vector of the simulation box local to this process
1485!> \param tile_vertices coordinates of the vertices of the input tile
1486!> \param sigma the width of the transition region between 0 and 1 at the edges of the tile
1487!> \param tile_pw the output pw data
1488!> \par History
1489!> 11.2015 created [Hossein Bani-Hashemian]
1490!> \author Mohammad Hossein Bani-Hashemian
1491! **************************************************************************************************
1492 SUBROUTINE arbitrary_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, tile_vertices, &
1493 sigma, tile_pw)
1494
1495 REAL(dp), DIMENSION(2, 3), INTENT(IN) :: cell_xtnts
1496 REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: x_locl, y_locl, z_locl
1497 REAL(dp), DIMENSION(3, 8), INTENT(IN) :: tile_vertices
1498 REAL(dp), INTENT(IN) :: sigma
1499 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: tile_pw
1500
1501 CHARACTER(LEN=*), PARAMETER :: routinen = 'arbitrary_tile_pw_compute'
1502
1503 INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
1504 ub2, ub3
1505 INTEGER, DIMENSION(2, 3) :: bounds_local
1506 REAL(dp) :: cm_x, cm_y, cm_z, fx, fy, fz, xi0, yj0, &
1507 zk0
1508 REAL(dp), ALLOCATABLE, DIMENSION(:) :: xlocll, xloclr, ylocll, yloclr, zlocll, &
1509 zloclr
1510 REAL(dp), DIMENSION(2) :: transfdx_xtnt, transfdy_xtnt, &
1511 transfdz_xtnt, x_xtnt, y_xtnt, z_xtnt
1512 REAL(dp), DIMENSION(3) :: per, pnt0, tpnt
1513 REAL(dp), DIMENSION(3, 3) :: prefactor, rmat
1514 REAL(dp), DIMENSION(3, 8) :: tile_transfd_vertices
1515
1516 CALL timeset(routinen, handle)
1517
1518 bounds_local = tile_pw%pw_grid%bounds_local
1519 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1)
1520 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2)
1521 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3)
1522
1523 ALLOCATE (xlocll(lb1:ub1), ylocll(lb2:ub2), zlocll(lb3:ub3))
1524 ALLOCATE (xloclr(lb1:ub1), yloclr(lb2:ub2), zloclr(lb3:ub3))
1525
1526 per = cell_xtnts(2, :) - cell_xtnts(1, :)
1527
1528 CALL rotate_translate_cuboid(tile_vertices, rmat, tpnt, tile_transfd_vertices)
1529
1530 transfdx_xtnt(1) = minval(tile_transfd_vertices(1, :))
1531 transfdx_xtnt(2) = maxval(tile_transfd_vertices(1, :))
1532 transfdy_xtnt(1) = minval(tile_transfd_vertices(2, :))
1533 transfdy_xtnt(2) = maxval(tile_transfd_vertices(2, :))
1534 transfdz_xtnt(1) = minval(tile_transfd_vertices(3, :))
1535 transfdz_xtnt(2) = maxval(tile_transfd_vertices(3, :))
1536
1537 cm_x = 0.5_dp*(transfdx_xtnt(2) + transfdx_xtnt(1))
1538 cm_y = 0.5_dp*(transfdy_xtnt(2) + transfdy_xtnt(1))
1539 cm_z = 0.5_dp*(transfdz_xtnt(2) + transfdz_xtnt(1))
1540
1541 x_xtnt = transfdx_xtnt - cm_x
1542 y_xtnt = transfdy_xtnt - cm_y
1543 z_xtnt = transfdz_xtnt - cm_z
1544
1545 prefactor = 0.5_dp
1546
1547 DO k = lb3, ub3
1548 DO j = lb2, ub2
1549 DO i = lb1, ub1
1550 pnt0 = rotate_translate_vector(rmat, tpnt, bwrot, (/x_locl(i), y_locl(j), z_locl(k)/))
1551
1552 xi0 = pnt0(1); yj0 = pnt0(2); zk0 = pnt0(3)
1553
1554 fx = 0.5_dp*(erf((xi0 - x_xtnt(1))/sigma) - erf((xi0 - x_xtnt(2))/sigma))
1555 fy = 0.5_dp*(erf((yj0 - y_xtnt(1))/sigma) - erf((yj0 - y_xtnt(2))/sigma))
1556 fz = 0.5_dp*(erf((zk0 - z_xtnt(1))/sigma) - erf((zk0 - z_xtnt(2))/sigma))
1557
1558 tile_pw%array(i, j, k) = fx*fy*fz
1559 END DO
1560 END DO
1561 END DO
1562
1563 CALL timestop(handle)
1564
1565 END SUBROUTINE arbitrary_tile_pw_compute
1566
1567END MODULE dirichlet_bc_methods
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_zero
Definition pw_methods.F:82
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
dirichlet_bc_methods
subroutines for defining and creating Dirichlet type subdomains
Definition dirichlet_bc_methods.F:15
dirichlet_bc_methods::dirichlet_boundary_region_setup
subroutine, public dirichlet_boundary_region_setup(pw_pool, poisson_params, dbcs)
Sets up the Dirichlet boundary condition.
Definition dirichlet_bc_methods.F:65
dirichlet_bc_types
Dirichlet boundary condition data types.
Definition dirichlet_bc_types.F:15
dirichlet_bc_types::dbc_parameters_dealloc
subroutine, public dbc_parameters_dealloc(dbc_params)
deallocates dirichlet_bc_parameters type
Definition dirichlet_bc_types.F:234
dirichlet_bc_types::z_axis
integer, parameter, public z_axis
Definition dirichlet_bc_types.F:140
dirichlet_bc_types::inscribed
integer, parameter, public inscribed
Definition dirichlet_bc_types.F:133
dirichlet_bc_types::cylindrical
integer, parameter, public cylindrical
Definition dirichlet_bc_types.F:128
dirichlet_bc_types::aa_cuboidal
integer, parameter, public aa_cuboidal
Definition dirichlet_bc_types.F:128
dirichlet_bc_types::y_axis
integer, parameter, public y_axis
Definition dirichlet_bc_types.F:140
dirichlet_bc_types::circumscribed
integer, parameter, public circumscribed
Definition dirichlet_bc_types.F:133
dirichlet_bc_types::x_axis
integer, parameter, public x_axis
Definition dirichlet_bc_types.F:140
dirichlet_bc_types::aa_planar
integer, parameter, public aa_planar
Definition dirichlet_bc_types.F:128
dirichlet_bc_types::xz_plane
integer, parameter, public xz_plane
Definition dirichlet_bc_types.F:136
dirichlet_bc_types::planar
integer, parameter, public planar
Definition dirichlet_bc_types.F:128
dirichlet_bc_types::yz_plane
integer, parameter, public yz_plane
Definition dirichlet_bc_types.F:136
dirichlet_bc_types::xy_plane
integer, parameter, public xy_plane
Definition dirichlet_bc_types.F:136
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::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::inv_3x3
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition mathlib.F:516
mathlib::vector_product
pure real(kind=dp) function, dimension(3), public vector_product(a, b)
Calculation of the vector product c = a x b.
Definition mathlib.F:1274
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
ps_implicit_types
Types containing essential information for running implicit (iterative) Poisson solver.
Definition ps_implicit_types.F:15
ps_implicit_types::neumann_bc
integer, parameter, public neumann_bc
Definition ps_implicit_types.F:32
ps_implicit_types::mixed_bc
integer, parameter, public mixed_bc
Definition ps_implicit_types.F:32
ps_implicit_types::mixed_periodic_bc
integer, parameter, public mixed_periodic_bc
Definition ps_implicit_types.F:32
ps_implicit_types::periodic_bc
integer, parameter, public periodic_bc
Definition ps_implicit_types.F:32
pw_grid_types
Definition pw_grid_types.F:13
pw_methods
Definition pw_methods.F:25
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
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_types
Definition pw_types.F:24
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
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
dirichlet_bc_types::dirichlet_bc_p_type
Definition dirichlet_bc_types.F:104
dirichlet_bc_types::dirichlet_bc_type
Definition dirichlet_bc_types.F:90
dirichlet_bc_types::tile_p_type
Definition dirichlet_bc_types.F:115
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:62
pw_poisson_types::pw_poisson_parameter_type
parameters for the poisson solver independet of input_section
Definition pw_poisson_types.F:96
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_r3d_rs_type
Definition pw_types.F:62