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