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., &
49 SUBROUTINE pw_to_cube(pw, unit_nr, title, particles_r, particles_z, stride, zero_tails, &
52 INTEGER,
INTENT(IN) :: unit_nr
53 CHARACTER(*),
INTENT(IN),
OPTIONAL :: title
54 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
55 OPTIONAL :: particles_r
56 INTEGER,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_z
57 INTEGER,
DIMENSION(:),
OPTIONAL,
POINTER :: stride
58 LOGICAL,
INTENT(IN),
OPTIONAL :: zero_tails, silent, mpi_io
60 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_to_cube'
61 INTEGER,
PARAMETER :: entry_len = 13, num_entries_line = 6
63 INTEGER :: checksum, dest, handle, i, i1, i2, i3, iat, ip, l1, l2, l3, msglen, my_rank, &
64 my_stride(3), np, num_linebreak, num_pe, rank(2), size_of_z, source, tag, u1, u2, u3
65 LOGICAL :: be_silent, my_zero_tails, parallel_write
66 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buf
70 CALL timeset(routinen, handle)
72 my_zero_tails = .false.
74 parallel_write = .false.
75 IF (
PRESENT(zero_tails)) my_zero_tails = zero_tails
76 IF (
PRESENT(silent)) be_silent = silent
77 IF (
PRESENT(mpi_io)) parallel_write = mpi_io
79 IF (
PRESENT(stride))
THEN
80 IF (
SIZE(stride) /= 1 .AND.
SIZE(stride) /= 3) &
81 CALL cp_abort(__location__,
"STRIDE keyword can accept only 1 "// &
82 "(the same for X,Y,Z) or 3 values. Correct your input file.")
83 IF (
SIZE(stride) == 1)
THEN
85 my_stride(i) = stride(1)
88 my_stride = stride(1:3)
90 cpassert(my_stride(1) > 0)
91 cpassert(my_stride(2) > 0)
92 cpassert(my_stride(3) > 0)
95 IF (.NOT. parallel_write)
THEN
99 WRITE (unit_nr,
'(a11)')
"-Quickstep-"
100 IF (
PRESENT(title))
THEN
101 WRITE (unit_nr, *) trim(title)
103 WRITE (unit_nr, *)
"No Title"
106 cpassert(
PRESENT(particles_z) .EQV.
PRESENT(particles_r))
108 IF (
PRESENT(particles_z))
THEN
109 cpassert(
SIZE(particles_z) ==
SIZE(particles_r, dim=2))
112 np = min(99999,
SIZE(particles_z))
115 WRITE (unit_nr,
'(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp
117 WRITE (unit_nr,
'(I5,3f12.6)') (pw%pw_grid%npts(1) + my_stride(1) - 1)/my_stride(1), &
118 pw%pw_grid%dh(1, 1)*real(my_stride(1),
dp), pw%pw_grid%dh(2, 1)*real(my_stride(1),
dp), &
119 pw%pw_grid%dh(3, 1)*real(my_stride(1),
dp)
120 WRITE (unit_nr,
'(I5,3f12.6)') (pw%pw_grid%npts(2) + my_stride(2) - 1)/my_stride(2), &
121 pw%pw_grid%dh(1, 2)*real(my_stride(2),
dp), pw%pw_grid%dh(2, 2)*real(my_stride(2),
dp), &
122 pw%pw_grid%dh(3, 2)*real(my_stride(2),
dp)
123 WRITE (unit_nr,
'(I5,3f12.6)') (pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(3), &
124 pw%pw_grid%dh(1, 3)*real(my_stride(3),
dp), pw%pw_grid%dh(2, 3)*real(my_stride(3),
dp), &
125 pw%pw_grid%dh(3, 3)*real(my_stride(3),
dp)
127 IF (
PRESENT(particles_z))
THEN
129 WRITE (unit_nr,
'(I5,4f12.6)') particles_z(iat), 0._dp, particles_r(:, iat)
135 l1 = pw%pw_grid%bounds(1, 1)
136 l2 = pw%pw_grid%bounds(1, 2)
137 l3 = pw%pw_grid%bounds(1, 3)
138 u1 = pw%pw_grid%bounds(2, 1)
139 u2 = pw%pw_grid%bounds(2, 2)
140 u3 = pw%pw_grid%bounds(2, 3)
142 ALLOCATE (buf(l3:u3))
144 my_rank = pw%pw_grid%para%group%mepos
145 gid = pw%pw_grid%para%group
146 num_pe = pw%pw_grid%para%group%num_pe
152 IF (unit_nr > 0) checksum = 1
154 CALL gid%sum(checksum)
155 cpassert(checksum == 1)
157 CALL gid%maxloc(rank)
158 cpassert(rank(1) > 0)
161 DO i1 = l1, u1, my_stride(1)
162 DO i2 = l2, u2, my_stride(2)
166 DO ip = 0, num_pe - 1
167 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. &
168 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
176 IF (source == dest)
THEN
177 IF (my_rank == source)
THEN
178 buf(:) = pw%array(i1, i2, :)
181 IF (my_rank == source)
THEN
182 buf(:) = pw%array(i1, i2, :)
183 CALL gid%send(buf, dest, tag)
185 IF (my_rank == dest)
THEN
186 CALL gid%recv(buf, source, tag)
190 IF (unit_nr > 0)
THEN
191 IF (my_zero_tails)
THEN
193 IF (buf(i3) < 1.e-7_dp) buf(i3) = 0.0_dp
196 WRITE (unit_nr,
'(6E13.5)') (buf(i3), i3=l3, u3, my_stride(3))
212 size_of_z = ceiling(real(pw%pw_grid%bounds(2, 3) - pw%pw_grid%bounds(1, 3) + 1,
dp)/real(my_stride(3),
dp))
213 num_linebreak = size_of_z/num_entries_line
214 IF (
modulo(size_of_z, num_entries_line) /= 0) &
215 num_linebreak = num_linebreak + 1
217 CALL mp_unit%set_handle(unit_nr)
218 CALL pw_to_cube_parallel(pw, mp_unit, title, particles_r, particles_z, my_stride, my_zero_tails, msglen)
221 CALL timestop(handle)
237 SUBROUTINE cube_to_pw(grid, filename, scaling, parallel_read, silent)
240 CHARACTER(len=*),
INTENT(in) :: filename
241 REAL(kind=
dp),
INTENT(in) :: scaling
242 LOGICAL,
INTENT(in) :: parallel_read
243 LOGICAL,
INTENT(in),
OPTIONAL :: silent
245 CHARACTER(len=*),
PARAMETER :: routinen =
'cube_to_pw'
246 INTEGER,
PARAMETER :: entry_len = 13, num_entries_line = 6
248 INTEGER :: extunit, handle, i, j, k, msglen, &
249 my_rank, nat, ndum, num_linebreak, &
250 num_pe, output_unit, size_of_z, tag
251 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, npoints, &
252 npoints_local, ubounds, ubounds_local
254 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buffer
255 REAL(kind=
dp),
DIMENSION(3) :: dr, rdum
260 CALL timeset(routinen, handle)
263 IF (
PRESENT(silent))
THEN
267 gid = grid%pw_grid%para%group
268 my_rank = grid%pw_grid%para%group%mepos
269 num_pe = grid%pw_grid%para%group%num_pe
272 lbounds_local = grid%pw_grid%bounds_local(1, :)
273 ubounds_local = grid%pw_grid%bounds_local(2, :)
274 size_of_z = ubounds_local(3) - lbounds_local(3) + 1
276 IF (.NOT. parallel_read)
THEN
277 npoints = grid%pw_grid%npts
278 lbounds = grid%pw_grid%bounds(1, :)
279 ubounds = grid%pw_grid%bounds(2, :)
282 dr(i) = grid%pw_grid%dh(i, i)
285 npoints_local = grid%pw_grid%npts_local
287 ALLOCATE (buffer(lbounds(3):ubounds(3)))
289 IF (my_rank == 0)
THEN
290 IF (output_unit > 0 .AND. .NOT. be_silent)
THEN
291 WRITE (output_unit, fmt=
"(/,T2,A,/,/,T2,A,/)")
"Reading the cube file: ", trim(filename)
296 file_form=
"FORMATTED", &
297 file_action=
"READ", &
304 READ (extunit, *) nat, rdum
306 READ (extunit, *) ndum, rdum
307 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
308 output_unit > 0)
THEN
309 WRITE (output_unit, *)
"Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
310 WRITE (output_unit, *)
"Restart from density | ", ndum,
" DIFFERS FROM ", npoints(i)
311 WRITE (output_unit, *)
"Restart from density | ", rdum,
" DIFFERS FROM ", dr(i)
321 DO i = lbounds(1), ubounds(1)
322 DO j = lbounds(2), ubounds(2)
323 IF (my_rank .EQ. 0)
THEN
324 READ (extunit, *) (buffer(k), k=lbounds(3), ubounds(3))
326 CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0)
329 IF ((lbounds_local(1) .LE. i) .AND. (i .LE. ubounds_local(1)) .AND. (lbounds_local(2) .LE. j) &
330 .AND. (j .LE. ubounds_local(2)))
THEN
332 grid%array(i, j, lbounds(3):ubounds(3)) = buffer(lbounds(3):ubounds(3))*scaling
338 IF (my_rank == 0)
CALL close_file(unit_number=extunit)
347 num_linebreak = size_of_z/num_entries_line
348 IF (
modulo(size_of_z, num_entries_line) /= 0) &
349 num_linebreak = num_linebreak + 1
351 CALL cube_to_pw_parallel(grid, filename, scaling, msglen, silent=silent)
354 CALL timestop(handle)
369 SUBROUTINE cube_to_pw_parallel(grid, filename, scaling, msglen, silent)
372 CHARACTER(len=*),
INTENT(in) :: filename
373 REAL(kind=
dp),
INTENT(in) :: scaling
374 INTEGER,
INTENT(in) :: msglen
375 LOGICAL,
INTENT(in),
OPTIONAL :: silent
377 INTEGER,
PARAMETER :: entry_len = 13, num_entries_line = 6
379 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, npoints, &
380 npoints_local, ubounds, ubounds_local
381 INTEGER,
ALLOCATABLE,
DIMENSION(:),
TARGET :: blocklengths
382 INTEGER(kind=file_offset),
ALLOCATABLE, &
383 DIMENSION(:),
TARGET :: displacements
384 INTEGER(kind=file_offset) :: bof
385 INTEGER :: extunit_handle, i, ientry, islice, j, k, my_rank, nat, ndum, nslices, num_pe, &
386 offset_global, output_unit, pos, readstat, size_of_z, tag
387 CHARACTER(LEN=msglen),
ALLOCATABLE,
DIMENSION(:) ::
readbuffer
388 CHARACTER(LEN=msglen) :: tmp
389 LOGICAL :: be_silent, should_read(2)
390 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buffer
391 REAL(kind=
dp),
DIMENSION(3) :: dr, rdum
399 IF (
PRESENT(silent))
THEN
404 gid = grid%pw_grid%para%group
405 my_rank = grid%pw_grid%para%group%mepos
406 num_pe = grid%pw_grid%para%group%num_pe
410 dr(i) = grid%pw_grid%dh(i, i)
413 npoints = grid%pw_grid%npts
414 lbounds = grid%pw_grid%bounds(1, :)
415 ubounds = grid%pw_grid%bounds(2, :)
417 npoints_local = grid%pw_grid%npts_local
418 lbounds_local = grid%pw_grid%bounds_local(1, :)
419 ubounds_local = grid%pw_grid%bounds_local(2, :)
420 size_of_z = ubounds_local(3) - lbounds_local(3) + 1
421 nslices = (ubounds_local(1) - lbounds_local(1) + 1)*(ubounds_local(2) - lbounds_local(2) + 1)
425 IF (my_rank == 0)
THEN
426 IF (output_unit > 0 .AND. .NOT. be_silent)
THEN
427 WRITE (output_unit, fmt=
"(/,T2,A,/,/,T2,A,/)")
"Reading the cube file: ", trim(filename)
432 file_form=
"FORMATTED", &
433 file_action=
"READ", &
434 file_access=
"STREAM", &
435 unit_number=extunit_handle)
439 READ (extunit_handle, *)
441 READ (extunit_handle, *) nat, rdum
443 READ (extunit_handle, *) ndum, rdum
444 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
445 output_unit > 0)
THEN
446 WRITE (output_unit, *)
"Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
447 WRITE (output_unit, *)
"Restart from density | ", ndum,
" DIFFERS FROM ", npoints(i)
448 WRITE (output_unit, *)
"Restart from density | ", rdum,
" DIFFERS FROM ", dr(i)
453 READ (extunit_handle, *)
456 INQUIRE (extunit_handle, pos=offset_global)
460 CALL gid%bcast(offset_global, grid%pw_grid%para%group%source)
462 CALL extunit%open(groupid=gid, filepath=filename, amode_status=
file_amode_rdonly)
464 ALLOCATE (displacements(nslices))
466 DO i = lbounds(1), ubounds(1)
467 should_read(:) = .true.
468 IF (i .LT. lbounds_local(1))
THEN
469 should_read(1) = .false.
470 ELSE IF (i .GT. ubounds_local(1))
THEN
473 DO j = lbounds(2), ubounds(2)
474 should_read(2) = .true.
475 IF (j .LT. lbounds_local(2) .OR. j .GT. ubounds_local(2))
THEN
476 should_read(2) = .false.
478 IF (all(should_read .EQV. .true.))
THEN
479 IF (islice > nslices) cpabort(
"Index out of bounds.")
480 displacements(islice) = bof
488 ALLOCATE (blocklengths(nslices))
489 blocklengths(:) = msglen
497 CALL extunit%read_all(msglen, nslices,
readbuffer, mp_file_desc)
503 ALLOCATE (buffer(lbounds(3):ubounds(3)))
507 READ (
readbuffer(1), *, iostat=readstat) (buffer(k), k=lbounds(3), ubounds(3))
508 IF (readstat == 0)
THEN
509 parses_linebreaks = .true.
511 parses_linebreaks = .false.
516 DO islice = 1, nslices
517 IF (parses_linebreaks)
THEN
520 READ (
readbuffer(islice), *) (buffer(k), k=lbounds(3), ubounds(3))
526 DO k = lbounds_local(3), ubounds_local(3)
527 IF (
modulo(ientry, num_entries_line) == 0 .OR. k == ubounds_local(3))
THEN
529 READ (tmp(pos:pos + (entry_len - 2)),
'(E12.5)') buffer(k)
530 pos = pos + (entry_len + 1)
532 READ (tmp(pos:pos + (entry_len - 1)),
'(E13.5)') buffer(k)
533 pos = pos + entry_len
539 grid%array(i, j, lbounds(3):ubounds(3)) = scaling*buffer(lbounds(3):ubounds(3))
541 IF (j > ubounds_local(2))
THEN
547 DEALLOCATE (blocklengths, displacements)
548 IF (debug_this_module)
THEN
551 IF (my_rank == 0)
THEN
552 IF (output_unit > 0 .AND. .NOT. be_silent)
THEN
553 WRITE (output_unit, fmt=
"(/,T2,A,/,/,T2,A)")
"Reading the cube file: ", filename
558 file_form=
"FORMATTED", &
559 file_action=
"READ", &
560 unit_number=extunit_handle)
564 READ (extunit_handle, *)
566 READ (extunit_handle, *) nat, rdum
568 READ (extunit_handle, *) ndum, rdum
569 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
570 output_unit > 0)
THEN
571 WRITE (output_unit, *)
"Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
572 WRITE (output_unit, *)
"Restart from density | ", ndum,
" DIFFERS FROM ", npoints(i)
573 WRITE (output_unit, *)
"Restart from density | ", rdum,
" DIFFERS FROM ", dr(i)
578 READ (extunit_handle, *)
583 DO i = lbounds(1), ubounds(1)
584 DO j = lbounds(2), ubounds(2)
585 IF (my_rank .EQ. 0)
THEN
586 READ (extunit_handle, *) (buffer(k), k=lbounds(3), ubounds(3))
588 CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0)
591 IF ((lbounds_local(1) .LE. i) .AND. (i .LE. ubounds_local(1)) .AND. (lbounds_local(2) .LE. j) &
592 .AND. (j .LE. ubounds_local(2)))
THEN
594 IF (any(grid%array(i, j, lbounds(3):ubounds(3)) /= buffer(lbounds(3):ubounds(3))*scaling)) &
595 CALL cp_abort(__location__, &
596 "Error in parallel read of input cube file.")
602 IF (my_rank == 0)
CALL close_file(unit_number=extunit_handle)
608 END SUBROUTINE cube_to_pw_parallel
623 SUBROUTINE pw_to_cube_parallel(grid, unit_nr, title, particles_r, particles_z, stride, zero_tails, &
628 CHARACTER(*),
INTENT(IN),
OPTIONAL :: title
629 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN), &
630 OPTIONAL :: particles_r
631 INTEGER,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: particles_z
632 INTEGER,
INTENT(IN) :: stride(3)
633 LOGICAL,
INTENT(IN) :: zero_tails
634 INTEGER,
INTENT(IN) :: msglen
636 INTEGER,
PARAMETER :: entry_len = 13, header_len = 41, &
637 header_len_z = 53, num_entries_line = 6
639 CHARACTER(LEN=entry_len) :: value
640 CHARACTER(LEN=header_len) ::
header
641 CHARACTER(LEN=header_len_z) :: header_z
642 INTEGER,
DIMENSION(3) :: lbounds, lbounds_local, ubounds, &
644 INTEGER,
ALLOCATABLE,
DIMENSION(:),
TARGET :: blocklengths
645 INTEGER(kind=file_offset),
ALLOCATABLE, &
646 DIMENSION(:),
TARGET :: displacements
647 INTEGER(kind=file_offset) :: bof
648 INTEGER :: counter, i, islice, j, k, last_z, &
649 my_rank, np, nslices, size_of_z
650 CHARACTER(LEN=msglen),
ALLOCATABLE,
DIMENSION(:) ::
writebuffer
651 CHARACTER(LEN=msglen) :: tmp
652 LOGICAL :: should_write(2)
657 gid = grid%pw_grid%para%group
658 my_rank = grid%pw_grid%para%group%mepos
661 lbounds = grid%pw_grid%bounds(1, :)
662 ubounds = grid%pw_grid%bounds(2, :)
663 lbounds_local = grid%pw_grid%bounds_local(1, :)
664 ubounds_local = grid%pw_grid%bounds_local(2, :)
666 size_of_z = ceiling(real(ubounds_local(3) - lbounds_local(3) + 1,
dp)/real(stride(3),
dp))
668 DO i = lbounds(1), ubounds(1), stride(1)
669 should_write(:) = .true.
670 IF (i .LT. lbounds_local(1))
THEN
671 should_write(1) = .false.
672 ELSE IF (i .GT. ubounds_local(1))
THEN
675 DO j = lbounds(2), ubounds(2), stride(2)
676 should_write(2) = .true.
677 IF (j .LT. lbounds_local(2) .OR. j .GT. ubounds_local(2))
THEN
678 should_write(2) = .false.
680 IF (all(should_write .EQV. .true.))
THEN
686 DO k = lbounds(3), ubounds(3), stride(3)
687 IF (k + stride(3) > ubounds(3)) last_z = k
691 CALL unit_nr%get_position(bof)
693 IF (my_rank == 0)
THEN
696 CALL unit_nr%write_at(bof,
"-Quickstep-"//new_line(
"C"))
698 IF (
PRESENT(title))
THEN
699 CALL unit_nr%write_at(bof, trim(title)//new_line(
"C"))
702 CALL unit_nr%write_at(bof,
"No Title"//new_line(
"C"))
706 cpassert(
PRESENT(particles_z) .EQV.
PRESENT(particles_r))
708 IF (
PRESENT(particles_z))
THEN
709 cpassert(
SIZE(particles_z) ==
SIZE(particles_r, dim=2))
712 np = min(99999,
SIZE(particles_z))
715 WRITE (
header,
'(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp
716 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
719 WRITE (
header,
'(I5,3f12.6)') (grid%pw_grid%npts(1) + stride(1) - 1)/stride(1), &
720 grid%pw_grid%dh(1, 1)*real(stride(1),
dp), grid%pw_grid%dh(2, 1)*real(stride(1),
dp), &
721 grid%pw_grid%dh(3, 1)*real(stride(1),
dp)
722 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
725 WRITE (
header,
'(I5,3f12.6)') (grid%pw_grid%npts(2) + stride(2) - 1)/stride(2), &
726 grid%pw_grid%dh(1, 2)*real(stride(2),
dp), grid%pw_grid%dh(2, 2)*real(stride(2),
dp), &
727 grid%pw_grid%dh(3, 2)*real(stride(2),
dp)
728 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
731 WRITE (
header,
'(I5,3f12.6)') (grid%pw_grid%npts(3) + stride(3) - 1)/stride(3), &
732 grid%pw_grid%dh(1, 3)*real(stride(3),
dp), grid%pw_grid%dh(2, 3)*real(stride(3),
dp), &
733 grid%pw_grid%dh(3, 3)*real(stride(3),
dp)
734 CALL unit_nr%write_at(bof,
header//new_line(
"C"))
737 IF (
PRESENT(particles_z))
THEN
739 WRITE (header_z,
'(I5,4f12.6)') particles_z(i), 0._dp, particles_r(:, i)
740 CALL unit_nr%write_at(bof, header_z//new_line(
"C"))
746 CALL gid%bcast(bof, grid%pw_grid%para%group%source)
749 ALLOCATE (displacements(nslices))
753 DO i = lbounds(1), ubounds(1), stride(1)
754 should_write(:) = .true.
755 IF (i .LT. lbounds_local(1))
THEN
756 should_write(1) = .false.
757 ELSE IF (i .GT. ubounds_local(1))
THEN
760 DO j = lbounds(2), ubounds(2), stride(2)
761 should_write(2) = .true.
762 IF (j .LT. lbounds_local(2) .OR. j .GT. ubounds_local(2))
THEN
763 should_write(2) = .false.
765 IF (all(should_write .EQV. .true.))
THEN
766 IF (islice > nslices) cpabort(
"Index out of bounds.")
767 displacements(islice) = bof
770 DO k = lbounds(3), ubounds(3), stride(3)
771 IF (zero_tails .AND. grid%array(i, j, k) < 1.e-7_dp)
THEN
772 WRITE (
value,
'(E13.5)') 0.0_dp
774 WRITE (
value,
'(E13.5)') grid%array(i, j, k)
776 tmp = trim(tmp)//trim(
value)
777 counter = counter + 1
778 IF (
modulo(counter, num_entries_line) == 0 .OR. k == last_z) &
779 tmp = trim(tmp)//new_line(
'C')
790 ALLOCATE (blocklengths(nslices))
791 blocklengths(:) = msglen
800 CALL unit_nr%write_all(msglen, nslices,
writebuffer, mp_desc)
804 DEALLOCATE (blocklengths, displacements)
806 END SUBROUTINE pw_to_cube_parallel
820 INTEGER,
INTENT(IN) :: unit_nr
821 INTEGER,
DIMENSION(:),
INTENT(IN),
OPTIONAL :: stride
824 CHARACTER(len=*),
PARAMETER :: routinen =
'pw_to_simple_volumetric'
826 INTEGER :: checksum, dest, handle, i, i1, i2, i3, &
827 ip, l1, l2, l3, my_rank, my_stride(3), &
828 ngrids, npoints, num_pe, rank(2), &
829 source, tag, u1, u2, u3
831 REAL(kind=
dp) :: x, y, z
832 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: buf, buf2
835 CALL timeset(routinen, handle)
839 IF (
PRESENT(pw2)) double = .true.
842 IF (
PRESENT(stride))
THEN
843 IF (
SIZE(stride) /= 1 .AND.
SIZE(stride) /= 3) &
844 CALL cp_abort(__location__,
"STRIDE keyword can accept only 1 "// &
845 "(the same for X,Y,Z) or 3 values. Correct your input file.")
846 IF (
SIZE(stride) == 1)
THEN
848 my_stride(i) = stride(1)
851 my_stride = stride(1:3)
853 cpassert(my_stride(1) > 0)
854 cpassert(my_stride(2) > 0)
855 cpassert(my_stride(3) > 0)
859 l1 = pw%pw_grid%bounds(1, 1)
860 l2 = pw%pw_grid%bounds(1, 2)
861 l3 = pw%pw_grid%bounds(1, 3)
862 u1 = pw%pw_grid%bounds(2, 1)
863 u2 = pw%pw_grid%bounds(2, 2)
864 u3 = pw%pw_grid%bounds(2, 3)
868 IF (double) ngrids = 2
869 npoints = ((pw%pw_grid%npts(1) + my_stride(1) - 1)/my_stride(1))* &
870 ((pw%pw_grid%npts(2) + my_stride(2) - 1)/my_stride(1))* &
871 ((pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(1))
872 IF (unit_nr > 1)
WRITE (unit_nr,
'(I7,I5)') npoints, ngrids
874 ALLOCATE (buf(l3:u3))
875 IF (double)
ALLOCATE (buf2(l3:u3))
877 my_rank = pw%pw_grid%para%group%mepos
878 gid = pw%pw_grid%para%group
879 num_pe = pw%pw_grid%para%group%num_pe
885 IF (unit_nr > 0) checksum = 1
887 CALL gid%sum(checksum)
888 cpassert(checksum == 1)
890 CALL gid%maxloc(rank)
891 cpassert(rank(1) > 0)
894 DO i1 = l1, u1, my_stride(1)
895 DO i2 = l2, u2, my_stride(2)
899 DO ip = 0, num_pe - 1
900 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. &
901 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
909 IF (source == dest)
THEN
910 IF (my_rank == source)
THEN
911 buf(:) = pw%array(i1, i2, :)
912 IF (double) buf2(:) = pw2%array(i1, i2, :)
915 IF (my_rank == source)
THEN
916 buf(:) = pw%array(i1, i2, :)
917 CALL gid%send(buf, dest, tag)
919 buf2(:) = pw2%array(i1, i2, :)
920 CALL gid%send(buf2, dest, tag)
923 IF (my_rank == dest)
THEN
924 CALL gid%recv(buf, source, tag)
925 IF (double)
CALL gid%recv(buf2, source, tag)
929 IF (.NOT. double)
THEN
930 DO i3 = l3, u3, my_stride(3)
931 x = pw%pw_grid%dh(1, 1)*i1 + &
932 pw%pw_grid%dh(2, 1)*i2 + &
933 pw%pw_grid%dh(3, 1)*i3
935 y = pw%pw_grid%dh(1, 2)*i1 + &
936 pw%pw_grid%dh(2, 2)*i2 + &
937 pw%pw_grid%dh(3, 2)*i3
939 z = pw%pw_grid%dh(1, 3)*i1 + &
940 pw%pw_grid%dh(2, 3)*i2 + &
941 pw%pw_grid%dh(3, 3)*i3
943 IF (unit_nr > 0)
THEN
944 WRITE (unit_nr,
'(6E13.5, 6E13.5, 6E13.5, 6E13.5)') x, y, z, buf(i3)
950 DO i3 = l3, u3, my_stride(3)
951 x = pw%pw_grid%dh(1, 1)*i1 + &
952 pw%pw_grid%dh(2, 1)*i2 + &
953 pw%pw_grid%dh(3, 1)*i3
955 y = pw%pw_grid%dh(1, 2)*i1 + &
956 pw%pw_grid%dh(2, 2)*i2 + &
957 pw%pw_grid%dh(3, 2)*i3
959 z = pw%pw_grid%dh(1, 3)*i1 + &
960 pw%pw_grid%dh(2, 3)*i2 + &
961 pw%pw_grid%dh(3, 3)*i3
963 IF (unit_nr > 0)
THEN
964 WRITE (unit_nr,
'(6E13.5, 6E13.5, 6E13.5, 6E13.5)') x, y, z, buf(i3), buf2(i3)
982 IF (double)
DEALLOCATE (buf2)
984 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, stride, 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.