41 #include "../base/base_uses.f90"
45 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'dirichlet_bc_methods'
49 REAL(dp),
PARAMETER,
PRIVATE :: small_value = 1.0e-8_dp
50 INTEGER,
PARAMETER,
PRIVATE :: FWROT = 1, &
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
71 CHARACTER(LEN=*),
PARAMETER :: routinen =
'dirichlet_boundary_region_setup'
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, &
81 REAL(dp),
ALLOCATABLE,
DIMENSION(:) :: x_glbl, x_locl, y_glbl, y_locl, z_glbl, &
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
89 CALL timeset(routinen, handle)
92 IF (logger%para_env%is_source())
THEN
98 dr = pw_pool%pw_grid%dr
99 npts = pw_pool%pw_grid%npts
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)
105 verbose = poisson_params%dbc_params%verbose_output
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)
116 SELECT CASE (poisson_params%ps_implicit_params%boundary_condition)
118 IF (n_dbcs .EQ. 0)
THEN
119 cpabort(
"No Dirichlet region is defined.")
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 ..."
128 DO j = 1, n_aa_planar
129 ALLOCATE (dbcs(j)%dirichlet_bc)
130 n_prtn = poisson_params%dbc_params%aa_planar_nprtn(:, 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)
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]"
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)
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)
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)
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)
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]"
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)
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)
185 l = n_aa_planar + n_aa_cuboidal
186 DO j = l + 1, l + n_planar
187 ALLOCATE (dbcs(j)%dirichlet_bc)
189 n_prtn(1:2) = poisson_params%dbc_params%planar_nprtn(:, 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)
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]"
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)
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)
215 l = n_aa_planar + n_aa_cuboidal + n_planar
216 DO j = 1, n_aa_cylindrical
218 ind_end = l + aa_cylindrical_nsides(j)
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)
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]"
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))
246 l = l + aa_cylindrical_nsides(j)
256 CALL timestop(handle)
277 SUBROUTINE dirichlet_bc_partition(sigma, n_prtn, pw_pool, cell_xtnts, &
278 x_locl, y_locl, z_locl, is_periodic, verbose, dirichlet_bc)
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
288 CHARACTER(LEN=*),
PARAMETER :: routinen =
'dirichlet_bc_partition'
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
297 CALL timeset(routinen, handle)
300 IF (logger%para_env%is_source())
THEN
306 SELECT CASE (dirichlet_bc%dbc_geom)
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
316 ALLOCATE (tile_maxvals(n_tiles))
317 tile_maxvals = 0.0_dp
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)
323 IF ((unit_nr .GT. 0) .AND. verbose)
THEN
324 WRITE (unit_nr,
'(T7,A,I5)')
"tile", k
326 WRITE (unit_nr,
'(T10,A,I1,3F10.3)') &
327 " vertex ", i,
angstrom*dirichlet_bc%tiles(k)%tile%vertices(:, i)
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
338 IF ((unit_nr .GT. 0) .AND. verbose)
WRITE (unit_nr,
'(T3,A)') repeat(
'=', 78)
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
351 NULLIFY (cylinder_pw)
352 ALLOCATE (cylinder_pw)
353 CALL pw_pool%create_pw(cylinder_pw)
354 CALL pw_zero(cylinder_pw)
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)
361 IF ((unit_nr .GT. 0) .AND. verbose)
THEN
362 WRITE (unit_nr,
'(T7,A,I5)')
"tile", k
364 WRITE (unit_nr,
'(T10,A,I1,3F10.3)') &
365 " vertex ", i,
angstrom*dirichlet_bc%tiles(k)%tile%vertices(:, i)
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
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).")
386 CALL pw_pool%give_back_pw(cylinder_pw)
387 DEALLOCATE (cylinder_pw)
389 IF ((unit_nr .GT. 0) .AND. verbose)
WRITE (unit_nr,
'(T3,A)') repeat(
'=', 78)
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
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)
407 IF ((unit_nr .GT. 0) .AND. verbose)
THEN
408 WRITE (unit_nr,
'(T7,A,I5)')
"tile", k
410 WRITE (unit_nr,
'(T10,A,I1,3F10.3)') &
411 " vertex ", i,
angstrom*dirichlet_bc%tiles(k)%tile%vertices(:, i)
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
422 IF ((unit_nr .GT. 0) .AND. verbose)
WRITE (unit_nr,
'(T3,A)') repeat(
'=', 78)
428 CALL timestop(handle)
430 END SUBROUTINE dirichlet_bc_partition
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)
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, &
460 INTEGER,
INTENT(IN) :: dbc_id
461 LOGICAL,
INTENT(IN) :: verbose
462 TYPE(dirichlet_bc_type),
INTENT(INOUT),
POINTER :: dirichlet_bc
464 CHARACTER(LEN=*),
PARAMETER :: routinen =
'aa_planar_dbc_setup'
465 LOGICAL,
DIMENSION(6),
PARAMETER :: test_forb_xtnts = .true.
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
473 CALL timeset(routinen, handle)
476 IF (logger%para_env%is_source())
THEN
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)
486 SELECT CASE (parallel_plane)
488 zlb = zlb - thickness*0.5_dp
489 zub = zub + thickness*0.5_dp
491 ylb = ylb - thickness*0.5_dp
492 yub = yub + thickness*0.5_dp
494 xlb = xlb - thickness*0.5_dp
495 xub = xub + thickness*0.5_dp
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.")
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
518 dirichlet_bc%smoothing_width = sigma
519 dirichlet_bc%n_tiles = 1
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/)
530 IF ((unit_nr .GT. 0) .AND. verbose)
THEN
531 WRITE (unit_nr,
'(T3,A,A)')
"======== verbose ", repeat(
'=', 61)
533 WRITE (unit_nr,
'(T10,A,I1,3F10.3)')
" vertex ", i,
angstrom*dirichlet_bc%vertices(:, i)
537 CALL timestop(handle)
539 END SUBROUTINE aa_planar_dbc_setup
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)
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
569 CHARACTER(LEN=*),
PARAMETER :: routinen =
'aa_cuboidal_dbc_setup'
570 LOGICAL,
DIMENSION(6),
PARAMETER :: test_forb_xtnts = .true.
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
577 CALL timeset(routinen, handle)
580 IF (logger%para_env%is_source())
THEN
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.")
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
606 dirichlet_bc%smoothing_width = sigma
607 dirichlet_bc%n_tiles = 1
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)/)
618 IF ((unit_nr .GT. 0) .AND. verbose)
THEN
619 WRITE (unit_nr,
'(T3,A,A)')
"======== verbose ", repeat(
'=', 61)
621 WRITE (unit_nr,
'(T10,A,I1,3F10.3)')
" vertex ", i,
angstrom*dirichlet_bc%vertices(:, i)
625 CALL timestop(handle)
627 END SUBROUTINE aa_cuboidal_dbc_setup
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)
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, &
656 INTEGER,
INTENT(IN) :: dbc_id
657 LOGICAL,
INTENT(IN) :: verbose
658 TYPE(dirichlet_bc_type),
INTENT(INOUT),
POINTER :: dirichlet_bc
660 CHARACTER(LEN=*),
PARAMETER :: routinen =
'arbitrary_planar_dbc_setup'
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, &
669 TYPE(cp_logger_type),
POINTER :: logger
671 CALL timeset(routinen, handle)
674 IF (logger%para_env%is_source())
THEN
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.")
689 unit_normal = normal_vector/sqrt(sum(normal_vector**2))
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.")
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.")
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.")
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
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
745 IF ((unit_nr .GT. 0) .AND. verbose)
THEN
746 WRITE (unit_nr,
'(T3,A,A)')
"======== verbose ", repeat(
'=', 61)
748 WRITE (unit_nr,
'(T10,A,I1,3F10.3)')
" vertex ", i,
angstrom*dirichlet_bc%vertices(:, i)
752 CALL timestop(handle)
754 END SUBROUTINE arbitrary_planar_dbc_setup
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)
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
807 CHARACTER(LEN=*),
PARAMETER :: routinen =
'aa_cylindrical_dbc_setup'
809 INTEGER :: handle, i, intern_apx_type, j, n_dbcs, &
811 LOGICAL :: forb_xtnt1, forb_xtnt2, forb_xtnt3, &
813 REAL(dp) :: cntrlaxis_lb, cntrlaxis_ub, h, lx, ly, &
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, &
820 CALL timeset(routinen, handle)
823 IF (logger%para_env%is_source())
THEN
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)
833 SELECT CASE (parallel_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.")
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.")
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.")
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.")
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.")
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.")
876 IF (
PRESENT(apx_type)) intern_apx_type = apx_type
880 cntrlaxis_lb = xtnt(1)
881 cntrlaxis_ub = xtnt(2)
888 h = base_radius/cos(0.5*theta)
890 cpabort(
"Unknown approximation type for cylinder.")
892 SELECT CASE (parallel_axis)
894 IF (h .GT. minval((/ly, lz/)))
THEN
895 cpabort(
"Reduce the base radius!")
898 IF (h .GT. minval((/lx, lz/)))
THEN
899 cpabort(
"Reduce the base radius!")
902 IF (h .GT. minval((/lx, ly/)))
THEN
903 cpabort(
"Reduce the base radius!")
907 alpha = linspace(0.0_dp, 2*
pi, n_dbcs + 1)
908 alpha_rotd = alpha + delta_alpha;
909 SELECT CASE (parallel_axis)
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))
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)
922 ALLOCATE (dbcs(j)%dirichlet_bc)
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
929 dbcs(j)%dirichlet_bc%smoothing_width = sigma
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)/)
936 unit_normal = normal_vec/sqrt(sum(normal_vec**2))
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
947 dbcs(j)%dirichlet_bc%n_tiles = 1
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"
954 WRITE (unit_nr,
'(T10,A,I1,3F10.3)')
" vertex ", i, &
955 angstrom*dbcs(j)%dirichlet_bc%vertices(:, i)
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)
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))
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)
974 ALLOCATE (dbcs(j)%dirichlet_bc)
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
981 dbcs(j)%dirichlet_bc%smoothing_width = sigma
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)/)
988 unit_normal = normal_vec/sqrt(sum(normal_vec**2))
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
999 dbcs(j)%dirichlet_bc%n_tiles = 1
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"
1006 WRITE (unit_nr,
'(T10,A,I1,3F10.3)')
" vertex ", i, &
1007 angstrom*dbcs(j)%dirichlet_bc%vertices(:, i)
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)
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))
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)
1026 ALLOCATE (dbcs(j)%dirichlet_bc)
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
1033 dbcs(j)%dirichlet_bc%smoothing_width = sigma
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/)
1040 unit_normal = normal_vec/sqrt(sum(normal_vec**2))
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
1051 dbcs(j)%dirichlet_bc%n_tiles = 1
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"
1058 WRITE (unit_nr,
'(T10,A,I1,3F10.3)')
" vertex ", i, &
1059 angstrom*dbcs(j)%dirichlet_bc%vertices(:, i)
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)
1068 CALL timestop(handle)
1070 END SUBROUTINE aa_cylindrical_dbc_setup
1082 FUNCTION linspace(a, b, N)
RESULT(arr)
1084 REAL(dp),
INTENT(IN) :: a, b
1085 INTEGER,
INTENT(IN) :: n
1086 REAL(dp),
DIMENSION(N) :: arr
1091 dx = (b - a)/(n - 1)
1092 arr = (/(a + (i - 1)*dx, i=1, n)/)
1094 END FUNCTION linspace
1106 SUBROUTINE rotate_translate_cuboid(cuboid_vtx, Rmat, Tpnt, cuboid_transfd_vtx)
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
1113 CHARACTER(LEN=*),
PARAMETER :: routinen =
'rotate_translate_cuboid'
1116 REAL(dp),
DIMENSION(3) :: a, a2, aa2, ab, ad, b, b2, c, c2, d, d2, &
1117 e1_locl, e2_locl, e3_locl
1119 CALL timeset(routinen, handle)
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)
1128 tpnt = (a + c2)/2.0_dp
1135 e1_locl = ab/sqrt(sum(ab**2))
1136 e2_locl = ad/sqrt(sum(ad**2))
1137 e3_locl = aa2/sqrt(sum(aa2**2))
1139 rmat(1:3, 1) = e1_locl
1140 rmat(1:3, 2) = e2_locl
1141 rmat(1:3, 3) = e3_locl
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)
1152 CALL timestop(handle)
1154 END SUBROUTINE rotate_translate_cuboid
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
1175 REAL(dp),
DIMENSION(3) :: tpoint, vec_tmp
1176 REAL(dp),
DIMENSION(3, 3) :: rmat_inv
1180 IF (direction .EQ. fwrot)
THEN
1181 vec_tmp = matmul(rmat, vec)
1182 vec_transfd = vec_tmp + tp
1183 ELSEIF (direction .EQ. bwrot)
THEN
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
1192 END FUNCTION rotate_translate_vector
1204 SUBROUTINE aa_dbc_partition(dbc_vertices, n_prtn, tiles)
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), &
1211 CHARACTER(LEN=*),
PARAMETER :: routinen =
'aa_dbc_partition'
1213 INTEGER :: handle, i, ii, jj, k, kk
1215 REAL(dp),
ALLOCATABLE,
DIMENSION(:) :: xprtn_lb, xprtn_ub, yprtn_lb, yprtn_ub, &
1217 REAL(dp),
DIMENSION(2) :: x_xtnt, y_xtnt, z_xtnt
1219 CALL timeset(routinen, handle)
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)))
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, :))
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
1241 xprtn_ub(:) = x_xtnt(1) + (/(i, i=1, n_prtn(1))/)*step
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
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
1251 ALLOCATE (tiles(n_prtn(1)*n_prtn(2)*n_prtn(3)))
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
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)/)
1268 tiles(k)%tile%volume = 0.0_dp
1274 CALL timestop(handle)
1276 END SUBROUTINE aa_dbc_partition
1287 SUBROUTINE arbitrary_dbc_partition(dbc_vertices, n_prtn, tiles)
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), &
1294 CHARACTER(LEN=*),
PARAMETER :: routinen =
'arbitrary_dbc_partition'
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, &
1303 CALL timeset(routinen, handle)
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)
1315 unit_normal = normal_vector/sqrt(sum(normal_vector**2))
1316 thickness = sqrt(sum((d - d2)**2))
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)
1325 np1 = min(n_prtn(1) + 1, n_prtn(2) + 1)
1326 np2 = max(n_prtn(1) + 1, n_prtn(2) + 1)
1329 ALLOCATE (x(np1, np2, 3))
1330 ALLOCATE (end_points(2, np1, 3))
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
1342 point1(:) = end_points(1, l, :)
1343 point2(:) = end_points(2, l, :)
1345 step = real((i - 1), kind=dp)/real((np2 - 1), kind=dp)
1346 x(l, i, :) = point1*(1.0_dp - step) + point2*step
1350 ALLOCATE (tiles((np1 - 1)*(np2 - 1)))
1354 ALLOCATE (tiles(k)%tile)
1355 tiles(k)%tile%tile_id = k
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
1366 tiles(k)%tile%volume = 0.0_dp
1371 CALL timestop(handle)
1373 END SUBROUTINE arbitrary_dbc_partition
1389 SUBROUTINE aa_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, tile_vertices, &
1390 sigma, is_periodic, tile_pw)
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
1399 CHARACTER(LEN=*),
PARAMETER :: routinen =
'aa_tile_pw_compute'
1401 INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
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, &
1408 REAL(dp),
DIMENSION(2) :: x_xtnt, y_xtnt, z_xtnt
1409 REAL(dp),
DIMENSION(3) :: per
1410 REAL(dp),
DIMENSION(3, 3) :: prefactor
1412 CALL timeset(routinen, handle)
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)
1419 ALLOCATE (xlocll(lb1:ub1), ylocll(lb2:ub2), zlocll(lb3:ub3))
1420 ALLOCATE (xloclr(lb1:ub1), yloclr(lb2:ub2), zloclr(lb3:ub3))
1423 per = cell_xtnts(2, :) - cell_xtnts(1, :)
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, :))
1432 IF (is_periodic)
THEN
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)
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))
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))
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))
1455 tile_pw%array(i, j, k) = (fx + gx + hx)*(fy + gy + hy)*(fz + gz + hz)
1463 xi0 = x_locl(i); yj0 = y_locl(j); zk0 = z_locl(k)
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))
1469 tile_pw%array(i, j, k) = fx*fy*fz
1475 CALL timestop(handle)
1477 END SUBROUTINE aa_tile_pw_compute
1492 SUBROUTINE arbitrary_tile_pw_compute(cell_xtnts, x_locl, y_locl, z_locl, tile_vertices, &
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
1501 CHARACTER(LEN=*),
PARAMETER :: routinen =
'arbitrary_tile_pw_compute'
1503 INTEGER :: handle, i, j, k, lb1, lb2, lb3, ub1, &
1505 INTEGER,
DIMENSION(2, 3) :: bounds_local
1506 REAL(dp) :: cm_x, cm_y, cm_z, fx, fy, fz, xi0, yj0, &
1508 REAL(dp),
ALLOCATABLE,
DIMENSION(:) :: xlocll, xloclr, ylocll, yloclr, zlocll, &
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
1516 CALL timeset(routinen, handle)
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)
1523 ALLOCATE (xlocll(lb1:ub1), ylocll(lb2:ub2), zlocll(lb3:ub3))
1524 ALLOCATE (xloclr(lb1:ub1), yloclr(lb2:ub2), zloclr(lb3:ub3))
1526 per = cell_xtnts(2, :) - cell_xtnts(1, :)
1528 CALL rotate_translate_cuboid(tile_vertices, rmat, tpnt, tile_transfd_vertices)
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, :))
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))
1541 x_xtnt = transfdx_xtnt - cm_x
1542 y_xtnt = transfdy_xtnt - cm_y
1543 z_xtnt = transfdz_xtnt - cm_z
1550 pnt0 = rotate_translate_vector(rmat, tpnt, bwrot, (/x_locl(i), y_locl(j), z_locl(k)/))
1552 xi0 = pnt0(1); yj0 = pnt0(2); zk0 = pnt0(3)
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))
1558 tile_pw%array(i, j, k) = fx*fy*fz
1563 CALL timestop(handle)
1565 END SUBROUTINE arbitrary_tile_pw_compute
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.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public twopi
Collection of simple mathematical functions and subroutines.
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
pure real(kind=dp) function, dimension(3), public vector_product(a, b)
Calculation of the vector product c = a x b.
Definition of physical constants:
real(kind=dp), parameter, public angstrom
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 ...
numerical operations on real-space grid
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