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