(git:0de0cc2)
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 ! **************************************************************************************************
16 
19  cp_logger_type
20  USE dirichlet_bc_types, ONLY: &
22  dbc_parameters_dealloc, dirichlet_bc_p_type, dirichlet_bc_type, tile_p_type, x_axis, &
24  USE kinds, ONLY: dp
25  USE mathconstants, ONLY: pi,&
26  twopi
27  USE mathlib, ONLY: inv_3x3,&
29  USE physcon, ONLY: angstrom
30  USE ps_implicit_types, ONLY: mixed_bc,&
32  neumann_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 
48 
49  REAL(dp), PARAMETER, PRIVATE :: small_value = 1.0e-8_dp
50  INTEGER, PARAMETER, PRIVATE :: FWROT = 1, &
51  bwrot = -1
52 
53 CONTAINS
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)
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 
1567 END MODULE dirichlet_bc_methods
various routines to log and control the output. The idea is that decisions about where to log should ...
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
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
subroutines for defining and creating Dirichlet type subdomains
subroutine, public dirichlet_boundary_region_setup(pw_pool, poisson_params, dbcs)
Sets up the Dirichlet boundary condition.
Dirichlet boundary condition data types.
subroutine, public dbc_parameters_dealloc(dbc_params)
deallocates dirichlet_bc_parameters type
integer, parameter, public z_axis
integer, parameter, public inscribed
integer, parameter, public cylindrical
integer, parameter, public aa_cuboidal
integer, parameter, public y_axis
integer, parameter, public circumscribed
integer, parameter, public x_axis
integer, parameter, public aa_planar
integer, parameter, public xz_plane
integer, parameter, public planar
integer, parameter, public yz_plane
integer, parameter, public xy_plane
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public twopi
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
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:520
pure real(kind=dp) function, dimension(3), public vector_product(a, b)
Calculation of the vector product c = a x b.
Definition: mathlib.F:1282
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
Types containing essential information for running implicit (iterative) Poisson solver.
integer, parameter, public neumann_bc
integer, parameter, public mixed_bc
integer, parameter, public mixed_periodic_bc
integer, parameter, public periodic_bc
functions related to the poisson solver on regular grids
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
numerical operations on real-space grid
Definition: rs_methods.F:14
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