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