22#include "../base/base_uses.f90"
30 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'realspace_grid_cube'
31 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
32 LOGICAL,
PRIVATE :: parses_linebreaks = .false., &
51 SUBROUTINE pw_to_cube(pw, unit_nr, title, particles_r, particles_z, particles_zeff, &
52 stride, max_file_size_mb, zero_tails, silent, mpi_io)
54 INTEGER,
INTENT(IN) :: unit_nr
55 CHARACTER(*),
INTENT(IN),
OPTIONAL :: title
56 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
57 OPTIONAL :: particles_r
58 INTEGER,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_z
59 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_zeff
60 INTEGER,
DIMENSION(:),
OPTIONAL,
POINTER :: stride
61 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: max_file_size_mb
62 LOGICAL,
INTENT(IN),
OPTIONAL :: zero_tails, silent, mpi_io
64 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_to_cube'
65 INTEGER,
PARAMETER :: entry_len = 13, num_entries_line = 6
67 INTEGER :: checksum, dest, handle, i, i1, i2, i3, iat, ip, l1, l2, l3, msglen, my_rank, &
68 my_stride(3), np, num_linebreak, num_pe, rank(2), size_of_z, source, tag, u1, u2, u3
69 LOGICAL :: be_silent, my_zero_tails, parallel_write
70 REAL(kind=
dp) :: compression_factor, my_max_file_size_mb
71 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buf
75 CALL timeset(routinen, handle)
77 my_zero_tails = .false.
79 parallel_write = .false.
80 my_max_file_size_mb = 0.0_dp
81 IF (
PRESENT(zero_tails)) my_zero_tails = zero_tails
82 IF (
PRESENT(silent)) be_silent = silent
83 IF (
PRESENT(mpi_io)) parallel_write = mpi_io
84 IF (
PRESENT(max_file_size_mb)) my_max_file_size_mb = max_file_size_mb
85 cpassert(my_max_file_size_mb >= 0)
88 IF (
PRESENT(stride))
THEN
89 IF (
SIZE(stride) /= 1 .AND.
SIZE(stride) /= 3) &
90 CALL cp_abort(__location__,
"STRIDE keyword can accept only 1 "// &
91 "(the same for X,Y,Z) or 3 values. Correct your input file.")
92 IF (
SIZE(stride) == 1)
THEN
94 my_stride(i) = stride(1)
97 my_stride = stride(1:3)
101 IF (my_max_file_size_mb > 0)
THEN
103 compression_factor = 1.3e-05_dp*product(real(pw%pw_grid%npts,
dp))/max_file_size_mb
104 my_stride(:) = int(compression_factor**(1.0/3.0)) + 1
107 cpassert(my_stride(1) > 0)
108 cpassert(my_stride(2) > 0)
109 cpassert(my_stride(3) > 0)
111 IF (.NOT. parallel_write)
THEN
112 IF (unit_nr > 0)
THEN
115 WRITE (unit_nr,
'(a11)')
"-Quickstep-"
116 IF (
PRESENT(title))
THEN
117 WRITE (unit_nr, *) trim(title)
119 WRITE (unit_nr, *)
"No Title"
122 cpassert(
PRESENT(particles_z) .EQV.
PRESENT(particles_r))
124 IF (
PRESENT(particles_z))
THEN
125 cpassert(
SIZE(particles_z) ==
SIZE(particles_r, dim=2))
128 np = min(99999,
SIZE(particles_z))
131 WRITE (unit_nr,
'(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp
133 WRITE (unit_nr,
'(I5,3f12.6)') (pw%pw_grid%npts(1) + my_stride(1) - 1)/my_stride(1), &
134 pw%pw_grid%dh(1, 1)*real(my_stride(1),
dp), pw%pw_grid%dh(2, 1)*real(my_stride(1),
dp), &
135 pw%pw_grid%dh(3, 1)*real(my_stride(1),
dp)
136 WRITE (unit_nr,
'(I5,3f12.6)') (pw%pw_grid%npts(2) + my_stride(2) - 1)/my_stride(2), &
137 pw%pw_grid%dh(1, 2)*real(my_stride(2),
dp), pw%pw_grid%dh(2, 2)*real(my_stride(2),
dp), &
138 pw%pw_grid%dh(3, 2)*real(my_stride(2),
dp)
139 WRITE (unit_nr,
'(I5,3f12.6)') (pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(3), &
140 pw%pw_grid%dh(1, 3)*real(my_stride(3),
dp), pw%pw_grid%dh(2, 3)*real(my_stride(3),
dp), &
141 pw%pw_grid%dh(3, 3)*real(my_stride(3),
dp)
143 IF (
PRESENT(particles_z))
THEN
144 IF (
PRESENT(particles_zeff))
THEN
146 WRITE (unit_nr,
'(I5,4f12.6)') particles_z(iat), particles_zeff(iat), particles_r(:, iat)
150 WRITE (unit_nr,
'(I5,4f12.6)') particles_z(iat), 0._dp, particles_r(:, iat)
157 l1 = pw%pw_grid%bounds(1, 1)
158 l2 = pw%pw_grid%bounds(1, 2)
159 l3 = pw%pw_grid%bounds(1, 3)
160 u1 = pw%pw_grid%bounds(2, 1)
161 u2 = pw%pw_grid%bounds(2, 2)
162 u3 = pw%pw_grid%bounds(2, 3)
164 ALLOCATE (buf(l3:u3))
166 my_rank = pw%pw_grid%para%group%mepos
167 gid = pw%pw_grid%para%group
168 num_pe = pw%pw_grid%para%group%num_pe
174 IF (unit_nr > 0) checksum = 1
176 CALL gid%sum(checksum)
177 cpassert(checksum == 1)
179 CALL gid%maxloc(rank)
180 cpassert(rank(1) > 0)
183 DO i1 = l1, u1, my_stride(1)
184 DO i2 = l2, u2, my_stride(2)
188 DO ip = 0, num_pe - 1
189 IF (pw%pw_grid%para%bo(1, 1, ip, 1) <= i1 - l1 + 1 .AND. pw%pw_grid%para%bo(2, 1, ip, 1) >= i1 - l1 + 1 .AND. &
190 pw%pw_grid%para%bo(1, 2, ip, 1) <= i2 - l2 + 1 .AND. pw%pw_grid%para%bo(2, 2, ip, 1) >= i2 - l2 + 1)
THEN
198 IF (source == dest)
THEN
199 IF (my_rank == source)
THEN
200 buf(:) = pw%array(i1, i2, :)
203 IF (my_rank == source)
THEN
204 buf(:) = pw%array(i1, i2, :)
205 CALL gid%send(buf, dest, tag)
207 IF (my_rank == dest)
THEN
208 CALL gid%recv(buf, source, tag)
212 IF (unit_nr > 0)
THEN
213 IF (my_zero_tails)
THEN
215 IF (buf(i3) < 1.e-7_dp) buf(i3) = 0.0_dp
218 WRITE (unit_nr,
'(6E13.5)') (buf(i3), i3=l3, u3, my_stride(3))
234 size_of_z = ceiling(real(pw%pw_grid%bounds(2, 3) - pw%pw_grid%bounds(1, 3) + 1,
dp)/real(my_stride(3),
dp))
235 num_linebreak = size_of_z/num_entries_line
236 IF (
modulo(size_of_z, num_entries_line) /= 0) &
237 num_linebreak = num_linebreak + 1
239 CALL mp_unit%set_handle(unit_nr)
240 CALL pw_to_cube_parallel(pw, mp_unit, title, particles_r, particles_z, particles_zeff, &
241 my_stride, my_zero_tails, msglen)
244 CALL timestop(handle)
260 SUBROUTINE cube_to_pw(grid, filename, scaling, parallel_read, silent)
263 CHARACTER(len=*),
INTENT(in) :: filename
264 REAL(kind=
dp),
INTENT(in) :: scaling
265 LOGICAL,
INTENT(in) :: parallel_read
266 LOGICAL,
INTENT(in),
OPTIONAL :: silent
268 CHARACTER(len=*),
PARAMETER :: routinen =
'cube_to_pw'
269 INTEGER,
PARAMETER :: entry_len = 13, num_entries_line = 6
271 INTEGER :: extunit, handle, i, j, k, msglen, &
272 my_rank, nat, ndum, num_linebreak, &
273 num_pe, output_unit, size_of_z, tag
274 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, npoints, &
275 npoints_local, ubounds, ubounds_local
277 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buffer
278 REAL(kind=
dp),
DIMENSION(3) :: dr, rdum
283 CALL timeset(routinen, handle)
286 IF (
PRESENT(silent))
THEN
290 gid = grid%pw_grid%para%group
291 my_rank = grid%pw_grid%para%group%mepos
292 num_pe = grid%pw_grid%para%group%num_pe
295 lbounds_local = grid%pw_grid%bounds_local(1, :)
296 ubounds_local = grid%pw_grid%bounds_local(2, :)
297 size_of_z = ubounds_local(3) - lbounds_local(3) + 1
299 IF (.NOT. parallel_read)
THEN
300 npoints = grid%pw_grid%npts
301 lbounds = grid%pw_grid%bounds(1, :)
302 ubounds = grid%pw_grid%bounds(2, :)
305 dr(i) = grid%pw_grid%dh(i, i)
308 npoints_local = grid%pw_grid%npts_local
310 ALLOCATE (buffer(lbounds(3):ubounds(3)))
312 IF (my_rank == 0)
THEN
313 IF (output_unit > 0 .AND. .NOT. be_silent)
THEN
314 WRITE (output_unit, fmt=
"(/,T2,A,/,/,T2,A,/)")
"Reading the cube file: ", trim(filename)
319 file_form=
"FORMATTED", &
320 file_action=
"READ", &
327 READ (extunit, *) nat, rdum
329 READ (extunit, *) ndum, rdum
330 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
331 output_unit > 0)
THEN
332 WRITE (output_unit, *)
"Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
333 WRITE (output_unit, *)
"Restart from density | ", ndum,
" DIFFERS FROM ", npoints(i)
334 WRITE (output_unit, *)
"Restart from density | ", rdum,
" DIFFERS FROM ", dr(i)
344 DO i = lbounds(1), ubounds(1)
345 DO j = lbounds(2), ubounds(2)
346 IF (my_rank == 0)
THEN
347 READ (extunit, *) (buffer(k), k=lbounds(3), ubounds(3))
349 CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0)
352 IF ((lbounds_local(1) <= i) .AND. (i <= ubounds_local(1)) .AND. (lbounds_local(2) <= j) &
353 .AND. (j <= ubounds_local(2)))
THEN
355 grid%array(i, j, lbounds(3):ubounds(3)) = buffer(lbounds(3):ubounds(3))*scaling
361 IF (my_rank == 0)
CALL close_file(unit_number=extunit)
370 num_linebreak = size_of_z/num_entries_line
371 IF (
modulo(size_of_z, num_entries_line) /= 0) &
372 num_linebreak = num_linebreak + 1
374 CALL cube_to_pw_parallel(grid, filename, scaling, msglen, silent=silent)
377 CALL timestop(handle)
392 SUBROUTINE cube_to_pw_parallel(grid, filename, scaling, msglen, silent)
395 CHARACTER(len=*),
INTENT(in) :: filename
396 REAL(kind=
dp),
INTENT(in) :: scaling
397 INTEGER,
INTENT(in) :: msglen
398 LOGICAL,
INTENT(in),
OPTIONAL :: silent
400 INTEGER,
PARAMETER :: entry_len = 13, num_entries_line = 6
402 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, npoints, &
403 npoints_local, ubounds, ubounds_local
404 INTEGER,
ALLOCATABLE,
DIMENSION(:),
TARGET :: blocklengths
405 INTEGER(kind=file_offset),
ALLOCATABLE, &
406 DIMENSION(:),
TARGET :: displacements
407 INTEGER(kind=file_offset) :: bof
408 INTEGER :: extunit_handle, i, ientry, islice, j, k, my_rank, nat, ndum, nslices, num_pe, &
409 offset_global, output_unit, pos, readstat, size_of_z, tag
410 CHARACTER(LEN=msglen),
ALLOCATABLE,
DIMENSION(:) ::
readbuffer
411 CHARACTER(LEN=msglen) :: tmp
412 LOGICAL :: be_silent, should_read(2)
413 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buffer
414 REAL(kind=
dp),
DIMENSION(3) :: dr, rdum
422 IF (
PRESENT(silent))
THEN
427 gid = grid%pw_grid%para%group
428 my_rank = grid%pw_grid%para%group%mepos
429 num_pe = grid%pw_grid%para%group%num_pe
433 dr(i) = grid%pw_grid%dh(i, i)
436 npoints = grid%pw_grid%npts
437 lbounds = grid%pw_grid%bounds(1, :)
438 ubounds = grid%pw_grid%bounds(2, :)
440 npoints_local = grid%pw_grid%npts_local
441 lbounds_local = grid%pw_grid%bounds_local(1, :)
442 ubounds_local = grid%pw_grid%bounds_local(2, :)
443 size_of_z = ubounds_local(3) - lbounds_local(3) + 1
444 nslices = (ubounds_local(1) - lbounds_local(1) + 1)*(ubounds_local(2) - lbounds_local(2) + 1)
448 IF (my_rank == 0)
THEN
449 IF (output_unit > 0 .AND. .NOT. be_silent)
THEN
450 WRITE (output_unit, fmt=
"(/,T2,A,/,/,T2,A,/)")
"Reading the cube file: ", trim(filename)
455 file_form=
"FORMATTED", &
456 file_action=
"READ", &
457 file_access=
"STREAM", &
458 unit_number=extunit_handle)
462 READ (extunit_handle, *)
464 READ (extunit_handle, *) nat, rdum
466 READ (extunit_handle, *) ndum, rdum
467 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
468 output_unit > 0)
THEN
469 WRITE (output_unit, *)
"Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
470 WRITE (output_unit, *)
"Restart from density | ", ndum,
" DIFFERS FROM ", npoints(i)
471 WRITE (output_unit, *)
"Restart from density | ", rdum,
" DIFFERS FROM ", dr(i)
476 READ (extunit_handle, *)
479 INQUIRE (extunit_handle, pos=offset_global)
483 CALL gid%bcast(offset_global, grid%pw_grid%para%group%source)
485 CALL extunit%open(groupid=gid, filepath=filename, amode_status=
file_amode_rdonly)
487 ALLOCATE (displacements(nslices))
489 DO i = lbounds(1), ubounds(1)
490 should_read(:) = .true.
491 IF (i < lbounds_local(1))
THEN
492 should_read(1) = .false.
493 ELSE IF (i > ubounds_local(1))
THEN
496 DO j = lbounds(2), ubounds(2)
497 should_read(2) = .true.
498 IF (j < lbounds_local(2) .OR. j > ubounds_local(2))
THEN
499 should_read(2) = .false.
501 IF (all(should_read .EQV. .true.))
THEN
502 IF (islice > nslices) cpabort(
"Index out of bounds.")
503 displacements(islice) = bof
511 ALLOCATE (blocklengths(nslices))
512 blocklengths(:) = msglen
520 CALL extunit%read_all(msglen, nslices,
readbuffer, mp_file_desc)
526 ALLOCATE (buffer(lbounds(3):ubounds(3)))
530 READ (
readbuffer(1), *, iostat=readstat) (buffer(k), k=lbounds(3), ubounds(3))
531 IF (readstat == 0)
THEN
532 parses_linebreaks = .true.
534 parses_linebreaks = .false.
539 DO islice = 1, nslices
540 IF (parses_linebreaks)
THEN
543 READ (
readbuffer(islice), *) (buffer(k), k=lbounds(3), ubounds(3))
549 DO k = lbounds_local(3), ubounds_local(3)
550 IF (
modulo(ientry, num_entries_line) == 0 .OR. k == ubounds_local(3))
THEN
552 READ (tmp(pos:pos + (entry_len - 2)),
'(E12.5)') buffer(k)
553 pos = pos + (entry_len + 1)
555 READ (tmp(pos:pos + (entry_len - 1)),
'(E13.5)') buffer(k)
556 pos = pos + entry_len
562 grid%array(i, j, lbounds(3):ubounds(3)) = scaling*buffer(lbounds(3):ubounds(3))
564 IF (j > ubounds_local(2))
THEN
570 DEALLOCATE (blocklengths, displacements)
571 IF (debug_this_module)
THEN
574 IF (my_rank == 0)
THEN
575 IF (output_unit > 0 .AND. .NOT. be_silent)
THEN
576 WRITE (output_unit, fmt=
"(/,T2,A,/,/,T2,A)")
"Reading the cube file: ", filename
581 file_form=
"FORMATTED", &
582 file_action=
"READ", &
583 unit_number=extunit_handle)
587 READ (extunit_handle, *)
589 READ (extunit_handle, *) nat, rdum
591 READ (extunit_handle, *) ndum, rdum
592 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
593 output_unit > 0)
THEN
594 WRITE (output_unit, *)
"Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
595 WRITE (output_unit, *)
"Restart from density | ", ndum,
" DIFFERS FROM ", npoints(i)
596 WRITE (output_unit, *)
"Restart from density | ", rdum,
" DIFFERS FROM ", dr(i)
601 READ (extunit_handle, *)
606 DO i = lbounds(1), ubounds(1)
607 DO j = lbounds(2), ubounds(2)
608 IF (my_rank == 0)
THEN
609 READ (extunit_handle, *) (buffer(k), k=lbounds(3), ubounds(3))
611 CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0)
614 IF ((lbounds_local(1) <= i) .AND. (i <= ubounds_local(1)) .AND. (lbounds_local(2) <= j) &
615 .AND. (j <= ubounds_local(2)))
THEN
617 IF (any(grid%array(i, j, lbounds(3):ubounds(3)) /= buffer(lbounds(3):ubounds(3))*scaling)) &
618 CALL cp_abort(__location__, &
619 "Error in parallel read of input cube file.")
625 IF (my_rank == 0)
CALL close_file(unit_number=extunit_handle)
631 END SUBROUTINE cube_to_pw_parallel
647 SUBROUTINE pw_to_cube_parallel(grid, unit_nr, title, particles_r, particles_z, particles_zeff, &
648 stride, zero_tails, msglen)
652 CHARACTER(*),
INTENT(IN),
OPTIONAL :: title
653 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
654 OPTIONAL :: particles_r
655 INTEGER,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_z
656 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_zeff
657 INTEGER,
INTENT(IN) :: stride(3)
658 LOGICAL,
INTENT(IN) :: zero_tails
659 INTEGER,
INTENT(IN) :: msglen
661 INTEGER,
PARAMETER :: entry_len = 13, header_len = 41, &
662 header_len_z = 53, num_entries_line = 6
664 CHARACTER(LEN=entry_len) :: value
665 CHARACTER(LEN=header_len) ::
header
666 CHARACTER(LEN=header_len_z) :: header_z
667 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, ubounds, &
669 INTEGER,
ALLOCATABLE,
DIMENSION(:),
TARGET :: blocklengths
670 INTEGER(kind=file_offset),
ALLOCATABLE, &
671 DIMENSION(:),
TARGET :: displacements
672 INTEGER(kind=file_offset) :: bof
673 INTEGER :: counter, i, islice, j, k, last_z, &
674 my_rank, np, nslices, size_of_z
675 CHARACTER(LEN=msglen),
ALLOCATABLE,
DIMENSION(:) ::
writebuffer
676 CHARACTER(LEN=msglen) :: tmp
677 LOGICAL :: should_write(2)
682 gid = grid%pw_grid%para%group
683 my_rank = grid%pw_grid%para%group%mepos
686 lbounds = grid%pw_grid%bounds(1, :)
687 ubounds = grid%pw_grid%bounds(2, :)
688 lbounds_local = grid%pw_grid%bounds_local(1, :)
689 ubounds_local = grid%pw_grid%bounds_local(2, :)
691 size_of_z = ceiling(real(ubounds_local(3) - lbounds_local(3) + 1,
dp)/real(stride(3),
dp))
693 DO i = lbounds(1), ubounds(1), stride(1)
694 should_write(:) = .true.
695 IF (i < lbounds_local(1))
THEN
696 should_write(1) = .false.
697 ELSE IF (i > ubounds_local(1))
THEN
700 DO j = lbounds(2), ubounds(2), stride(2)
701 should_write(2) = .true.
702 IF (j < lbounds_local(2) .OR. j > ubounds_local(2))
THEN
703 should_write(2) = .false.
705 IF (all(should_write .EQV. .true.))
THEN
711 DO k = lbounds(3), ubounds(3), stride(3)
712 IF (k + stride(3) > ubounds(3)) last_z = k
716 CALL unit_nr%get_position(bof)
718 IF (my_rank == 0)
THEN
721 CALL unit_nr%write_at(bof,
"-Quickstep-"//new_line(
"C"))
723 IF (
PRESENT(title))
THEN
724 CALL unit_nr%write_at(bof, trim(title)//new_line(
"C"))
727 CALL unit_nr%write_at(bof,
"No Title"//new_line(
"C"))
731 cpassert(
PRESENT(particles_z) .EQV.
PRESENT(particles_r))
733 IF (
PRESENT(particles_z))
THEN
734 cpassert(
SIZE(particles_z) ==
SIZE(particles_r, dim=2))
737 np = min(99999,
SIZE(particles_z))
740 WRITE (
header,
'(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp
741 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
744 WRITE (
header,
'(I5,3f12.6)') (grid%pw_grid%npts(1) + stride(1) - 1)/stride(1), &
745 grid%pw_grid%dh(1, 1)*real(stride(1),
dp), grid%pw_grid%dh(2, 1)*real(stride(1),
dp), &
746 grid%pw_grid%dh(3, 1)*real(stride(1),
dp)
747 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
750 WRITE (
header,
'(I5,3f12.6)') (grid%pw_grid%npts(2) + stride(2) - 1)/stride(2), &
751 grid%pw_grid%dh(1, 2)*real(stride(2),
dp), grid%pw_grid%dh(2, 2)*real(stride(2),
dp), &
752 grid%pw_grid%dh(3, 2)*real(stride(2),
dp)
753 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
756 WRITE (
header,
'(I5,3f12.6)') (grid%pw_grid%npts(3) + stride(3) - 1)/stride(3), &
757 grid%pw_grid%dh(1, 3)*real(stride(3),
dp), grid%pw_grid%dh(2, 3)*real(stride(3),
dp), &
758 grid%pw_grid%dh(3, 3)*real(stride(3),
dp)
759 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
762 IF (
PRESENT(particles_z))
THEN
763 IF (
PRESENT(particles_zeff))
THEN
765 WRITE (header_z,
'(I5,4f12.6)') particles_z(i), particles_zeff(i), particles_r(:, i)
766 CALL unit_nr%write_at(bof, header_z//new_line(
"C"))
771 WRITE (header_z,
'(I5,4f12.6)') particles_z(i), 0._dp, particles_r(:, i)
772 CALL unit_nr%write_at(bof, header_z//new_line(
"C"))
779 CALL gid%bcast(bof, grid%pw_grid%para%group%source)
782 ALLOCATE (displacements(nslices))
786 DO i = lbounds(1), ubounds(1), stride(1)
787 should_write(:) = .true.
788 IF (i < lbounds_local(1))
THEN
789 should_write(1) = .false.
790 ELSE IF (i > ubounds_local(1))
THEN
793 DO j = lbounds(2), ubounds(2), stride(2)
794 should_write(2) = .true.
795 IF (j < lbounds_local(2) .OR. j > ubounds_local(2))
THEN
796 should_write(2) = .false.
798 IF (all(should_write .EQV. .true.))
THEN
799 IF (islice > nslices) cpabort(
"Index out of bounds.")
800 displacements(islice) = bof
803 DO k = lbounds(3), ubounds(3), stride(3)
804 IF (zero_tails .AND. grid%array(i, j, k) < 1.e-7_dp)
THEN
805 WRITE (
value,
'(E13.5)') 0.0_dp
807 WRITE (
value,
'(E13.5)') grid%array(i, j, k)
809 tmp = trim(tmp)//trim(
value)
810 counter = counter + 1
811 IF (
modulo(counter, num_entries_line) == 0 .OR. k == last_z) &
812 tmp = trim(tmp)//new_line(
'C')
823 ALLOCATE (blocklengths(nslices))
824 blocklengths(:) = msglen
833 CALL unit_nr%write_all(msglen, nslices,
writebuffer, mp_desc)
837 DEALLOCATE (blocklengths, displacements)
839 END SUBROUTINE pw_to_cube_parallel
853 INTEGER,
INTENT(IN) :: unit_nr
854 INTEGER,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: stride
857 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_to_simple_volumetric'
859 INTEGER :: checksum, dest, handle, i, i1, i2, i3, &
860 ip, l1, l2, l3, my_rank, my_stride(3), &
861 ngrids, npoints, num_pe, rank(2), &
862 source, tag, u1, u2, u3
864 REAL(kind=
dp) :: x, y, z
865 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buf, buf2
868 CALL timeset(routinen, handle)
872 IF (
PRESENT(pw2)) double = .true.
875 IF (
PRESENT(stride))
THEN
876 IF (
SIZE(stride) /= 1 .AND.
SIZE(stride) /= 3) &
877 CALL cp_abort(__location__,
"STRIDE keyword can accept only 1 "// &
878 "(the same for X,Y,Z) or 3 values. Correct your input file.")
879 IF (
SIZE(stride) == 1)
THEN
881 my_stride(i) = stride(1)
884 my_stride = stride(1:3)
886 cpassert(my_stride(1) > 0)
887 cpassert(my_stride(2) > 0)
888 cpassert(my_stride(3) > 0)
892 l1 = pw%pw_grid%bounds(1, 1)
893 l2 = pw%pw_grid%bounds(1, 2)
894 l3 = pw%pw_grid%bounds(1, 3)
895 u1 = pw%pw_grid%bounds(2, 1)
896 u2 = pw%pw_grid%bounds(2, 2)
897 u3 = pw%pw_grid%bounds(2, 3)
901 IF (double) ngrids = 2
902 npoints = ((pw%pw_grid%npts(1) + my_stride(1) - 1)/my_stride(1))* &
903 ((pw%pw_grid%npts(2) + my_stride(2) - 1)/my_stride(1))* &
904 ((pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(1))
905 IF (unit_nr > 1)
WRITE (unit_nr,
'(I7,I5)') npoints, ngrids
907 ALLOCATE (buf(l3:u3))
908 IF (double)
ALLOCATE (buf2(l3:u3))
910 my_rank = pw%pw_grid%para%group%mepos
911 gid = pw%pw_grid%para%group
912 num_pe = pw%pw_grid%para%group%num_pe
918 IF (unit_nr > 0) checksum = 1
920 CALL gid%sum(checksum)
921 cpassert(checksum == 1)
923 CALL gid%maxloc(rank)
924 cpassert(rank(1) > 0)
927 DO i1 = l1, u1, my_stride(1)
928 DO i2 = l2, u2, my_stride(2)
932 DO ip = 0, num_pe - 1
933 IF (pw%pw_grid%para%bo(1, 1, ip, 1) <= i1 - l1 + 1 .AND. pw%pw_grid%para%bo(2, 1, ip, 1) >= i1 - l1 + 1 .AND. &
934 pw%pw_grid%para%bo(1, 2, ip, 1) <= i2 - l2 + 1 .AND. pw%pw_grid%para%bo(2, 2, ip, 1) >= i2 - l2 + 1)
THEN
942 IF (source == dest)
THEN
943 IF (my_rank == source)
THEN
944 buf(:) = pw%array(i1, i2, :)
945 IF (double) buf2(:) = pw2%array(i1, i2, :)
948 IF (my_rank == source)
THEN
949 buf(:) = pw%array(i1, i2, :)
950 CALL gid%send(buf, dest, tag)
952 buf2(:) = pw2%array(i1, i2, :)
953 CALL gid%send(buf2, dest, tag)
956 IF (my_rank == dest)
THEN
957 CALL gid%recv(buf, source, tag)
958 IF (double)
CALL gid%recv(buf2, source, tag)
962 IF (.NOT. double)
THEN
963 DO i3 = l3, u3, my_stride(3)
964 x = pw%pw_grid%dh(1, 1)*i1 + &
965 pw%pw_grid%dh(2, 1)*i2 + &
966 pw%pw_grid%dh(3, 1)*i3
968 y = pw%pw_grid%dh(1, 2)*i1 + &
969 pw%pw_grid%dh(2, 2)*i2 + &
970 pw%pw_grid%dh(3, 2)*i3
972 z = pw%pw_grid%dh(1, 3)*i1 + &
973 pw%pw_grid%dh(2, 3)*i2 + &
974 pw%pw_grid%dh(3, 3)*i3
976 IF (unit_nr > 0)
THEN
977 WRITE (unit_nr,
'(6E13.5, 6E13.5, 6E13.5, 6E13.5)') x, y, z, buf(i3)
983 DO i3 = l3, u3, my_stride(3)
984 x = pw%pw_grid%dh(1, 1)*i1 + &
985 pw%pw_grid%dh(2, 1)*i2 + &
986 pw%pw_grid%dh(3, 1)*i3
988 y = pw%pw_grid%dh(1, 2)*i1 + &
989 pw%pw_grid%dh(2, 2)*i2 + &
990 pw%pw_grid%dh(3, 2)*i3
992 z = pw%pw_grid%dh(1, 3)*i1 + &
993 pw%pw_grid%dh(2, 3)*i2 + &
994 pw%pw_grid%dh(3, 3)*i3
996 IF (unit_nr > 0)
THEN
997 WRITE (unit_nr,
'(6E13.5, 6E13.5, 6E13.5, 6E13.5)') x, y, z, buf(i3), buf2(i3)
1015 IF (double)
DEALLOCATE (buf2)
1017 CALL timestop(handle)
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Utility routines to open and close files. Tracking of preconnections.
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Defines the basic variable types.
integer, parameter, public dp
Interface to the message passing library MPI.
type(mp_file_descriptor_type) function, public mp_file_type_hindexed_make_chv(count, lengths, displs)
Creates an indexed MPI type for arrays of strings using bytes for spacing (hindexed type)
subroutine, public mp_file_type_free(type_descriptor)
Releases the type used for MPI I/O.
integer, parameter, public mpi_character_size
integer, parameter, public file_offset
integer, parameter, public file_amode_rdonly
subroutine, public mp_file_type_set_view_chv(fh, offset, type_descriptor)
Uses a previously created indexed MPI character type to tell the MPI processes how to partition (set_...
integer, parameter, public pw_mode_local
Generate Gaussian cube files.
subroutine, public cube_to_pw(grid, filename, scaling, parallel_read, silent)
Computes the external density on the grid hacked from external_read_density.
subroutine, public pw_to_simple_volumetric(pw, unit_nr, stride, pw2)
Prints a simple grid file: X Y Z value.
subroutine, public pw_to_cube(pw, unit_nr, title, particles_r, particles_z, particles_zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
void writebuffer(int *psockfd, char *data, int *plen)
Writes to a socket.
void readbuffer(int *psockfd, char *data, int *plen)
Reads from a socket.