Octopus
pes_flux.F90
Go to the documentation of this file.
1!! Copyright (C) 2015 P. Wopperer and U. De Giovannini
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module pes_flux_oct_m
24 use box_oct_m
28 use comm_oct_m
29 use debug_oct_m
33 use global_oct_m
37 use io_oct_m
38 use, intrinsic :: iso_fortran_env
39 use lasers_oct_m
40 use loct_oct_m
42 use math_oct_m
44 use mesh_oct_m
47 use mpi_oct_m
49 use parser_oct_m
52 use space_oct_m
55 use string_oct_m
56 use unit_oct_m
58 use utils_oct_m
60
61 implicit none
62
63 private
64
65 public :: &
66 pes_flux_t, &
78
79 type pes_flux_t
80 private
81
84
85 integer :: nkpnts
86 integer :: nkpnts_start, nkpnts_end
87 integer :: nk
88 integer :: nk_start, nk_end
89
90 ! spherical momentum grid
91 integer :: nstepsthetak, nstepsphik
92 real(real64) :: thetak_rng(1:2)
93 real(real64) :: phik_rng(1:2)
94 integer :: nstepsomegak
95 integer :: nstepsomegak_start, nstepsomegak_end
96 real(real64) :: ktransf(1:3,1:3)
97
98 real(real64), allocatable :: klinear(:,:)
99 ! !< polar: klinear(nk,1) defines the radial part of the k-mesh
100 ! !< cartesian: klinear(nk,1:3) defines the the k-mesh along each axes
101
102 real(real64) :: dk
103 real(real64), allocatable :: kcoords_cub(:,:,:)
104 real(real64), allocatable :: kcoords_sph(:,:,:)
105 integer :: kgrid
106
107 ! Surface relates quantities
108 integer, public :: surf_shape
109 integer :: nsrfcpnts
110 integer :: nsrfcpnts_start, nsrfcpnts_end
111 real(real64), allocatable :: srfcnrml(:,:)
112 real(real64), allocatable :: rcoords(:,:)
113 integer, allocatable :: srfcpnt(:)
114 integer, allocatable :: rankmin(:)
115 integer :: lmax
116 complex(real64), allocatable :: ylm_r(:,:,:)
117 complex(real64), allocatable :: ylm_k(:,:,:)
118 real(real64), allocatable :: j_l(:,:)
119 real(real64) :: radius
120 integer, allocatable :: face_idx_range(:,:)
121 ! !! 1 (2) start (end) idx of face nface in rcoords(:,:)
122
123
124 complex(real64), allocatable :: bvk_phase(:,:)
125 complex(real64), allocatable :: expkr(:,:,:,:)
126 complex(real64), allocatable :: expkr_perp(:,:)
127 ! !! direction perpendicular to the face
128 real(real64), allocatable :: LLr(:,:)
129 integer, allocatable :: NN(:,:)
130 integer, allocatable :: Lkpuvz_inv(:,:,:)
131 ! !! (u,v parametric on face and z perpendicular)
132
133 ! Surface and time integration
134 integer :: tdsteps
135 ! !< = mesh%mpi_grp%size (PES_SPHERICAL)
136 integer :: max_iter
137 integer :: save_iter
138 integer :: itstep
139
140 complex(real64), allocatable :: wf(:,:,:,:,:)
141 complex(real64), allocatable :: gwf(:,:,:,:,:,:)
142 real(real64), allocatable :: veca(:,:)
143 complex(real64), allocatable :: conjgphase_prev(:,:)
144 complex(real64), allocatable :: spctramp_cub(:,:,:,:)
145 complex(real64), allocatable :: spctramp_sph(:,:,:,:,:)
146
147 integer, public :: ll(3)
148 ! !< mesh. Used when working with semi-periodic systems
149
150 type(mesh_interpolation_t) :: interp
151
152 logical :: parallel_in_momentum
153 logical :: arpes_grid
154 logical :: surf_interp
155 logical :: use_symmetry
156 logical :: runtime_output
157 logical :: anisotrpy_correction
158
159 integer :: par_strategy
160 integer :: dim
161 integer :: pdim
162
163 end type pes_flux_t
164
165 integer, parameter :: &
166 PES_CUBIC = 1, &
167 pes_spherical = 2, &
168 pes_plane = 3
169
170 integer, parameter :: &
171 PES_POLAR = 1, &
173
174
175contains
176
177
178 ! ---------------------------------------------------------
179 subroutine pes_flux_init(this, namespace, space, mesh, st, ext_partners, kpoints, abb, save_iter, max_iter)
180 type(pes_flux_t), intent(inout) :: this
181 type(namespace_t), intent(in) :: namespace
182 class(space_t), intent(in) :: space
183 type(mesh_t), intent(in) :: mesh
184 type(states_elec_t), intent(in) :: st
185 type(partner_list_t), intent(in) :: ext_partners
186 type(kpoints_t), intent(in) :: kpoints
187 type(absorbing_boundaries_t), intent(in) :: abb
188 integer, intent(in) :: save_iter
189 integer, intent(in) :: max_iter
190
191 type(block_t) :: blk
192 real(real64) :: border(space%dim) ! distance of surface from border
193 real(real64) :: offset(space%dim) ! offset for border
194 integer :: stst, stend, kptst, kptend, sdim, mdim, pdim
195 integer :: imdim
196 integer :: isp
197 integer :: il
199 integer :: nstepsphir, nstepsthetar
200 integer :: ll, mm
201 integer :: default_shape
202 real(real64) :: fc_ptdens
204 integer :: par_strategy
205 logical :: use_symmetry
207 type(lasers_t), pointer :: lasers
209 push_sub_with_profile(pes_flux_init)
211 stst = st%st_start
212 stend = st%st_end
213 kptst = st%d%kpt%start
214 kptend = st%d%kpt%end
215 sdim = st%d%dim
216 mdim = space%dim
217 pdim = space%periodic_dim
219 this%surf_interp = .false.
220
222 lasers => list_get_lasers(ext_partners)
223 if(associated(lasers)) then
224 do il = 1, lasers%no_lasers
225 if (laser_kind(lasers%lasers(il)) /= e_field_vector_potential) then
226 message(1) = 't-surff only works in velocity gauge.'
227 call messages_fatal(1, namespace=namespace)
228 end if
229 end do
230 end if
232 message(1) = 'Info: Calculating PES using t-surff technique.'
233 call messages_info(1, namespace=namespace)
235 ! -----------------------------------------------------------------
236 ! Setting up r-mesh (the surface)
237 ! -----------------------------------------------------------------
238 !%Variable PES_Flux_Shape
239 !%Type integer
240 !%Section Time-Dependent::PhotoElectronSpectrum
241 !%Description
242 !% The shape of the surface.
243 !%Option cub 1
244 !% Uses a parallelepiped as surface. All surface points are grid points.
245 !% Choose the location of the surface with variable <tt>PES_Flux_Lsize</tt>
246 !% (default for 1D and 2D).
247 !%Option sph 2
248 !% Constructs a sphere with radius <tt>PES_Flux_Radius</tt>. Points on the sphere
249 !% are interpolated by trilinear interpolation (default for 3D).
250 !%Option pln 3
251 !% This option is for periodic systems.
252 !% Constructs a plane perpendicular to the non-periodic dimension
253 !% at <tt>PES_Flux_Lsize</tt>.
254 !%End
255
256 default_shape = pes_spherical
257 if (space%is_periodic()) then
258 default_shape = pes_plane
259 else if (mdim <= 2) then
260 default_shape = pes_cubic
261 else
262 select type (box => mesh%box)
264 default_shape = pes_cubic
265 end select
266 end if
267
268 call parse_variable(namespace, 'PES_Flux_Shape', default_shape, this%surf_shape)
269 if (.not. varinfo_valid_option('PES_Flux_Shape', this%surf_shape, is_flag = .true.)) then
270 call messages_input_error(namespace,'PES_Flux_Shape')
271 end if
272 if (this%surf_shape == pes_spherical .and. mdim /= 3) then
273 message(1) = 'Spherical grid works only in 3d.'
274 call messages_fatal(1, namespace=namespace)
275 end if
276 call messages_print_var_option('PES_Flux_Shape', this%surf_shape, namespace=namespace)
277
278
279 !%Variable PES_Flux_AnisotropyCorrection
280 !%Type logical
281 !%Section Time-Dependent::PhotoElectronSpectrum
282 !%Description
283 !% Apply anisotropy correction.
284 !%
285 !%End
286 if (this%surf_shape == pes_cubic) then
287 call parse_variable(namespace, 'PES_Flux_AnisotropyCorrection', .false., this%anisotrpy_correction)
288 call messages_print_var_value('PES_Flux_AnisotropyCorrection', this%anisotrpy_correction, namespace=namespace)
289 else
290 this%anisotrpy_correction = .false.
291 end if
292
293 !%Variable PES_Flux_Offset
294 !%Type block
295 !%Section Time-Dependent::PhotoElectronSpectrum
296 !%Description
297 !% Shifts the surface for <tt>PES_Flux_Shape = cub</tt>. The syntax is:
298 !%
299 !% <tt>%PES_Flux_Offset
300 !% <br>&nbsp;&nbsp;xshift | yshift | zshift
301 !% <br>%
302 !% </tt>
303 !%End
304 offset = m_zero
305 if (parse_block(namespace, 'PES_Flux_Offset', blk) == 0) then
306 do imdim = 1, mdim
307 call parse_block_float(blk, 0, imdim - 1, offset(imdim))
308 end do
309 call parse_block_end(blk)
310 end if
311
312 !%Variable PES_Flux_Lmax
313 !%Type integer
314 !%Default 80
315 !%Section Time-Dependent::PhotoElectronSpectrum
316 !%Description
317 !% Maximum order of the spherical harmonic to be integrated on an equidistant spherical
318 !% grid (to be changed to Gauss-Legendre quadrature).
319 !%End
320 if (this%surf_shape == pes_spherical) then
321 call parse_variable(namespace, 'PES_Flux_Lmax', 80, this%lmax)
322 if (this%lmax < 1) call messages_input_error(namespace,'PES_Flux_Lmax', 'must be > 0')
323 call messages_print_var_value('PES_Flux_Lmax', this%lmax, namespace=namespace)
324 end if
325
326
327 if (this%surf_shape == pes_cubic .or. this%surf_shape == pes_plane) then
328
329
330 !%Variable PES_Flux_Lsize
331 !%Type block
332 !%Section Time-Dependent::PhotoElectronSpectrum
333 !%Description
334 !% For <tt>PES_Flux_Shape = cub</tt> sets the dimensions along each direction. The syntax is:
335 !%
336 !% <tt>%PES_Flux_Lsize
337 !% <br>&nbsp;&nbsp;xsize | ysize | zsize
338 !% <br>%
339 !% </tt>
340 !%
341 !% where <tt>xsize</tt>, <tt>ysize</tt>, and <tt>zsize</tt> are with respect to the origin. The surface can
342 !% be shifted with <tt>PES_Flux_Offset</tt>.
343 !% If <tt>PES_Flux_Shape = pln</tt>, specifies the position of two planes perpendicular to
344 !% the non-periodic dimension symmetrically placed at <tt>PES_Flux_Lsize</tt> distance from
345 !% the origin.
346 !%End
347 if (parse_block(namespace, 'PES_Flux_Lsize', blk) == 0) then
348 do imdim = 1, mdim
349 call parse_block_float(blk, 0, imdim - 1, border(imdim))
350 end do
351 ! Snap the face to the closest grid point
352 border(1:mdim) = int(border(1:mdim) / mesh%spacing(1:mdim)) * mesh%spacing(1:mdim)
353 call parse_block_end(blk)
354
355 else if (parse_is_defined(namespace, 'PES_Flux_Lsize')) then
356 border(mdim) = mesh%box%bounding_box_l(mdim) * m_half
357 if (space%is_periodic()) then
358 ! the cube sides along the periodic directions are out of the simulation box
359 border(1:pdim)= mesh%box%bounding_box_l(1:pdim) * m_two
360 call parse_variable(namespace, 'PES_Flux_Lsize', border(mdim), border(mdim))
361 ! Snap the plane to the closest grid point
362 border(mdim) = floor(border(mdim) / mesh%spacing(mdim)) * mesh%spacing(mdim)
363 else
364 call parse_variable(namespace, 'PES_Flux_Lsize', border(mdim), border(mdim))
365 border(1:mdim - 1) = border(mdim)
366 border(1:mdim) = floor(border(1:mdim) / mesh%spacing(1:mdim)) * mesh%spacing(1:mdim)
367 end if
368 else
369 select type (box => mesh%box)
370 type is (box_sphere_t)
371 border(1:mdim) = box%radius / sqrt(m_two) * m_half
372 type is (box_parallelepiped_t)
373 border(1:mdim) = box%half_length(1:mdim) * m_half
374 class default
375 call messages_write('PES_Flux_Lsize not specified. No default values available for this box shape.')
376 call messages_new_line()
377 call messages_write('Specify the location of the parallelepiped with block PES_Flux_Lsize.')
378 call messages_fatal(namespace=namespace)
379 end select
380 call messages_write('PES_Flux_Lsize not specified. Using default values.')
381 call messages_info(namespace=namespace)
382 end if
383
384 call messages_print_var_value('PES_Flux_Lsize', border(1:mdim), namespace=namespace)
385
386 !%Variable PES_Flux_Face_Dens
387 !%Type block
388 !%Section Time-Dependent::PhotoElectronSpectrum
389 !%Description
390 !% Define the number of points density per unit of area (in au) on the
391 !% face of the 'cub' surface.
392 !%End
393 if (parse_is_defined(namespace, 'PES_Flux_Face_Dens')) then
394 this%surf_interp = .true.
395 call parse_variable(namespace, 'PES_Flux_Face_Dens', m_one, fc_ptdens)
396 call messages_print_var_value('PES_Flux_Face_Dens', fc_ptdens, namespace=namespace)
397 end if
398
399 else
400
401 this%surf_interp = .true.
402
403 !%Variable PES_Flux_Radius
404 !%Type float
405 !%Section Time-Dependent::PhotoElectronSpectrum
406 !%Description
407 !% The radius of the sphere, if <tt>PES_Flux_Shape == sph</tt>.
408 !%End
409 if (parse_is_defined(namespace, 'PES_Flux_Radius')) then
410 call parse_variable(namespace, 'PES_Flux_Radius', m_zero, this%radius)
411 if (this%radius <= m_zero) call messages_input_error(namespace, 'PES_Flux_Radius')
412 call messages_print_var_value('PES_Flux_Radius', this%radius, namespace=namespace)
413 else
414 select type (box => mesh%box)
415 type is (box_sphere_t)
416 this%radius = box%radius
417 type is (box_parallelepiped_t)
418 this%radius = minval(box%half_length(1:mdim))
419 class default
420 message(1) = 'PES_Flux_Radius not specified. No default values available for this box shape.'
421 message(2) = 'Specify the radius of the sphere with variable PES_Flux_Radius.'
422 call messages_fatal(2, namespace=namespace)
423 end select
424 message(1) = 'PES_Flux_Radius not specified. Using default values.'
425 call messages_info(1, namespace=namespace)
426 call messages_print_var_value('PES_Flux_Radius', this%radius, namespace=namespace)
427 end if
428
429 !%Variable PES_Flux_StepsThetaR
430 !%Type integer
431 !%Default: 2 <tt>PES_Flux_Lmax</tt> + 1
432 !%Section Time-Dependent::PhotoElectronSpectrum
433 !%Description
434 !% Number of steps in <math>\theta</math> (<math>0 \le \theta \le \pi</math>) for the spherical surface.
435 !%End
436 call parse_variable(namespace, 'PES_Flux_StepsThetaR', 2 * this%lmax + 1, nstepsthetar)
437 if (nstepsthetar < 0) call messages_input_error(namespace, 'PES_Flux_StepsThetaR')
438 call messages_print_var_value("PES_Flux_StepsThetaR", nstepsthetar, namespace=namespace)
439
440 !%Variable PES_Flux_StepsPhiR
441 !%Type integer
442 !%Default: 2 <tt>PES_Flux_Lmax</tt> + 1
443 !%Section Time-Dependent::PhotoElectronSpectrum
444 !%Description
445 !% Number of steps in <math>\phi</math> (<math>0 \le \phi \le 2 \pi</math>) for the spherical surface.
446 !%End
447 call parse_variable(namespace, 'PES_Flux_StepsPhiR', 2 * this%lmax + 1, nstepsphir)
448 if (nstepsphir < 0) call messages_input_error(namespace, 'PES_Flux_StepsPhiR')
449 call messages_print_var_value("PES_Flux_StepsPhiR", nstepsphir, namespace=namespace)
450
451 end if
452
453
454 if (this%surf_interp) call mesh_interpolation_init(this%interp, mesh)
455
456 ! -----------------------------------------------------------------
457 ! Get the surface points
458 ! -----------------------------------------------------------------
459 if (this%surf_shape == pes_cubic .or. this%surf_shape == pes_plane) then
460 call pes_flux_getcube(this, mesh, abb, border, offset, fc_ptdens)
461 else
462 ! equispaced grid in theta & phi (Gauss-Legendre would optimize to nstepsthetar = this%lmax & nstepsphir = 2*this%lmax + 1):
463 ! nstepsthetar = M_TWO * this%lmax + 1
464 ! nstepsphir = M_TWO * this%lmax + 1
465 call pes_flux_getsphere(this, mesh, nstepsthetar, nstepsphir)
466 end if
467
468
469 !%Variable PES_Flux_Parallelization
470 !%Type flag
471 !%Section Time-Dependent::PhotoElectronSpectrum
472 !%Description
473 !% The parallelization strategy to be used to calculate the PES spectrum
474 !% using the resources available in the domain parallelization pool.
475 !% This option is not available without domain parallelization.
476 !% Parallelization over k-points and states is always enabled when available.
477 !%Option pf_none bit(1)
478 !% No parallelization.
479 !%Option pf_time bit(2)
480 !% Parallelize time integration. This requires to store some quantities over a
481 !% number of time steps equal to the number of cores available.
482 !%Option pf_momentum bit(3)
483 !% Parallelize over the final momentum grid. This strategy has a much lower
484 !% memory footprint than the one above (time) but seems to provide a smaller
485 !% speedup.
486 !%Option pf_surface bit(4)
487 !% Parallelize over surface points.
488 !%
489 !%
490 !% Option pf_time and pf_surface can be combined: pf_time + pf_surface.
491 !%
492 !%End
493
494 this%par_strategy = option__pes_flux_parallelization__pf_none
495 if (mesh%parallel_in_domains) then
496
497 if (this%surf_shape == pes_spherical) then
498 this%par_strategy = option__pes_flux_parallelization__pf_time &
499 + option__pes_flux_parallelization__pf_surface
500 else
501 this%par_strategy = option__pes_flux_parallelization__pf_surface
502 if (space%dim == 1) this%par_strategy = option__pes_flux_parallelization__pf_time
503 end if
504 par_strategy = this%par_strategy
505 call parse_variable(namespace, 'PES_Flux_Parallelization', par_strategy, this%par_strategy)
506 if (.not. varinfo_valid_option('PES_Flux_Parallelization', this%par_strategy, is_flag = .true.)) then
507 call messages_input_error(namespace,'PES_Flux_Parallelization')
508 end if
509
510 end if
511
512 !Sanity check
513 if (bitand(this%par_strategy, option__pes_flux_parallelization__pf_surface) /= 0 .and. &
514 bitand(this%par_strategy, option__pes_flux_parallelization__pf_momentum) /= 0) then
515 call messages_input_error(namespace,'PES_Flux_Parallelization', "Cannot combine pf_surface and pf_momentum")
516 end if
517
518 write(message(1),'(a)') 'Input: [PES_Flux_Parallelization = '
519 if (bitand(this%par_strategy, option__pes_flux_parallelization__pf_none) /= 0) then
520 write(message(1),'(a,1x,a)') trim(message(1)), 'pf_none '
521 end if
522 if (bitand(this%par_strategy, option__pes_flux_parallelization__pf_time) /= 0) then
523 write(message(1),'(a,1x,a)') trim(message(1)), 'pf_time '
524 end if
525 if (bitand(this%par_strategy, option__pes_flux_parallelization__pf_momentum) /= 0) then
526 write(message(1),'(a,1x,a)') trim(message(1)), 'pf_momentum'
527 end if
528 if (bitand(this%par_strategy, option__pes_flux_parallelization__pf_surface) /= 0) then
529 write(message(1),'(a,1x,a)') trim(message(1)), 'pf_surface'
530 end if
531 write(message(1),'(2a)') trim(message(1)), ']'
532 call messages_info(1, namespace=namespace)
533
534
535
536 ! distribute the surface points on nodes,
537 ! since mesh domains may have different numbers of surface points.
538 if (bitand(this%par_strategy, option__pes_flux_parallelization__pf_surface) /= 0) then
539
540 call pes_flux_distribute(1, this%nsrfcpnts, this%nsrfcpnts_start, this%nsrfcpnts_end, mesh%mpi_grp%comm)
541
542 if (this%surf_shape /= pes_spherical) call pes_flux_distribute_facepnts_cub(this, mesh)
543
544! Keep this because is useful for debug but not enough to bother having it polished out
545! if (debug%info) then
546! call mpi_world%barrier()
547! write(*,*) &
548! 'Debug: surface points on node ', mpi_world%rank, ' : ', this%nsrfcpnts_start, this%nsrfcpnts_end
549! call mpi_world%barrier()
550! end if
551
552 else
553 this%nsrfcpnts_start = 1
554 this%nsrfcpnts_end = this%nsrfcpnts
555
556 end if
557
558! Keep this because is useful for debug but not enough to bother having it polished out
559! if (debug%info .and. mpi_grp_is_root(mpi_world)) then
560! do isp = 1, this%nsrfcpnts
561! write(223,*) isp, this%rcoords(:, isp), this%srfcnrml(:,isp)
562! end do
563! if (this%nsrfcpnts > 0) flush(223)
564! end if
565
566 ! Generate the momentum space mesh grid
567 call pes_flux_reciprocal_mesh_gen(this, namespace, space, st, kpoints, mesh%mpi_grp%comm)
568
569
570 !%Variable PES_Flux_UseSymmetries
571 !%Type logical
572 !%Section Time-Dependent::PhotoElectronSpectrum
573 !%Description
574 !% Use surface and momentum grid symmetries to speed up calculation and
575 !% lower memory footprint.
576 !% By default available only when the surface shape matches the grid symmetry i.e.:
577 !% PES_Flux_Shape = m_cub or m_pln and PES_Flux_Momenutum_Grid = m_cartesian
578 !% or
579 !% PES_Flux_Shape = m_sph and PES_Flux_Momenutum_Grid = m_polar
580 !%End
581 use_symmetry = .false.
582 if ((this%surf_shape == pes_cubic .or. this%surf_shape == pes_plane) &
583 .and. this%kgrid == pes_cartesian .and. mdim == 3) use_symmetry = .true.
584 if (this%surf_shape == pes_spherical .and. this%kgrid == pes_polar) use_symmetry = .true.
585 call parse_variable(namespace, 'PES_Flux_UseSymmetries', use_symmetry, this%use_symmetry)
586 call messages_print_var_value('PES_Flux_UseSymmetries', this%use_symmetry, namespace=namespace)
587
588
589
590 ! get the values of the spherical harmonics for the surface points for PES_SPHERICAL
591 if (this%surf_shape == pes_spherical) then
592 safe_allocate(this%ylm_r(this%nsrfcpnts_start:this%nsrfcpnts_end, 0:this%lmax, -this%lmax:this%lmax))
593 this%ylm_r = m_z0
594
595 if (this%nsrfcpnts_start > 0) then
596 do ll = 0, this%lmax
597 do mm = -ll, ll
598 do isp = this%nsrfcpnts_start, this%nsrfcpnts_end
599 call ylmr_cmplx(this%rcoords(1:3, isp), ll, mm, this%ylm_r(isp, ll, mm))
600 this%ylm_r(isp, ll, mm) = conjg(this%ylm_r(isp, ll, mm))
601 end do
602 end do
603 end do
604 end if
605
606 else
607
608 call pes_flux_integrate_cub_tabulate(this, space, mesh, st, kpoints)
609
610 end if
611
612
613 !%Variable PES_Flux_RuntimeOutput
614 !%Type logical
615 !%Section Time-Dependent::PhotoElectronSpectrum
616 !%Description
617 !% Write output in ascii format at runtime.
618 !%
619 !%End
620 call parse_variable(namespace, 'PES_Flux_RuntimeOutput', .false., this%runtime_output)
621 call messages_print_var_value('PES_Flux_RuntimeOutput', this%runtime_output, namespace=namespace)
622
623
624
625 ! -----------------------------------------------------------------
626 ! Options for time integration
627 ! -----------------------------------------------------------------
628 this%max_iter = max_iter
629 this%save_iter = save_iter
630 this%itstep = 0
631 this%tdsteps = 1
632
633 if (mesh%parallel_in_domains .and. bitand(this%par_strategy, option__pes_flux_parallelization__pf_time) /= 0) then
634 this%tdsteps = mesh%mpi_grp%size
635 end if
636
637 ! -----------------------------------------------------------------
638 ! Allocations
639 ! -----------------------------------------------------------------
640
641 safe_allocate(this%wf(0:this%nsrfcpnts, stst:stend, 1:sdim, kptst:kptend, 1:this%tdsteps))
642 this%wf = m_z0
643
644 safe_allocate(this%gwf(0:this%nsrfcpnts, stst:stend, 1:sdim, kptst:kptend, 1:this%tdsteps, 1:mdim))
645 this%gwf = m_z0
646
647 safe_allocate(this%veca(1:mdim, 1:this%tdsteps))
648 this%veca = m_zero
649
650 if (this%surf_shape == pes_spherical) then
651 safe_allocate(this%spctramp_sph(1:this%nstepsomegak, stst:stend, 1:sdim, kptst:kptend, 1:this%nk))
652 this%spctramp_sph = m_z0
653
654 safe_allocate(this%conjgphase_prev(1:this%nstepsomegak, 1:this%nk))
655
656 else
657 safe_allocate(this%spctramp_cub(stst:stend, 1:sdim, kptst:kptend, 1:this%nkpnts))
658 this%spctramp_cub = m_z0
659
660 safe_allocate(this%conjgphase_prev(1:this%nkpnts, kptst:kptend))
661
662 end if
663
664 this%conjgphase_prev = m_z1
665
666 call messages_write('Info: Total number of surface points = ')
667 call messages_write(this%nsrfcpnts)
668 call messages_info(namespace=namespace)
669
670 call messages_write('Info: Total number of momentum points = ')
671 call messages_write(this%nkpnts)
672 call messages_info(namespace=namespace)
673
674 pop_sub_with_profile(pes_flux_init)
675 end subroutine pes_flux_init
676
677 ! ---------------------------------------------------------
678 subroutine pes_flux_end(this)
679 type(pes_flux_t), intent(inout) :: this
680
681 push_sub(pes_flux_end)
682
683 if (this%surf_shape == pes_spherical) then
684 safe_deallocate_a(this%kcoords_sph)
685 safe_deallocate_a(this%ylm_k)
686 safe_deallocate_a(this%j_l)
687 safe_deallocate_a(this%ylm_r)
688 safe_deallocate_a(this%conjgphase_prev)
689 safe_deallocate_a(this%spctramp_sph)
690 else
691 safe_deallocate_a(this%kcoords_cub)
692 safe_deallocate_a(this%conjgphase_prev)
693 safe_deallocate_a(this%spctramp_cub)
694
695 if (.not. this%surf_interp) then
696 safe_deallocate_a(this%srfcpnt)
697 end if
698 safe_deallocate_a(this%rankmin)
699
700 safe_deallocate_a(this%face_idx_range)
701 safe_deallocate_a(this%LLr)
702 safe_deallocate_a(this%NN)
703
704 safe_deallocate_a(this%expkr)
705 safe_deallocate_a(this%expkr_perp)
706 safe_deallocate_a(this%bvk_phase)
707
708 safe_deallocate_a(this%Lkpuvz_inv)
709
710 end if
711
712 safe_deallocate_a(this%klinear)
713
714 safe_deallocate_a(this%srfcnrml)
715 safe_deallocate_a(this%rcoords)
716
717 safe_deallocate_a(this%wf)
718 safe_deallocate_a(this%gwf)
719 safe_deallocate_a(this%veca)
720
721 pop_sub(pes_flux_end)
722 end subroutine pes_flux_end
723
724 ! ---------------------------------------------------------
725 subroutine pes_flux_reciprocal_mesh_gen(this, namespace, space, st, kpoints, comm, post)
726 type(pes_flux_t), intent(inout) :: this
727 type(namespace_t), intent(in) :: namespace
728 class(space_t), intent(in) :: space
729 type(states_elec_t), intent(in) :: st
730 type(kpoints_t), intent(in) :: kpoints
731 type(mpi_comm), intent(in) :: comm
732 logical, optional, intent(in) :: post
733
734 integer :: mdim, pdim
735 integer :: kptst, kptend
736 integer :: ikp, ikpt
737 integer :: ll, mm, idim, idir, ispin, nspin
738 integer :: ikk, ith, iph, iomk,ie, ik1, ik2, ik3, kgrid_block_dim
739 real(real64) :: kmax(space%dim), kmin(space%dim), kact, thetak, phik
740 type(block_t) :: blk
741
742 real(real64) :: emin, emax, de , kvec(1:3), xx(3)
743 integer :: nkp_out, nkmin, nkmax
744
745 integer :: kgrid
746 real(real64), allocatable :: gpoints(:,:), gpoints_reduced(:,:)
747 real(real64) :: dk(1:3), kpoint(1:3), dthetak, dphik
748 logical :: use_enegy_grid, arpes_grid
749
751
752 kptst = st%d%kpt%start
753 kptend = st%d%kpt%end
754 mdim = space%dim
755 pdim = space%periodic_dim
756
757 this%dim = mdim
758 this%pdim = pdim
759
760 !%Variable PES_Flux_Momenutum_Grid
761 !%Type integer
762 !%Section Time-Dependent::PhotoElectronSpectrum
763 !%Description
764 !% Decides how the grid in momentum space is generated.
765 !%Option polar 1
766 !% The grid is in polar coordinates with the zenith axis is along z.
767 !% The grid parameters are defined by PES_Flux_Kmax, PES_Flux_DeltaK,
768 !% PES_Flux_StepsThetaK, PES_Flux_StepsPhiK.
769 !% This is the default choice for PES_Flux_Shape = sph or cub.
770 !%Option cartesian 2
771 !% The grid is in cartesian coordinates with parameters defined by
772 !% PES_Flux_ARPES_grid, PES_Flux_EnergyGrid.
773 !% This is the default choice for PES_Flux_Shape = sph or cub.
774 !%End
775
776 ! default values
777 select case (this%surf_shape)
778 case (pes_spherical)
779 kgrid = pes_polar
780 case (pes_plane)
781 kgrid = pes_cartesian
782 case (pes_cubic)
783 kgrid = pes_cartesian
784 if (mdim == 1) kgrid = pes_polar
785
786 case default
787 assert(.false.)
788
789 end select
790
791 call parse_variable(namespace, 'PES_Flux_Momenutum_Grid', kgrid, this%kgrid)
792 if (.not. varinfo_valid_option('PES_Flux_Momenutum_Grid', this%kgrid, is_flag = .true.)) then
793 call messages_input_error(namespace,'PES_Flux_Momenutum_Grid')
794 end if
795 call messages_print_var_option('PES_Flux_Momenutum_Grid', this%kgrid, namespace=namespace)
796
797
798
799 !Check availability of the calculation requested
800 if (this%surf_shape == pes_spherical) then
801 if (this%kgrid == pes_cartesian) then
802 call messages_not_implemented('Cartesian momentum grid with a spherical surface', namespace=namespace)
803 end if
804 if (space%is_periodic()) then
805 call messages_not_implemented('Spherical surface flux for periodic systems', namespace=namespace)
806 end if
807 if (mdim == 1) then
808 call messages_not_implemented('Spherical surface flux for one-dimensional systems', namespace=namespace)
809 end if
810 end if
811
812 if (this%surf_shape == pes_cubic) then
813 if (space%is_periodic()) then
814 call messages_not_implemented('Use of cubic surface for periodic systems (use pln)', namespace=namespace)
815 end if
816 end if
817
819
820 !%Variable PES_Flux_EnergyGrid
821 !%Type block
822 !%Section Time-Dependent::PhotoElectronSpectrum
823 !%Description
824 !% The block <tt>PES_Flux_EnergyGrid</tt> specifies the energy grid
825 !% in momentum space.
826 !% <tt><br>%PES_Flux_EnergyGrid
827 !% <br>&nbsp;&nbsp; Emin | Emax | DeltaE
828 !% <br>%</tt>
829 !%End
830
831 emin = 0
832 emax = 10
833 de = 0.1_real64
834 use_enegy_grid = .false.
835
836 if (parse_block(namespace, 'PES_Flux_EnergyGrid', blk) == 0) then
837
838 call parse_block_float(blk, 0, 0, emin)
839 call parse_block_float(blk, 0, 1, emax)
840 call parse_block_float(blk, 0, 2, de)
841
842 call parse_block_end(blk)
843 use_enegy_grid = .true.
844
845
846 call messages_write("Energy grid (Emin, Emax, DE) [")
847 call messages_write(trim(units_abbrev(units_out%energy)))
848 call messages_write("]: (")
849 call messages_write(units_from_atomic(units_out%energy, emin),fmt = '(f8.3)')
850 call messages_write(", ")
851 call messages_write(units_from_atomic(units_out%energy, emax), fmt = '(f8.3)')
852 call messages_write(", ")
853 call messages_write(units_from_atomic(units_out%energy, de), fmt = '(e9.2)')
854 call messages_write(")")
855 call messages_info(namespace=namespace)
856
857 kmax(1:mdim) = sqrt(m_two*emax)
858 kmin(1:mdim) = sqrt(m_two*emin)
859 this%dk = sqrt(m_two*de)
860
861 end if
862
863 kgrid_block_dim = 1
864 !ugly hack (since this input variable is properly handled below) but effective
865 call parse_variable(namespace, 'PES_Flux_ARPES_grid', .false., this%arpes_grid)
866 if (.not. use_enegy_grid .or. this%arpes_grid) then
867
868 !%Variable PES_Flux_Kmax
869 !%Type float
870 !%Default 1.0
871 !%Section Time-Dependent::PhotoElectronSpectrum
872 !%Description
873 !% The maximum value of the photoelectron momentum.
874 !% For cartesian momentum grids one can specify a value different
875 !% for cartesian direction using a block input.
876 !%End
877 if (parse_block(namespace, 'PES_Flux_Kmax', blk) == 0) then
878 if (this%kgrid == pes_cartesian) then
879 do idim = 1, mdim
880 call parse_block_float(blk, 0, idim - 1, kmax(idim))
881 end do
882 kgrid_block_dim = mdim
883 else
884 message(1) = 'Wrong block format for PES_Flux_Kmax and non-cartesian grid'
885 call messages_fatal(1, namespace=namespace)
886 end if
887 else
888 call parse_variable(namespace, 'PES_Flux_Kmax', m_one, kmax(1))
889 kmax(:)=kmax(1)
890 end if
891
892 !%Variable PES_Flux_Kmin
893 !%Type float
894 !%Default 0.0
895 !%Section Time-Dependent::PhotoElectronSpectrum
896 !%Description
897 !% The minimum value of the photoelectron momentum.
898 !% For cartesian momentum grids one can specify a value different
899 !% for cartesian direction using a block input.
900 !%End
901 if (parse_block(namespace, 'PES_Flux_Kmin', blk) == 0) then
902 if (this%kgrid == pes_cartesian) then
903 do idim = 1, mdim
904 call parse_block_float(blk, 0, idim - 1, kmin(idim))
905 end do
906 kgrid_block_dim = mdim
907 else
908 message(1) = 'Wrong block format for PES_Flux_Kmin and non-cartesian grid'
909 call messages_fatal(1, namespace=namespace)
910 end if
911 else
912 call parse_variable(namespace, 'PES_Flux_Kmin', m_zero, kmin(1))
913 kmin(:)=kmin(1)
914 end if
915
916
917 !%Variable PES_Flux_DeltaK
918 !%Type float
919 !%Default 0.02
920 !%Section Time-Dependent::PhotoElectronSpectrum
921 !%Description
922 !% Spacing of the the photoelectron momentum grid.
923 !%End
924 call parse_variable(namespace, 'PES_Flux_DeltaK', 0.02_real64, this%dk)
925 if (this%dk <= m_zero) call messages_input_error(namespace,'PES_Flux_DeltaK')
926
927 end if
928
929 do idim = 1, kgrid_block_dim
930 if (kgrid_block_dim == 1) then
931 call messages_write("Momentum linear grid (Pmin, Pmax, DP) [")
932 else
933 call messages_write("Momentum linear grid (Pmin, Pmax, DP) "//index2axis(idim)//"-axis [")
934 end if
935 call messages_write(trim(units_abbrev(units_out%mass*units_out%velocity)))
936 call messages_write("]: (")
937 call messages_write(kmin(idim),fmt = '(f8.3)')
938 call messages_write(", ")
939 call messages_write(kmax(idim), fmt = '(f8.3)')
940 call messages_write(", ")
941 call messages_write(this%dk, fmt = '(e9.2)')
942 call messages_write(")")
943 call messages_info(namespace=namespace)
944 end do
945
946
947
948! if (this%surf_shape == PES_SPHERICAL .or. this%surf_shape == PES_CUBIC) then
949 if (this%kgrid == pes_polar) then
950 ! -----------------------------------------------------------------
951 ! Setting up k-mesh
952 ! 1D = 2 points, 2D = polar coordinates, 3D = spherical coordinates
953 ! -----------------------------------------------------------------
954
955
956 !%Variable PES_Flux_ThetaK
957 !%Type block
958 !%Section Time-Dependent::PhotoElectronSpectrum
959 !%Description
960 !% Define the grid points on theta (<math>0 \le \theta \le \pi</math>) when
961 !% using a spherical grid in momentum.
962 !% The block defines the maximum and minimum values of theta and the number of
963 !% of points for the discretization.
964 !%
965 !% <tt>%PES_Flux_ThetaK
966 !% <br>&nbsp;&nbsp; theta_min | theta_max | npoints
967 !% <br>%
968 !% </tt>
969 !%
970 !% By default theta_min=0, theta_max = pi, npoints = 45.
971 !%End
972 this%nstepsthetak = 45
973 this%thetak_rng(1:2) = (/m_zero, m_pi/)
974 if (parse_block(namespace, 'PES_Flux_ThetaK', blk) == 0) then
975 call parse_block_float(blk, 0, 0, this%thetak_rng(1))
976 call parse_block_float(blk, 0, 1, this%thetak_rng(2))
977 call parse_block_integer(blk, 0, 2, this%nstepsthetak)
978 call parse_block_end(blk)
979 do idim = 1,2
980 if (this%thetak_rng(idim) < m_zero .or. this%thetak_rng(idim) > m_pi) then
981 call messages_input_error(namespace,'PES_Flux_ThetaK')
982 end if
983 end do
984 if (this%nstepsthetak < 0) call messages_input_error(namespace,'PES_Flux_ThetaK')
985 else
986
987 !%Variable PES_Flux_StepsThetaK
988 !%Type integer
989 !%Default 45
990 !%Section Time-Dependent::PhotoElectronSpectrum
991 !%Description
992 !% Number of steps in <math>\theta</math> (<math>0 \le \theta \le \pi</math>) for the spherical grid in k.
993 !%End
994 call parse_variable(namespace, 'PES_Flux_StepsThetaK', 45, this%nstepsthetak)
995 if (this%nstepsthetak < 0) call messages_input_error(namespace,'PES_Flux_StepsThetaK')
996
997 ! should do this at some point
998 ! call messages_obsolete_variable(namespace, 'PES_Flux_StepsThetaK')
999 end if
1000
1001
1002
1003 !%Variable PES_Flux_PhiK
1004 !%Type block
1005 !%Section Time-Dependent::PhotoElectronSpectrum
1006 !%Description
1007 !% Define the grid points on theta (<math>0 \le \theta \le 2\pi</math>) when
1008 !% using a spherical grid in momentum.
1009 !% The block defines the maximum and minimum values of theta and the number of
1010 !% of points for the discretization.
1011 !%
1012 !% <tt>%PES_Flux_PhiK
1013 !% <br>&nbsp;&nbsp; theta_min | theta_max | npoints
1014 !% <br>%
1015 !% </tt>
1016 !%
1017 !% By default theta_min=0, theta_max = pi, npoints = 90.
1018 !%End
1019 this%nstepsphik = 90
1020 this%phik_rng(1:2) = (/m_zero, 2 * m_pi/)
1021 if (mdim == 1) this%phik_rng(1:2) = (/m_pi, m_zero/)
1022 if (parse_block(namespace, 'PES_Flux_PhiK', blk) == 0) then
1023 call parse_block_float(blk, 0, 0, this%phik_rng(1))
1024 call parse_block_float(blk, 0, 1, this%phik_rng(2))
1025 call parse_block_integer(blk, 0, 2, this%nstepsphik)
1026 call parse_block_end(blk)
1027 do idim = 1,2
1028 if (this%phik_rng(idim) < m_zero .or. this%phik_rng(idim) > m_two * m_pi) then
1029 call messages_input_error(namespace,'PES_Flux_PhiK')
1030 end if
1031 end do
1032 if (this%nstepsphik < 0) call messages_input_error(namespace,'PES_Flux_PhiK')
1033
1034 else
1035
1036 !%Variable PES_Flux_StepsPhiK
1037 !%Type integer
1038 !%Default 90
1039 !%Section Time-Dependent::PhotoElectronSpectrum
1040 !%Description
1041 !% Number of steps in <math>\phi</math> (<math>0 \le \phi \le 2 \pi</math>) for the spherical grid in k.
1042 !%End
1043 call parse_variable(namespace, 'PES_Flux_StepsPhiK', 90, this%nstepsphik)
1044 if (this%nstepsphik < 0) call messages_input_error(namespace,'PES_Flux_StepsPhiK')
1045 if (this%nstepsphik == 0) this%nstepsphik = 1
1046 end if
1047
1048
1049 dthetak = m_zero
1050 if (mdim == 3) dthetak = abs(this%thetak_rng(2) - this%thetak_rng(1)) / (this%nstepsthetak)
1051 dphik = abs(this%phik_rng(2) - this%phik_rng(1)) / (this%nstepsphik)
1052
1053
1054 select case (mdim)
1055 case (1)
1056 this%nstepsthetak = 0
1057 this%nstepsphik = 2
1058 this%nstepsomegak = this%nstepsphik
1059
1060 case (2)
1061 this%nstepsthetak = 0
1062 this%nstepsomegak = this%nstepsphik
1063
1064 case (3)
1065 if (this%nstepsthetak <= 1) then
1066 this%nstepsphik = 1
1067 this%nstepsthetak = 1
1068 end if
1069 ! count the omegak samples
1070 iomk = 0
1071 do ith = 0, this%nstepsthetak
1072 thetak = ith * dthetak + this%thetak_rng(1)
1073 do iph = 0, this%nstepsphik - 1
1074 phik = iph * dphik + this%phik_rng(1)
1075 iomk = iomk + 1
1076 if (thetak < m_epsilon .or. abs(thetak-m_pi) < m_epsilon) exit
1077 end do
1078 end do
1079 this%nstepsomegak = iomk
1080
1081 case default
1082 assert(.false.)
1083
1084 end select
1085
1086 write(message(1),'(a)') "Polar momentum grid:"
1087 call messages_info(1, namespace=namespace)
1088 if (mdim == 3) then
1089 write(message(1),'(a,f12.6,a,f12.6,a, i6)') &
1090 " Theta = (", this%thetak_rng(1), ",",this%thetak_rng(2), &
1091 "); n = ", this%nstepsthetak
1092 call messages_info(1, namespace=namespace)
1093 end if
1094 write(message(1),'(a,f12.6,a,f12.6,a, i6)') &
1095 " Phi = (", this%phik_rng(1), ",",this%phik_rng(2), &
1096 "); n = ", this%nstepsphik
1097 call messages_info(1, namespace=namespace)
1098
1099
1100
1101 if (use_enegy_grid) then
1102 this%nk = nint(abs(emax - emin) / de)
1103 else
1104 this%nk = nint(abs(kmax(1) - kmin(1)) / this%dk)
1105 end if
1106 this%nkpnts = this%nstepsomegak * this%nk
1107
1108 this%ll(1) = this%nk
1109 this%ll(2) = this%nstepsphik
1110 this%ll(3) = this%nstepsthetak !- 1
1111 this%ll(mdim+1:3) = 1
1112
1113 safe_allocate(this%klinear(1:this%nk, 1))
1114
1115
1116 else
1117 ! Cartesian
1118
1119 dk(1:mdim) = m_one/kpoints%nik_axis(1:mdim)
1120
1121 this%arpes_grid = .false.
1122 if (space%is_periodic()) then
1123 !%Variable PES_Flux_ARPES_grid
1124 !%Type logical
1125 !%Section Time-Dependent::PhotoElectronSpectrum
1126 !%Description
1127 !% Use a curvilinear momentum space grid that compensates the transformation
1128 !% used to obtain ARPES. With this choice ARPES data is laid out on a Cartesian
1129 !% regular grid.
1130 !% By default true when <tt>PES_Flux_Shape = pln</tt> and a <tt>KPointsPath</tt>
1131 !% is specified.
1132 !%End
1133 arpes_grid = kpoints%have_zero_weight_path()
1134 call parse_variable(namespace, 'PES_Flux_ARPES_grid', arpes_grid, this%arpes_grid)
1135 call messages_print_var_value("PES_Flux_ARPES_grid", this%arpes_grid, namespace=namespace)
1136 end if
1137
1138
1139
1140
1141 this%ll(:) = 1
1142
1143 if (kpoints%have_zero_weight_path()) then
1144
1145 if (this%arpes_grid) then
1146 nkmax = floor(emax / de)
1147 nkmin = floor(emin / de)
1148
1149 else
1150 nkmax = floor(kmax(mdim) / this%dk)
1151 nkmin = floor(kmin(mdim) / this%dk)
1152
1153 end if
1154
1155 this%ll(mdim) = abs(nkmax - nkmin) + 1
1156
1157 this%nk = this%ll(mdim)
1158
1159
1160 else
1161
1162 if (.not. this%arpes_grid) then
1163 this%ll(1:mdim) = floor(abs(kmax(1:mdim) - kmin(1:mdim))/this%dk) + 1
1164 this%nk = maxval(this%ll(1:mdim))
1165
1166 else
1167
1168 nkmax = floor(emax / de)
1169 nkmin = floor(emin / de)
1170 this%nk = abs(nkmax - nkmin) + 1
1171
1172 this%ll(1:pdim) = floor(abs(kmax(1:pdim) - kmin(1:pdim))/this%dk) + 1
1173 this%ll(mdim) = this%nk
1174
1175 end if
1176
1177 safe_allocate(this%klinear(1:maxval(this%ll(1:mdim)), 1:mdim))
1178 this%klinear = m_zero
1179
1180 end if
1181
1182 ! Total number of points
1183 this%nkpnts = product(this%ll(1:mdim))
1184
1185
1186 end if
1187
1188
1189
1190
1191
1192 this%parallel_in_momentum = .false.
1193
1194 ! Create the grid
1195 select case (this%kgrid)
1196
1197 case (pes_polar)
1198
1199 if (use_enegy_grid) then
1200 do ie = 1, this%nk
1201 this%klinear(ie, 1) = sqrt(m_two * (ie * de + emin))
1202 end do
1203 else
1204 do ikk = 1, this%nk
1205 this%klinear(ikk, 1) = ikk * this%dk + kmin(1)
1206 end do
1207 end if
1208
1209
1210 if (this%surf_shape == pes_spherical) then
1211
1212 if (optional_default(post, .false.)) then
1214 return
1215 end if
1216
1217 ! we split the k-mesh in radial & angular part
1218 call pes_flux_distribute(1, this%nk, this%nk_start, this%nk_end, comm)
1219 if ((this%nk_end - this%nk_start + 1) < this%nk) this%parallel_in_momentum = .true.
1220 call pes_flux_distribute(1, this%nstepsomegak, this%nstepsomegak_start, this%nstepsomegak_end, comm)
1221
1222! Keep this because is useful for debug but not enough to bother having it polished out
1223! if (debug%info) then
1224! call mpi_world%barrier()
1225! write(*,*) &
1226! 'Debug: momentum points on node ', mpi_world%rank, ' : ', this%nk_start, this%nk_end
1227! call mpi_world%barrier()
1228! write(*,*) &
1229! 'Debug: momentum directions on node ', mpi_world%rank, ' : ', this%nstepsomegak_start, this%nstepsomegak_end
1230! call mpi_world%barrier()
1231! end if
1232 safe_allocate(this%j_l(this%nk_start:this%nk_end, 0:this%lmax))
1233 this%j_l = m_zero
1234
1235 safe_allocate(this%kcoords_sph(1:3, this%nk_start:this%nk_end, 1:this%nstepsomegak))
1236 this%kcoords_sph = m_zero
1237
1238 safe_allocate(this%ylm_k(this%nstepsomegak_start:this%nstepsomegak_end, -this%lmax:this%lmax, 0:this%lmax))
1239 this%ylm_k = m_z0
1240
1241 ! spherical harmonics & kcoords_sph
1242 iomk = 0
1243 do ith = 0, this%nstepsthetak
1244 thetak = ith * dthetak + this%thetak_rng(1)
1245 do iph = 0, this%nstepsphik - 1
1246 phik = iph * dphik + this%phik_rng(1)
1247 iomk = iomk + 1
1248 if (iomk >= this%nstepsomegak_start .and. iomk <= this%nstepsomegak_end) then
1249 xx(1) = cos(phik) * sin(thetak)
1250 xx(2) = sin(phik) * sin(thetak)
1251 xx(3) = cos(thetak)
1252 !$omp parallel do private(ll, mm)
1253 do ll = 0, this%lmax
1254 do mm = -ll, ll
1255 call ylmr_cmplx(xx, ll, mm, this%ylm_k(iomk, mm, ll))
1256 end do
1257 end do
1258 end if
1259 if (this%nk_start > 0) then
1260 this%kcoords_sph(1, this%nk_start:this%nk_end, iomk) = cos(phik) * sin(thetak)
1261 this%kcoords_sph(2, this%nk_start:this%nk_end, iomk) = sin(phik) * sin(thetak)
1262 this%kcoords_sph(3, this%nk_start:this%nk_end, iomk) = cos(thetak)
1263 end if
1264 if (thetak < m_epsilon .or. abs(thetak-m_pi) < m_epsilon) exit
1265 end do
1266 end do
1267
1268 if (this%nk_start > 0) then
1269 ! Bessel functions & kcoords_sph
1270 do ikk = this%nk_start, this%nk_end
1271 kact = this%klinear(ikk,1)
1272 do ll = 0, this%lmax
1273 this%j_l(ikk, ll) = loct_sph_bessel(ll, kact * this%radius) * &
1274 m_two * m_pi / (m_two * m_pi)**m_three / m_two
1275 end do
1276 this%kcoords_sph(:, ikk, :) = kact * this%kcoords_sph(:, ikk, :)
1277 end do
1278 end if
1279
1280 else
1281 !planar or cubic surface
1282
1283 ! no distribution
1284 this%nkpnts_start = 1
1285 this%nkpnts_end = this%nkpnts
1286
1287
1288 safe_allocate(this%kcoords_cub(1:mdim, this%nkpnts_start:this%nkpnts_end, 1))
1289
1290 this%kcoords_cub = m_zero
1291
1292 if (mdim == 1) then
1293
1294 ikp = 0
1295 do ikk = -this%nk, this%nk
1296 if (ikk == 0) cycle
1297 ikp = ikp + 1
1298 kact = sign(1,ikk) * this%klinear(abs(ikk), 1)
1299 this%kcoords_cub(1, ikp, 1) = kact
1300 end do
1301
1302 else !mdim = 2,3
1303 thetak = m_pi / m_two
1304
1305 ikp = 0
1306 do ikk = 1, this%nk
1307 do ith = 0, this%nstepsthetak
1308 if (mdim == 3) thetak = ith * dthetak + this%thetak_rng(1)
1309 do iph = 0, this%nstepsphik - 1
1310 ikp = ikp + 1
1311 phik = iph * dphik + this%phik_rng(1)
1312 kact = this%klinear(ikk,1)
1313 this%kcoords_cub(1, ikp,1) = kact * cos(phik) * sin(thetak)
1314 this%kcoords_cub(2, ikp,1) = kact * sin(phik) * sin(thetak)
1315 if (mdim == 3) this%kcoords_cub(3, ikp,1) = kact * cos(thetak)
1316 if (mdim == 3 .and. (thetak < m_epsilon .or. abs(thetak-m_pi) < m_epsilon)) exit
1317 end do
1318 end do
1319 end do
1320
1321 end if
1322 end if
1323
1324 case (pes_cartesian)
1325
1326 this%nkpnts_start = 1
1327 this%nkpnts_end = this%nkpnts
1328
1329 if (kpoints%have_zero_weight_path()) then
1330 ! its a special case because we are generating a different (1D) grid for
1331 ! each kpoint and then we combine it in postprocessing
1332
1333 ! Note that kptst:kptend includes spin
1334 safe_allocate(this%kcoords_cub(1:mdim, this%nkpnts_start:this%nkpnts_end, kptst:kptend))
1335 this%kcoords_cub = m_zero
1336
1337 nkp_out = 0
1338
1339 ! For the below code to work, we need to reinitialize ikp for each spin
1340 nspin = 1
1341 if (st%d%ispin == spin_polarized) nspin = 2
1342
1343 do ispin = 1, nspin
1344 do ikpt = (kptst-1)+ispin, kptend, nspin
1345 ikp = 0
1346 do ikk = nkmin, nkmax
1347 kvec(1:mdim) = - kpoints%get_point(st%d%get_kpoint_index(ikpt))
1349 end do
1350 end do
1351 end do
1352
1353 else
1354
1355 safe_allocate(this%kcoords_cub(1:mdim, this%nkpnts_start:this%nkpnts_end, 1))
1356
1357 safe_allocate(this%Lkpuvz_inv(1:this%ll(1), 1:this%ll(2), 1:this%ll(3)))
1358
1359 do idir = 1, mdim
1360 do ikk = 1, this%ll(idir)
1361 this%klinear(ikk, idir) = ikk * this%dk + kmin(idir)
1362 end do
1363 end do
1364
1365
1366 if (.not. this%arpes_grid) then
1367 ! Normal velocity map
1368
1369
1370 ikpt = 1
1371 ikp = 0
1372 kvec(:) = m_zero
1373 do ik3 = 1, this%ll(3)
1374 if (mdim>2) kvec(3) = this%klinear(ik3, 3)
1375 do ik2 = 1, this%ll(2)
1376 if (mdim>1) kvec(2) = this%klinear(ik2, 2)
1377 do ik1 = 1, this%ll(1)
1378 ikp = ikp + 1
1379 kvec(1) = this%klinear(ik1, 1)
1380 this%kcoords_cub(1:mdim, ikp, ikpt) = kvec(1:mdim)
1381 this%Lkpuvz_inv(ik1, ik2, ik3) = ikp
1382
1383 end do
1384 end do
1385 end do
1386
1387 else ! we want an ARPES-friendly grid layout
1388
1389 nkp_out = 0
1390 ikpt = 1
1391 ikp = 0
1392 kvec(:) = m_zero
1393 do ikk = nkmin, nkmax ! this is going to be turned into energy
1394
1395 ! loop over periodic directions
1396 select case (pdim)
1397 case (1)
1398
1399 do ik1 = 1, this%ll(1)
1400 kvec(1) = this%klinear(ik1,1)
1401 kvec(1:pdim) = matmul(kpoints%latt%klattice_primitive(1:pdim, 1:pdim), kvec(1:pdim))
1403 end do
1404
1405 case (2)
1406
1407 do ik2 = 1, this%ll(2)
1408 do ik1 = 1, this%ll(1)
1409 kvec(1:2) = (/this%klinear(ik1,1), this%klinear(ik2,2)/)
1410 kvec(1:pdim) = matmul(kpoints%latt%klattice_primitive(1:pdim, 1:pdim), kvec(1:pdim))
1412
1413 this%Lkpuvz_inv(ik1, ik2, ikk - nkmin + 1) = ikp
1414
1415
1416 end do
1417 end do
1418
1419 case default
1420 assert(.false.)
1421
1422 end select
1423
1424 end do
1425
1426
1427 end if
1428
1429 end if
1430
1431! Keep this because is useful for debug but not enough to bother having it polished out
1432! if (debug%info .and. mpi_grp_is_root(mpi_world)) then
1433! ! this does not work for parallel in kpoint
1434! ! you need to gather kcoords_pln
1435! if (kpoints_have_zero_weight_path(sb%kpoints)) then
1436! write(229,*) "# ikpt (kpoint index), ikp (momentum index), this%kcoords_cub(1:mdim, ikp, ikpt)"
1437! do ikpt = kptst, kptend
1438! do ikp = 1, this%nkpnts
1439! write(229,*) ikpt, ikp, this%kcoords_cub(1:mdim, ikp, ikpt)
1440! end do
1441! end do
1442! else
1443! write(229,*) "# ikp (momentum index), this%kcoords_cub(1:mdim, ikp, ikpt)"
1444! do ikp = 1, this%nkpnts
1445! write(229,*) ikp, this%kcoords_cub(1:mdim, ikp, 1)
1446! end do
1447! end if
1448! flush(229)
1449!
1450!
1451! if (mdim == 3) then
1452! write(230,*) "# ik1, ik2, ik3, this%Lkpuvz_inv(ik1,ik2,ik3) "
1453! do ik1 = 1, this%ll(1)
1454! do ik2 = 1, this%ll(2)
1455! do ik3 = 1, this%ll(3)
1456! write(230,*) ik1, ik2, ik3, this%Lkpuvz_inv(ik1,ik2,ik3)
1457! end do
1458! end do
1459! end do
1460!
1461! flush(230)
1462! end if
1463!
1464! end if
1465
1466 if (this%arpes_grid) then
1467 call messages_write("Number of points with E<p//^2/2 = ")
1468 call messages_write(nkp_out)
1469 call messages_write(" [of ")
1470 call messages_write(this%nkpnts*kpoints_number(kpoints))
1471 call messages_write("]")
1472 call messages_info(namespace=namespace)
1473 end if
1474
1475 safe_deallocate_a(gpoints)
1476 safe_deallocate_a(gpoints_reduced)
1477
1478 case default
1479 assert(.false.)
1480
1481 end select
1482
1483
1484 !%Variable PES_Flux_GridTransformMatrix
1485 !%Type block
1486 !%Section Time-Dependent::PhotoElectronSpectrum
1487 !%Description
1488 !% Define an optional transformation matrix for the momentum grid.
1489 !%
1490 !% <tt>%PES_Flux_GridTransformMatrix
1491 !% <br>&nbsp;&nbsp; M_11 | M_12 | M_13
1492 !% <br>&nbsp;&nbsp; M_21 | M_22 | M_23
1493 !% <br>&nbsp;&nbsp; M_31 | M_32 | M_33
1494 !% <br>%
1495 !% </tt>
1496 !%End
1497 this%ktransf(:,:) = m_zero
1498 do idim = 1,mdim
1499 this%ktransf(idim,idim) = m_one
1500 end do
1501
1502 if (parse_block(namespace, 'PES_Flux_GridTransformMatrix', blk) == 0) then
1503 do idim = 1,mdim
1504 do idir = 1, mdim
1505 call parse_block_float(blk, idir - 1, idim - 1, this%ktransf(idim, idir))
1506 end do
1507 end do
1508 call parse_block_end(blk)
1509
1510 write(message(1),'(a)') 'Momentum grid transformation matrix :'
1511 do idir = 1, space%dim
1512 write(message(1 + idir),'(9f12.6)') ( this%ktransf(idim, idir), idim = 1, mdim)
1513 end do
1514 call messages_info(1 + mdim, namespace=namespace)
1515
1516
1517 !Apply the transformation
1518 if (this%surf_shape == pes_spherical) then
1519
1520 iomk = 0
1521 do ith = 0, this%nstepsthetak
1522 do iph = 0, this%nstepsphik - 1
1523 iomk = iomk + 1
1524 do ikk = this%nk_start, this%nk_end
1525 kvec(1:mdim) = this%kcoords_sph(1:mdim, ikk, iomk)
1526 this%kcoords_sph(1:mdim, ikk, iomk) = matmul(this%ktransf(1:mdim, 1:mdim), kvec(1:mdim))
1527 end do
1528 if (ith == 0 .or. ith == this%nstepsthetak) exit
1529 end do
1530 end do
1531
1532 else !planar or cubic surface
1533
1534 do ikpt = kptst, kptend + 1
1535 if (ikpt == kptend + 1) then
1536 kpoint(1:space%dim) = m_zero
1537 else
1538 kpoint(1:space%dim) = kpoints%get_point(st%d%get_kpoint_index(ikpt))
1539 end if
1540
1541 do ikp = 1, this%nkpnts
1542
1543 kvec(1:mdim) = this%kcoords_cub(1:mdim, ikp, ikpt) - kpoint(1:mdim)
1544 this%kcoords_cub(1:mdim, ikp, ikpt) = matmul(this%ktransf(1:mdim, 1:mdim), kvec(1:mdim)) &
1545 + kpoint(1:mdim)
1546 end do
1547 end do
1548
1549 end if
1550
1551
1552 end if
1553
1554
1555
1556
1558
1559 contains
1560
1561 ! Fill the non-periodic direction
1562 subroutine fill_non_periodic_dimension(this)
1563 type(pes_flux_t), intent(inout) :: this
1564
1565 integer :: sign
1566 real(real64) :: kpar(1:pdim), val
1567
1568 ikp = ikp + 1
1569
1570 sign = 1
1571 if (ikk /= 0) sign= ikk / abs(ikk)
1572
1573 kpar(1:pdim) = kvec(1:pdim)
1574 val = abs(ikk) * de * m_two - sum(kpar(1:pdim)**2)
1575 if (val >= 0) then
1576 kvec(mdim) = sign * sqrt(val)
1577 else ! if E < p//^2/2
1578 !FIXME: Should handle the exception doing something smarter than this
1579 kvec(mdim) = sqrt(val) ! set to NaN
1580 nkp_out = nkp_out + 1
1581 end if
1582
1583 this%kcoords_cub(1:mdim, ikp, ikpt) = kvec(1:mdim)
1584
1585
1586 end subroutine fill_non_periodic_dimension
1587
1588
1589
1590 end subroutine pes_flux_reciprocal_mesh_gen
1591
1592 ! ---------------------------------------------------------
1593 subroutine pes_flux_calc(this, space, mesh, st, der, ext_partners, kpoints, iter, dt, stopping)
1594 type(pes_flux_t), intent(inout) :: this
1595 class(space_t), intent(in) :: space
1596 type(mesh_t), intent(in) :: mesh
1597 type(states_elec_t), intent(inout) :: st
1598 type(derivatives_t), intent(in) :: der
1599 type(partner_list_t), intent(in) :: ext_partners
1600 type(kpoints_t), intent(in) :: kpoints
1601 integer, intent(in) :: iter
1602 real(real64), intent(in) :: dt
1603 logical, intent(in) :: stopping
1604
1605 integer :: stst, stend, kptst, kptend, sdim, mdim
1606 integer :: ist, ik, isdim, imdim
1607 integer :: isp,ip
1608 integer :: il, ip_local
1609 complex(real64), allocatable :: gpsi(:,:), psi(:)
1610 complex(real64), allocatable :: interp_values(:)
1611 logical :: contains_isp
1612 real(real64) :: kpoint(1:3)
1613 complex(real64), allocatable :: tmp_array(:)
1614
1615 type(lasers_t), pointer :: lasers
1616
1617 push_sub_with_profile(pes_flux_calc)
1618
1619 if (iter > 0) then
1620
1621 stst = st%st_start
1622 stend = st%st_end
1623 kptst = st%d%kpt%start
1624 kptend = st%d%kpt%end
1625 sdim = st%d%dim
1626 mdim = space%dim
1627
1628 safe_allocate(psi(1:mesh%np_part))
1629 safe_allocate(gpsi(1:mesh%np_part, 1:mdim))
1630
1631 if (this%surf_interp) then
1632 safe_allocate(interp_values(1:this%nsrfcpnts))
1633 end if
1634
1635 this%itstep = this%itstep + 1
1636
1637 ! get and save current laser field
1638 lasers => list_get_lasers(ext_partners)
1639 if(associated(lasers)) then
1640 do il = 1, lasers%no_lasers
1641 call laser_field(lasers%lasers(il), this%veca(1:mdim, this%itstep), iter*dt)
1642 end do
1643 end if
1644 this%veca(:, this%itstep) = - this%veca(:, this%itstep)
1645
1646! Ideally one could directly access uniform_vector_potential
1647! this%veca(:, this%itstep) = hm%hm_base%uniform_vector_potential(:)
1648
1649 if(mesh%parallel_in_domains) then
1650 safe_allocate(tmp_array(1:this%nsrfcpnts))
1651 end if
1652
1653 ! save wavefunctions & gradients
1654 do ik = kptst, kptend
1655 do ist = stst, stend
1656 do isdim = 1, sdim
1657 call states_elec_get_state(st, mesh, isdim, ist, ik, psi)
1658
1659 ! We apply the PBC first before applying the phase
1660 call boundaries_set(der%boundaries, mesh, psi)
1661
1662 if (this%surf_shape == pes_plane) then
1663 ! Apply the phase containing kpoint only
1664 kpoint(1:mdim) = kpoints%get_point(st%d%get_kpoint_index(ik))
1665
1666 !$omp parallel do schedule(static)
1667 do ip = 1, mesh%np_part
1668 psi(ip) = exp(- m_zi * sum(mesh%x(ip, 1:mdim) * kpoint(1:mdim))) * psi(ip)
1669 end do
1670 !$omp end parallel do
1671 end if
1672
1673 call zderivatives_grad(der, psi, gpsi, set_bc = .false.)
1674
1675 if (this%surf_interp) then
1676 call mesh_interpolation_evaluate(this%interp, this%nsrfcpnts, psi(1:mesh%np_part), &
1677 this%rcoords(1:mdim, 1:this%nsrfcpnts), interp_values(1:this%nsrfcpnts))
1678 this%wf(1:this%nsrfcpnts, ist, isdim, ik, this%itstep) = st%occ(ist, ik) * interp_values(1:this%nsrfcpnts)
1679 do imdim = 1, mdim
1680 call mesh_interpolation_evaluate(this%interp, this%nsrfcpnts, gpsi(1:mesh%np_part, imdim), &
1681 this%rcoords(1:mdim, 1:this%nsrfcpnts), interp_values(1:this%nsrfcpnts))
1682 this%gwf(1:this%nsrfcpnts, ist, isdim, ik, this%itstep, imdim) = st%occ(ist, ik) * interp_values(1:this%nsrfcpnts)
1683 end do
1684
1685 else
1687 contains_isp = .true.
1688 do isp = 1, this%nsrfcpnts
1689 if (mesh%parallel_in_domains) then
1690 if (mesh%mpi_grp%rank == this%rankmin(isp)) then
1691 contains_isp = .true.
1692 else
1693 contains_isp = .false.
1694 end if
1695 end if
1696
1697 if (contains_isp) then
1698 ip_local = this%srfcpnt(isp)
1699 this%wf(isp, ist, isdim, ik, this%itstep) = st%occ(ist, ik) * psi(ip_local)
1700 this%gwf(isp, ist, isdim, ik, this%itstep, 1:mdim) = &
1701 st%occ(ist, ik) * gpsi(ip_local, 1:mdim)
1702 end if
1703 end do
1704 if (mesh%parallel_in_domains) then
1705
1706 tmp_array = this%wf(1:this%nsrfcpnts, ist, isdim, ik, this%itstep)
1707 call mesh%allreduce(tmp_array)
1708 this%wf(1:this%nsrfcpnts, ist, isdim, ik, this%itstep) = tmp_array
1709
1710 do imdim = 1, mdim
1711 tmp_array = this%gwf(1:this%nsrfcpnts, ist, isdim, ik, this%itstep, imdim)
1712 call mesh%allreduce(tmp_array)
1713 this%gwf(1:this%nsrfcpnts, ist, isdim, ik, this%itstep, imdim) = tmp_array
1714 end do
1715 end if
1716 end if
1717 end do
1718 end do
1719 end do
1720
1721 if(mesh%parallel_in_domains) then
1722 safe_deallocate_a(tmp_array)
1723 end if
1724
1725 safe_deallocate_a(psi)
1726 safe_deallocate_a(gpsi)
1727
1728 if (this%itstep == this%tdsteps .or. mod(iter, this%save_iter) == 0 .or. iter == this%max_iter .or. stopping) then
1729 if (this%surf_shape == pes_cubic .or. this%surf_shape == pes_plane) then
1730 call pes_flux_integrate_cub(this, space, mesh, st, kpoints, dt)
1731 else
1732 call pes_flux_integrate_sph(this, mesh, st, dt)
1733 end if
1734
1735 ! clean up fields
1736 this%wf = m_z0
1737 this%gwf = m_z0
1738 this%veca = m_zero
1739 this%itstep = 0
1740 end if
1741
1742 end if
1743
1744 pop_sub_with_profile(pes_flux_calc)
1745 end subroutine pes_flux_calc
1746
1747
1748
1749 ! ---------------------------------------------------------
1750 subroutine pes_flux_integrate_cub_tabulate(this, space, mesh, st, kpoints)
1751 type(pes_flux_t), intent(inout) :: this
1752 class(space_t), intent(in) :: space
1753 type(mesh_t), intent(in) :: mesh
1754 type(states_elec_t), intent(in) :: st
1755 type(kpoints_t), intent(in) :: kpoints
1756
1757 integer :: kptst, kptend, mdim
1758 integer :: ik, j1, j2, jvec(1:2)
1759 integer :: isp, ikp, ikp_start, ikp_end
1760 integer :: ik_map
1761
1762 integer :: idir, pdim, nfaces, ifc, n_dir
1763 complex(real64) :: tmp
1764 real(real64) :: kpoint(1:3), vec(1:3), lvec(1:3)
1765
1767
1768 if (kpoints%have_zero_weight_path()) then
1769 kptst = st%d%kpt%start
1770 kptend = st%d%kpt%end
1771 else
1772 kptst = 1
1773 kptend = 1
1774 end if
1775
1776 mdim = space%dim
1777 pdim = space%periodic_dim
1778
1779 ikp_start = this%nkpnts_start
1780 ikp_end = this%nkpnts_end
1781
1782
1783 if (this%surf_shape == pes_plane) then
1784 ! This is not general but should work in the specific case where it is relevant
1785 !i.e. when the system is semiperiodic in <=2 dimensions
1786 this%srfcnrml(:, 1:this%nsrfcpnts) = this%srfcnrml(:, 1:this%nsrfcpnts) * mesh%coord_system%surface_element(3)
1787 end if
1788
1789 if (.not. this%use_symmetry .or. kpoints%have_zero_weight_path()) then
1790
1791 safe_allocate(this%expkr(1:this%nsrfcpnts,ikp_start:ikp_end,kptst:kptend,1))
1792
1793 do ik = kptst, kptend
1794 do ikp = ikp_start, ikp_end
1795 do isp = 1, this%nsrfcpnts
1796 this%expkr(isp,ikp,ik,1) = exp(-m_zi * sum(this%rcoords(1:mdim,isp) &
1797 * this%kcoords_cub(1:mdim, ikp, ik))) &
1798 * (m_two * m_pi)**(-mdim / m_two)
1799
1800 end do
1801 end do
1802 end do
1803
1804
1805
1806 else !do something smart to exploit symmetries
1807
1808 nfaces = mdim*2
1809 if (this%surf_shape == pes_plane) nfaces = 1
1810
1811
1812 safe_allocate(this%expkr(1:this%nsrfcpnts,maxval(this%ll(1:mdim)),kptst:kptend,1:mdim))
1813 this%expkr(:,:,:,:) = m_z1
1814
1815 do ik = kptst, kptend !should only be ik=1
1816 do idir = 1, mdim
1817 do ikp = 1, this%ll(idir)
1818 do isp = 1, this%nsrfcpnts
1819 this%expkr(isp,ikp,ik,idir) = exp(-m_zi * this%rcoords(idir,isp) &
1820 * this%klinear(ikp, idir)) &
1821 * (m_two * m_pi)**(-m_one / m_two)
1822
1823 end do
1824 end do
1825 end do
1826 end do
1827
1828 safe_allocate(this%expkr_perp(1:maxval(this%ll(1:mdim)), 1:nfaces))
1829 this%expkr_perp(:,:) = m_z1
1830
1831 do ifc = 1, nfaces
1832 if (this%face_idx_range(ifc, 1) < 0) cycle ! this face have no local surface point
1833 isp = this%face_idx_range(ifc, 1)
1834 do idir = 1, mdim
1835 if (abs(this%srfcnrml(idir, isp)) >= m_epsilon) n_dir = idir
1836 end do
1837
1838 do ikp = 1, this%ll(n_dir)
1839 this%expkr_perp(ikp, ifc) = exp(- m_zi * this%rcoords(n_dir,isp) &
1840 * this%klinear(ikp, n_dir)) &
1841 * (m_two * m_pi)**(-m_one / m_two)
1842
1844 end do
1845 end do
1846
1847 end if
1848
1849
1850 if (space%is_periodic()) then
1851 !Tabulate the Born-von Karman phase
1852 safe_allocate(this%bvk_phase(ikp_start:ikp_end,st%d%kpt%start:st%d%kpt%end))
1853
1854 this%bvk_phase(:,:) = m_z0
1855 vec(:) = m_zero
1856
1857 do ik = st%d%kpt%start, st%d%kpt%end
1858 kpoint(1:mdim) = kpoints%get_point(st%d%get_kpoint_index(ik))
1859 if (kpoints%have_zero_weight_path()) then
1860 ik_map = ik
1861 else
1862 ik_map = 1
1863 end if
1864 do ikp = ikp_start, ikp_end
1865 vec(1:pdim) = this%kcoords_cub(1:pdim, ikp, ik_map) + kpoint(1:pdim)
1866 do j1 = 0, kpoints%nik_axis(1) - 1
1867 do j2 = 0, kpoints%nik_axis(2) - 1
1868 jvec(1:2)=(/j1, j2/)
1869 lvec(1:pdim) = matmul(kpoints%latt%rlattice(1:pdim, 1:2), jvec(1:2))
1870 tmp = sum(lvec(1:pdim) * vec(1:pdim))
1871 this%bvk_phase(ikp, ik) = this%bvk_phase(ikp,ik) &
1872 + exp(m_zi * tmp)
1873
1874 end do
1875 end do
1876
1877 end do
1878 end do
1879 this%bvk_phase(:,:) = this%bvk_phase(:,:) * m_one / product(kpoints%nik_axis(1:pdim))
1880
1881
1882! Keep this because is useful for debug but not enough to bother having it polished out
1883! if (debug%info .and. mpi_grp_is_root(mpi_world)) then
1884! write(225,*) "ik, ikp, this%bvk_phase(ikp,ik)"
1885! do ik = st%d%kpt%start, st%d%kpt%end
1886! do ikp = ikp_start, ikp_end
1887! write(225,*) ik, ikp, this%bvk_phase(ikp,ik)
1888! end do
1889! end do
1890! flush(225)
1891! end if
1892
1893 end if
1894
1896
1897 end subroutine pes_flux_integrate_cub_tabulate
1898
1899
1900 ! ---------------------------------------------------------
1901 pure function get_ikp(this, ikpu, ikpv, ikpz, n_dir) result(ikp)
1902 type(pes_flux_t), intent(in) :: this
1903 integer, intent(in) :: ikpu
1904 integer, intent(in) :: ikpv
1905 integer, intent(in) :: ikpz
1906 integer, intent(in) :: n_dir
1907 integer :: ikp
1908
1909 select case (n_dir)
1910 case (1)
1911 ikp = this%Lkpuvz_inv(ikpz, ikpu, ikpv)
1912 case (2)
1913 ikp = this%Lkpuvz_inv(ikpu, ikpz, ikpv)
1914 case (3)
1915 ikp = this%Lkpuvz_inv(ikpu, ikpv, ikpz)
1916
1917 case default
1918 ! should die here but cannot use assert in a pure function
1919
1920 end select
1921
1922 end function get_ikp
1923
1924 ! ---------------------------------------------------------
1925 subroutine pes_flux_distribute_facepnts_cub(this, mesh)
1926 type(pes_flux_t), intent(inout) :: this
1927 type(mesh_t), intent(in) :: mesh
1928
1929
1930 integer :: mdim, nfaces, ifc, ifp_start, ifp_end
1931
1933
1934 mdim = mesh%box%dim
1935
1936 nfaces = mdim*2
1937 if (this%surf_shape == pes_plane) nfaces = 1
1938
1939 do ifc = 1, nfaces
1940
1941 ifp_start = this%face_idx_range(ifc, 1)
1942 ifp_end = this%face_idx_range(ifc, 2)
1943
1944
1945 if (this%nsrfcpnts_start <= ifp_end) then ! the local domain could be in this face
1946
1947 if (this%nsrfcpnts_start >= ifp_start) then
1948 if (this%nsrfcpnts_start <= ifp_end) then
1949 this%face_idx_range(ifc, 1) = this%nsrfcpnts_start
1950 else
1951 this%face_idx_range(ifc, 1:2) = -1 ! the local domain is not in this face
1952 end if
1953 end if
1954
1955 if (this%nsrfcpnts_end <= ifp_end) then
1956 if (this%nsrfcpnts_end >= ifp_start) then
1957 this%face_idx_range(ifc, 2) = this%nsrfcpnts_end
1958 else
1959 this%face_idx_range(ifc, 1:2) = -1 ! the local domain is not in this face
1960 end if
1961 end if
1962
1963 else
1964
1965 this%face_idx_range(ifc, 1:2) = -1 ! the local domain is not in this face
1966
1967 end if
1968
1969! Keep this because is useful for debug but not enough to bother having it polished out
1970! if (debug%info) then
1971! print *, mpi_world%rank, ifc, ifp_start, ifp_end, this%face_idx_range(ifc, 1:2), this%nsrfcpnts_start, this%nsrfcpnts_end
1972! end if
1973
1974 end do
1975
1976
1979
1980
1981 ! ---------------------------------------------------------
1982 subroutine pes_flux_integrate_cub(this, space, mesh, st, kpoints, dt)
1983 type(pes_flux_t), intent(inout) :: this
1984 class(space_t), intent(in) :: space
1985 type(mesh_t), intent(in) :: mesh
1986 type(states_elec_t), intent(inout) :: st
1987 type(kpoints_t), intent(in) :: kpoints
1988 real(real64), intent(in) :: dt
1989
1990 integer :: stst, stend, kptst, kptend, sdim, mdim
1991 integer :: ist, ik, isdim, imdim, ik_map
1992 integer :: ikp, itstep
1993 integer :: idir, n_dir, nfaces
1994 complex(real64), allocatable :: jk_cub(:,:,:,:), spctramp_cub(:,:,:,:)
1995 integer :: ikp_start, ikp_end, isp_start, isp_end
1996 real(real64) :: vec, kpoint(1:3)
1997 integer :: ifc
1998 complex(real64), allocatable :: wfpw(:), gwfpw(:)
1999 complex(real64), allocatable :: phase(:,:),vphase(:,:)
2000
2001 integer :: tdstep_on_node
2002 integer :: nfp
2003
2004 !Symmetry helpers
2005 integer :: ikpu, ikpv, ikpz, dir_on_face(1:2)
2006 complex(real64) :: face_int_gwf, face_int_wf
2007
2008
2009 push_sub(pes_flux_integrate_cub)
2010
2011 ! this routine is parallelized over time steps and surface points
2012
2013 call profiling_in('pes_flux_integrate_cub')
2014
2015 call messages_write("Debug: calculating pes_flux cub surface integral (accelerated direct expression)")
2016 call messages_info(debug_only=.true.)
2017
2019 tdstep_on_node = 1
2020 if (bitand(this%par_strategy, option__pes_flux_parallelization__pf_time) /= 0) then
2021 if (mesh%parallel_in_domains) tdstep_on_node = mesh%mpi_grp%rank + 1
2022 end if
2023
2024 stst = st%st_start
2025 stend = st%st_end
2026 kptst = st%d%kpt%start
2027 kptend = st%d%kpt%end
2028 sdim = st%d%dim
2029 mdim = mesh%box%dim
2030
2031 ikp_start = this%nkpnts_start
2032 ikp_end = this%nkpnts_end
2033
2034
2035
2036 safe_allocate(jk_cub(stst:stend, 1:sdim, kptst:kptend, ikp_start:ikp_end))
2037 safe_allocate(spctramp_cub(stst:stend, 1:sdim, kptst:kptend, ikp_start:ikp_end))
2038 spctramp_cub = m_z0
2039
2040
2041 safe_allocate(vphase(ikp_start:ikp_end, kptst:kptend))
2042 safe_allocate(phase(ikp_start:ikp_end, kptst:kptend))
2043 vphase(:,:) = m_z0
2044 phase(:,:) = m_z0
2045
2046
2047 nfaces = mdim * 2
2048 if (this%surf_shape == pes_plane) nfaces = 1
2049
2050 safe_allocate( wfpw(ikp_start:ikp_end))
2051 safe_allocate(gwfpw(ikp_start:ikp_end))
2052
2053
2054 do ifc = 1, nfaces
2055
2056 if (this%face_idx_range(ifc, 1)<0) cycle ! this face have no local surface point
2057
2058
2059 isp_start = this%face_idx_range(ifc, 1)
2060 isp_end = this%face_idx_range(ifc, 2)
2061
2062
2063 nfp = isp_end - isp_start + 1 ! faces can have a different number of points
2064
2065 wfpw = m_z0
2066 gwfpw = m_z0
2067
2068 ! get the directions normal to the surface and parallel to it
2069 imdim = 1
2070 do idir = 1, mdim
2071 if (abs(this%srfcnrml(idir, isp_start)) >= m_epsilon) then
2072 n_dir = idir
2073 else
2074 dir_on_face(imdim)=idir
2075 imdim = imdim + 1
2076 end if
2077 end do
2078
2079
2080 ! get the previous Volkov phase
2081 vphase(ikp_start:ikp_end,:) = this%conjgphase_prev(ikp_start:ikp_end,:)
2082
2083 jk_cub(:, :, :, :) = m_z0
2084
2085 do itstep = 1, this%itstep
2086
2087 do ik = kptst, kptend
2088
2089 if (kpoints%have_zero_weight_path()) then
2090 ik_map = ik
2091 else
2092 ik_map = 1
2093 end if
2094
2095 kpoint(1:mdim) = kpoints%get_point(st%d%get_kpoint_index(ik))
2096 do ikp = ikp_start, ikp_end
2097 vec = sum((this%kcoords_cub(1:mdim, ikp, ik_map) - this%veca(1:mdim, itstep) / p_c)**2)
2098 vphase(ikp, ik) = vphase(ikp, ik) * exp(m_zi * vec * dt / m_two)
2099
2100
2101 if (space%is_periodic()) then
2102 phase(ikp, ik) = vphase(ikp, ik) * this%bvk_phase(ikp,ik)
2103 else
2104 phase(ikp, ik) = vphase(ikp, ik)
2105 end if
2106
2107 end do
2108
2109 if (itstep /= tdstep_on_node) cycle
2110
2111 do ist = stst, stend
2112 do isdim = 1, sdim
2113
2114
2115 ! calculate the surface integrals
2116 if (.not. this%use_symmetry .or. kpoints%have_zero_weight_path()) then
2117
2118 do ikp = ikp_start , ikp_end
2119
2120 gwfpw(ikp) = &
2121 sum(this%gwf(isp_start:isp_end, ist, isdim, ik, itstep, n_dir) &
2122 * this%expkr(isp_start:isp_end, ikp, ik_map, 1) &
2123 * this%srfcnrml(n_dir, isp_start:isp_end))
2124
2125
2126 wfpw(ikp) = &
2127 sum(this%wf(isp_start:isp_end, ist, isdim, ik, itstep) &
2128 * this%expkr(isp_start:isp_end, ikp,ik_map, 1) &
2129 * this%srfcnrml(n_dir, isp_start:isp_end))
2130 end do
2131
2132 else
2133 !$omp parallel do private (ikpv,ikpz,face_int_gwf,face_int_wf) shared(gwfpw, wfpw)
2134 do ikpu = 1, this%ll(dir_on_face(1))
2135 do ikpv = 1, this%ll(dir_on_face(2))
2136
2137
2138 face_int_gwf = sum(this%gwf(isp_start:isp_end, ist, isdim, ik, itstep, n_dir) &
2139 * this%expkr(isp_start:isp_end, ikpu, ik_map, dir_on_face(1)) &
2140 * this%expkr(isp_start:isp_end, ikpv, ik_map, dir_on_face(2)) &
2141 * this%srfcnrml(n_dir, isp_start:isp_end))
2142
2143 face_int_wf = sum(this%wf(isp_start:isp_end, ist, isdim, ik, itstep) &
2144 * this%expkr(isp_start:isp_end, ikpu, ik_map, dir_on_face(1)) &
2145 * this%expkr(isp_start:isp_end, ikpv, ik_map, dir_on_face(2)) &
2146 * this%srfcnrml(n_dir, isp_start:isp_end))
2147
2148 do ikpz = 1, this%ll(n_dir)
2149
2150 gwfpw(get_ikp(this, ikpu, ikpv, ikpz, n_dir)) = face_int_gwf &
2151 * this%expkr_perp(ikpz, ifc)
2152 wfpw( get_ikp(this, ikpu, ikpv, ikpz, n_dir)) = face_int_wf &
2153 * this%expkr_perp(ikpz, ifc)
2154
2155 end do
2156
2157 end do
2158 end do
2159
2160
2161 end if
2162
2163
2164 ! Sum it up
2165 jk_cub(ist, isdim, ik, ikp_start:ikp_end) = jk_cub(ist, isdim, ik, ikp_start:ikp_end) + &
2166 phase(ikp_start:ikp_end, ik) * (wfpw(ikp_start:ikp_end) * &
2167 (m_two * this%veca(n_dir, itstep) / p_c - this%kcoords_cub(n_dir, ikp_start:ikp_end, ik_map)) + &
2168 m_zi * gwfpw(ikp_start:ikp_end))
2169
2170
2171 end do ! isdim
2172 end do ! ist
2173 end do ! is
2174
2175 end do !istep
2176
2177
2178 spctramp_cub(:,:,:,:) = spctramp_cub(:,:,:,:) + jk_cub(:,:,:,:) * m_half
2179
2180
2181 end do ! face loop
2182
2183
2184 this%conjgphase_prev(ikp_start:ikp_end, :) = vphase(ikp_start:ikp_end, :)
2185
2186
2187 if (mesh%parallel_in_domains .and.(bitand(this%par_strategy, option__pes_flux_parallelization__pf_time) /= 0 &
2188 .or. bitand(this%par_strategy, option__pes_flux_parallelization__pf_surface) /= 0)) then
2189 call mesh%allreduce(spctramp_cub)
2190 end if
2191
2192
2193
2194 this%spctramp_cub = this%spctramp_cub + spctramp_cub * dt
2195
2196 safe_deallocate_a(gwfpw)
2197 safe_deallocate_a(wfpw)
2198
2199 safe_deallocate_a(jk_cub)
2200 safe_deallocate_a(spctramp_cub)
2201 safe_deallocate_a(phase)
2202 safe_deallocate_a(vphase)
2203
2204 call profiling_out('pes_flux_integrate_cub')
2205
2206 pop_sub(pes_flux_integrate_cub)
2207 end subroutine pes_flux_integrate_cub
2208
2209
2210
2211 ! ---------------------------------------------------------
2212 subroutine pes_flux_integrate_sph(this, mesh, st, dt)
2213 type(pes_flux_t), intent(inout) :: this
2214 type(mesh_t), intent(in) :: mesh
2215 type(states_elec_t), intent(inout) :: st
2216 real(real64), intent(in) :: dt
2217
2218 integer :: stst, stend, kptst, kptend, sdim
2219 integer :: ist, ik, isdim, imdim
2220 integer :: itstep, lmax
2221 integer :: ll, mm, ip
2222 complex(real64), allocatable :: s1_node(:,:,:,:,:,:), s1_act(:,:,:)
2223 complex(real64), allocatable :: s2_node(:,:,:,:,:), s2_act(:,:)
2224 complex(real64), allocatable :: integ11_t(:,:,:), integ12_t(:,:,:)
2225 complex(real64), allocatable :: integ21_t(:,:), integ22_t(:,:)
2226 complex(real64), allocatable :: spctramp_sph(:,:,:,:,:)
2227 integer :: ikk, ikk_start, ikk_end
2228 integer :: isp_start, isp_end
2229 integer :: iomk, iomk_start, iomk_end
2230 complex(real64), allocatable :: phase_act(:,:)
2231 real(real64) :: vec
2232 complex(real64) :: weight_x, weight_y, weight_z
2233 integer :: tdstep_on_node
2234
2235 push_sub_with_profile(pes_flux_integrate_sph)
2236
2237 ! this routine is parallelized over time steps and surface points
2238
2239 call messages_write("Debug: calculating pes_flux sph surface integral")
2240 call messages_info(debug_only=.true.)
2241
2242 tdstep_on_node = 1
2243 if (mesh%parallel_in_domains) tdstep_on_node = mesh%mpi_grp%rank + 1
2244
2245 stst = st%st_start
2246 stend = st%st_end
2247 kptst = st%d%kpt%start
2248 kptend = st%d%kpt%end
2249 sdim = st%d%dim
2250
2251 lmax = this%lmax
2252 ikk_start = this%nk_start
2253 ikk_end = this%nk_end
2254 isp_start = this%nsrfcpnts_start
2255 isp_end = this%nsrfcpnts_end
2256 iomk_start = this%nstepsomegak_start
2257 iomk_end = this%nstepsomegak_end
2258
2259 !TODO: The storage (-lmax:lmax, 0:lmax) is not efficient for large l
2260
2261 ! surface integral S_lm (for time step on node)
2262 safe_allocate(s1_node(1:3, -lmax:lmax, 0:lmax, stst:stend, 1:sdim, kptst:kptend))
2263 safe_allocate(s2_node(-lmax:lmax, 0:lmax, stst:stend, 1:sdim, kptst:kptend))
2264
2265 safe_allocate(s1_act(1:3, -lmax:lmax, 0:lmax))
2266 safe_allocate(s2_act(-lmax:lmax, 0:lmax))
2267
2268 call profiling_in("PES_FLU_SPH_INT")
2269
2270 do itstep = 1, this%itstep
2271
2272 do ik = kptst, kptend
2273 do ist = stst, stend
2274 do isdim = 1, sdim
2275
2276 ! lmax can be as large as 80, so this is better to parallelize over ll than paying the price
2277 ! of (really) many reductions
2278 !$omp parallel do private(mm, ip, weight_x, weight_y, weight_z)
2279 do ll = 0, lmax
2280 do mm = -ll, ll
2281 s1_act(1:3, mm, ll) = m_zero
2282 s2_act(mm, ll) = m_zero
2283 ! surface integrals
2284 do ip = isp_start, isp_end
2285 weight_x = this%ylm_r(ip, ll, mm) * this%srfcnrml(1, ip)
2286 weight_y = this%ylm_r(ip, ll, mm) * this%srfcnrml(2, ip)
2287 weight_z = this%ylm_r(ip, ll, mm) * this%srfcnrml(3, ip)
2288 s1_act(1, mm, ll) = s1_act(1, mm, ll) + weight_x * this%wf(ip, ist, isdim, ik, itstep)
2289 s1_act(2, mm, ll) = s1_act(2, mm, ll) + weight_y * this%wf(ip, ist, isdim, ik, itstep)
2290 s1_act(3, mm, ll) = s1_act(3, mm, ll) + weight_z * this%wf(ip, ist, isdim, ik, itstep)
2291
2292 s2_act(mm, ll) = s2_act(mm, ll) + weight_x * this%gwf(ip, ist, isdim, ik, itstep, 1) &
2293 + weight_y * this%gwf(ip, ist, isdim, ik, itstep, 2) &
2294 + weight_z * this%gwf(ip, ist, isdim, ik, itstep, 3)
2295 end do
2296 end do
2297 end do
2298 !$omp end parallel do
2299
2300 if (mesh%parallel_in_domains .and. bitand(this%par_strategy, option__pes_flux_parallelization__pf_surface) /= 0) then
2301 call mesh%allreduce(s1_act)
2302 call mesh%allreduce(s2_act)
2303 end if
2304
2305 if (itstep == tdstep_on_node) then
2306 s1_node(:, :, :, ist, isdim, ik) = s1_act(:, :, :)
2307 s2_node(:, :, ist, isdim, ik) = s2_act(:, :)
2308 end if
2309
2310 end do
2311 end do
2312 end do
2313 end do
2314
2315 call profiling_out("PES_FLU_SPH_INT")
2316
2317 ! spectral amplitude for k-points on node
2318 safe_allocate(integ11_t(1:this%nstepsomegak, 1:3, 0:lmax))
2319 safe_allocate(integ21_t(1:this%nstepsomegak, 0:lmax))
2320 safe_allocate(integ12_t(1:this%nstepsomegak, 1:3, ikk_start:ikk_end))
2321 safe_allocate(integ22_t(1:this%nstepsomegak, ikk_start:ikk_end))
2322
2323 safe_allocate(spctramp_sph(1:this%nstepsomegak, stst:stend, 1:sdim, kptst:kptend, 1:this%nk))
2324 spctramp_sph = m_z0
2325
2326 ! get the previous Volkov phase
2327 safe_allocate(phase_act(1:this%nstepsomegak, ikk_start:ikk_end))
2328 if (ikk_start > 0) then
2329 phase_act(:, ikk_start:ikk_end) = this%conjgphase_prev(:, ikk_start:ikk_end)
2330 else
2331 phase_act(:, ikk_start:ikk_end) = m_z0
2332 end if
2333
2334 do itstep = 1, this%itstep
2335 ! calculate the current Volkov phase
2336 do ikk = ikk_start, ikk_end
2337 do iomk = 1, this%nstepsomegak
2338 vec = sum((this%kcoords_sph(1:3, ikk, iomk) - this%veca(1:3, itstep) / p_c)**2)
2339 phase_act(iomk, ikk) = phase_act(iomk, ikk) * exp(m_zi * vec * dt * m_half)
2340 end do
2341 end do
2342
2343 do ik = kptst, kptend
2344 do ist = stst, stend
2345 do isdim = 1, sdim
2346
2347 ! communicate the current surface integrals
2348 if (itstep == tdstep_on_node) then
2349 s1_act(:,:,:) = s1_node(:, :, :, ist, isdim, ik)
2350 s2_act(:,:) = s2_node(:, :, ist, isdim, ik)
2351 end if
2352 if (mesh%parallel_in_domains) then
2353 call mesh%mpi_grp%bcast(s1_act, (lmax + 1) * (2 * lmax + 1) * 3, mpi_double_complex, itstep - 1)
2354 call mesh%mpi_grp%bcast(s2_act, (lmax + 1) * (2 * lmax + 1), mpi_double_complex, itstep - 1)
2355 end if
2356
2357 integ12_t = m_z0
2358 integ22_t = m_z0
2359 integ11_t = m_z0
2360 integ21_t = m_z0
2361
2362 call profiling_in("PES_FLU_SPH_INTEGXX")
2363
2364 !$omp parallel do private(ll, mm, imdim, ikk)
2365 do ll = 0, lmax
2366 do mm = -ll, ll
2367 ! multiply with spherical harmonics & sum over all mm
2368 do imdim = 1, 3
2369 !$omp simd
2370 do ikk = iomk_start, iomk_end
2371 integ11_t(ikk, imdim, ll) = integ11_t(ikk, imdim, ll) + s1_act(imdim, mm, ll) * this%ylm_k(ikk, mm, ll)
2372 end do
2373 end do
2374 !$omp simd
2375 do ikk = iomk_start, iomk_end
2376 integ21_t(ikk, ll) = integ21_t(ikk, ll) + s2_act(mm, ll) * this%ylm_k(ikk, mm, ll)
2377 end do
2378 end do
2379 end do
2380 !$omp end parallel do
2381
2382 call mesh%allreduce(integ11_t)
2383 call mesh%allreduce(integ21_t)
2384
2385 do ll = 0, lmax
2386 weight_x = ( - m_zi)**ll
2387 ! multiply with Bessel function & sum over all ll
2388 !$omp parallel do private(imdim)
2389 do ikk = ikk_start, ikk_end
2390 do imdim = 1, 3
2391 integ12_t(:, imdim, ikk) = integ12_t(:, imdim, ikk) + this%j_l(ikk, ll) * weight_x * integ11_t(:, imdim, ll)
2392 end do
2393 integ22_t(:, ikk) = integ22_t(:, ikk) + this%j_l(ikk, ll) * weight_x * integ21_t(:, ll)
2394 end do
2395 end do
2396
2397 call profiling_out("PES_FLU_SPH_INTEGXX")
2398
2399 do ikk = ikk_start, ikk_end
2400 ! sum over dimensions
2401 do imdim = 1, 3
2402 spctramp_sph(:, ist, isdim, ik, ikk) = &
2403 spctramp_sph(:, ist, isdim, ik, ikk) + &
2404 phase_act(:, ikk) * (integ12_t(:, imdim, ikk) * &
2405 (m_two * this%veca(imdim, itstep) / p_c - this%kcoords_sph(imdim, ikk,:)))
2406 end do
2407 spctramp_sph(:, ist, isdim, ik, ikk) = &
2408 spctramp_sph(:, ist, isdim, ik, ikk) + phase_act(:, ikk) * integ22_t(:, ikk) * m_zi
2409 end do
2410 end do
2411 end do
2412 end do
2413 end do
2414
2415 safe_deallocate_a(s1_node)
2416 safe_deallocate_a(s2_node)
2417 safe_deallocate_a(s1_act)
2418 safe_deallocate_a(s2_act)
2419
2420 safe_deallocate_a(integ11_t)
2421 safe_deallocate_a(integ12_t)
2422 safe_deallocate_a(integ21_t)
2423 safe_deallocate_a(integ22_t)
2424
2425 ! save the Volkov phase and the spectral amplitude
2426 this%conjgphase_prev = m_z0
2427
2428 if (ikk_start > 0) then
2429 this%conjgphase_prev(:, ikk_start:ikk_end) = phase_act(:, ikk_start:ikk_end)
2430 end if
2431 safe_deallocate_a(phase_act)
2432
2433 if (mesh%parallel_in_domains .and. bitand(this%par_strategy, option__pes_flux_parallelization__pf_time) /= 0) then
2434 call mesh%allreduce(this%conjgphase_prev)
2435 call mesh%allreduce(spctramp_sph)
2436 end if
2437
2438 !TODO: Use axpy here
2439 do ikk = 1, this%nk
2440 do ik = kptst, kptend
2441 call lalg_axpy(this%nstepsomegak, stend-stst+1, sdim, dt, spctramp_sph(1:this%nstepsomegak, stst:stend, 1:sdim, ik, ikk), &
2442 this%spctramp_sph(1:this%nstepsomegak, stst:stend, 1:sdim, ik, ikk))
2443 end do
2444 end do
2445 safe_deallocate_a(spctramp_sph)
2446
2447 pop_sub_with_profile(pes_flux_integrate_sph)
2448 end subroutine pes_flux_integrate_sph
2449
2450 ! ---------------------------------------------------------
2451 subroutine pes_flux_getcube(this, mesh, abb, border, offset, fc_ptdens)
2452 type(mesh_t), intent(in) :: mesh
2453 type(pes_flux_t), intent(inout) :: this
2454 type(absorbing_boundaries_t), intent(in) :: abb
2455 real(real64), intent(in) :: border(:)
2456 real(real64), intent(in) :: offset(:)
2457 real(real64), intent(in) :: fc_ptdens
2458
2459 integer, allocatable :: which_surface(:)
2460 real(real64) :: xx(mesh%box%dim), dd, area, dS(mesh%box%dim, 2), factor
2461 integer :: mdim, imdim, idir, isp, pm,nface, idim, ndir, iu,iv, iuv(1:2)
2462 integer(int64) :: ip_global
2463 integer :: npface
2464 integer :: rankmin, nsurfaces
2465 logical :: in_ab
2466 integer :: ip_local, NN(mesh%box%dim, 2), idx(mesh%box%dim, 2)
2467 integer :: isp_end, isp_start, ifc, n_dir, nfaces, mindim
2468 real(real64) :: RSmax(2), RSmin(2), RS(2), dRS(2), weight
2469
2470
2471 push_sub(pes_flux_getcube)
2472
2473 ! this routine works on absolute coordinates
2474
2475
2476 mdim = mesh%box%dim
2477
2478
2479
2480 safe_allocate(this%face_idx_range(1:mdim * 2, 1:2))
2481 safe_allocate(this%LLr(mdim, 1:2))
2482 safe_allocate(this%NN(mdim, 1:2))
2483
2484 this%face_idx_range(:, :) = 0
2485 this%LLr(:, :) = m_zero
2486 this%NN(:, :) = 1
2487 nn(:, :) = 1
2488
2489 if (this%surf_interp) then
2490 ! Create a surface with points not on the mesh
2491
2492 idx(:,:) = 0
2493
2494 mindim = 1
2495 factor = m_two
2496 if (this%surf_shape == pes_plane) then
2497 mindim = mdim ! We only have two planes along the non periodic dimension
2498 ! this is due to the fact that for semiperiodic systems we multiply border
2499 ! by two to prenvet the creation of surfaces at the edges
2500 factor = m_one
2501 end if
2502
2503
2504 this%nsrfcpnts = 0
2505
2506 do ndir = mdim, mindim, -1
2507 area = m_one
2508 do idir=1, mdim
2509 if (idir == ndir) cycle
2510 area = area * border(idir) * factor
2511 end do
2512
2513 npface = int(fc_ptdens * area) ! number of points on the face
2514
2515 idim = 1
2516 do idir=1, mdim
2517 if (idir == ndir) cycle
2518 nn(ndir, idim) = int(sqrt(npface * (border(idir) * factor)**2 / area))
2519 ds(ndir, idim) = border(idir) * factor / nn(ndir, idim)
2520 this%LLr(ndir, idim) = nn(ndir, idim) * ds(ndir, idim)
2521 idx(ndir, idim) = idir
2522 idim = idim + 1
2523 end do
2524 this%nsrfcpnts = this%nsrfcpnts + 2 * product(nn(ndir, 1:mdim - 1))
2525 end do
2526
2527 this%NN(1:mdim, :) = nn(1:mdim, :)
2528
2529
2530 assert(this%nsrfcpnts > 0) !at this point should be satisfied otherwise the point density is way too small
2531
2532
2533 safe_allocate(this%srfcnrml(1:mdim, 0:this%nsrfcpnts))
2534 safe_allocate(this%rcoords(1:mdim, 0:this%nsrfcpnts))
2535
2536 this%srfcnrml(:, :) = m_zero
2537 this%rcoords(:, :) = m_zero
2538
2539
2540 isp = 1
2541 ifc = 1
2542 do ndir = mdim, mindim, -1
2543
2544 !Up face
2545 this%face_idx_range(ifc,1) = isp
2546 do iu = 1, nn(ndir,1)
2547 do iv = 1, nn(ndir,2)
2548 this%rcoords(ndir, isp) = border(ndir)
2549 this%srfcnrml(ndir, isp) = product(ds(ndir,1:mdim-1))
2550 iuv =(/iu, iv/)
2551 do idim = 1, mdim-1
2552 this%rcoords(idx(ndir, idim), isp) = (-nn(ndir, idim) * m_half - m_half + iuv(idim)) * ds(ndir, idim)
2553 end do
2554 isp = isp + 1
2555 end do
2556 end do
2557 this%face_idx_range(ifc, 2) = isp-1
2558
2559 ifc = ifc + 1
2560
2561 !Down face
2562 this%face_idx_range(ifc, 1) = isp
2563 do iu = 1, nn(ndir, 1)
2564 do iv = 1, nn(ndir, 2)
2565 this%rcoords(ndir, isp) = -border(ndir)
2566 this%srfcnrml(ndir, isp) = -product(ds(ndir, 1:mdim - 1))
2567 iuv =(/iu, iv/)
2568 do idim = 1, mdim-1
2569 this%rcoords(idx(ndir, idim), isp) = (-nn(ndir, idim) * m_half -m_half + iuv(idim)) * ds(ndir, idim)
2570 end do
2571 isp = isp + 1
2572 end do
2573 end do
2574 this%face_idx_range(ifc, 2) = isp-1
2575
2576 ifc = ifc + 1
2577 end do
2578
2579 do isp = 1, this%nsrfcpnts
2580 this%rcoords(1:mdim, isp) = mesh%coord_system%to_cartesian(this%rcoords(1:mdim, isp))
2581 end do
2582
2583 else
2584 ! Surface points are on the mesh
2585
2586 nfaces = mdim * 2
2587 if (this%surf_shape == pes_plane) nfaces = 2
2588
2589 in_ab = .false.
2590
2591 safe_allocate(which_surface(1:mesh%np_global))
2592 which_surface = 0
2593
2594 ! get the surface points
2595 this%nsrfcpnts = 0
2596 do ip_local = 1, mesh%np
2597 ip_global = mesh_local2global(mesh, ip_local)
2598
2599 nsurfaces = 0
2600
2601 xx(1:mdim) = mesh%x(ip_local, 1:mdim) - offset(1:mdim)
2602
2603 ! eventually check whether we are in absorbing zone
2604 select case (abb%abtype)
2605 case (mask_absorbing)
2606 in_ab = .not. is_close(abb%mf(ip_local), m_one)
2607 case (imaginary_absorbing)
2608 in_ab = abs(abb%mf(ip_local)) > m_epsilon
2609 case default
2610 assert(.false.)
2611 end select
2612
2613 ! check whether the point is inside the cube
2614 if (all(abs(xx(1:mdim)) <= border(1:mdim)) .and. .not. in_ab) then
2615 ! check whether the point is close to any border
2616 do imdim = 1, mdim
2617 if (this%surf_shape == pes_plane) then
2618 dd = border(imdim) - xx(imdim)
2619 else
2620 dd = border(imdim) - abs(xx(imdim))
2621 end if
2622 if (dd < mesh%spacing(imdim) / m_two) then
2623 nsurfaces = nsurfaces + 1
2624 which_surface(ip_global) = int(sign(m_one, xx(imdim))) * imdim ! +-x=+-1, +-y=+-2, +-z=+-3
2625 end if
2626 end do
2627
2628 ! check whether the point is close to one border only (not an edge)
2629 if (nsurfaces == 1) then
2630 this%nsrfcpnts = this%nsrfcpnts + 1
2631 else
2632 which_surface(ip_global) = 0
2633 end if
2634 end if
2635
2636 end do
2637
2638 if (mesh%parallel_in_domains) then
2639 assert(mesh%np_global < huge(0_int32))
2640 call mesh%allreduce(this%nsrfcpnts)
2641 call mesh%allreduce(which_surface)
2642 end if
2643
2644 safe_allocate(this%srfcpnt(1:this%nsrfcpnts))
2645 safe_allocate(this%srfcnrml(1:mdim, 0:this%nsrfcpnts))
2646 safe_allocate(this%rcoords(1:mdim, 0:this%nsrfcpnts))
2647 safe_allocate(this%rankmin(1:this%nsrfcpnts))
2648
2649 this%srfcnrml = m_zero
2650 this%rcoords = m_zero
2651
2652 this%face_idx_range(:,:) = 0
2653
2654 isp = 0
2655 nface = 0
2656 do idir = mdim, 1, -1 ! Start counting from z to simplify implementation for semiperiodic systems (pln)
2657 do pm = 1, -1, -2
2658 nface = nface + 1
2659 this%face_idx_range(nface, 1) = isp + 1
2660
2661 do ip_global = 1, mesh%np_global
2662 if (abs(which_surface(ip_global)) == idir .and. sign(1, which_surface(ip_global)) == pm) then
2663 isp = isp + 1
2664 ! coordinate of surface point
2665 xx(1:mdim) = mesh_x_global(mesh, ip_global)
2666 this%rcoords(1:mdim, isp) = xx(1:mdim)
2667 ! local ip & node which has the surface point
2668 this%srfcpnt(isp) = mesh_nearest_point(mesh, this%rcoords(1:mdim, isp), dd, rankmin)
2669 this%rankmin(isp) = rankmin
2670 ! surface normal
2671 this%srfcnrml(idir, isp) = sign(1, which_surface(ip_global))
2672 ! add the surface element (of directions orthogonal to the normal vector)
2673 do imdim = 1, mdim
2674 if (imdim == idir) cycle
2675 this%srfcnrml(idir, isp) = this%srfcnrml(idir, isp) * mesh%spacing(imdim)
2676 end do
2677 end if
2678 end do
2679
2680 this%face_idx_range(nface,2) = isp
2681
2682 end do
2683 end do
2684
2685
2686 !Get dimensions and spacing on each (pair of) face
2687 do ifc = 1, nint((nfaces + 0.5) / 2)
2688 isp_start = this%face_idx_range(ifc, 1)
2689 isp_end = this%face_idx_range(ifc, 2)
2690
2691 !retrieve face normal direction
2692 n_dir = 0
2693 do idir = 1, mdim
2694 if (abs(this%srfcnrml(idir, isp_start)) >= m_epsilon) n_dir = idir
2695 end do
2696
2697 !retrieve the spacing on the face
2698 drs(:) = m_zero
2699 idim = 1
2700 do idir = 1, mdim
2701 if (idir == n_dir) cycle
2702 drs(idim)= mesh%spacing(idir)
2703 idim = idim + 1
2704 end do
2705
2706
2707 !retrieve the dimensions of the face
2708 rsmin = m_zero
2709 rsmax = m_zero
2710 do isp = isp_start, isp_end
2711 !measure in reduced coordinates
2712 xx = mesh%coord_system%from_cartesian(this%rcoords(1:mdim, isp))
2713 idim = 1
2714 do idir = 1, mdim
2715 if (idir == n_dir) cycle
2716 rs(idim)=xx(idir)
2717 if (rs(idim) < rsmin(idim)) rsmin(idim) = rs(idim)
2718 if (rs(idim) > rsmax(idim)) rsmax(idim) = rs(idim)
2719 idim = idim + 1
2720 end do
2721 end do
2722
2723
2724 do idir = 1, mdim - 1
2725 this%LLr(n_dir, idir) = rsmax(idir) - rsmin(idir) + drs(idir)
2726 if (drs(idir) > m_zero) this%NN(n_dir, idir) = int(this%LLr(n_dir, idir) / drs(idir))
2727 end do
2728
2729 end do
2730
2731
2732 end if
2733
2734 if (this%anisotrpy_correction) then
2735 !Compensate the fact that the angular distribution of points is not uniform
2736 !and have peaks in correspondence to the edges and corners of the parallelepiped
2737 do isp = 1, this%nsrfcpnts
2738 weight = sum(this%rcoords(1:mdim, isp) * this%srfcnrml(1:mdim, isp))
2739 weight = weight/sum(this%rcoords(1:mdim, isp)**2) / sum(this%srfcnrml(1:mdim, isp)**2)
2740 this%srfcnrml(1:mdim, isp) = this%srfcnrml(1:mdim, isp) * abs(weight)
2741 end do
2742
2743 end if
2744
2745
2746 safe_deallocate_a(which_surface)
2747
2748 pop_sub(pes_flux_getcube)
2749 end subroutine pes_flux_getcube
2750
2751 ! ---------------------------------------------------------
2752 subroutine pes_flux_getsphere(this, mesh, nstepsthetar, nstepsphir)
2753 type(pes_flux_t), intent(inout) :: this
2754 type(mesh_t), intent(in) :: mesh
2755 integer, intent(inout) :: nstepsthetar, nstepsphir
2756
2757 integer :: mdim
2758 integer :: isp, ith, iph
2759 real(real64) :: thetar, phir
2760 real(real64) :: weight, dthetar
2761
2762 push_sub(pes_flux_getsphere)
2763
2764 mdim = mesh%box%dim
2765
2766 if (nstepsphir == 0) nstepsphir = 1
2767
2768 if (nstepsthetar <= 1) then
2769 nstepsphir = 1
2770 nstepsthetar = 1
2771 end if
2772 dthetar = m_pi / nstepsthetar
2773 this%nsrfcpnts = nstepsphir * (nstepsthetar - 1) + 2
2774
2775 safe_allocate(this%srfcnrml(1:mdim, 0:this%nsrfcpnts))
2776 safe_allocate(this%rcoords(1:mdim, 0:this%nsrfcpnts))
2777
2778 this%srfcnrml = m_zero
2779 this%rcoords = m_zero
2780
2781 ! initializing spherical grid
2782 isp = 0
2783 do ith = 0, nstepsthetar
2784 thetar = ith * dthetar
2785
2786 if (ith == 0 .or. ith == nstepsthetar) then
2787 weight = (m_one - cos(dthetar / m_two)) * m_two * m_pi
2788 else
2789 weight = abs(cos(thetar - dthetar / m_two) - cos(thetar + dthetar / m_two)) &
2790 * m_two * m_pi / nstepsphir
2791 end if
2792
2793 do iph = 0, nstepsphir - 1 ! 2*pi is the same than zero
2794 isp = isp + 1
2795 phir = iph * m_two * m_pi / nstepsphir
2796 this%srfcnrml(1, isp) = cos(phir) * sin(thetar)
2797 if (mdim >= 2) this%srfcnrml(2, isp) = sin(phir) * sin(thetar)
2798 if (mdim == 3) this%srfcnrml(3, isp) = cos(thetar)
2799 this%rcoords(1:mdim, isp) = this%radius * this%srfcnrml(1:mdim, isp)
2800 ! here we also include the surface elements
2801 this%srfcnrml(1:mdim, isp) = this%radius**m_two * weight * this%srfcnrml(1:mdim, isp)
2802 if (ith == 0 .or. ith == nstepsthetar) exit
2803 end do
2804 end do
2805
2806 pop_sub(pes_flux_getsphere)
2807 end subroutine pes_flux_getsphere
2808
2809
2810
2811
2812 ! ---------------------------------------------------------
2813 subroutine pes_flux_distribute(istart_global, iend_global, istart, iend, comm)
2814 integer, intent(in) :: istart_global
2815 integer, intent(in) :: iend_global
2816 integer, intent(inout) :: istart
2817 integer, intent(inout) :: iend
2818 type(mpi_comm), intent(in) :: comm
2819
2820#if defined(HAVE_MPI)
2821 integer, allocatable :: dimrank(:)
2822 integer :: mpisize, mpirank, irest, irank
2823 integer :: inumber
2824#endif
2825
2826 push_sub(pes_flux_distribute)
2827
2828#if defined(HAVE_MPI)
2829 call mpi_comm_size(comm, mpisize)
2830 call mpi_comm_rank(comm, mpirank)
2831
2832 safe_allocate(dimrank(0:mpisize-1))
2833
2834 inumber = iend_global - istart_global + 1
2835 irest = mod(inumber, mpisize)
2836 do irank = 0, mpisize - 1
2837 if (irest > 0) then
2838 dimrank(irank) = int(inumber / mpisize) + 1
2839 irest = irest - 1
2840 else
2841 dimrank(irank) = int(inumber / mpisize)
2842 end if
2843 end do
2844
2845 iend = istart_global + sum(dimrank(:mpirank)) - 1
2846 istart = iend - dimrank(mpirank) + 1
2847
2848 if (istart > iend) then
2849 iend = 0
2850 istart = 0
2851 end if
2852
2853 safe_deallocate_a(dimrank)
2854
2855#else
2856 istart = istart_global
2857 iend = iend_global
2858#endif
2859
2860 pop_sub(pes_flux_distribute)
2861 end subroutine pes_flux_distribute
2862
2863
2864#include "pes_flux_out_inc.F90"
2865
2866end module pes_flux_oct_m
2867
2868!! Local Variables:
2869!! mode: f90
2870!! coding: utf-8
2871!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:171
Prints out to iunit a message in the form: ["InputVariable" = value] where "InputVariable" is given b...
Definition: messages.F90:180
double exp(double __x) __attribute__((__nothrow__
double sin(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
double floor(double __x) __attribute__((__nothrow__
integer, parameter, public imaginary_absorbing
integer, parameter, public mask_absorbing
Module implementing boundary conditions in Octopus.
Definition: boundaries.F90:122
This module handles the calculation mode.
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public zderivatives_grad(der, ff, op_ff, ghost_update, set_bc, to_cartesian)
apply the gradient to a mesh function
integer, parameter, public spin_polarized
type(lasers_t) function, pointer, public list_get_lasers(partners)
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
complex(real64), parameter, public m_z0
Definition: global.F90:198
complex(real64), parameter, public m_zi
Definition: global.F90:202
real(real64), parameter, public m_epsilon
Definition: global.F90:204
complex(real64), parameter, public m_z1
Definition: global.F90:199
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public p_c
Electron gyromagnetic ratio, see Phys. Rev. Lett. 130, 071801 (2023)
Definition: global.F90:224
real(real64), parameter, public m_one
Definition: global.F90:189
real(real64), parameter, public m_three
Definition: global.F90:191
This module defines classes and functions for interaction partners.
Definition: io.F90:114
integer pure function, public kpoints_number(this)
Definition: kpoints.F90:1099
integer, parameter, public e_field_vector_potential
Definition: lasers.F90:177
integer pure elemental function, public laser_kind(laser)
Definition: lasers.F90:715
subroutine, public laser_field(laser, field, time)
Retrieves the value of either the electric or the magnetic field. If the laser is given by a scalar p...
Definition: lasers.F90:1115
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public ylmr_cmplx(xx, li, mi, ylm)
Computes spherical harmonics ylm at position (x, y, z)
Definition: math.F90:305
subroutine, public mesh_interpolation_init(this, mesh)
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
integer function, public mesh_nearest_point(mesh, pos, dmin, rankmin)
Returns the index of the point which is nearest to a given vector position pos.
Definition: mesh.F90:382
integer(int64) function, public mesh_local2global(mesh, ip)
This function returns the global mesh index for a given local index.
Definition: mesh.F90:967
real(real64) function, dimension(1:mesh%box%dim), public mesh_x_global(mesh, ipg)
Definition: mesh.F90:812
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
subroutine, public messages_new_line()
Definition: messages.F90:1134
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:502
integer function, public parse_block(namespace, name, blk, check_varinfo_)
Definition: parser.F90:618
subroutine, public pes_flux_map_from_states(this, restart, st, ll, pesP, krng, Lp, istin)
Definition: pes_flux.F90:3020
subroutine pes_flux_getcube(this, mesh, abb, border, offset, fc_ptdens)
Definition: pes_flux.F90:2545
subroutine, public pes_flux_calc(this, space, mesh, st, der, ext_partners, kpoints, iter, dt, stopping)
Definition: pes_flux.F90:1687
subroutine, public pes_flux_out_energy(this, pesK, file, namespace, ll, Ekin, dim)
Definition: pes_flux.F90:3677
subroutine pes_flux_getsphere(this, mesh, nstepsthetar, nstepsphir)
Definition: pes_flux.F90:2846
subroutine, public pes_flux_output(this, box, st, namespace)
Definition: pes_flux.F90:3883
subroutine, public pes_flux_pmesh(this, namespace, dim, kpoints, ll, pmesh, idxZero, krng, nspin, Lp, Ekin)
Definition: pes_flux.F90:2979
subroutine, public pes_flux_reciprocal_mesh_gen(this, namespace, space, st, kpoints, comm, post)
Definition: pes_flux.F90:819
subroutine, public pes_flux_load(restart, this, st, ierr)
Definition: pes_flux.F90:4506
subroutine, public pes_flux_dump(restart, this, mesh, st, ierr)
Definition: pes_flux.F90:4413
subroutine pes_flux_distribute_facepnts_cub(this, mesh)
Definition: pes_flux.F90:2019
subroutine pes_flux_integrate_sph(this, mesh, st, dt)
Definition: pes_flux.F90:2306
integer, parameter pes_cartesian
Definition: pes_flux.F90:263
integer, parameter pes_plane
Definition: pes_flux.F90:258
subroutine pes_flux_distribute(istart_global, iend_global, istart, iend, comm)
Definition: pes_flux.F90:2907
subroutine pes_flux_integrate_cub(this, space, mesh, st, kpoints, dt)
Definition: pes_flux.F90:2076
subroutine, public pes_flux_end(this)
Definition: pes_flux.F90:772
pure integer function get_ikp(this, ikpu, ikpv, ikpz, n_dir)
Definition: pes_flux.F90:1995
subroutine pes_flux_integrate_cub_tabulate(this, space, mesh, st, kpoints)
Definition: pes_flux.F90:1844
subroutine, public pes_flux_init(this, namespace, space, mesh, st, ext_partners, kpoints, abb, save_iter, max_iter)
Definition: pes_flux.F90:273
integer, parameter pes_spherical
Definition: pes_flux.F90:258
subroutine, public pes_flux_out_vmap(this, pesK, file, namespace, ll, pmesh, dim)
Definition: pes_flux.F90:3767
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
This module handles spin dimensions of the states and the k-point distribution.
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:223
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:118
character pure function, public index2axis(idir)
Definition: utils.F90:202
subroutine fill_non_periodic_dimension(this)
Definition: pes_flux.F90:1656
Class implementing a parallelepiped box. Currently this is restricted to a rectangular cuboid (all th...
Class implementing a spherical box.
Definition: box_sphere.F90:132
class representing derivatives
Describes mesh distribution to nodes.
Definition: mesh.F90:186
The states_elec_t class contains all electronic wave functions.
int true(void)