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