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