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