(git:24d69ee)
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! **************************************************************************************************
11MODULE 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
42CONTAINS
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 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 READ (values, *, iostat=readstat) buffer
57 IF (readstat == 0) RETURN
58
59 pos = 1
60 DO i = 1, SIZE(buffer)
61 IF (pos + cube_entry_len - 1 > len(values)) cpabort("Unexpected end of cube data.")
62 value = values(pos:pos + cube_entry_len - 1)
63 READ (value, '(E13.5)', iostat=readstat) buffer(i)
64 IF (readstat /= 0) cpabort("Bad value while reading cube data.")
65 pos = pos + cube_entry_len
66 IF (modulo(i, cube_num_entries_line) == 0) THEN
67 IF (pos <= len(values)) THEN
68 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buf
110 TYPE(mp_comm_type) :: gid
111 TYPE(mp_file_type) :: mp_unit
112
113 CALL timeset(routinen, handle)
114
115 my_zero_tails = .false.
116 be_silent = .false.
117 parallel_write = .false.
118 my_max_file_size_mb = 0.0_dp
119 IF (PRESENT(zero_tails)) my_zero_tails = zero_tails
120 IF (PRESENT(silent)) be_silent = silent
121 IF (PRESENT(mpi_io)) parallel_write = mpi_io
122 IF (PRESENT(max_file_size_mb)) my_max_file_size_mb = max_file_size_mb
123 cpassert(my_max_file_size_mb >= 0)
124 IF (.NOT. parallel_write .AND. PRESENT(particles_z) .AND. &
125 .NOT. PRESENT(particles_zeff)) THEN
126 cpwarn(missing_zeff_warning)
127 END IF
128
129 my_stride = 1
130 IF (PRESENT(stride)) THEN
131 IF (SIZE(stride) /= 1 .AND. SIZE(stride) /= 3) THEN
132 CALL cp_abort(__location__, "STRIDE keyword can accept only 1 "// &
133 "(the same for X,Y,Z) or 3 values. Correct your input file.")
134 END IF
135 IF (SIZE(stride) == 1) THEN
136 DO i = 1, 3
137 my_stride(i) = stride(1)
138 END DO
139 ELSE
140 my_stride = stride(1:3)
141 END IF
142 END IF
143
144 IF (my_max_file_size_mb > 0) THEN
145 ! A single grid point takes up 13 bytes, which is 1.3e-05 MB.
146 compression_factor = 1.3e-05_dp*product(real(pw%pw_grid%npts, dp))/max_file_size_mb
147 my_stride(:) = int(compression_factor**(1.0/3.0)) + 1
148 END IF
149
150 cpassert(my_stride(1) > 0)
151 cpassert(my_stride(2) > 0)
152 cpassert(my_stride(3) > 0)
153
154 IF (.NOT. parallel_write) THEN
155 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 WRITE (unit_nr, '(a11)') "-Quickstep-"
159 IF (PRESENT(title)) THEN
160 WRITE (unit_nr, *) trim(title)
161 ELSE
162 WRITE (unit_nr, *) "No Title"
163 END IF
164
165 cpassert(PRESENT(particles_z) .EQV. PRESENT(particles_r))
166 np = 0
167 IF (PRESENT(particles_z)) THEN
168 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 np = min(99999, SIZE(particles_z))
172 END IF
173
174 WRITE (unit_nr, '(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp !start of cube
175
176 WRITE (unit_nr, '(I5,3f12.6)') (pw%pw_grid%npts(1) + my_stride(1) - 1)/my_stride(1), &
177 pw%pw_grid%dh(1, 1)*real(my_stride(1), dp), pw%pw_grid%dh(2, 1)*real(my_stride(1), dp), &
178 pw%pw_grid%dh(3, 1)*real(my_stride(1), dp)
179 WRITE (unit_nr, '(I5,3f12.6)') (pw%pw_grid%npts(2) + my_stride(2) - 1)/my_stride(2), &
180 pw%pw_grid%dh(1, 2)*real(my_stride(2), dp), pw%pw_grid%dh(2, 2)*real(my_stride(2), dp), &
181 pw%pw_grid%dh(3, 2)*real(my_stride(2), dp)
182 WRITE (unit_nr, '(I5,3f12.6)') (pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(3), &
183 pw%pw_grid%dh(1, 3)*real(my_stride(3), dp), pw%pw_grid%dh(2, 3)*real(my_stride(3), dp), &
184 pw%pw_grid%dh(3, 3)*real(my_stride(3), dp)
185
186 IF (PRESENT(particles_z)) THEN
187 IF (PRESENT(particles_zeff)) THEN
188 DO iat = 1, np
189 WRITE (unit_nr, '(I5,4f12.6)') particles_z(iat), particles_zeff(iat), particles_r(:, iat)
190 END DO
191 ELSE
192 DO iat = 1, np
193 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 l1 = pw%pw_grid%bounds(1, 1)
201 l2 = pw%pw_grid%bounds(1, 2)
202 l3 = pw%pw_grid%bounds(1, 3)
203 u1 = pw%pw_grid%bounds(2, 1)
204 u2 = pw%pw_grid%bounds(2, 2)
205 u3 = pw%pw_grid%bounds(2, 3)
206
207 ALLOCATE (buf(l3:u3))
208
209 my_rank = pw%pw_grid%para%group%mepos
210 gid = pw%pw_grid%para%group
211 num_pe = pw%pw_grid%para%group%num_pe
212 tag = 1
213
214 rank(1) = unit_nr
215 rank(2) = my_rank
216 checksum = 0
217 IF (unit_nr > 0) checksum = 1
218
219 CALL gid%sum(checksum)
220 cpassert(checksum == 1)
221
222 CALL gid%maxloc(rank)
223 cpassert(rank(1) > 0)
224
225 dest = rank(2)
226 DO i1 = l1, u1, my_stride(1)
227 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 IF (pw%pw_grid%para%mode /= pw_mode_local) THEN
231 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 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 source = ip
235 END IF
236 END DO
237 ELSE
238 source = dest
239 END IF
240
241 IF (source == dest) THEN
242 IF (my_rank == source) THEN
243 buf(:) = pw%array(i1, i2, :)
244 END IF
245 ELSE
246 IF (my_rank == source) THEN
247 buf(:) = pw%array(i1, i2, :)
248 CALL gid%send(buf, dest, tag)
249 END IF
250 IF (my_rank == dest) THEN
251 CALL gid%recv(buf, source, tag)
252 END IF
253 END IF
254
255 IF (unit_nr > 0) THEN
256 IF (my_zero_tails) THEN
257 DO i3 = l3, u3
258 IF (buf(i3) < 1.e-7_dp) buf(i3) = 0.0_dp
259 END DO
260 END IF
261 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 CALL gid%sync()
271
272 END DO
273 END DO
274
275 DEALLOCATE (buf)
276 ELSE
277 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 num_linebreak = size_of_z/num_entries_line
279 IF (modulo(size_of_z, num_entries_line) /= 0) THEN
280 num_linebreak = num_linebreak + 1
281 END IF
282 msglen = (size_of_z*entry_len + num_linebreak)*mpi_character_size
283 CALL mp_unit%set_handle(unit_nr)
284 CALL pw_to_cube_parallel(pw, mp_unit, title, particles_r, particles_z, particles_zeff, &
285 my_stride, my_zero_tails, msglen)
286 END IF
287
288 CALL timestop(handle)
289
290 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buffer
324 REAL(kind=dp), DIMENSION(3) :: dr, rdum
325 TYPE(mp_comm_type) :: gid
326
327 output_unit = cp_logger_get_default_io_unit()
328
329 CALL timeset(routinen, handle)
330
331 be_silent = .false.
332 IF (PRESENT(silent)) THEN
333 be_silent = silent
334 END IF
335 !get rs grids and parallel environment
336 gid = grid%pw_grid%para%group
337 my_rank = grid%pw_grid%para%group%mepos
338 num_pe = grid%pw_grid%para%group%num_pe
339 tag = 1
340
341 lbounds_local = grid%pw_grid%bounds_local(1, :)
342 ubounds_local = grid%pw_grid%bounds_local(2, :)
343 size_of_z = ubounds_local(3) - lbounds_local(3) + 1
344
345 IF (.NOT. parallel_read) THEN
346 npoints = grid%pw_grid%npts
347 lbounds = grid%pw_grid%bounds(1, :)
348 ubounds = grid%pw_grid%bounds(2, :)
349
350 DO i = 1, 3
351 dr(i) = grid%pw_grid%dh(i, i)
352 END DO
353
354 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 ALLOCATE (buffer(lbounds(3):ubounds(3)))
357
358 IF (my_rank == 0) THEN
359 IF (output_unit > 0 .AND. .NOT. be_silent) THEN
360 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 unit_number=extunit)
368
369 !skip header comments
370 DO i = 1, 2
371 READ (extunit, *)
372 END DO
373 READ (extunit, *) nat, rdum
374 DO i = 1, 3
375 READ (extunit, *) ndum, rdum
376 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
377 output_unit > 0) THEN
378 WRITE (output_unit, *) "Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
379 WRITE (output_unit, *) "Restart from density | ", ndum, " DIFFERS FROM ", npoints(i)
380 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 DO i = 1, nat
385 READ (extunit, *)
386 END DO
387 END IF
388
389 !master sends all data to everyone
390 DO i = lbounds(1), ubounds(1)
391 DO j = lbounds(2), ubounds(2)
392 IF (my_rank == 0) THEN
393 DO k = lbounds(3), ubounds(3), cube_num_entries_line
394 last_z = min(k + cube_num_entries_line - 1, ubounds(3))
395 READ (extunit, '(A)') value_line
396 CALL cube_read_values(value_line, buffer(k:last_z))
397 END DO
398 END IF
399 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 .AND. (j <= ubounds_local(2))) THEN
404 !allow scaling of external potential values by factor 'scaling' (SCALING_FACTOR in input file)
405 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 IF (my_rank == 0) CALL close_file(unit_number=extunit)
412
413 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 num_linebreak = size_of_z/num_entries_line
421 IF (modulo(size_of_z, num_entries_line) /= 0) THEN
422 num_linebreak = num_linebreak + 1
423 END IF
424 msglen = (size_of_z*entry_len + num_linebreak)*mpi_character_size
425 CALL cube_to_pw_parallel(grid, filename, scaling, msglen, silent=silent)
426 END IF
427
428 CALL timestop(handle)
429
430 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 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 INTEGER, ALLOCATABLE, DIMENSION(:), TARGET :: blocklengths
455 INTEGER(kind=file_offset), ALLOCATABLE, &
456 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 CHARACTER(LEN=msglen), ALLOCATABLE, DIMENSION(:) :: readbuffer
463 LOGICAL :: be_silent, should_read(2)
464 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 output_unit = cp_logger_get_default_io_unit()
471
472 be_silent = .false.
473 IF (PRESENT(silent)) THEN
474 be_silent = silent
475 END IF
476
477 !get rs grids and parallel envnment
478 gid = grid%pw_grid%para%group
479 my_rank = grid%pw_grid%para%group%mepos
480 num_pe = grid%pw_grid%para%group%num_pe
481 tag = 1
482
483 DO i = 1, 3
484 dr(i) = grid%pw_grid%dh(i, i)
485 END DO
486
487 npoints = grid%pw_grid%npts
488 lbounds = grid%pw_grid%bounds(1, :)
489 ubounds = grid%pw_grid%bounds(2, :)
490
491 npoints_local = grid%pw_grid%npts_local
492 lbounds_local = grid%pw_grid%bounds_local(1, :)
493 ubounds_local = grid%pw_grid%bounds_local(2, :)
494 size_of_z = ubounds_local(3) - lbounds_local(3) + 1
495 nslices = (ubounds_local(1) - lbounds_local(1) + 1)*(ubounds_local(2) - lbounds_local(2) + 1)
496 islice = 1
497
498 ! Read header information and determine byte offset of cube data on master process
499 IF (my_rank == 0) THEN
500 IF (output_unit > 0 .AND. .NOT. be_silent) THEN
501 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 unit_number=extunit_handle)
510
511 !skip header comments
512 DO i = 1, 2
513 READ (extunit_handle, *)
514 END DO
515 READ (extunit_handle, *) nat, rdum
516 DO i = 1, 3
517 READ (extunit_handle, *) ndum, rdum
518 IF ((ndum /= npoints(i) .OR. (abs(rdum(i) - dr(i)) > 1e-4)) .AND. &
519 output_unit > 0) THEN
520 WRITE (output_unit, *) "Restart from density | ERROR! | CUBE FILE NOT COINCIDENT WITH INTERNAL GRID ", i
521 WRITE (output_unit, *) "Restart from density | ", ndum, " DIFFERS FROM ", npoints(i)
522 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 DO i = 1, nat
527 READ (extunit_handle, *)
528 END DO
529 ! Get byte offset
530 INQUIRE (extunit_handle, pos=offset_global)
531 CALL close_file(unit_number=extunit_handle)
532 END IF
533 ! Sync offset and start parallel read
534 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 bof = offset_global - 1
538 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 ALLOCATE (displacements(nslices))
541 displacements = 0
542 DO i = lbounds(1), ubounds(1)
543 should_read(:) = .true.
544 IF (i < lbounds_local(1)) THEN
545 should_read(1) = .false.
546 ELSE IF (i > ubounds_local(1)) THEN
547 EXIT
548 END IF
549 DO j = lbounds(2), ubounds(2)
550 should_read(2) = .true.
551 IF (j < lbounds_local(2) .OR. j > ubounds_local(2)) THEN
552 should_read(2) = .false.
553 END IF
554 IF (all(should_read .EQV. .true.)) THEN
555 IF (islice > nslices) cpabort("Index out of bounds.")
556 displacements(islice) = bof
557 islice = islice + 1
558 END IF
559 ! Update global byte offset
560 bof = bof + msglen
561 END DO
562 END DO
563 ! Size of each z-slice is msglen
564 ALLOCATE (blocklengths(nslices))
565 blocklengths(:) = msglen
566 ! Create indexed MPI type using calculated byte offsets as displacements and use it as a file view
567 mp_file_desc = mp_file_type_hindexed_make_chv(nslices, blocklengths, displacements)
568 bof = 0
569 CALL mp_file_type_set_view_chv(extunit, bof, mp_file_desc)
570 ! Collective read of cube
571 ALLOCATE (readbuffer(nslices))
572 readbuffer(:) = ''
573 CALL extunit%read_all(msglen, nslices, readbuffer, mp_file_desc)
574 CALL mp_file_type_free(mp_file_desc)
575 CALL extunit%close()
576 ! Convert cube values string -> real
577 i = lbounds_local(1)
578 j = lbounds_local(2)
579 ALLOCATE (buffer(lbounds(3):ubounds(3)))
580 buffer = 0.0_dp
581 DO islice = 1, nslices
582 CALL cube_read_values(readbuffer(islice), buffer(lbounds(3):ubounds(3)))
583 ! Optionally scale cube file values
584 grid%array(i, j, lbounds(3):ubounds(3)) = scaling*buffer(lbounds(3):ubounds(3))
585 j = j + 1
586 IF (j > ubounds_local(2)) THEN
587 j = lbounds_local(2)
588 i = i + 1
589 END IF
590 END DO
591 DEALLOCATE (readbuffer)
592 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 DEALLOCATE (buffer)
657
658 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 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 INTEGER, ALLOCATABLE, DIMENSION(:), TARGET :: blocklengths
697 INTEGER(kind=file_offset), ALLOCATABLE, &
698 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 CHARACTER(LEN=msglen), ALLOCATABLE, DIMENSION(:) :: writebuffer
703 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 gid = grid%pw_grid%para%group
710 my_rank = grid%pw_grid%para%group%mepos
711 IF (PRESENT(particles_z) .AND. .NOT. PRESENT(particles_zeff)) THEN
712 cpwarn(missing_zeff_warning)
713 END IF
714
715 ! Shortcut
716 lbounds = grid%pw_grid%bounds(1, :)
717 ubounds = grid%pw_grid%bounds(2, :)
718 lbounds_local = grid%pw_grid%bounds_local(1, :)
719 ubounds_local = grid%pw_grid%bounds_local(2, :)
720 ! Determine the total number of z-slices and the number of values per slice
721 size_of_z = ceiling(real(ubounds_local(3) - lbounds_local(3) + 1, dp)/real(stride(3), dp))
722 islice = 1
723 DO i = lbounds(1), ubounds(1), stride(1)
724 should_write(:) = .true.
725 IF (i < lbounds_local(1)) THEN
726 should_write(1) = .false.
727 ELSE IF (i > ubounds_local(1)) THEN
728 EXIT
729 END IF
730 DO j = lbounds(2), ubounds(2), stride(2)
731 should_write(2) = .true.
732 IF (j < lbounds_local(2) .OR. j > ubounds_local(2)) THEN
733 should_write(2) = .false.
734 END IF
735 IF (all(should_write .EQV. .true.)) THEN
736 islice = islice + 1
737 END IF
738 END DO
739 END DO
740 nslices = islice - 1
741 DO k = lbounds(3), ubounds(3), stride(3)
742 IF (k + stride(3) > ubounds(3)) last_z = k
743 END DO
744 islice = 1
745 ! Determine initial byte offset (0 or EOF if data is appended)
746 CALL unit_nr%get_position(bof)
747 ! Write header information on master process and update byte offset accordingly
748 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 CALL unit_nr%write_at(bof, "-Quickstep-"//new_line("C"))
752 bof = bof + len("-Quickstep-"//new_line("C"))*mpi_character_size
753 IF (PRESENT(title)) THEN
754 CALL unit_nr%write_at(bof, trim(title)//new_line("C"))
755 bof = bof + len(trim(title)//new_line("C"))*mpi_character_size
756 ELSE
757 CALL unit_nr%write_at(bof, "No Title"//new_line("C"))
758 bof = bof + len("No Title"//new_line("C"))*mpi_character_size
759 END IF
760
761 cpassert(PRESENT(particles_z) .EQV. PRESENT(particles_r))
762 np = 0
763 IF (PRESENT(particles_z)) THEN
764 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 np = min(99999, SIZE(particles_z))
768 END IF
769
770 WRITE (header, '(I5,3f12.6)') np, 0.0_dp, 0._dp, 0._dp !start of cube
771 CALL unit_nr%write_at(bof, header//new_line("C"))
772 bof = bof + len(header//new_line("C"))*mpi_character_size
773
774 WRITE (header, '(I5,3f12.6)') (grid%pw_grid%npts(1) + stride(1) - 1)/stride(1), &
775 grid%pw_grid%dh(1, 1)*real(stride(1), dp), grid%pw_grid%dh(2, 1)*real(stride(1), dp), &
776 grid%pw_grid%dh(3, 1)*real(stride(1), dp)
777 CALL unit_nr%write_at(bof, header//new_line("C"))
778 bof = bof + len(header//new_line("C"))*mpi_character_size
779
780 WRITE (header, '(I5,3f12.6)') (grid%pw_grid%npts(2) + stride(2) - 1)/stride(2), &
781 grid%pw_grid%dh(1, 2)*real(stride(2), dp), grid%pw_grid%dh(2, 2)*real(stride(2), dp), &
782 grid%pw_grid%dh(3, 2)*real(stride(2), dp)
783 CALL unit_nr%write_at(bof, header//new_line("C"))
784 bof = bof + len(header//new_line("C"))*mpi_character_size
785
786 WRITE (header, '(I5,3f12.6)') (grid%pw_grid%npts(3) + stride(3) - 1)/stride(3), &
787 grid%pw_grid%dh(1, 3)*real(stride(3), dp), grid%pw_grid%dh(2, 3)*real(stride(3), dp), &
788 grid%pw_grid%dh(3, 3)*real(stride(3), dp)
789 CALL unit_nr%write_at(bof, header//new_line("C"))
790 bof = bof + len(header//new_line("C"))*mpi_character_size
791
792 IF (PRESENT(particles_z)) THEN
793 IF (PRESENT(particles_zeff)) THEN
794 DO i = 1, np
795 WRITE (header_z, '(I5,4f12.6)') particles_z(i), particles_zeff(i), particles_r(:, i)
796 CALL unit_nr%write_at(bof, header_z//new_line("C"))
797 bof = bof + len(header_z//new_line("C"))*mpi_character_size
798 END DO
799 ELSE
800 DO i = 1, np
801 WRITE (header_z, '(I5,4f12.6)') particles_z(i), 0._dp, particles_r(:, i)
802 CALL unit_nr%write_at(bof, header_z//new_line("C"))
803 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 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 ALLOCATE (displacements(nslices))
813 displacements = 0
814 ALLOCATE (writebuffer(nslices))
815 writebuffer(:) = ''
816 DO i = lbounds(1), ubounds(1), stride(1)
817 should_write(:) = .true.
818 IF (i < lbounds_local(1)) THEN
819 should_write(1) = .false.
820 ELSE IF (i > ubounds_local(1)) THEN
821 EXIT
822 END IF
823 DO j = lbounds(2), ubounds(2), stride(2)
824 should_write(2) = .true.
825 IF (j < lbounds_local(2) .OR. j > ubounds_local(2)) THEN
826 should_write(2) = .false.
827 END IF
828 IF (all(should_write .EQV. .true.)) THEN
829 IF (islice > nslices) cpabort("Index out of bounds.")
830 displacements(islice) = bof
831 tmp = ''
832 counter = 0
833 DO k = lbounds(3), ubounds(3), stride(3)
834 IF (zero_tails .AND. grid%array(i, j, k) < 1.e-7_dp) THEN
835 WRITE (value, cube_value_format) 0.0_dp
836 ELSE
837 WRITE (value, cube_value_format) grid%array(i, j, k)
838 END IF
839 tmp = trim(tmp)//trim(value)
840 counter = counter + 1
841 IF (modulo(counter, num_entries_line) == 0 .OR. k == last_z) THEN
842 tmp = trim(tmp)//new_line('C')
843 END IF
844 END DO
845 writebuffer(islice) = tmp
846 islice = islice + 1
847 END IF
848 ! Update global byte offset
849 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 ALLOCATE (blocklengths(nslices))
855 blocklengths(:) = msglen
856 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 bof = 0
862 CALL mp_file_type_set_view_chv(unit_nr, bof, mp_desc)
863 ! Collective write of cube
864 CALL unit_nr%write_all(msglen, nslices, writebuffer, mp_desc)
865 ! Clean up
866 CALL mp_file_type_free(mp_desc)
867 DEALLOCATE (writebuffer)
868 DEALLOCATE (blocklengths, displacements)
869
870 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buf, buf2
897 TYPE(mp_comm_type) :: gid
898
899 CALL timeset(routinen, handle)
900
901 ! Check if we write two grids
902 double = .false.
903 IF (PRESENT(pw2)) double = .true.
904
905 my_stride = 1
906 IF (PRESENT(stride)) THEN
907 IF (SIZE(stride) /= 1 .AND. SIZE(stride) /= 3) THEN
908 CALL cp_abort(__location__, "STRIDE keyword can accept only 1 "// &
909 "(the same for X,Y,Z) or 3 values. Correct your input file.")
910 END IF
911 IF (SIZE(stride) == 1) THEN
912 DO i = 1, 3
913 my_stride(i) = stride(1)
914 END DO
915 ELSE
916 my_stride = stride(1:3)
917 END IF
918 cpassert(my_stride(1) > 0)
919 cpassert(my_stride(2) > 0)
920 cpassert(my_stride(3) > 0)
921 END IF
922
923 ! shortcut
924 l1 = pw%pw_grid%bounds(1, 1)
925 l2 = pw%pw_grid%bounds(1, 2)
926 l3 = pw%pw_grid%bounds(1, 3)
927 u1 = pw%pw_grid%bounds(2, 1)
928 u2 = pw%pw_grid%bounds(2, 2)
929 u3 = pw%pw_grid%bounds(2, 3)
930
931 ! Write the header: number of points and number of spins
932 ngrids = 1
933 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 ((pw%pw_grid%npts(3) + my_stride(3) - 1)/my_stride(1))
937 IF (unit_nr > 1) WRITE (unit_nr, '(I7,I5)') npoints, ngrids
938
939 ALLOCATE (buf(l3:u3))
940 IF (double) ALLOCATE (buf2(l3:u3))
941
942 my_rank = pw%pw_grid%para%group%mepos
943 gid = pw%pw_grid%para%group
944 num_pe = pw%pw_grid%para%group%num_pe
945 tag = 1
946
947 rank(1) = unit_nr
948 rank(2) = my_rank
949 checksum = 0
950 IF (unit_nr > 0) checksum = 1
951
952 CALL gid%sum(checksum)
953 cpassert(checksum == 1)
954
955 CALL gid%maxloc(rank)
956 cpassert(rank(1) > 0)
957
958 dest = rank(2)
959 DO i1 = l1, u1, my_stride(1)
960 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 IF (pw%pw_grid%para%mode /= pw_mode_local) THEN
964 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 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 source = ip
968 END IF
969 END DO
970 ELSE
971 source = dest
972 END IF
973
974 IF (source == dest) THEN
975 IF (my_rank == source) THEN
976 buf(:) = pw%array(i1, i2, :)
977 IF (double) buf2(:) = pw2%array(i1, i2, :)
978 END IF
979 ELSE
980 IF (my_rank == source) THEN
981 buf(:) = pw%array(i1, i2, :)
982 CALL gid%send(buf, dest, tag)
983 IF (double) THEN
984 buf2(:) = pw2%array(i1, i2, :)
985 CALL gid%send(buf2, dest, tag)
986 END IF
987 END IF
988 IF (my_rank == dest) THEN
989 CALL gid%recv(buf, source, tag)
990 IF (double) CALL gid%recv(buf2, source, tag)
991 END IF
992 END IF
993
994 IF (.NOT. double) THEN
995 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 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 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 pw%pw_grid%dh(3, 3)*i3
1007
1008 IF (unit_nr > 0) THEN
1009 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 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 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 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 pw%pw_grid%dh(3, 3)*i3
1027
1028 IF (unit_nr > 0) THEN
1029 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 CALL gid%sync()
1042
1043 END DO
1044 END DO
1045
1046 DEALLOCATE (buf)
1047 IF (double) DEALLOCATE (buf2)
1048
1049 CALL timestop(handle)
1050
1051 END SUBROUTINE pw_to_simple_volumetric
1052
1053END MODULE realspace_grid_cube
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
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
cp_files::close_file
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
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
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...
Definition cp_log_handling.F:523
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
message_passing::mp_file_type_hindexed_make_chv
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)
Definition message_passing.F:4453
message_passing::mp_file_type_free
subroutine, public mp_file_type_free(type_descriptor)
Releases the type used for MPI I/O.
Definition message_passing.F:4642
message_passing::mpi_character_size
integer, parameter, public mpi_character_size
Definition message_passing.F:142
message_passing::file_offset
integer, parameter, public file_offset
Definition message_passing.F:130
message_passing::file_amode_rdonly
integer, parameter, public file_amode_rdonly
Definition message_passing.F:133
message_passing::mp_file_type_set_view_chv
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_...
Definition message_passing.F:4496
pw_grid_types
Definition pw_grid_types.F:13
pw_grid_types::pw_mode_local
integer, parameter, public pw_mode_local
Definition pw_grid_types.F:27
pw_types
Definition pw_types.F:24
realspace_grid_cube
Generate Gaussian cube files.
Definition realspace_grid_cube.F:11
realspace_grid_cube::cube_to_pw
subroutine, public cube_to_pw(grid, filename, scaling, parallel_read, silent)
Computes the external density on the grid hacked from external_read_density.
Definition realspace_grid_cube.F:305
realspace_grid_cube::cube_values_format
character(len= *), parameter, public cube_values_format
Definition realspace_grid_cube.F:35
realspace_grid_cube::pw_to_simple_volumetric
subroutine, public pw_to_simple_volumetric(pw, unit_nr, stride, pw2)
Prints a simple grid file: X Y Z value.
Definition realspace_grid_cube.F:883
realspace_grid_cube::cube_value_format
character(len= *), parameter, public cube_value_format
Definition realspace_grid_cube.F:35
realspace_grid_cube::pw_to_cube
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)
...
Definition realspace_grid_cube.F:91
realspace_grid_cube::cube_read_values
subroutine, public cube_read_values(values, buffer)
Read cube values from a character buffer.
Definition realspace_grid_cube.F:50
writebuffer
void writebuffer(int *psockfd, char *data, int *plen)
Writes to a socket.
Definition sockets.c:203
readbuffer
void readbuffer(int *psockfd, char *data, int *plen)
Reads from a socket.
Definition sockets.c:221
message_passing::mp_comm_type
Definition message_passing.F:178
message_passing::mp_file_descriptor_type
Definition message_passing.F:971
message_passing::mp_file_type
Definition message_passing.F:673
pw_types::pw_r3d_rs_type
Definition pw_types.F:62