22#include "../base/base_uses.f90"
31 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'realspace_grid_cube'
32 INTEGER,
PARAMETER,
PRIVATE :: cube_entry_len = 13, &
33 cube_num_entries_line = 6
34 INTEGER,
PARAMETER,
PRIVATE :: cube_line_len = cube_entry_len*cube_num_entries_line
37 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .false.
47 CHARACTER(LEN=*),
INTENT(IN) :: values
48 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: buffer
50 CHARACTER(LEN=cube_entry_len) :: value
51 INTEGER :: i, pos, readstat
53 READ (values, *, iostat=readstat) buffer
54 IF (readstat == 0)
RETURN
57 DO i = 1,
SIZE(buffer)
58 IF (pos + cube_entry_len - 1 > len(values)) cpabort(
"Unexpected end of cube data.")
59 value = values(pos:pos + cube_entry_len - 1)
60 READ (
value,
'(E13.5)', iostat=readstat) buffer(i)
61 IF (readstat /= 0) cpabort(
"Bad value while reading cube data.")
62 pos = pos + cube_entry_len
63 IF (
modulo(i, cube_num_entries_line) == 0)
THEN
64 IF (pos <= len(values))
THEN
65 IF (values(pos:pos) == new_line(
'C')) pos = pos + 1
86 SUBROUTINE pw_to_cube(pw, unit_nr, title, particles_r, particles_z, particles_zeff, &
87 stride, max_file_size_mb, zero_tails, silent, mpi_io)
89 INTEGER,
INTENT(IN) :: unit_nr
90 CHARACTER(*),
INTENT(IN),
OPTIONAL :: title
91 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
92 OPTIONAL :: particles_r
93 INTEGER,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_z
94 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_zeff
95 INTEGER,
DIMENSION(:),
OPTIONAL,
POINTER :: stride
96 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: max_file_size_mb
97 LOGICAL,
INTENT(IN),
OPTIONAL :: zero_tails, silent, mpi_io
99 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_to_cube'
100 INTEGER,
PARAMETER :: entry_len = 13, num_entries_line = 6
102 INTEGER :: checksum, dest, handle, i, i1, i2, i3, iat, ip, l1, l2, l3, msglen, my_rank, &
103 my_stride(3), np, num_linebreak, num_pe, rank(2), size_of_z, source, tag, u1, u2, u3
104 LOGICAL :: be_silent, my_zero_tails, parallel_write
105 REAL(kind=
dp) :: compression_factor, my_max_file_size_mb
106 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buf
110 CALL timeset(routinen, handle)
112 my_zero_tails = .false.
114 parallel_write = .false.
115 my_max_file_size_mb = 0.0_dp
116 IF (
PRESENT(zero_tails)) my_zero_tails = zero_tails
117 IF (
PRESENT(silent)) be_silent = silent
118 IF (
PRESENT(mpi_io)) parallel_write = mpi_io
119 IF (
PRESENT(max_file_size_mb)) my_max_file_size_mb = max_file_size_mb
120 cpassert(my_max_file_size_mb >= 0)
123 IF (
PRESENT(stride))
THEN
124 IF (
SIZE(stride) /= 1 .AND.
SIZE(stride) /= 3) &
125 CALL cp_abort(__location__,
"STRIDE keyword can accept only 1 "// &
126 "(the same for X,Y,Z) or 3 values. Correct your input file.")
127 IF (
SIZE(stride) == 1)
THEN
129 my_stride(i) = stride(1)
132 my_stride = stride(1:3)
136 IF (my_max_file_size_mb > 0)
THEN
138 compression_factor = 1.3e-05_dp*product(real(pw%pw_grid%npts,
dp))/max_file_size_mb
139 my_stride(:) = int(compression_factor**(1.0/3.0)) + 1
142 cpassert(my_stride(1) > 0)
143 cpassert(my_stride(2) > 0)
144 cpassert(my_stride(3) > 0)
146 IF (.NOT. parallel_write)
THEN
147 IF (unit_nr > 0)
THEN
150 WRITE (unit_nr,
'(a11)')
"-Quickstep-"
151 IF (
PRESENT(title))
THEN
152 WRITE (unit_nr, *) trim(title)
154 WRITE (unit_nr, *)
"No Title"
157 cpassert(
PRESENT(particles_z) .EQV.
PRESENT(particles_r))
159 IF (
PRESENT(particles_z))
THEN
160 cpassert(
SIZE(particles_z) ==
SIZE(particles_r, dim=2))
163 np = min(99999,
SIZE(particles_z))
166 WRITE (unit_nr,
'(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp
168 WRITE (unit_nr,
'(I5,3f12.6)') (pw%pw_grid%npts(1) + my_stride(1) - 1)/my_stride(1), &
169 pw%pw_grid%dh(1, 1)*real(my_stride(1),
dp), pw%pw_grid%dh(2, 1)*real(my_stride(1),
dp), &
170 pw%pw_grid%dh(3, 1)*real(my_stride(1),
dp)
171 WRITE (unit_nr,
'(I5,3f12.6)') (pw%pw_grid%npts(2) + my_stride(2) - 1)/my_stride(2), &
172 pw%pw_grid%dh(1, 2)*real(my_stride(2),
dp), pw%pw_grid%dh(2, 2)*real(my_stride(2),
dp), &
173 pw%pw_grid%dh(3, 2)*real(my_stride(2),
dp)
174 WRITE (unit_nr,
'(I5,3f12.6)') (pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(3), &
175 pw%pw_grid%dh(1, 3)*real(my_stride(3),
dp), pw%pw_grid%dh(2, 3)*real(my_stride(3),
dp), &
176 pw%pw_grid%dh(3, 3)*real(my_stride(3),
dp)
178 IF (
PRESENT(particles_z))
THEN
179 IF (
PRESENT(particles_zeff))
THEN
181 WRITE (unit_nr,
'(I5,4f12.6)') particles_z(iat), particles_zeff(iat), particles_r(:, iat)
185 WRITE (unit_nr,
'(I5,4f12.6)') particles_z(iat), 0._dp, particles_r(:, iat)
192 l1 = pw%pw_grid%bounds(1, 1)
193 l2 = pw%pw_grid%bounds(1, 2)
194 l3 = pw%pw_grid%bounds(1, 3)
195 u1 = pw%pw_grid%bounds(2, 1)
196 u2 = pw%pw_grid%bounds(2, 2)
197 u3 = pw%pw_grid%bounds(2, 3)
199 ALLOCATE (buf(l3:u3))
201 my_rank = pw%pw_grid%para%group%mepos
202 gid = pw%pw_grid%para%group
203 num_pe = pw%pw_grid%para%group%num_pe
209 IF (unit_nr > 0) checksum = 1
211 CALL gid%sum(checksum)
212 cpassert(checksum == 1)
214 CALL gid%maxloc(rank)
215 cpassert(rank(1) > 0)
218 DO i1 = l1, u1, my_stride(1)
219 DO i2 = l2, u2, my_stride(2)
223 DO ip = 0, num_pe - 1
224 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. &
225 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
233 IF (source == dest)
THEN
234 IF (my_rank == source)
THEN
235 buf(:) = pw%array(i1, i2, :)
238 IF (my_rank == source)
THEN
239 buf(:) = pw%array(i1, i2, :)
240 CALL gid%send(buf, dest, tag)
242 IF (my_rank == dest)
THEN
243 CALL gid%recv(buf, source, tag)
247 IF (unit_nr > 0)
THEN
248 IF (my_zero_tails)
THEN
250 IF (buf(i3) < 1.e-7_dp) buf(i3) = 0.0_dp
269 size_of_z = ceiling(real(pw%pw_grid%bounds(2, 3) - pw%pw_grid%bounds(1, 3) + 1,
dp)/real(my_stride(3),
dp))
270 num_linebreak = size_of_z/num_entries_line
271 IF (
modulo(size_of_z, num_entries_line) /= 0) &
272 num_linebreak = num_linebreak + 1
274 CALL mp_unit%set_handle(unit_nr)
275 CALL pw_to_cube_parallel(pw, mp_unit, title, particles_r, particles_z, particles_zeff, &
276 my_stride, my_zero_tails, msglen)
279 CALL timestop(handle)
295 SUBROUTINE cube_to_pw(grid, filename, scaling, parallel_read, silent)
298 CHARACTER(len=*),
INTENT(in) :: filename
299 REAL(kind=
dp),
INTENT(in) :: scaling
300 LOGICAL,
INTENT(in) :: parallel_read
301 LOGICAL,
INTENT(in),
OPTIONAL :: silent
303 CHARACTER(len=*),
PARAMETER :: routinen =
'cube_to_pw'
304 INTEGER,
PARAMETER :: entry_len = 13, num_entries_line = 6
306 CHARACTER(LEN=cube_line_len) :: value_line
307 INTEGER :: extunit, handle, i, j, k, last_z, &
308 msglen, my_rank, nat, ndum, &
309 num_linebreak, num_pe, output_unit, &
311 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, npoints, &
312 npoints_local, ubounds, ubounds_local
314 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buffer
315 REAL(kind=
dp),
DIMENSION(3) :: dr, rdum
320 CALL timeset(routinen, handle)
323 IF (
PRESENT(silent))
THEN
327 gid = grid%pw_grid%para%group
328 my_rank = grid%pw_grid%para%group%mepos
329 num_pe = grid%pw_grid%para%group%num_pe
332 lbounds_local = grid%pw_grid%bounds_local(1, :)
333 ubounds_local = grid%pw_grid%bounds_local(2, :)
334 size_of_z = ubounds_local(3) - lbounds_local(3) + 1
336 IF (.NOT. parallel_read)
THEN
337 npoints = grid%pw_grid%npts
338 lbounds = grid%pw_grid%bounds(1, :)
339 ubounds = grid%pw_grid%bounds(2, :)
342 dr(i) = grid%pw_grid%dh(i, i)
345 npoints_local = grid%pw_grid%npts_local
347 ALLOCATE (buffer(lbounds(3):ubounds(3)))
349 IF (my_rank == 0)
THEN
350 IF (output_unit > 0 .AND. .NOT. be_silent)
THEN
351 WRITE (output_unit, fmt=
"(/,T2,A,/,/,T2,A,/)")
"Reading the cube file: ", trim(filename)
356 file_form=
"FORMATTED", &
357 file_action=
"READ", &
364 READ (extunit, *) nat, rdum
366 READ (extunit, *) ndum, rdum
367 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
368 output_unit > 0)
THEN
369 WRITE (output_unit, *)
"Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
370 WRITE (output_unit, *)
"Restart from density | ", ndum,
" DIFFERS FROM ", npoints(i)
371 WRITE (output_unit, *)
"Restart from density | ", rdum,
" DIFFERS FROM ", dr(i)
381 DO i = lbounds(1), ubounds(1)
382 DO j = lbounds(2), ubounds(2)
383 IF (my_rank == 0)
THEN
384 DO k = lbounds(3), ubounds(3), cube_num_entries_line
385 last_z = min(k + cube_num_entries_line - 1, ubounds(3))
386 READ (extunit,
'(A)') value_line
390 CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0)
393 IF ((lbounds_local(1) <= i) .AND. (i <= ubounds_local(1)) .AND. (lbounds_local(2) <= j) &
394 .AND. (j <= ubounds_local(2)))
THEN
396 grid%array(i, j, lbounds(3):ubounds(3)) = buffer(lbounds(3):ubounds(3))*scaling
402 IF (my_rank == 0)
CALL close_file(unit_number=extunit)
411 num_linebreak = size_of_z/num_entries_line
412 IF (
modulo(size_of_z, num_entries_line) /= 0) &
413 num_linebreak = num_linebreak + 1
415 CALL cube_to_pw_parallel(grid, filename, scaling, msglen, silent=silent)
418 CALL timestop(handle)
433 SUBROUTINE cube_to_pw_parallel(grid, filename, scaling, msglen, silent)
436 CHARACTER(len=*),
INTENT(in) :: filename
437 REAL(kind=
dp),
INTENT(in) :: scaling
438 INTEGER,
INTENT(in) :: msglen
439 LOGICAL,
INTENT(in),
OPTIONAL :: silent
441 CHARACTER(LEN=cube_line_len) :: value_line
442 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, npoints, &
443 npoints_local, ubounds, ubounds_local
444 INTEGER,
ALLOCATABLE,
DIMENSION(:),
TARGET :: blocklengths
445 INTEGER(kind=file_offset),
ALLOCATABLE, &
446 DIMENSION(:),
TARGET :: displacements
447 INTEGER(kind=file_offset) :: bof
448 INTEGER :: extunit_handle, i, islice, j, k, last_z, &
449 my_rank, nat, ndum, nslices, num_pe, &
450 offset_global, output_unit, size_of_z, &
452 CHARACTER(LEN=msglen),
ALLOCATABLE,
DIMENSION(:) ::
readbuffer
453 LOGICAL :: be_silent, should_read(2)
454 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buffer
455 REAL(kind=
dp),
DIMENSION(3) :: dr, rdum
463 IF (
PRESENT(silent))
THEN
468 gid = grid%pw_grid%para%group
469 my_rank = grid%pw_grid%para%group%mepos
470 num_pe = grid%pw_grid%para%group%num_pe
474 dr(i) = grid%pw_grid%dh(i, i)
477 npoints = grid%pw_grid%npts
478 lbounds = grid%pw_grid%bounds(1, :)
479 ubounds = grid%pw_grid%bounds(2, :)
481 npoints_local = grid%pw_grid%npts_local
482 lbounds_local = grid%pw_grid%bounds_local(1, :)
483 ubounds_local = grid%pw_grid%bounds_local(2, :)
484 size_of_z = ubounds_local(3) - lbounds_local(3) + 1
485 nslices = (ubounds_local(1) - lbounds_local(1) + 1)*(ubounds_local(2) - lbounds_local(2) + 1)
489 IF (my_rank == 0)
THEN
490 IF (output_unit > 0 .AND. .NOT. be_silent)
THEN
491 WRITE (output_unit, fmt=
"(/,T2,A,/,/,T2,A,/)")
"Reading the cube file: ", trim(filename)
496 file_form=
"FORMATTED", &
497 file_action=
"READ", &
498 file_access=
"STREAM", &
499 unit_number=extunit_handle)
503 READ (extunit_handle, *)
505 READ (extunit_handle, *) nat, rdum
507 READ (extunit_handle, *) ndum, rdum
508 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
509 output_unit > 0)
THEN
510 WRITE (output_unit, *)
"Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
511 WRITE (output_unit, *)
"Restart from density | ", ndum,
" DIFFERS FROM ", npoints(i)
512 WRITE (output_unit, *)
"Restart from density | ", rdum,
" DIFFERS FROM ", dr(i)
517 READ (extunit_handle, *)
520 INQUIRE (extunit_handle, pos=offset_global)
524 CALL gid%bcast(offset_global, grid%pw_grid%para%group%source)
527 bof = offset_global - 1
528 CALL extunit%open(groupid=gid, filepath=filename, amode_status=
file_amode_rdonly)
530 ALLOCATE (displacements(nslices))
532 DO i = lbounds(1), ubounds(1)
533 should_read(:) = .true.
534 IF (i < lbounds_local(1))
THEN
535 should_read(1) = .false.
536 ELSE IF (i > ubounds_local(1))
THEN
539 DO j = lbounds(2), ubounds(2)
540 should_read(2) = .true.
541 IF (j < lbounds_local(2) .OR. j > ubounds_local(2))
THEN
542 should_read(2) = .false.
544 IF (all(should_read .EQV. .true.))
THEN
545 IF (islice > nslices) cpabort(
"Index out of bounds.")
546 displacements(islice) = bof
554 ALLOCATE (blocklengths(nslices))
555 blocklengths(:) = msglen
563 CALL extunit%read_all(msglen, nslices,
readbuffer, mp_file_desc)
569 ALLOCATE (buffer(lbounds(3):ubounds(3)))
571 DO islice = 1, nslices
574 grid%array(i, j, lbounds(3):ubounds(3)) = scaling*buffer(lbounds(3):ubounds(3))
576 IF (j > ubounds_local(2))
THEN
582 DEALLOCATE (blocklengths, displacements)
583 IF (debug_this_module)
THEN
586 IF (my_rank == 0)
THEN
587 IF (output_unit > 0 .AND. .NOT. be_silent)
THEN
588 WRITE (output_unit, fmt=
"(/,T2,A,/,/,T2,A)")
"Reading the cube file: ", filename
593 file_form=
"FORMATTED", &
594 file_action=
"READ", &
595 unit_number=extunit_handle)
599 READ (extunit_handle, *)
601 READ (extunit_handle, *) nat, rdum
603 READ (extunit_handle, *) ndum, rdum
604 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
605 output_unit > 0)
THEN
606 WRITE (output_unit, *)
"Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
607 WRITE (output_unit, *)
"Restart from density | ", ndum,
" DIFFERS FROM ", npoints(i)
608 WRITE (output_unit, *)
"Restart from density | ", rdum,
" DIFFERS FROM ", dr(i)
613 READ (extunit_handle, *)
618 DO i = lbounds(1), ubounds(1)
619 DO j = lbounds(2), ubounds(2)
620 IF (my_rank == 0)
THEN
621 DO k = lbounds(3), ubounds(3), cube_num_entries_line
622 last_z = min(k + cube_num_entries_line - 1, ubounds(3))
623 READ (extunit_handle,
'(A)') value_line
627 CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0)
630 IF ((lbounds_local(1) <= i) .AND. (i <= ubounds_local(1)) .AND. (lbounds_local(2) <= j) &
631 .AND. (j <= ubounds_local(2)))
THEN
633 IF (any(grid%array(i, j, lbounds(3):ubounds(3)) /= buffer(lbounds(3):ubounds(3))*scaling)) &
634 CALL cp_abort(__location__, &
635 "Error in parallel read of input cube file.")
641 IF (my_rank == 0)
CALL close_file(unit_number=extunit_handle)
647 END SUBROUTINE cube_to_pw_parallel
663 SUBROUTINE pw_to_cube_parallel(grid, unit_nr, title, particles_r, particles_z, particles_zeff, &
664 stride, zero_tails, msglen)
668 CHARACTER(*),
INTENT(IN),
OPTIONAL :: title
669 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
670 OPTIONAL :: particles_r
671 INTEGER,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_z
672 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_zeff
673 INTEGER,
INTENT(IN) :: stride(3)
674 LOGICAL,
INTENT(IN) :: zero_tails
675 INTEGER,
INTENT(IN) :: msglen
677 INTEGER,
PARAMETER :: entry_len = 13, header_len = 41, &
678 header_len_z = 53, num_entries_line = 6
680 CHARACTER(LEN=entry_len) :: value
681 CHARACTER(LEN=header_len) ::
header
682 CHARACTER(LEN=header_len_z) :: header_z
683 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, ubounds, &
685 INTEGER,
ALLOCATABLE,
DIMENSION(:),
TARGET :: blocklengths
686 INTEGER(kind=file_offset),
ALLOCATABLE, &
687 DIMENSION(:),
TARGET :: displacements
688 INTEGER(kind=file_offset) :: bof
689 INTEGER :: counter, i, islice, j, k, last_z, &
690 my_rank, np, nslices, size_of_z
691 CHARACTER(LEN=msglen),
ALLOCATABLE,
DIMENSION(:) ::
writebuffer
692 CHARACTER(LEN=msglen) :: tmp
693 LOGICAL :: should_write(2)
698 gid = grid%pw_grid%para%group
699 my_rank = grid%pw_grid%para%group%mepos
702 lbounds = grid%pw_grid%bounds(1, :)
703 ubounds = grid%pw_grid%bounds(2, :)
704 lbounds_local = grid%pw_grid%bounds_local(1, :)
705 ubounds_local = grid%pw_grid%bounds_local(2, :)
707 size_of_z = ceiling(real(ubounds_local(3) - lbounds_local(3) + 1,
dp)/real(stride(3),
dp))
709 DO i = lbounds(1), ubounds(1), stride(1)
710 should_write(:) = .true.
711 IF (i < lbounds_local(1))
THEN
712 should_write(1) = .false.
713 ELSE IF (i > ubounds_local(1))
THEN
716 DO j = lbounds(2), ubounds(2), stride(2)
717 should_write(2) = .true.
718 IF (j < lbounds_local(2) .OR. j > ubounds_local(2))
THEN
719 should_write(2) = .false.
721 IF (all(should_write .EQV. .true.))
THEN
727 DO k = lbounds(3), ubounds(3), stride(3)
728 IF (k + stride(3) > ubounds(3)) last_z = k
732 CALL unit_nr%get_position(bof)
734 IF (my_rank == 0)
THEN
737 CALL unit_nr%write_at(bof,
"-Quickstep-"//new_line(
"C"))
739 IF (
PRESENT(title))
THEN
740 CALL unit_nr%write_at(bof, trim(title)//new_line(
"C"))
743 CALL unit_nr%write_at(bof,
"No Title"//new_line(
"C"))
747 cpassert(
PRESENT(particles_z) .EQV.
PRESENT(particles_r))
749 IF (
PRESENT(particles_z))
THEN
750 cpassert(
SIZE(particles_z) ==
SIZE(particles_r, dim=2))
753 np = min(99999,
SIZE(particles_z))
756 WRITE (
header,
'(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp
757 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
760 WRITE (
header,
'(I5,3f12.6)') (grid%pw_grid%npts(1) + stride(1) - 1)/stride(1), &
761 grid%pw_grid%dh(1, 1)*real(stride(1),
dp), grid%pw_grid%dh(2, 1)*real(stride(1),
dp), &
762 grid%pw_grid%dh(3, 1)*real(stride(1),
dp)
763 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
766 WRITE (
header,
'(I5,3f12.6)') (grid%pw_grid%npts(2) + stride(2) - 1)/stride(2), &
767 grid%pw_grid%dh(1, 2)*real(stride(2),
dp), grid%pw_grid%dh(2, 2)*real(stride(2),
dp), &
768 grid%pw_grid%dh(3, 2)*real(stride(2),
dp)
769 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
772 WRITE (
header,
'(I5,3f12.6)') (grid%pw_grid%npts(3) + stride(3) - 1)/stride(3), &
773 grid%pw_grid%dh(1, 3)*real(stride(3),
dp), grid%pw_grid%dh(2, 3)*real(stride(3),
dp), &
774 grid%pw_grid%dh(3, 3)*real(stride(3),
dp)
775 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
778 IF (
PRESENT(particles_z))
THEN
779 IF (
PRESENT(particles_zeff))
THEN
781 WRITE (header_z,
'(I5,4f12.6)') particles_z(i), particles_zeff(i), particles_r(:, i)
782 CALL unit_nr%write_at(bof, header_z//new_line(
"C"))
787 WRITE (header_z,
'(I5,4f12.6)') particles_z(i), 0._dp, particles_r(:, i)
788 CALL unit_nr%write_at(bof, header_z//new_line(
"C"))
795 CALL gid%bcast(bof, grid%pw_grid%para%group%source)
798 ALLOCATE (displacements(nslices))
802 DO i = lbounds(1), ubounds(1), stride(1)
803 should_write(:) = .true.
804 IF (i < lbounds_local(1))
THEN
805 should_write(1) = .false.
806 ELSE IF (i > ubounds_local(1))
THEN
809 DO j = lbounds(2), ubounds(2), stride(2)
810 should_write(2) = .true.
811 IF (j < lbounds_local(2) .OR. j > ubounds_local(2))
THEN
812 should_write(2) = .false.
814 IF (all(should_write .EQV. .true.))
THEN
815 IF (islice > nslices) cpabort(
"Index out of bounds.")
816 displacements(islice) = bof
819 DO k = lbounds(3), ubounds(3), stride(3)
820 IF (zero_tails .AND. grid%array(i, j, k) < 1.e-7_dp)
THEN
825 tmp = trim(tmp)//trim(
value)
826 counter = counter + 1
827 IF (
modulo(counter, num_entries_line) == 0 .OR. k == last_z) &
828 tmp = trim(tmp)//new_line(
'C')
839 ALLOCATE (blocklengths(nslices))
840 blocklengths(:) = msglen
849 CALL unit_nr%write_all(msglen, nslices,
writebuffer, mp_desc)
853 DEALLOCATE (blocklengths, displacements)
855 END SUBROUTINE pw_to_cube_parallel
869 INTEGER,
INTENT(IN) :: unit_nr
870 INTEGER,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: stride
873 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_to_simple_volumetric'
875 INTEGER :: checksum, dest, handle, i, i1, i2, i3, &
876 ip, l1, l2, l3, my_rank, my_stride(3), &
877 ngrids, npoints, num_pe, rank(2), &
878 source, tag, u1, u2, u3
880 REAL(kind=
dp) :: x, y, z
881 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buf, buf2
884 CALL timeset(routinen, handle)
888 IF (
PRESENT(pw2)) double = .true.
891 IF (
PRESENT(stride))
THEN
892 IF (
SIZE(stride) /= 1 .AND.
SIZE(stride) /= 3) &
893 CALL cp_abort(__location__,
"STRIDE keyword can accept only 1 "// &
894 "(the same for X,Y,Z) or 3 values. Correct your input file.")
895 IF (
SIZE(stride) == 1)
THEN
897 my_stride(i) = stride(1)
900 my_stride = stride(1:3)
902 cpassert(my_stride(1) > 0)
903 cpassert(my_stride(2) > 0)
904 cpassert(my_stride(3) > 0)
908 l1 = pw%pw_grid%bounds(1, 1)
909 l2 = pw%pw_grid%bounds(1, 2)
910 l3 = pw%pw_grid%bounds(1, 3)
911 u1 = pw%pw_grid%bounds(2, 1)
912 u2 = pw%pw_grid%bounds(2, 2)
913 u3 = pw%pw_grid%bounds(2, 3)
917 IF (double) ngrids = 2
918 npoints = ((pw%pw_grid%npts(1) + my_stride(1) - 1)/my_stride(1))* &
919 ((pw%pw_grid%npts(2) + my_stride(2) - 1)/my_stride(1))* &
920 ((pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(1))
921 IF (unit_nr > 1)
WRITE (unit_nr,
'(I7,I5)') npoints, ngrids
923 ALLOCATE (buf(l3:u3))
924 IF (double)
ALLOCATE (buf2(l3:u3))
926 my_rank = pw%pw_grid%para%group%mepos
927 gid = pw%pw_grid%para%group
928 num_pe = pw%pw_grid%para%group%num_pe
934 IF (unit_nr > 0) checksum = 1
936 CALL gid%sum(checksum)
937 cpassert(checksum == 1)
939 CALL gid%maxloc(rank)
940 cpassert(rank(1) > 0)
943 DO i1 = l1, u1, my_stride(1)
944 DO i2 = l2, u2, my_stride(2)
948 DO ip = 0, num_pe - 1
949 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. &
950 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
958 IF (source == dest)
THEN
959 IF (my_rank == source)
THEN
960 buf(:) = pw%array(i1, i2, :)
961 IF (double) buf2(:) = pw2%array(i1, i2, :)
964 IF (my_rank == source)
THEN
965 buf(:) = pw%array(i1, i2, :)
966 CALL gid%send(buf, dest, tag)
968 buf2(:) = pw2%array(i1, i2, :)
969 CALL gid%send(buf2, dest, tag)
972 IF (my_rank == dest)
THEN
973 CALL gid%recv(buf, source, tag)
974 IF (double)
CALL gid%recv(buf2, source, tag)
978 IF (.NOT. double)
THEN
979 DO i3 = l3, u3, my_stride(3)
980 x = pw%pw_grid%dh(1, 1)*i1 + &
981 pw%pw_grid%dh(2, 1)*i2 + &
982 pw%pw_grid%dh(3, 1)*i3
984 y = pw%pw_grid%dh(1, 2)*i1 + &
985 pw%pw_grid%dh(2, 2)*i2 + &
986 pw%pw_grid%dh(3, 2)*i3
988 z = pw%pw_grid%dh(1, 3)*i1 + &
989 pw%pw_grid%dh(2, 3)*i2 + &
990 pw%pw_grid%dh(3, 3)*i3
992 IF (unit_nr > 0)
THEN
993 WRITE (unit_nr,
'(6E13.5E3, 6E13.5E3, 6E13.5E3, 6E13.5E3)') x, y, z, buf(i3)
999 DO i3 = l3, u3, my_stride(3)
1000 x = pw%pw_grid%dh(1, 1)*i1 + &
1001 pw%pw_grid%dh(2, 1)*i2 + &
1002 pw%pw_grid%dh(3, 1)*i3
1004 y = pw%pw_grid%dh(1, 2)*i1 + &
1005 pw%pw_grid%dh(2, 2)*i2 + &
1006 pw%pw_grid%dh(3, 2)*i3
1008 z = pw%pw_grid%dh(1, 3)*i1 + &
1009 pw%pw_grid%dh(2, 3)*i2 + &
1010 pw%pw_grid%dh(3, 3)*i3
1012 IF (unit_nr > 0)
THEN
1013 WRITE (unit_nr,
'(6E13.5E3, 6E13.5E3, 6E13.5E3, 6E13.5E3)') x, y, z, buf(i3), buf2(i3)
1031 IF (double)
DEALLOCATE (buf2)
1033 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.
character(len= *), parameter, public cube_values_format
subroutine, public pw_to_simple_volumetric(pw, unit_nr, stride, pw2)
Prints a simple grid file: X Y Z value.
character(len= *), parameter, public cube_value_format
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)
...
subroutine, public cube_read_values(values, buffer)
Read cube values from a character buffer.
void writebuffer(int *psockfd, char *data, int *plen)
Writes to a socket.
void readbuffer(int *psockfd, char *data, int *plen)
Reads from a socket.