Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Generate Gaussian cube files
10 : ! **************************************************************************************************
11 : MODULE realspace_grid_cube
12 : USE cp_files, ONLY: close_file,&
13 : open_file
14 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
15 : USE kinds, ONLY: dp
16 : USE message_passing, ONLY: &
17 : file_amode_rdonly, file_offset, mp_comm_type, mp_file_descriptor_type, mp_file_type, &
18 : mp_file_type_free, mp_file_type_hindexed_make_chv, mp_file_type_set_view_chv, &
19 : mpi_character_size
20 : USE pw_grid_types, ONLY: PW_MODE_LOCAL
21 : USE pw_types, ONLY: pw_r3d_rs_type
22 : #include "../base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 :
28 : PUBLIC :: pw_to_cube, cube_to_pw, pw_to_simple_volumetric
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'realspace_grid_cube'
31 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
32 : LOGICAL, PRIVATE :: parses_linebreaks = .FALSE., &
33 : parse_test = .TRUE.
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief ...
39 : !> \param pw ...
40 : !> \param unit_nr ...
41 : !> \param title ...
42 : !> \param particles_r ...
43 : !> \param particles_z ...
44 : !> \param stride ...
45 : !> \param zero_tails ...
46 : !> \param silent ...
47 : !> \param mpi_io ...
48 : ! **************************************************************************************************
49 1920 : SUBROUTINE pw_to_cube(pw, unit_nr, title, particles_r, particles_z, stride, zero_tails, &
50 : silent, mpi_io)
51 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
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
59 :
60 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_to_cube'
61 : INTEGER, PARAMETER :: entry_len = 13, num_entries_line = 6
62 :
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 1920 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buf
67 : TYPE(mp_comm_type) :: gid
68 : TYPE(mp_file_type) :: mp_unit
69 :
70 1920 : CALL timeset(routineN, handle)
71 :
72 1920 : my_zero_tails = .FALSE.
73 1920 : be_silent = .FALSE.
74 1920 : parallel_write = .FALSE.
75 1920 : IF (PRESENT(zero_tails)) my_zero_tails = zero_tails
76 1920 : IF (PRESENT(silent)) be_silent = silent
77 1920 : IF (PRESENT(mpi_io)) parallel_write = mpi_io
78 7680 : my_stride = 1
79 1920 : IF (PRESENT(stride)) THEN
80 1872 : IF (SIZE(stride) /= 1 .AND. SIZE(stride) /= 3) &
81 : CALL cp_abort(__LOCATION__, "STRIDE keyword can accept only 1 "// &
82 1872 : "(the same for X,Y,Z) or 3 values. Correct your input file.")
83 1872 : IF (SIZE(stride) == 1) THEN
84 88 : DO i = 1, 3
85 88 : my_stride(i) = stride(1)
86 : END DO
87 : ELSE
88 7400 : my_stride = stride(1:3)
89 : END IF
90 1872 : CPASSERT(my_stride(1) > 0)
91 1872 : CPASSERT(my_stride(2) > 0)
92 1872 : CPASSERT(my_stride(3) > 0)
93 : END IF
94 :
95 1920 : IF (.NOT. parallel_write) THEN
96 70 : IF (unit_nr > 0) THEN
97 : ! this format seems to work for e.g. molekel and gOpenmol
98 : ! latest version of VMD can read non orthorhombic cells
99 35 : WRITE (unit_nr, '(a11)') "-Quickstep-"
100 35 : IF (PRESENT(title)) THEN
101 35 : WRITE (unit_nr, *) TRIM(title)
102 : ELSE
103 0 : WRITE (unit_nr, *) "No Title"
104 : END IF
105 :
106 35 : CPASSERT(PRESENT(particles_z) .EQV. PRESENT(particles_r))
107 35 : np = 0
108 35 : IF (PRESENT(particles_z)) THEN
109 35 : CPASSERT(SIZE(particles_z) == SIZE(particles_r, dim=2))
110 : ! cube files can only be written for 99999 particles due to a format limitation (I5)
111 : ! so we limit the number of particles written.
112 35 : np = MIN(99999, SIZE(particles_z))
113 : END IF
114 :
115 35 : WRITE (unit_nr, '(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp !start of cube
116 :
117 35 : WRITE (unit_nr, '(I5,3f12.6)') (pw%pw_grid%npts(1) + my_stride(1) - 1)/my_stride(1), &
118 35 : pw%pw_grid%dh(1, 1)*REAL(my_stride(1), dp), pw%pw_grid%dh(2, 1)*REAL(my_stride(1), dp), &
119 70 : pw%pw_grid%dh(3, 1)*REAL(my_stride(1), dp)
120 35 : WRITE (unit_nr, '(I5,3f12.6)') (pw%pw_grid%npts(2) + my_stride(2) - 1)/my_stride(2), &
121 35 : pw%pw_grid%dh(1, 2)*REAL(my_stride(2), dp), pw%pw_grid%dh(2, 2)*REAL(my_stride(2), dp), &
122 70 : pw%pw_grid%dh(3, 2)*REAL(my_stride(2), dp)
123 35 : WRITE (unit_nr, '(I5,3f12.6)') (pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(3), &
124 35 : pw%pw_grid%dh(1, 3)*REAL(my_stride(3), dp), pw%pw_grid%dh(2, 3)*REAL(my_stride(3), dp), &
125 70 : pw%pw_grid%dh(3, 3)*REAL(my_stride(3), dp)
126 :
127 35 : IF (PRESENT(particles_z)) THEN
128 165 : DO iat = 1, np
129 165 : WRITE (unit_nr, '(I5,4f12.6)') particles_z(iat), 0._dp, particles_r(:, iat)
130 : END DO
131 : END IF
132 : END IF
133 :
134 : ! shortcut
135 70 : L1 = pw%pw_grid%bounds(1, 1)
136 70 : L2 = pw%pw_grid%bounds(1, 2)
137 70 : L3 = pw%pw_grid%bounds(1, 3)
138 70 : U1 = pw%pw_grid%bounds(2, 1)
139 70 : U2 = pw%pw_grid%bounds(2, 2)
140 70 : U3 = pw%pw_grid%bounds(2, 3)
141 :
142 210 : ALLOCATE (buf(L3:U3))
143 :
144 70 : my_rank = pw%pw_grid%para%group%mepos
145 70 : gid = pw%pw_grid%para%group
146 70 : num_pe = pw%pw_grid%para%group%num_pe
147 70 : tag = 1
148 :
149 70 : rank (1) = unit_nr
150 70 : rank (2) = my_rank
151 70 : checksum = 0
152 70 : IF (unit_nr > 0) checksum = 1
153 :
154 70 : CALL gid%sum(checksum)
155 70 : CPASSERT(checksum == 1)
156 :
157 70 : CALL gid%maxloc(rank)
158 70 : CPASSERT(rank(1) > 0)
159 :
160 70 : dest = rank(2)
161 2110 : DO I1 = L1, U1, my_stride(1)
162 66146 : DO I2 = L2, U2, my_stride(2)
163 :
164 : ! cycling through the CPUs, check if the current ray (I1,I2) is local to that CPU
165 64036 : IF (pw%pw_grid%para%mode .NE. PW_MODE_LOCAL) THEN
166 192108 : 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 192108 : 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
169 64036 : source = ip
170 : END IF
171 : END DO
172 : ELSE
173 0 : source = dest
174 : END IF
175 :
176 64036 : IF (source == dest) THEN
177 32216 : IF (my_rank == source) THEN
178 717159 : buf(:) = pw%array(I1, I2, :)
179 : END IF
180 : ELSE
181 31820 : IF (my_rank == source) THEN
182 703236 : buf(:) = pw%array(I1, I2, :)
183 15910 : CALL gid%send(buf, dest, tag)
184 : END IF
185 31820 : IF (my_rank == dest) THEN
186 15910 : CALL gid%recv(buf, source, tag)
187 : END IF
188 : END IF
189 :
190 64036 : IF (unit_nr > 0) THEN
191 32018 : IF (my_zero_tails) THEN
192 0 : DO I3 = L3, U3
193 0 : IF (buf(I3) < 1.E-7_dp) buf(I3) = 0.0_dp
194 : END DO
195 : END IF
196 32018 : WRITE (unit_nr, '(6E13.5)') (buf(I3), I3=L3, U3, my_stride(3))
197 : END IF
198 :
199 : ! this double loop generates so many messages that it can overload
200 : ! the message passing system, e.g. on XT3
201 : ! we therefore put a barrier here that limits the amount of message
202 : ! that flies around at any given time.
203 : ! if ever this routine becomes a bottleneck, we should go for a
204 : ! more complicated rewrite
205 66076 : CALL gid%sync()
206 :
207 : END DO
208 : END DO
209 :
210 70 : DEALLOCATE (buf)
211 : ELSE
212 1850 : 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 1850 : num_linebreak = size_of_z/num_entries_line
214 1850 : IF (MODULO(size_of_z, num_entries_line) /= 0) &
215 1544 : num_linebreak = num_linebreak + 1
216 1850 : msglen = (size_of_z*entry_len + num_linebreak)*mpi_character_size
217 1850 : CALL mp_unit%set_handle(unit_nr)
218 1850 : CALL pw_to_cube_parallel(pw, mp_unit, title, particles_r, particles_z, my_stride, my_zero_tails, msglen)
219 : END IF
220 :
221 1920 : CALL timestop(handle)
222 :
223 3840 : END SUBROUTINE pw_to_cube
224 :
225 : ! **************************************************************************************************
226 : !> \brief Computes the external density on the grid
227 : !> hacked from external_read_density
228 : !> \param grid pw to read from cube file
229 : !> \param filename name of cube file
230 : !> \param scaling scale values before storing
231 : !> \param parallel_read ...
232 : !> \param silent ...
233 : !> \par History
234 : !> Created [M.Watkins] (01.2014)
235 : !> Use blocking, collective MPI read for parallel simulations [Nico Holmberg] (05.2017)
236 : ! **************************************************************************************************
237 38 : SUBROUTINE cube_to_pw(grid, filename, scaling, parallel_read, silent)
238 :
239 : TYPE(pw_r3d_rs_type), INTENT(IN) :: grid
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
244 :
245 : CHARACTER(len=*), PARAMETER :: routineN = 'cube_to_pw'
246 : INTEGER, PARAMETER :: entry_len = 13, num_entries_line = 6
247 :
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
253 : LOGICAL :: be_silent
254 38 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buffer
255 : REAL(kind=dp), DIMENSION(3) :: dr, rdum
256 : TYPE(mp_comm_type) :: gid
257 :
258 76 : output_unit = cp_logger_get_default_io_unit()
259 :
260 38 : CALL timeset(routineN, handle)
261 :
262 38 : be_silent = .FALSE.
263 38 : IF (PRESENT(silent)) THEN
264 0 : be_silent = silent
265 : END IF
266 : !get rs grids and parallel environment
267 38 : gid = grid%pw_grid%para%group
268 38 : my_rank = grid%pw_grid%para%group%mepos
269 38 : num_pe = grid%pw_grid%para%group%num_pe
270 38 : tag = 1
271 :
272 152 : lbounds_local = grid%pw_grid%bounds_local(1, :)
273 152 : ubounds_local = grid%pw_grid%bounds_local(2, :)
274 38 : size_of_z = ubounds_local(3) - lbounds_local(3) + 1
275 :
276 38 : IF (.NOT. parallel_read) THEN
277 0 : npoints = grid%pw_grid%npts
278 0 : lbounds = grid%pw_grid%bounds(1, :)
279 0 : ubounds = grid%pw_grid%bounds(2, :)
280 :
281 0 : DO i = 1, 3
282 0 : dr(i) = grid%pw_grid%dh(i, i)
283 : END DO
284 :
285 0 : npoints_local = grid%pw_grid%npts_local
286 : !pw grids at most pencils - all processors have a full set of z data for x,y
287 0 : ALLOCATE (buffer(lbounds(3):ubounds(3)))
288 :
289 0 : IF (my_rank == 0) THEN
290 0 : IF (output_unit > 0 .AND. .NOT. be_silent) THEN
291 0 : WRITE (output_unit, FMT="(/,T2,A,/,/,T2,A,/)") "Reading the cube file: ", TRIM(filename)
292 : END IF
293 :
294 : CALL open_file(file_name=filename, &
295 : file_status="OLD", &
296 : file_form="FORMATTED", &
297 : file_action="READ", &
298 0 : unit_number=extunit)
299 :
300 : !skip header comments
301 0 : DO i = 1, 2
302 0 : READ (extunit, *)
303 : END DO
304 0 : READ (extunit, *) nat, rdum
305 0 : DO i = 1, 3
306 0 : READ (extunit, *) ndum, rdum
307 0 : IF ((ndum /= npoints(i) .OR. (ABS(rdum(i) - dr(i)) > 1e-4)) .AND. &
308 0 : output_unit > 0) THEN
309 0 : WRITE (output_unit, *) "Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
310 0 : WRITE (output_unit, *) "Restart from density | ", ndum, " DIFFERS FROM ", npoints(i)
311 0 : WRITE (output_unit, *) "Restart from density | ", rdum, " DIFFERS FROM ", dr(i)
312 : END IF
313 : END DO
314 : !ignore atomic position data - read from coord or topology instead
315 0 : DO i = 1, nat
316 0 : READ (extunit, *)
317 : END DO
318 : END IF
319 :
320 : !master sends all data to everyone
321 0 : DO i = lbounds(1), ubounds(1)
322 0 : DO j = lbounds(2), ubounds(2)
323 0 : IF (my_rank .EQ. 0) THEN
324 0 : READ (extunit, *) (buffer(k), k=lbounds(3), ubounds(3))
325 : END IF
326 0 : CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0)
327 :
328 : !only use data that is local to me - i.e. in slice of pencil I own
329 : IF ((lbounds_local(1) .LE. i) .AND. (i .LE. ubounds_local(1)) .AND. (lbounds_local(2) .LE. j) &
330 0 : .AND. (j .LE. ubounds_local(2))) THEN
331 : !allow scaling of external potential values by factor 'scaling' (SCALING_FACTOR in input file)
332 0 : grid%array(i, j, lbounds(3):ubounds(3)) = buffer(lbounds(3):ubounds(3))*scaling
333 : END IF
334 :
335 : END DO
336 : END DO
337 :
338 0 : IF (my_rank == 0) CALL close_file(unit_number=extunit)
339 :
340 0 : CALL gid%sync()
341 : ELSE
342 : ! Parallel routine needs as input the byte size of each grid z-slice
343 : ! This is a hack to prevent compilation errors with gcc -Wall (up to versions 6.3)
344 : ! related to allocatable-length string declaration CHARACTER(LEN=:), ALLOCATABLE, DIMENSION(:) :: string
345 : ! Each data line of a Gaussian cube contains max 6 entries with last line potentially containing less if nz % 6 /= 0
346 : ! Thus, this size is simply the number of entries multiplied by the entry size + the number of line breaks
347 38 : num_linebreak = size_of_z/num_entries_line
348 38 : IF (MODULO(size_of_z, num_entries_line) /= 0) &
349 36 : num_linebreak = num_linebreak + 1
350 38 : msglen = (size_of_z*entry_len + num_linebreak)*mpi_character_size
351 38 : CALL cube_to_pw_parallel(grid, filename, scaling, msglen, silent=silent)
352 : END IF
353 :
354 38 : CALL timestop(handle)
355 :
356 76 : END SUBROUTINE cube_to_pw
357 :
358 : ! **************************************************************************************************
359 : !> \brief Reads a realspace potential/density from a cube file using collective MPI I/O and
360 : !> stores it in grid.
361 : !> \param grid pw to read from cube file
362 : !> \param filename name of cube file
363 : !> \param scaling scale values before storing
364 : !> \param msglen the size of each grid slice along z-axis in bytes
365 : !> \param silent ...
366 : !> \par History
367 : !> Created [Nico Holmberg] (05.2017)
368 : ! **************************************************************************************************
369 38 : SUBROUTINE cube_to_pw_parallel(grid, filename, scaling, msglen, silent)
370 :
371 : TYPE(pw_r3d_rs_type), INTENT(IN) :: grid
372 : CHARACTER(len=*), INTENT(in) :: filename
373 : REAL(kind=dp), INTENT(in) :: scaling
374 : INTEGER, INTENT(in) :: msglen
375 : LOGICAL, INTENT(in), OPTIONAL :: silent
376 :
377 : INTEGER, PARAMETER :: entry_len = 13, num_entries_line = 6
378 :
379 : INTEGER, DIMENSION(3) :: lbounds, lbounds_local, npoints, &
380 : npoints_local, ubounds, ubounds_local
381 38 : INTEGER, ALLOCATABLE, DIMENSION(:), TARGET :: blocklengths
382 : INTEGER(kind=file_offset), ALLOCATABLE, &
383 38 : 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 38 : CHARACTER(LEN=msglen), ALLOCATABLE, DIMENSION(:) :: readbuffer
388 38 : CHARACTER(LEN=msglen) :: tmp
389 : LOGICAL :: be_silent, should_read(2)
390 38 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buffer
391 : REAL(kind=dp), DIMENSION(3) :: dr, rdum
392 : TYPE(mp_comm_type) :: gid
393 : TYPE(mp_file_descriptor_type) :: mp_file_desc
394 : TYPE(mp_file_type) :: extunit
395 :
396 76 : output_unit = cp_logger_get_default_io_unit()
397 :
398 38 : be_silent = .FALSE.
399 38 : IF (PRESENT(silent)) THEN
400 0 : be_silent = silent
401 : END IF
402 :
403 : !get rs grids and parallel envnment
404 38 : gid = grid%pw_grid%para%group
405 38 : my_rank = grid%pw_grid%para%group%mepos
406 38 : num_pe = grid%pw_grid%para%group%num_pe
407 38 : tag = 1
408 :
409 152 : DO i = 1, 3
410 152 : dr(i) = grid%pw_grid%dh(i, i)
411 : END DO
412 :
413 152 : npoints = grid%pw_grid%npts
414 152 : lbounds = grid%pw_grid%bounds(1, :)
415 152 : ubounds = grid%pw_grid%bounds(2, :)
416 :
417 152 : npoints_local = grid%pw_grid%npts_local
418 152 : lbounds_local = grid%pw_grid%bounds_local(1, :)
419 152 : ubounds_local = grid%pw_grid%bounds_local(2, :)
420 38 : size_of_z = ubounds_local(3) - lbounds_local(3) + 1
421 38 : nslices = (ubounds_local(1) - lbounds_local(1) + 1)*(ubounds_local(2) - lbounds_local(2) + 1)
422 38 : islice = 1
423 :
424 : ! Read header information and determine byte offset of cube data on master process
425 38 : IF (my_rank == 0) THEN
426 19 : IF (output_unit > 0 .AND. .NOT. be_silent) THEN
427 19 : WRITE (output_unit, FMT="(/,T2,A,/,/,T2,A,/)") "Reading the cube file: ", TRIM(filename)
428 : END IF
429 :
430 : CALL open_file(file_name=filename, &
431 : file_status="OLD", &
432 : file_form="FORMATTED", &
433 : file_action="READ", &
434 : file_access="STREAM", &
435 19 : unit_number=extunit_handle)
436 :
437 : !skip header comments
438 57 : DO i = 1, 2
439 57 : READ (extunit_handle, *)
440 : END DO
441 19 : READ (extunit_handle, *) nat, rdum
442 76 : DO i = 1, 3
443 57 : READ (extunit_handle, *) ndum, rdum
444 57 : IF ((ndum /= npoints(i) .OR. (ABS(rdum(i) - dr(i)) > 1e-4)) .AND. &
445 19 : output_unit > 0) THEN
446 0 : WRITE (output_unit, *) "Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
447 0 : WRITE (output_unit, *) "Restart from density | ", ndum, " DIFFERS FROM ", npoints(i)
448 0 : WRITE (output_unit, *) "Restart from density | ", rdum, " DIFFERS FROM ", dr(i)
449 : END IF
450 : END DO
451 : !ignore atomic position data - read from coord or topology instead
452 44 : DO i = 1, nat
453 44 : READ (extunit_handle, *)
454 : END DO
455 : ! Get byte offset
456 19 : INQUIRE (extunit_handle, POS=offset_global)
457 19 : CALL close_file(unit_number=extunit_handle)
458 : END IF
459 : ! Sync offset and start parallel read
460 38 : CALL gid%bcast(offset_global, grid%pw_grid%para%group%source)
461 38 : BOF = offset_global
462 38 : CALL extunit%open(groupid=gid, filepath=filename, amode_status=file_amode_rdonly)
463 : ! Determine byte offsets for each grid z-slice which are local to a process
464 114 : ALLOCATE (displacements(nslices))
465 29450 : displacements = 0
466 1151 : DO i = lbounds(1), ubounds(1)
467 3396 : should_read(:) = .TRUE.
468 1132 : IF (i .LT. lbounds_local(1)) THEN
469 371 : should_read(1) = .FALSE.
470 761 : ELSE IF (i .GT. ubounds_local(1)) THEN
471 : EXIT
472 : END IF
473 45269 : DO j = lbounds(2), ubounds(2)
474 44118 : should_read(2) = .TRUE.
475 44118 : IF (j .LT. lbounds_local(2) .OR. j .GT. ubounds_local(2)) THEN
476 0 : should_read(2) = .FALSE.
477 : END IF
478 102942 : IF (ALL(should_read .EQV. .TRUE.)) THEN
479 29412 : IF (islice > nslices) CPABORT("Index out of bounds.")
480 29412 : displacements(islice) = BOF
481 29412 : islice = islice + 1
482 : END IF
483 : ! Update global byte offset
484 45231 : BOF = BOF + msglen
485 : END DO
486 : END DO
487 : ! Size of each z-slice is msglen
488 114 : ALLOCATE (blocklengths(nslices))
489 29450 : blocklengths(:) = msglen
490 : ! Create indexed MPI type using calculated byte offsets as displacements and use it as a file view
491 38 : mp_file_desc = mp_file_type_hindexed_make_chv(nslices, blocklengths, displacements)
492 38 : BOF = 0
493 38 : CALL mp_file_type_set_view_chv(extunit, BOF, mp_file_desc)
494 : ! Collective read of cube
495 114 : ALLOCATE (readbuffer(nslices))
496 29450 : readbuffer(:) = ''
497 38 : CALL extunit%read_all(msglen, nslices, readbuffer, mp_file_desc)
498 38 : CALL mp_file_type_free(mp_file_desc)
499 38 : CALL extunit%close()
500 : ! Convert cube values string -> real
501 38 : i = lbounds_local(1)
502 38 : j = lbounds_local(2)
503 114 : ALLOCATE (buffer(lbounds(3):ubounds(3)))
504 1522 : buffer = 0.0_dp
505 : ! Test if the compiler supports parsing lines with line breaks in them
506 38 : IF (parse_test) THEN
507 16 : READ (readbuffer(1), *, IOSTAT=readstat) (buffer(k), k=lbounds(3), ubounds(3))
508 16 : IF (readstat == 0) THEN
509 16 : parses_linebreaks = .TRUE.
510 : ELSE
511 0 : parses_linebreaks = .FALSE.
512 : END IF
513 16 : parse_test = .FALSE.
514 652 : buffer = 0.0_dp
515 : END IF
516 29450 : DO islice = 1, nslices
517 29412 : IF (parses_linebreaks) THEN
518 : ! Use list-directed conversion if supported
519 : ! Produces faster code, but maybe the latter method could be optimized
520 29412 : READ (readbuffer(islice), *) (buffer(k), k=lbounds(3), ubounds(3))
521 : ELSE
522 : ! Convert values by going through the z-slice one value at a time
523 0 : tmp = readbuffer(islice)
524 : pos = 1
525 : ientry = 1
526 0 : DO k = lbounds_local(3), ubounds_local(3)
527 0 : IF (MODULO(ientry, num_entries_line) == 0 .OR. k == ubounds_local(3)) THEN
528 : ! Last value on line, dont read line break
529 0 : READ (tmp(pos:pos + (entry_len - 2)), '(E12.5)') buffer(k)
530 0 : pos = pos + (entry_len + 1)
531 : ELSE
532 0 : READ (tmp(pos:pos + (entry_len - 1)), '(E13.5)') buffer(k)
533 0 : pos = pos + entry_len
534 : END IF
535 0 : ientry = ientry + 1
536 : END DO
537 : END IF
538 : ! Optionally scale cube file values
539 1213948 : grid%array(i, j, lbounds(3):ubounds(3)) = scaling*buffer(lbounds(3):ubounds(3))
540 29412 : j = j + 1
541 29450 : IF (j > ubounds_local(2)) THEN
542 742 : j = lbounds_local(2)
543 742 : i = i + 1
544 : END IF
545 : END DO
546 38 : DEALLOCATE (readbuffer)
547 38 : DEALLOCATE (blocklengths, displacements)
548 : IF (debug_this_module) THEN
549 : ! Check that cube was correctly read using intrinsic read on master who sends data to everyone
550 : buffer = 0.0_dp
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
554 : END IF
555 :
556 : CALL open_file(file_name=filename, &
557 : file_status="OLD", &
558 : file_form="FORMATTED", &
559 : file_action="READ", &
560 : unit_number=extunit_handle)
561 :
562 : !skip header comments
563 : DO i = 1, 2
564 : READ (extunit_handle, *)
565 : END DO
566 : READ (extunit_handle, *) nat, rdum
567 : DO i = 1, 3
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)
574 : END IF
575 : END DO
576 : !ignore atomic position data - read from coord or topology instead
577 : DO i = 1, nat
578 : READ (extunit_handle, *)
579 : END DO
580 : END IF
581 :
582 : !master sends all data to everyone
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))
587 : END IF
588 : CALL gid%bcast(buffer(lbounds(3):ubounds(3)), 0)
589 :
590 : !only use data that is local to me - i.e. in slice of pencil I own
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
593 : !allow scaling of external potential values by factor 'scaling' (SCALING_FACTOR in input file)
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.")
597 : END IF
598 :
599 : END DO
600 : END DO
601 :
602 : IF (my_rank == 0) CALL close_file(unit_number=extunit_handle)
603 :
604 : CALL gid%sync()
605 : END IF
606 38 : DEALLOCATE (buffer)
607 :
608 38 : END SUBROUTINE cube_to_pw_parallel
609 :
610 : ! **************************************************************************************************
611 : !> \brief Writes a realspace potential to a cube file using collective MPI I/O.
612 : !> \param grid the pw to output to the cube file
613 : !> \param unit_nr the handle associated with the cube file
614 : !> \param title title of the cube file
615 : !> \param particles_r Cartersian coordinates of the system
616 : !> \param particles_z atomic masses of atoms in the system
617 : !> \param stride every stride(i)th value of the potential is outputted (i=x,y,z)
618 : !> \param zero_tails flag that determines if small values of the potential should be zeroed
619 : !> \param msglen the size of each grid slice along z-axis in bytes
620 : !> \par History
621 : !> Created [Nico Holmberg] (11.2017)
622 : ! **************************************************************************************************
623 1850 : SUBROUTINE pw_to_cube_parallel(grid, unit_nr, title, particles_r, particles_z, stride, zero_tails, &
624 : msglen)
625 :
626 : TYPE(pw_r3d_rs_type), INTENT(IN) :: grid
627 : TYPE(mp_file_type), INTENT(IN) :: unit_nr
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
635 :
636 : INTEGER, PARAMETER :: entry_len = 13, header_len = 41, &
637 : header_len_z = 53, num_entries_line = 6
638 :
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, &
643 : ubounds_local
644 1850 : INTEGER, ALLOCATABLE, DIMENSION(:), TARGET :: blocklengths
645 : INTEGER(kind=file_offset), ALLOCATABLE, &
646 1850 : 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 1850 : CHARACTER(LEN=msglen), ALLOCATABLE, DIMENSION(:) :: writebuffer
651 1850 : CHARACTER(LEN=msglen) :: tmp
652 : LOGICAL :: should_write(2)
653 : TYPE(mp_comm_type) :: gid
654 : TYPE(mp_file_descriptor_type) :: mp_desc
655 :
656 : !get rs grids and parallel envnment
657 1850 : gid = grid%pw_grid%para%group
658 1850 : my_rank = grid%pw_grid%para%group%mepos
659 :
660 : ! Shortcut
661 7400 : lbounds = grid%pw_grid%bounds(1, :)
662 7400 : ubounds = grid%pw_grid%bounds(2, :)
663 7400 : lbounds_local = grid%pw_grid%bounds_local(1, :)
664 7400 : ubounds_local = grid%pw_grid%bounds_local(2, :)
665 : ! Determine the total number of z-slices and the number of values per slice
666 1850 : size_of_z = CEILING(REAL(ubounds_local(3) - lbounds_local(3) + 1, dp)/REAL(stride(3), dp))
667 1850 : islice = 1
668 28383 : DO i = lbounds(1), ubounds(1), stride(1)
669 82374 : should_write(:) = .TRUE.
670 27458 : IF (i .LT. lbounds_local(1)) THEN
671 8957 : should_write(1) = .FALSE.
672 18501 : ELSE IF (i .GT. ubounds_local(1)) THEN
673 : EXIT
674 : END IF
675 593956 : DO j = lbounds(2), ubounds(2), stride(2)
676 565573 : should_write(2) = .TRUE.
677 565573 : IF (j .LT. lbounds_local(2) .OR. j .GT. ubounds_local(2)) THEN
678 0 : should_write(2) = .FALSE.
679 : END IF
680 1342298 : IF (ALL(should_write .EQV. .TRUE.)) THEN
681 375096 : islice = islice + 1
682 : END IF
683 : END DO
684 : END DO
685 1850 : nslices = islice - 1
686 37612 : DO k = lbounds(3), ubounds(3), stride(3)
687 37612 : IF (k + stride(3) > ubounds(3)) last_z = k
688 : END DO
689 1850 : islice = 1
690 : ! Determine initial byte offset (0 or EOF if data is appended)
691 1850 : CALL unit_nr%get_position(BOF)
692 : ! Write header information on master process and update byte offset accordingly
693 1850 : IF (my_rank == 0) THEN
694 : ! this format seems to work for e.g. molekel and gOpenmol
695 : ! latest version of VMD can read non orthorhombic cells
696 925 : CALL unit_nr%write_at(BOF, "-Quickstep-"//NEW_LINE("C"))
697 925 : BOF = BOF + LEN("-Quickstep-"//NEW_LINE("C"))*mpi_character_size
698 925 : IF (PRESENT(title)) THEN
699 925 : CALL unit_nr%write_at(BOF, TRIM(title)//NEW_LINE("C"))
700 925 : BOF = BOF + LEN(TRIM(title)//NEW_LINE("C"))*mpi_character_size
701 : ELSE
702 0 : CALL unit_nr%write_at(BOF, "No Title"//NEW_LINE("C"))
703 0 : BOF = BOF + LEN("No Title"//NEW_LINE("C"))*mpi_character_size
704 : END IF
705 :
706 925 : CPASSERT(PRESENT(particles_z) .EQV. PRESENT(particles_r))
707 925 : np = 0
708 925 : IF (PRESENT(particles_z)) THEN
709 925 : CPASSERT(SIZE(particles_z) == SIZE(particles_r, dim=2))
710 : ! cube files can only be written for 99999 particles due to a format limitation (I5)
711 : ! so we limit the number of particles written.
712 925 : np = MIN(99999, SIZE(particles_z))
713 : END IF
714 :
715 925 : WRITE (header, '(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp !start of cube
716 925 : CALL unit_nr%write_at(BOF, header//NEW_LINE("C"))
717 925 : BOF = BOF + LEN(header//NEW_LINE("C"))*mpi_character_size
718 :
719 925 : WRITE (header, '(I5,3f12.6)') (grid%pw_grid%npts(1) + stride(1) - 1)/stride(1), &
720 925 : grid%pw_grid%dh(1, 1)*REAL(stride(1), dp), grid%pw_grid%dh(2, 1)*REAL(stride(1), dp), &
721 1850 : grid%pw_grid%dh(3, 1)*REAL(stride(1), dp)
722 925 : CALL unit_nr%write_at(BOF, header//NEW_LINE("C"))
723 925 : BOF = BOF + LEN(header//NEW_LINE("C"))*mpi_character_size
724 :
725 925 : WRITE (header, '(I5,3f12.6)') (grid%pw_grid%npts(2) + stride(2) - 1)/stride(2), &
726 925 : grid%pw_grid%dh(1, 2)*REAL(stride(2), dp), grid%pw_grid%dh(2, 2)*REAL(stride(2), dp), &
727 1850 : grid%pw_grid%dh(3, 2)*REAL(stride(2), dp)
728 925 : CALL unit_nr%write_at(BOF, header//NEW_LINE("C"))
729 925 : BOF = BOF + LEN(header//NEW_LINE("C"))*mpi_character_size
730 :
731 925 : WRITE (header, '(I5,3f12.6)') (grid%pw_grid%npts(3) + stride(3) - 1)/stride(3), &
732 925 : grid%pw_grid%dh(1, 3)*REAL(stride(3), dp), grid%pw_grid%dh(2, 3)*REAL(stride(3), dp), &
733 1850 : grid%pw_grid%dh(3, 3)*REAL(stride(3), dp)
734 925 : CALL unit_nr%write_at(BOF, header//NEW_LINE("C"))
735 925 : BOF = BOF + LEN(header//NEW_LINE("C"))*mpi_character_size
736 :
737 925 : IF (PRESENT(particles_z)) THEN
738 4067 : DO i = 1, np
739 3142 : WRITE (header_z, '(I5,4f12.6)') particles_z(i), 0._dp, particles_r(:, i)
740 3142 : CALL unit_nr%write_at(BOF, header_z//NEW_LINE("C"))
741 4067 : BOF = BOF + LEN(header_z//NEW_LINE("C"))*mpi_character_size
742 : END DO
743 : END IF
744 : END IF
745 : ! Sync offset
746 1850 : CALL gid%bcast(BOF, grid%pw_grid%para%group%source)
747 : ! Determine byte offsets for each grid z-slice which are local to a process
748 : ! and convert z-slices to cube format compatible strings
749 5550 : ALLOCATE (displacements(nslices))
750 376946 : displacements = 0
751 5550 : ALLOCATE (writebuffer(nslices))
752 376946 : writebuffer(:) = ''
753 28383 : DO i = lbounds(1), ubounds(1), stride(1)
754 82374 : should_write(:) = .TRUE.
755 27458 : IF (i .LT. lbounds_local(1)) THEN
756 8957 : should_write(1) = .FALSE.
757 18501 : ELSE IF (i .GT. ubounds_local(1)) THEN
758 : EXIT
759 : END IF
760 28383 : DO j = lbounds(2), ubounds(2), stride(2)
761 565573 : should_write(2) = .TRUE.
762 565573 : IF (j .LT. lbounds_local(2) .OR. j .GT. ubounds_local(2)) THEN
763 0 : should_write(2) = .FALSE.
764 : END IF
765 1315765 : IF (ALL(should_write .EQV. .TRUE.)) THEN
766 375096 : IF (islice > nslices) CPABORT("Index out of bounds.")
767 375096 : displacements(islice) = BOF
768 375096 : tmp = ''
769 375096 : counter = 0
770 9725443 : DO k = lbounds(3), ubounds(3), stride(3)
771 9350347 : IF (zero_tails .AND. grid%array(i, j, k) < 1.E-7_dp) THEN
772 54882 : WRITE (value, '(E13.5)') 0.0_dp
773 : ELSE
774 9295465 : WRITE (value, '(E13.5)') grid%array(i, j, k)
775 : END IF
776 9350347 : tmp = TRIM(tmp)//TRIM(value)
777 9350347 : counter = counter + 1
778 9350347 : IF (MODULO(counter, num_entries_line) == 0 .OR. k == last_z) &
779 2075781 : tmp = TRIM(tmp)//NEW_LINE('C')
780 : END DO
781 375096 : writebuffer(islice) = tmp
782 375096 : islice = islice + 1
783 : END IF
784 : ! Update global byte offset
785 565573 : BOF = BOF + msglen
786 : END DO
787 : END DO
788 : ! Create indexed MPI type using calculated byte offsets as displacements
789 : ! Size of each z-slice is msglen
790 5550 : ALLOCATE (blocklengths(nslices))
791 376946 : blocklengths(:) = msglen
792 1850 : mp_desc = mp_file_type_hindexed_make_chv(nslices, blocklengths, displacements)
793 : ! Use the created type as a file view
794 : ! NB. The vector 'displacements' contains the absolute offsets of each z-slice i.e.
795 : ! they are given relative to the beginning of the file. The global offset to
796 : ! set_view must therefore be set to 0
797 1850 : BOF = 0
798 1850 : CALL mp_file_type_set_view_chv(unit_nr, BOF, mp_desc)
799 : ! Collective write of cube
800 1850 : CALL unit_nr%write_all(msglen, nslices, writebuffer, mp_desc)
801 : ! Clean up
802 1850 : CALL mp_file_type_free(mp_desc)
803 1850 : DEALLOCATE (writebuffer)
804 1850 : DEALLOCATE (blocklengths, displacements)
805 :
806 1850 : END SUBROUTINE pw_to_cube_parallel
807 :
808 : ! **************************************************************************************************
809 : !> \brief Prints a simple grid file: X Y Z value
810 : !> \param pw ...
811 : !> \param unit_nr ...
812 : !> \param stride ...
813 : !> \param pw2 ...
814 : !> \par History
815 : !> Created [Vladimir Rybkin] (08.2018)
816 : !> \author Vladimir Rybkin
817 : ! **************************************************************************************************
818 16 : SUBROUTINE pw_to_simple_volumetric(pw, unit_nr, stride, pw2)
819 : TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
820 : INTEGER, INTENT(IN) :: unit_nr
821 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: stride
822 : TYPE(pw_r3d_rs_type), INTENT(IN), OPTIONAL :: pw2
823 :
824 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_to_simple_volumetric'
825 :
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
830 : LOGICAL :: DOUBLE
831 : REAL(KIND=dp) :: x, y, z
832 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buf, buf2
833 : TYPE(mp_comm_type) :: gid
834 :
835 16 : CALL timeset(routineN, handle)
836 :
837 : ! Check if we write two grids
838 16 : DOUBLE = .FALSE.
839 16 : IF (PRESENT(pw2)) DOUBLE = .TRUE.
840 :
841 64 : my_stride = 1
842 16 : IF (PRESENT(stride)) THEN
843 16 : IF (SIZE(stride) /= 1 .AND. SIZE(stride) /= 3) &
844 : CALL cp_abort(__LOCATION__, "STRIDE keyword can accept only 1 "// &
845 16 : "(the same for X,Y,Z) or 3 values. Correct your input file.")
846 16 : IF (SIZE(stride) == 1) THEN
847 0 : DO i = 1, 3
848 0 : my_stride(i) = stride(1)
849 : END DO
850 : ELSE
851 64 : my_stride = stride(1:3)
852 : END IF
853 16 : CPASSERT(my_stride(1) > 0)
854 16 : CPASSERT(my_stride(2) > 0)
855 16 : CPASSERT(my_stride(3) > 0)
856 : END IF
857 :
858 : ! shortcut
859 16 : L1 = pw%pw_grid%bounds(1, 1)
860 16 : L2 = pw%pw_grid%bounds(1, 2)
861 16 : L3 = pw%pw_grid%bounds(1, 3)
862 16 : U1 = pw%pw_grid%bounds(2, 1)
863 16 : U2 = pw%pw_grid%bounds(2, 2)
864 16 : U3 = pw%pw_grid%bounds(2, 3)
865 :
866 : ! Write the header: number of points and number of spins
867 16 : ngrids = 1
868 16 : 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 16 : ((pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(1))
872 16 : IF (unit_nr > 1) WRITE (unit_nr, '(I7,I5)') npoints, ngrids
873 :
874 48 : ALLOCATE (buf(L3:U3))
875 16 : IF (DOUBLE) ALLOCATE (buf2(L3:U3))
876 :
877 16 : my_rank = pw%pw_grid%para%group%mepos
878 16 : gid = pw%pw_grid%para%group
879 16 : num_pe = pw%pw_grid%para%group%num_pe
880 16 : tag = 1
881 :
882 16 : rank (1) = unit_nr
883 16 : rank (2) = my_rank
884 16 : checksum = 0
885 16 : IF (unit_nr > 0) checksum = 1
886 :
887 16 : CALL gid%sum(checksum)
888 16 : CPASSERT(checksum == 1)
889 :
890 16 : CALL gid%maxloc(rank)
891 16 : CPASSERT(rank(1) > 0)
892 :
893 16 : dest = rank(2)
894 500 : DO I1 = L1, U1, my_stride(1)
895 15288 : DO I2 = L2, U2, my_stride(2)
896 :
897 : ! cycling through the CPUs, check if the current ray (I1,I2) is local to that CPU
898 14788 : IF (pw%pw_grid%para%mode .NE. PW_MODE_LOCAL) THEN
899 44364 : 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 44364 : 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
902 14788 : source = ip
903 : END IF
904 : END DO
905 : ELSE
906 0 : source = dest
907 : END IF
908 :
909 14788 : IF (source == dest) THEN
910 7444 : IF (my_rank == source) THEN
911 118276 : buf(:) = pw%array(I1, I2, :)
912 3722 : IF (DOUBLE) buf2(:) = pw2%array(I1, I2, :)
913 : END IF
914 : ELSE
915 7344 : IF (my_rank == source) THEN
916 116976 : buf(:) = pw%array(I1, I2, :)
917 3672 : CALL gid%send(buf, dest, tag)
918 3672 : IF (DOUBLE) THEN
919 0 : buf2(:) = pw2%array(I1, I2, :)
920 0 : CALL gid%send(buf2, dest, tag)
921 : END IF
922 : END IF
923 7344 : IF (my_rank == dest) THEN
924 3672 : CALL gid%recv(buf, source, tag)
925 3672 : IF (DOUBLE) CALL gid%recv(buf2, source, tag)
926 : END IF
927 : END IF
928 :
929 7394 : IF (.NOT. DOUBLE) THEN
930 470504 : 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 455716 : pw%pw_grid%dh(3, 1)*I3
934 :
935 : y = pw%pw_grid%dh(1, 2)*I1 + &
936 : pw%pw_grid%dh(2, 2)*I2 + &
937 455716 : pw%pw_grid%dh(3, 2)*I3
938 :
939 : z = pw%pw_grid%dh(1, 3)*I1 + &
940 : pw%pw_grid%dh(2, 3)*I2 + &
941 455716 : pw%pw_grid%dh(3, 3)*I3
942 :
943 470504 : IF (unit_nr > 0) THEN
944 227858 : WRITE (unit_nr, '(6E13.5, 6E13.5, 6E13.5, 6E13.5)') x, y, z, buf(I3)
945 : END IF
946 : END DO
947 :
948 : ELSE
949 :
950 0 : 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 0 : pw%pw_grid%dh(3, 1)*I3
954 :
955 : y = pw%pw_grid%dh(1, 2)*I1 + &
956 : pw%pw_grid%dh(2, 2)*I2 + &
957 0 : pw%pw_grid%dh(3, 2)*I3
958 :
959 : z = pw%pw_grid%dh(1, 3)*I1 + &
960 : pw%pw_grid%dh(2, 3)*I2 + &
961 0 : pw%pw_grid%dh(3, 3)*I3
962 :
963 0 : IF (unit_nr > 0) THEN
964 0 : WRITE (unit_nr, '(6E13.5, 6E13.5, 6E13.5, 6E13.5)') x, y, z, buf(I3), buf2(I3)
965 : END IF
966 : END DO
967 :
968 : END IF ! Double
969 :
970 : ! this double loop generates so many messages that it can overload
971 : ! the message passing system, e.g. on XT3
972 : ! we therefore put a barrier here that limits the amount of message
973 : ! that flies around at any given time.
974 : ! if ever this routine becomes a bottleneck, we should go for a
975 : ! more complicated rewrite
976 15272 : CALL gid%sync()
977 :
978 : END DO
979 : END DO
980 :
981 16 : DEALLOCATE (buf)
982 16 : IF (DOUBLE) DEALLOCATE (buf2)
983 :
984 16 : CALL timestop(handle)
985 :
986 32 : END SUBROUTINE pw_to_simple_volumetric
987 :
988 : END MODULE realspace_grid_cube
|
