Octopus
pes_out.F90
Go to the documentation of this file.
1!! Copyright (C) 2016 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_out_oct_m
24 use cube_oct_m
25 use debug_oct_m
26 use global_oct_m
27 use io_oct_m
29 use loct_oct_m
30 use math_oct_m
33#if defined(HAVE_NETCDF)
34 use netcdf
35#endif
37 use qshep_oct_m
38 use space_oct_m
39 use sort_oct_m
40 use string_oct_m
41 use unit_oct_m
43 use vtk_oct_m
44
45 implicit none
46
47 private
48
49 public :: &
53
54 integer, parameter :: &
55 INTEGRATE_NONE = -1, &
56 integrate_phi = 1, &
57 integrate_theta = 2, &
58 integrate_r = 3, &
59 integrate_kx = 4, &
60 integrate_ky = 5, &
61 integrate_kz = 6
62
63
64
65contains
66
67 ! ---------------------------------------------------------
68 subroutine pes_out_velocity_map(pesK, file, namespace, space, Lk, ll, how, spacing, coord_system, pmesh)
69 real(real64), intent(in) :: pesK(:,:,:)
70 character(len=*), intent(in) :: file
71 type(namespace_t), intent(in) :: namespace
72 class(space_t), intent(in) :: space
73 real(real64), intent(in) :: Lk(:,:)
74 integer, intent(in) :: ll(:)
75 integer(int64), intent(in) :: how
76 real(real64), intent(in) :: spacing(3)
77 class(coordinate_system_t), intent(in) :: coord_system
78 real(real64), optional, intent(in) :: pmesh(:,:,:,:)
79
80 integer :: ierr
81 character(len=512) :: filename
82 type(cube_t) :: cube
83 type(cube_function_t) :: cf
84 real(real64) :: dk(3)
85
86 push_sub(pes_out_velocity_map)
87
88 call cube_init(cube, ll, namespace, space, spacing, coord_system)
89 call dcube_function_alloc_rs(cube, cf, force_alloc = .true.)
90 cf%dRS = pesk
91
92
93 if (.not. present(pmesh)) then
94 ! Ignore Lk and use pmesh
95 dk(:) = m_zero
96 dk(1:space%dim) = abs(lk(2, 1:space%dim) - lk(1,1:space%dim))
97 end if
98
99#if defined(HAVE_NETCDF)
100
101 if (bitand(how, option__outputformat__netcdf) /= 0) then
102 filename = trim(file)//".ncdf"
103 write(message(1), '(a)') 'Writing netcdf format file: '
104 call messages_info(1, namespace=namespace)
105
106 call dout_cf_netcdf(filename, ierr, cf, cube, space, dk(:), .false., unit_one/units_out%length, namespace)
107
108 end if
109
110#endif
111
112 if (bitand(how, option__outputformat__vtk) /= 0) then
113 filename = trim(file)//".vtk"
114 write(message(1), '(a)') 'Writing vtk format file: '
115 call messages_info(1, namespace=namespace)
116
117 if (present(pmesh)) then
118 call dvtk_out_cf_structured(filename, namespace, 'PES_vel_map', ierr, cf, cube,&
119 unit_one/units_out%length, pmesh, ascii = .false.)
120 else
121 call dvtk_out_cf(filename, namespace, 'PES_vel_map', ierr, cf, cube, dk(:),&
122 unit_one/units_out%length)
123 end if
124 end if
125
126 call cube_end(cube)
127 call dcube_function_free_rs(cube, cf)
128
129 pop_sub(pes_out_velocity_map)
130
131 end subroutine pes_out_velocity_map
132
133
134 ! ---------------------------------------------------------
135 subroutine pes_out_arpes_cut(namespace, arpes, file, dim, ll, pmesh, Ekin)
136 type(namespace_t), intent(in) :: namespace
137 real(real64), intent(in) :: arpes(:,:,:)
138 character(len=*), intent(in) :: file
139 integer, intent(in) :: dim
140 integer, intent(in) :: ll(:)
141 real(real64), intent(in) :: pmesh(:,:,:,:)
142 real(real64), intent(in) :: Ekin(:,:,:)
143
144 integer :: iunit, ip, ie
145 real(real64) :: dp, length, pp(1:2), pdiff(1:2)
146
148
149 iunit = io_open(file, namespace, action='write')
150 write(iunit, '(a)') '##################################################'
151 write(iunit, '(a1,a18,2x,a18,2x,a18,2x,a18,2x, a18,2x,a18)') '#', &
152 str_center("Ppath", 18), &
153 str_center("px", 18), str_center("py", 18), &
154 str_center("E", 18), str_center("P[Ppath,E]", 18)
155
156 write(iunit, '(a1,a18,2x,a18,2x,a18,2x,a18,2x, a18,2x,a18)') '#', &
157 str_center('[hbar/'//trim(units_abbrev(units_out%length)) // ']', 18), &
158 str_center('[hbar/'//trim(units_abbrev(units_out%length)) // ']', 18), &
159 str_center('[hbar/'//trim(units_abbrev(units_out%length)) // ']', 18), &
160 str_center('['//trim(units_abbrev(units_out%energy)) // ']', 18), &
161 str_center('[1/' //trim(units_abbrev(units_out%energy))//']', 18)
162 write(iunit, '(a)') '##################################################'
163
164
165 length = m_zero
166 pdiff(:) = m_zero
167
168 do ip = 1, ll(1)
169 pp(1:2) = pmesh(ip,1,1,1:2)
170 if (ip > 1) pdiff = pp(1:2) - pmesh(ip-1,1,1,1:2)
171 dp = norm2(pdiff(1:2))
172 length = length + dp
173
174 select case (dim)
175 case (2)
176 do ie = 1, ll(2)
177 write(iunit, '(es19.12,2x,es19.12,2x,es19.12,2x,es19.12,2x,es19.12,2x,es19.12)') &
178 units_from_atomic(unit_one/units_out%length, length),&
179 units_from_atomic(unit_one/units_out%length, pp(1)), &
180 units_from_atomic(unit_one/units_out%length, pp(2)), &
181 units_from_atomic(units_out%energy, ekin(ip,ie,1)), &
182 arpes(ip,ie,1)
183 end do
184
185 case (3)
186 do ie = 1, ll(3)
187 write(iunit, '(es19.12,2x,es19.12,2x,es19.12,2x,es19.12,2x,es19.12,2x,es19.12)') &
188 units_from_atomic(unit_one/units_out%length, length),&
189 units_from_atomic(unit_one/units_out%length, pp(1)), &
190 units_from_atomic(unit_one/units_out%length, pp(2)), &
191 units_from_atomic(units_out%energy, ekin(ip,1,ie)), &
192 arpes(ip,1,ie)
193 end do
194
195 end select
196
197
198 write(iunit, *)
199
200
201 end do
202
203
204 call io_close(iunit)
205
206 pop_sub(pes_out_arpes_cut)
207 end subroutine pes_out_arpes_cut
208
209
210 ! ---------------------------------------------------------
211 subroutine pes_out_velocity_map_cut(namespace, pesK, file, ll, dim, pol, dir, integrate, pos, Lk, pmesh)
212 type(namespace_t), intent(in) :: namespace
213 real(real64), intent(in) :: pesk(:,:,:)
214 character(len=*), intent(in) :: file
215 integer, intent(in) :: ll(:)
216 integer, intent(in) :: dim
217 real(real64), intent(in) :: pol(3)
218 integer, intent(in) :: dir
219 integer, intent(in) :: integrate
220 integer, optional, intent(in) :: pos(3)
221 real(real64), optional, intent(in) :: lk(:,:)
222 real(real64), optional, intent(in) :: pmesh(:,:,:,:)
223
224 integer :: ii, ix, iy, iunit,ldir(2), icut(3)
225 real(real64) :: kk(3),temp
226 integer, allocatable :: idx(:,:)
227 real(real64), allocatable :: lk_(:,:)
228 real(real64) :: rotation(1:dim,1:dim)
229 logical :: aligned_axis
230 ! integration
231 real(real64) :: k, kkk(3), theta, phi, dphi, dk(3)
232 integer :: iph, nphi
233 ! progress
234 integer :: idone, ntodo
235
236 real(real64), allocatable :: cube_f(:)
237 type(qshep_t) :: interp
238
240
241 iunit = io_open(file, namespace, action='write')
242
243
244
245 assert(size(pesk, 1) == ll(1))
246
247 if (.not. present(pmesh)) then
248 assert(present(lk))
249
250 safe_allocate(idx(1:maxval(ll(:)), 1:3))
251 safe_allocate(lk_(1:maxval(ll(:)), 1:3))
252
253 dk(:) = m_zero
254 dk(1:dim) = abs(lk(2,1:dim)-lk(1,1:dim))
255
256 do ii = 1, 3
257 lk_(:,ii) = lk(:,ii)
258 call sort(lk_(1:ll(ii), ii), idx(1:ll(ii), ii)) !We need to sort the k-vectors in order to dump in gnuplot format
259 end do
260 end if
261
262 aligned_axis = sum((pol-(/0 ,0 ,1/))**2) <= m_epsilon .or. &
263 sum((pol-(/0 ,1 ,0/))**2) <= m_epsilon .or. &
264 sum((pol-(/1 ,0 ,0/))**2) <= m_epsilon
265
266 if (present(pos)) then
267 icut(1:3) = pos(1:3)
268 else
269 icut(1:3) = ll(1:3)/2 + 1
270 end if
271
272 if (aligned_axis .and. integrate == integrate_none) then !no need to rotate and interpolate
273
274 select case (dir)
275 case (1)
276 ldir(:) =(/2,3/)
277 case (2)
278 ldir(:) =(/1,3/)
279 case (3)
280 ldir(:) =(/1,2/)
281
282 end select
283
284
285 if (present(pmesh)) then
286 do ix = 1, ll(ldir(1))
287 do iy = 1, ll(ldir(2))
288
289 select case (dir)
290 case (1)
291 temp = pesk(icut(dir), ix, iy)
292 kk(1) = pmesh(icut(dir), ix, iy, ldir(1))
293 kk(2) = pmesh(icut(dir), ix, iy, ldir(2))
294 case (2)
295 temp = pesk(ix, icut(dir), iy)
296 kk(1) = pmesh(ix, icut(dir), iy, ldir(1))
297 kk(2) = pmesh(ix, icut(dir), iy, ldir(2))
298 case (3)
299 temp = pesk(ix, iy, icut(dir))
300 kk(1) = pmesh(ix, iy, icut(dir), ldir(1))
301 kk(2) = pmesh(ix, iy, icut(dir), ldir(2))
302 end select
303
304 write(iunit, '(es19.12,2x,es19.12,2x,es19.12)') &
305 units_from_atomic(sqrt(units_out%energy), kk(1)),&
306 units_from_atomic(sqrt(units_out%energy), kk(2)),&
307 temp
308 end do
309 write(iunit, *)
310 end do
311
312 else
313 do ix = 1, ll(ldir(1))
314 kk(1) = lk_(ix, ldir(1))
315 do iy = 1, ll(ldir(2))
316 kk(2) = lk_(iy, ldir(2))
317
318 select case (dir)
319 case (1)
320 temp = pesk(idx(icut(dir), 1), idx(ix, 2), idx(iy, 3))
321 case (2)
322 temp = pesk(idx(ix, 1), idx(icut(dir), 2), idx(iy, 3))
323 case (3)
324 temp = pesk(idx(ix, 1), idx(iy, 2), idx(icut(dir), 3))
325 end select
326
327 write(iunit, '(es19.12,2x,es19.12,2x,es19.12)') &
328 units_from_atomic(sqrt(units_out%energy), kk(1)),&
329 units_from_atomic(sqrt(units_out%energy), kk(2)),&
330 temp
331 end do
332 write(iunit, *)
333 end do
334 end if
335
336 else
337 ! We set the z-axis along the pol vector
338 call generate_rotation_matrix(rotation, (/ m_zero, m_zero, m_one /), pol)
339
340 if (debug%info) then
341 print *,"Rotate z-axis over the zenith axis"
342 print *,rotation(1,:)
343 print *,rotation(2,:)
344 print *,rotation(3,:)
345 end if
346
347 call pes_out_interpolator_init(namespace, pesk, lk, ll, dim, cube_f, interp, pmesh)
348
349 ntodo = product(ll(1:2))
350 idone = 0
351 call loct_progress_bar(-1, ntodo)
352
353 do ix = 1, ll(1)
354 do iy = 1, ll(2)
355
356 !cut
357 select case (dir)
358 case (1)
359 kk(1) = m_zero
360 if (present(pmesh)) then
361 kk(2) = pmesh(ix, 1, 1, 1)
362 kk(3) = pmesh(1, iy, 1, 2)
363 else
364 kk(2) = lk_(ix, 1)
365 kk(3) = lk_(iy, 2)
366 end if
367 case (2)
368 kk(2) = m_zero
369 if (present(pmesh)) then
370 kk(1) = pmesh(ix, 1, 1, 1)
371 kk(3) = pmesh(1, iy, 1, 2)
372 else
373 kk(1) = lk_(ix, 1)
374 kk(3) = lk_(iy, 2)
375 end if
376
377 case (3)
378 kk(3) = m_zero
379 if (present(pmesh)) then
380 kk(1) = pmesh(ix, 1, 1, 1)
381 kk(2) = pmesh(1, iy, 1, 2)
382 else
383 kk(1) = lk_(ix, 1)
384 kk(2) = lk_(iy, 2)
385 end if
386
387 end select
388
389 temp = qshep_interpolate(interp, cube_f, matmul(rotation,kk(1:3)))
390
391 select case (integrate)
392 case (integrate_phi)
393 temp = m_zero
394 k = norm2(kk(1:3))
395
396 nphi = 360
397 dphi = m_two * m_pi / nphi
398
399 do iph = 0, nphi
400 phi = iph * dphi
401 theta = atan2(norm2(kk(1:2)), kk(3))
402
403 kkk(1) = k *sin(theta) * cos(phi)
404 kkk(2) = k *sin(theta) * sin(phi)
405 kkk(3) = k *cos(theta)
406
407 temp = temp + &
408 abs(qshep_interpolate(interp, cube_f, matmul(rotation,kkk(1:3))))
409 end do
410 temp = temp * dphi
411
412 case (integrate_kx)
413 temp = m_zero
414 do ii =1, ll(1)
415 kkk(:) = kk(:) + (/lk(ii,1), m_zero, m_zero/)
416 temp = temp + &
417 abs(qshep_interpolate(interp, cube_f, matmul(rotation,kkk(1:3))))
418 end do
419 temp = temp * dk(1)
420
421 case (integrate_ky)
422 temp = m_zero
423 do ii =1, ll(2)
424 kkk(:) = kk(:) + (/m_zero, lk(ii,2), m_zero/)
425 temp = temp + &
426 abs(qshep_interpolate(interp, cube_f, matmul(rotation,kkk(1:3))))
427 end do
428 temp = temp * dk(2)
429
430 case (integrate_kz)
431 temp = m_zero
432 do ii =1, ll(3)
433 kkk(:) = kk(:) + (/m_zero, m_zero, lk(ii,3)/)
434 temp = temp + &
435 abs(qshep_interpolate(interp, cube_f, matmul(rotation,kkk(1:3))))
436 end do
437 temp = temp * dk(3)
438
439 end select
440
441 write(iunit, '(es19.12,2x,es19.12,2x,es19.12)') &
442 units_from_atomic(sqrt(units_out%energy), kk(1)),&
443 units_from_atomic(sqrt(units_out%energy), kk( 2)),&
444 temp
445
446
447 idone = idone +1
448 call loct_progress_bar(idone, ntodo)
449
450 end do
451 write(iunit, *)
452 end do
453
454 write(stdout, '(1x)')
455
456
457 end if
458
459 call io_close(iunit)
460
461 call pes_out_interpolator_end(cube_f, interp)
462
463 safe_deallocate_a(idx)
464 safe_deallocate_a(lk_)
465
467 end subroutine pes_out_velocity_map_cut
468
469 ! --------------------------------------------------------
470 !
474 !
475 ! ---------------------------------------------------------
476 subroutine pes_out_interpolator_init(namespace, pesK, Lk, ll, dim, cube_f, interp, pmesh)
477 type(namespace_t), intent(in) :: namespace
478 real(real64), intent(in) :: pesk(:,:,:)
479 real(real64), intent(in) :: lk(:,:)
480 integer, intent(in) :: ll(:)
481 integer, intent(in) :: dim
482 real(real64), allocatable, intent(out) :: cube_f(:)
483 type(qshep_t), intent(out) :: interp
484 real(real64), optional, intent(in) :: pmesh(:,:,:,:)
485
486 integer :: np, ii, ix, iy, iz
487 real(real64) :: kk(3)
488 real(real64), allocatable :: kx(:),ky(:),kz(:)
489
491
492
493 call messages_write("Initializing Qshep interpolator. Be patient it may take a while... ")
494 call messages_info(namespace=namespace)
495
496 np = ll(1)*ll(2)*ll(3)
497
498 !check dim
499 if (dim < 2 .or. dim > 3) then
500 message(1) = "This interpolator works only for 2 <= dim <= 3."
501 call messages_fatal(1, namespace=namespace)
502 end if
503
504 safe_allocate(cube_f(1:np))
505
506 safe_allocate(kx(1:np))
507 safe_allocate(ky(1:np))
508 safe_allocate(kz(1:np))
509
510 cube_f = m_zero
511 kx = m_zero
512 ky = m_zero
513 kz = m_zero
514
515
516 ii=1
517 if (present(pmesh)) then
518 do ix = 1, ll(1)
519 do iy = 1, ll(2)
520 do iz = 1, ll(3)
521
522 cube_f(ii) = pesk(ix,iy,iz)
523
524 kx(ii) = pmesh(ix,iy,iz,1)
525 ky(ii) = pmesh(ix,iy,iz,2)
526 kz(ii) = pmesh(ix,iy,iz,3)
527
528 ii = ii +1
529 end do
530 end do
531 end do
532
533 else
534
535 do ix = 1, ll(1)
536 kk(1) = lk(ix, 1)
537 do iy = 1, ll(2)
538 kk(2) = lk(iy, 2)
539 do iz = 1, ll(3)
540 kk(3) = lk(iz, 3)
541
542 cube_f(ii) = pesk(ix,iy,iz)
543
544 kx(ii) = kk(1)
545 ky(ii) = kk(2)
546 kz(ii) = kk(3)
547
548 ii = ii +1
549 end do
550 end do
551 end do
552 end if
553
554
555 select case (dim)
556 case (2)
557 call qshep_init(interp, i4_to_i8(np), cube_f, kx, ky)
558 case (3)
559 call qshep_init(interp, i4_to_i8(np), cube_f, kx, ky, kz)
560 end select
561
562 safe_deallocate_a(kx)
563 safe_deallocate_a(ky)
564 safe_deallocate_a(kz)
565
566 call messages_write("done")
567 call messages_new_line()
568 call messages_info(namespace=namespace)
570
572 end subroutine pes_out_interpolator_init
573
574 ! ---------------------------------------------------------
576 ! ---------------------------------------------------------
577 subroutine pes_out_interpolator_end(cube_f, interp)
578 real(real64), allocatable, intent(inout) :: cube_f(:)
579 type(qshep_t), intent(inout) :: interp
580
582
583 call qshep_end(interp)
584
585 safe_deallocate_a(cube_f)
586
588 end subroutine pes_out_interpolator_end
589
590end module pes_out_oct_m
591
592
593!! Local Variables:
594!! mode: f90
595!! coding: utf-8
596!! End:
This is the common interface to a sorting routine. It performs the shell algorithm,...
Definition: sort.F90:149
double sin(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
double atan2(double __y, double __x) __attribute__((__nothrow__
subroutine, public dcube_function_alloc_rs(cube, cf, in_device, force_alloc)
Allocates locally the real space grid, if PFFT library is not used. Otherwise, it assigns the PFFT re...
subroutine, public dcube_function_free_rs(cube, cf)
Deallocates the real space grid.
subroutine, public cube_init(cube, nn, namespace, space, spacing, coord_system, fft_type, fft_library, dont_optimize, nn_out, mpi_grp, need_partition, tp_enlarge, blocksize)
Definition: cube.F90:202
subroutine, public cube_end(cube)
Definition: cube.F90:378
type(debug_t), save, public debug
Definition: debug.F90:154
real(real64), parameter, public m_two
Definition: global.F90:189
real(real64), parameter, public m_zero
Definition: global.F90:187
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:185
real(real64), parameter, public m_epsilon
Definition: global.F90:203
real(real64), parameter, public m_one
Definition: global.F90:188
subroutine, public dout_cf_netcdf(filename, ierr, cf, cube, space, spacing, transpose, unit, namespace)
Writes a cube_function in netcdf format.
Definition: io.F90:114
subroutine, public io_close(iunit, grp)
Definition: io.F90:468
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:395
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
subroutine, public generate_rotation_matrix(R, ff, tt)
Generates a rotation matrix R to rotate a vector f to t.
Definition: math.F90:1027
subroutine, public messages_info(no_lines, iunit, verbose_limit, stress, all_nodes, namespace)
Definition: messages.F90:624
subroutine, public messages_new_line()
Definition: messages.F90:1146
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:420
subroutine, public pes_out_velocity_map_cut(namespace, pesK, file, ll, dim, pol, dir, integrate, pos, Lk, pmesh)
Definition: pes_out.F90:305
subroutine pes_out_interpolator_init(namespace, pesK, Lk, ll, dim, cube_f, interp, pmesh)
Qshep interpolation helper function initialization. Generates the linearized version of pesK (cube_f)...
Definition: pes_out.F90:570
integer, parameter integrate_ky
Definition: pes_out.F90:147
subroutine, public pes_out_velocity_map(pesK, file, namespace, space, Lk, ll, how, spacing, coord_system, pmesh)
Definition: pes_out.F90:162
integer, parameter integrate_r
Definition: pes_out.F90:147
integer, parameter integrate_phi
Definition: pes_out.F90:147
integer, parameter integrate_theta
Definition: pes_out.F90:147
subroutine pes_out_interpolator_end(cube_f, interp)
Destroy the interpolation objects.
Definition: pes_out.F90:671
integer, parameter integrate_kx
Definition: pes_out.F90:147
subroutine, public pes_out_arpes_cut(namespace, arpes, file, dim, ll, pmesh, Ekin)
Definition: pes_out.F90:229
integer, parameter integrate_kz
Definition: pes_out.F90:147
subroutine, public qshep_end(interp)
Definition: qshep.F90:440
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:117
character(len=80) function, public str_center(s_in, l_in)
puts space around string, so that it is centered
Definition: string.F90:197
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
type(unit_t), public unit_one
some special units required for particular quantities
subroutine, public dvtk_out_cf_structured(filename, namespace, fieldname, ierr, cf_in, cube, unit, points, ascii)
Writes a mesh_function in VTK legacy format for structured grid data (i.e. curvilinear coordinates)....
Definition: vtk.F90:353
subroutine, public dvtk_out_cf(filename, namespace, fieldname, ierr, cf_in, cube, spacing, unit)
Writes a cube_function in VTK legacy format see http:
Definition: vtk.F90:215
int true(void)