(git:374b731)
Loading...
Searching...
No Matches
distribution_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Distribution methods for atoms, particles, or molecules
10!> \par History
11!> - 1d-distribution of molecules and particles (Sep. 2003, MK)
12!> - 2d-distribution for Quickstep updated with molecules (Oct. 2003, MK)
13!> \author MK (22.08.2003)
14! **************************************************************************************************
21 USE cell_types, ONLY: cell_type,&
22 pbc,&
31 USE cp_min_heap, ONLY: cp_heap_fill,&
37 USE cp_output_handling, ONLY: cp_p_file,&
41 USE dbcsr_api, ONLY: dbcsr_distribution_get_num_images
52 USE kinds, ONLY: dp,&
53 int_8
54 USE machine, ONLY: m_flush
55 USE mathconstants, ONLY: pi
56 USE mathlib, ONLY: gcd,&
57 lcm
62 USE parallel_rng_types, ONLY: uniform,&
65 USE qs_kind_types, ONLY: get_qs_kind,&
67 USE util, ONLY: sort
68#include "./base/base_uses.f90"
69
70 IMPLICIT NONE
71
72 PRIVATE
73
74! *** Global parameters (in this module) ***
75
76 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'distribution_methods'
77
78! *** Public subroutines ***
79
80 PUBLIC :: distribute_molecules_1d, &
82
83CONTAINS
84
85! **************************************************************************************************
86!> \brief Distribute molecules and particles
87!> \param atomic_kind_set particle (atomic) kind information
88!> \param particle_set particle information
89!> \param local_particles distribution of particles created by this routine
90!> \param molecule_kind_set molecule kind information
91!> \param molecule_set molecule information
92!> \param local_molecules distribution of molecules created by this routine
93!> \param force_env_section ...
94!> \param prev_molecule_kind_set previous molecule kind information, used with
95!> prev_local_molecules
96!> \param prev_local_molecules previous distribution of molecules, new one will
97!> be identical if all the prev_* arguments are present and associated
98!> \par History
99!> none
100!> \author MK (Jun. 2003)
101! **************************************************************************************************
102 SUBROUTINE distribute_molecules_1d(atomic_kind_set, particle_set, &
103 local_particles, &
104 molecule_kind_set, molecule_set, &
105 local_molecules, force_env_section, &
106 prev_molecule_kind_set, &
107 prev_local_molecules)
108
109 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
110 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
111 TYPE(distribution_1d_type), POINTER :: local_particles
112 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
113 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
114 TYPE(distribution_1d_type), POINTER :: local_molecules
115 TYPE(section_vals_type), POINTER :: force_env_section
116 TYPE(molecule_kind_type), DIMENSION(:), OPTIONAL, &
117 POINTER :: prev_molecule_kind_set
118 TYPE(distribution_1d_type), OPTIONAL, POINTER :: prev_local_molecules
119
120 CHARACTER(len=*), PARAMETER :: routinen = 'distribute_molecules_1d'
121
122 INTEGER :: atom_a, bin, handle, iatom, imolecule, imolecule_kind, imolecule_local, &
123 imolecule_prev_kind, iparticle_kind, ipe, iw, kind_a, molecule_a, n, natom, nbins, nload, &
124 nmolecule, nmolecule_kind, nparticle_kind, nsgf, output_unit
125 INTEGER(int_8) :: bin_price
126 INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: workload_count, workload_fill
127 INTEGER, ALLOCATABLE, DIMENSION(:) :: nmolecule_local, nparticle_local, work
128 INTEGER, DIMENSION(:), POINTER :: molecule_list
129 LOGICAL :: found, has_prev_subsys_info, is_local
130 TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:) :: local_molecule
131 TYPE(cp_heap_type) :: bin_heap_count, bin_heap_fill
132 TYPE(cp_logger_type), POINTER :: logger
133 TYPE(molecule_kind_type), POINTER :: molecule_kind
134
135 CALL timeset(routinen, handle)
136
137 has_prev_subsys_info = .false.
138 IF (PRESENT(prev_local_molecules) .AND. &
139 PRESENT(prev_molecule_kind_set)) THEN
140 IF (ASSOCIATED(prev_local_molecules) .AND. &
141 ASSOCIATED(prev_molecule_kind_set)) THEN
142 has_prev_subsys_info = .true.
143 END IF
144 END IF
145
146 logger => cp_get_default_logger()
147
148 associate(group => logger%para_env, mype => logger%para_env%mepos + 1, &
149 npe => logger%para_env%num_pe)
150
151 ALLOCATE (workload_count(npe))
152 workload_count(:) = 0
153
154 ALLOCATE (workload_fill(npe))
155 workload_fill(:) = 0
156
157 nmolecule_kind = SIZE(molecule_kind_set)
158
159 ALLOCATE (nmolecule_local(nmolecule_kind))
160 nmolecule_local(:) = 0
161
162 ALLOCATE (local_molecule(nmolecule_kind))
163
164 nparticle_kind = SIZE(atomic_kind_set)
165
166 ALLOCATE (nparticle_local(nparticle_kind))
167 nparticle_local(:) = 0
168
169 nbins = npe
170
171 CALL cp_heap_new(bin_heap_count, nbins)
172 CALL cp_heap_fill(bin_heap_count, workload_count)
173
174 CALL cp_heap_new(bin_heap_fill, nbins)
175 CALL cp_heap_fill(bin_heap_fill, workload_fill)
176
177 DO imolecule_kind = 1, nmolecule_kind
178
179 molecule_kind => molecule_kind_set(imolecule_kind)
180
181 NULLIFY (molecule_list)
182
183! *** Get the number of molecules and the number of ***
184! *** atoms in each molecule of that molecular kind ***
185
186 CALL get_molecule_kind(molecule_kind=molecule_kind, &
187 molecule_list=molecule_list, &
188 natom=natom, &
189 nsgf=nsgf)
190
191! *** Consider the number of atoms or basis ***
192! *** functions which depends on the method ***
193
194 nload = max(natom, nsgf)
195 nmolecule = SIZE(molecule_list)
196
197! *** Get the number of local molecules of the current molecule kind ***
198
199 DO imolecule = 1, nmolecule
200 IF (has_prev_subsys_info) THEN
201 DO imolecule_prev_kind = 1, SIZE(prev_molecule_kind_set)
202 IF (any(prev_local_molecules%list(imolecule_prev_kind)%array( &
203 1:prev_local_molecules%n_el(imolecule_prev_kind)) == molecule_list(imolecule))) THEN
204 ! molecule used to be local
205 nmolecule_local(imolecule_kind) = nmolecule_local(imolecule_kind) + 1
206 END IF
207 END DO
208 ELSE
209 CALL cp_heap_get_first(bin_heap_count, bin, bin_price, found)
210 IF (.NOT. found) &
211 cpabort("No topmost heap element found.")
212
213 ipe = bin
214 IF (bin_price /= workload_count(ipe)) &
215 cpabort("inconsistent heap")
216
217 workload_count(ipe) = workload_count(ipe) + nload
218 IF (ipe == mype) THEN
219 nmolecule_local(imolecule_kind) = nmolecule_local(imolecule_kind) + 1
220 END IF
221
222 bin_price = workload_count(ipe)
223 CALL cp_heap_reset_first(bin_heap_count, bin_price)
224 END IF
225 END DO
226
227! *** Distribute the molecules ***
228 n = nmolecule_local(imolecule_kind)
229
230 IF (n > 0) THEN
231 ALLOCATE (local_molecule(imolecule_kind)%array(n))
232 ELSE
233 NULLIFY (local_molecule(imolecule_kind)%array)
234 END IF
235
236 imolecule_local = 0
237 DO imolecule = 1, nmolecule
238 is_local = .false.
239 IF (has_prev_subsys_info) THEN
240 DO imolecule_prev_kind = 1, SIZE(prev_molecule_kind_set)
241 IF (any(prev_local_molecules%list(imolecule_prev_kind)%array( &
242 1:prev_local_molecules%n_el(imolecule_prev_kind)) == molecule_list(imolecule))) THEN
243 is_local = .true.
244 END IF
245 END DO
246 ELSE
247 CALL cp_heap_get_first(bin_heap_fill, bin, bin_price, found)
248 IF (.NOT. found) &
249 cpabort("No topmost heap element found.")
250
251 ipe = bin
252 IF (bin_price /= workload_fill(ipe)) &
253 cpabort("inconsistent heap")
254
255 workload_fill(ipe) = workload_fill(ipe) + nload
256 is_local = (ipe == mype)
257 END IF
258 IF (is_local) THEN
259 imolecule_local = imolecule_local + 1
260 molecule_a = molecule_list(imolecule)
261 local_molecule(imolecule_kind)%array(imolecule_local) = molecule_a
262 DO iatom = 1, natom
263 atom_a = molecule_set(molecule_a)%first_atom + iatom - 1
264
265 CALL get_atomic_kind(atomic_kind=particle_set(atom_a)%atomic_kind, &
266 kind_number=kind_a)
267 nparticle_local(kind_a) = nparticle_local(kind_a) + 1
268 END DO
269 END IF
270 IF (.NOT. has_prev_subsys_info) THEN
271 bin_price = workload_fill(ipe)
272 CALL cp_heap_reset_first(bin_heap_fill, bin_price)
273 END IF
274 END DO
275
276 END DO
277
278 IF (any(workload_fill .NE. workload_count)) &
279 cpabort("Inconsistent heaps encountered")
280
281 CALL cp_heap_release(bin_heap_count)
282 CALL cp_heap_release(bin_heap_fill)
283
284! *** Create the local molecule structure ***
285
286 CALL distribution_1d_create(local_molecules, &
287 n_el=nmolecule_local, &
288 para_env=logger%para_env)
289
290! *** Create the local particle structure ***
291
292 CALL distribution_1d_create(local_particles, &
293 n_el=nparticle_local, &
294 para_env=logger%para_env)
295
296! *** Store the generated local molecule and particle distributions ***
297
298 nparticle_local(:) = 0
299
300 DO imolecule_kind = 1, nmolecule_kind
301
302 IF (nmolecule_local(imolecule_kind) == 0) cycle
303
304 local_molecules%list(imolecule_kind)%array(:) = &
305 local_molecule(imolecule_kind)%array(:)
306
307 molecule_kind => molecule_kind_set(imolecule_kind)
308
309 CALL get_molecule_kind(molecule_kind=molecule_kind, &
310 natom=natom)
311
312 DO imolecule = 1, nmolecule_local(imolecule_kind)
313 molecule_a = local_molecule(imolecule_kind)%array(imolecule)
314 DO iatom = 1, natom
315 atom_a = molecule_set(molecule_a)%first_atom + iatom - 1
316 CALL get_atomic_kind(atomic_kind=particle_set(atom_a)%atomic_kind, &
317 kind_number=kind_a)
318 nparticle_local(kind_a) = nparticle_local(kind_a) + 1
319 local_particles%list(kind_a)%array(nparticle_local(kind_a)) = atom_a
320 END DO
321 END DO
322
323 END DO
324
325! *** Print distribution, if requested ***
326
327 IF (btest(cp_print_key_should_output(logger%iter_info, &
328 force_env_section, "PRINT%DISTRIBUTION1D"), cp_p_file)) THEN
329
330 output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%DISTRIBUTION1D", &
331 extension=".Log")
332
333 iw = output_unit
334 IF (output_unit < 0) iw = cp_logger_get_default_unit_nr(logger, local=.true.)
335
336! *** Print molecule distribution ***
337
338 ALLOCATE (work(npe))
339 work(:) = 0
340
341 work(mype) = sum(nmolecule_local)
342 CALL group%sum(work)
343
344 IF (output_unit > 0) THEN
345 WRITE (unit=output_unit, &
346 fmt="(/, T2, A, T51, A, /, (T52, I6, T73, I8))") &
347 "DISTRIBUTION OF THE MOLECULES", &
348 "Process Number of molecules", &
349 (ipe - 1, work(ipe), ipe=1, npe)
350 WRITE (unit=output_unit, fmt="(T55, A3, T73, I8)") &
351 "Sum", sum(work)
352 CALL m_flush(output_unit)
353 END IF
354
355 CALL group%sync()
356
357 DO ipe = 1, npe
358 IF (ipe == mype) THEN
359 WRITE (unit=iw, fmt="(/, T3, A)") &
360 "Process Kind Local molecules (global indices)"
361 DO imolecule_kind = 1, nmolecule_kind
362 IF (imolecule_kind == 1) THEN
363 WRITE (unit=iw, fmt="(T4, I6, 2X, I5, (T21, 10I6))") &
364 ipe - 1, imolecule_kind, &
365 (local_molecules%list(imolecule_kind)%array(imolecule), &
366 imolecule=1, nmolecule_local(imolecule_kind))
367 ELSE
368 WRITE (unit=iw, fmt="(T12, I5, (T21, 10I6))") &
369 imolecule_kind, &
370 (local_molecules%list(imolecule_kind)%array(imolecule), &
371 imolecule=1, nmolecule_local(imolecule_kind))
372 END IF
373 END DO
374 END IF
375 CALL m_flush(iw)
376 CALL group%sync()
377 END DO
378
379! *** Print particle distribution ***
380
381 work(:) = 0
382
383 work(mype) = sum(nparticle_local)
384 CALL group%sum(work)
385
386 IF (output_unit > 0) THEN
387 WRITE (unit=output_unit, &
388 fmt="(/, T2, A, T51, A, /, (T52, I6, T73, I8))") &
389 "DISTRIBUTION OF THE PARTICLES", &
390 "Process Number of particles", &
391 (ipe - 1, work(ipe), ipe=1, npe)
392 WRITE (unit=output_unit, fmt="(T55, A3, T73, I8)") &
393 "Sum", sum(work)
394 CALL m_flush(output_unit)
395 END IF
396
397 CALL group%sync()
398
399 DO ipe = 1, npe
400 IF (ipe == mype) THEN
401 WRITE (unit=iw, fmt="(/, T3, A)") &
402 "Process Kind Local particles (global indices)"
403 DO iparticle_kind = 1, nparticle_kind
404 IF (iparticle_kind == 1) THEN
405 WRITE (unit=iw, fmt="(T4, I6, 2X, I5, (T20, 10I6))") &
406 ipe - 1, iparticle_kind, &
407 (local_particles%list(iparticle_kind)%array(iatom), &
408 iatom=1, nparticle_local(iparticle_kind))
409 ELSE
410 WRITE (unit=iw, fmt="(T12, I5, (T20, 10I6))") &
411 iparticle_kind, &
412 (local_particles%list(iparticle_kind)%array(iatom), &
413 iatom=1, nparticle_local(iparticle_kind))
414 END IF
415 END DO
416 END IF
417 CALL m_flush(iw)
418 CALL group%sync()
419 END DO
420 DEALLOCATE (work)
421
422 CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
423 "PRINT%DISTRIBUTION1D")
424 END IF
425 END associate
426! *** Release work storage ***
427
428 DEALLOCATE (workload_count)
429
430 DEALLOCATE (workload_fill)
431
432 DEALLOCATE (nmolecule_local)
433
434 DEALLOCATE (nparticle_local)
435
436 DO imolecule_kind = 1, nmolecule_kind
437 IF (ASSOCIATED(local_molecule(imolecule_kind)%array)) THEN
438 DEALLOCATE (local_molecule(imolecule_kind)%array)
439 END IF
440 END DO
441 DEALLOCATE (local_molecule)
442
443 CALL timestop(handle)
444
445 END SUBROUTINE distribute_molecules_1d
446
447! **************************************************************************************************
448!> \brief Distributes the particle pairs creating a 2d distribution optimally
449!> suited for quickstep
450!> \param cell ...
451!> \param atomic_kind_set ...
452!> \param particle_set ...
453!> \param qs_kind_set ...
454!> \param molecule_kind_set ...
455!> \param molecule_set ...
456!> \param distribution_2d the distribution that will be created by this
457!> method
458!> \param blacs_env the parallel environment at the basis of the
459!> distribution
460!> \param force_env_section ...
461!> \par History
462!> - local_rows & cols blocksize optimizations (Aug. 2003, MK)
463!> - cleanup of distribution_2d (Sep. 2003, fawzi)
464!> - update for molecules (Oct. 2003, MK)
465!> \author fawzi (Feb. 2003)
466!> \note
467!> Intermediate generation of a 2d distribution of the molecules, but
468!> only the corresponding particle (atomic) distribution is currently
469!> used. The 2d distribution of the molecules is deleted, but may easily
470!> be recovered (MK).
471! **************************************************************************************************
472 SUBROUTINE distribute_molecules_2d(cell, atomic_kind_set, particle_set, &
473 qs_kind_set, molecule_kind_set, molecule_set, &
474 distribution_2d, blacs_env, force_env_section)
475 TYPE(cell_type), POINTER :: cell
476 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
477 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
478 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
479 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
480 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
481 TYPE(distribution_2d_type), POINTER :: distribution_2d
482 TYPE(cp_blacs_env_type), POINTER :: blacs_env
483 TYPE(section_vals_type), POINTER :: force_env_section
484
485 CHARACTER(len=*), PARAMETER :: routinen = 'distribute_molecules_2d'
486
487 INTEGER :: cluster_price, cost_model, handle, iatom, iatom_mol, iatom_one, ikind, imol, &
488 imolecule, imolecule_kind, iparticle_kind, ipcol, iprow, iw, kind_a, n, natom, natom_mol, &
489 nclusters, nmolecule, nmolecule_kind, nparticle_kind, nsgf, output_unit
490 INTEGER, ALLOCATABLE, DIMENSION(:) :: cluster_list, cluster_prices, &
491 nparticle_local_col, &
492 nparticle_local_row, work
493 INTEGER, DIMENSION(:), POINTER :: lmax_basis, molecule_list
494 INTEGER, DIMENSION(:, :), POINTER :: cluster_col_distribution, &
495 cluster_row_distribution, &
496 col_distribution, row_distribution
497 LOGICAL :: basic_cluster_optimization, basic_optimization, basic_spatial_optimization, &
498 molecular_distribution, skip_optimization
499 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coords, pbc_scaled_coords
500 REAL(kind=dp), DIMENSION(3) :: center
501 TYPE(cp_1d_i_p_type), DIMENSION(:), POINTER :: local_particle_col, local_particle_row
502 TYPE(cp_logger_type), POINTER :: logger
503 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
504 TYPE(molecule_kind_type), POINTER :: molecule_kind
505 TYPE(section_vals_type), POINTER :: distribution_section
506
507!...
508
509 CALL timeset(routinen, handle)
510
511 logger => cp_get_default_logger()
512
513 distribution_section => section_vals_get_subs_vals(force_env_section, "DFT%QS%DISTRIBUTION")
514
515 CALL section_vals_val_get(distribution_section, "2D_MOLECULAR_DISTRIBUTION", l_val=molecular_distribution)
516 CALL section_vals_val_get(distribution_section, "SKIP_OPTIMIZATION", l_val=skip_optimization)
517 CALL section_vals_val_get(distribution_section, "BASIC_OPTIMIZATION", l_val=basic_optimization)
518 CALL section_vals_val_get(distribution_section, "BASIC_SPATIAL_OPTIMIZATION", l_val=basic_spatial_optimization)
519 CALL section_vals_val_get(distribution_section, "BASIC_CLUSTER_OPTIMIZATION", l_val=basic_cluster_optimization)
520
521 CALL section_vals_val_get(distribution_section, "COST_MODEL", i_val=cost_model)
522 !
523
524 associate(group => blacs_env%para_env, myprow => blacs_env%mepos(1) + 1, mypcol => blacs_env%mepos(2) + 1, &
525 nprow => blacs_env%num_pe(1), npcol => blacs_env%num_pe(2))
526
527 nmolecule_kind = SIZE(molecule_kind_set)
528 CALL get_molecule_kind_set(molecule_kind_set, nmolecule=nmolecule)
529
530 nparticle_kind = SIZE(atomic_kind_set)
531 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
532
533 !
534 ! we need to generate two representations of the distribution, one as a straight array with global particles
535 ! one ordered wrt to kinds and only listing the local particles
536 !
537 ALLOCATE (row_distribution(natom, 2))
538 ALLOCATE (col_distribution(natom, 2))
539 ! Initialize the distributions to -1, as the second dimension only gets set with cluster optimization
540 ! but the information is needed by dbcsr
541 row_distribution = -1; col_distribution = -1
542
543 ALLOCATE (local_particle_col(nparticle_kind))
544 ALLOCATE (local_particle_row(nparticle_kind))
545 ALLOCATE (nparticle_local_row(nparticle_kind))
546 ALLOCATE (nparticle_local_col(nparticle_kind))
547
548 IF (basic_optimization .OR. basic_spatial_optimization .OR. basic_cluster_optimization) THEN
549
550 IF (molecular_distribution) THEN
551 nclusters = nmolecule
552 ELSE
553 nclusters = natom
554 END IF
555
556 ALLOCATE (cluster_list(nclusters))
557 ALLOCATE (cluster_prices(nclusters))
558 ALLOCATE (cluster_row_distribution(nclusters, 2))
559 ALLOCATE (cluster_col_distribution(nclusters, 2))
560 cluster_row_distribution = -1; cluster_col_distribution = -1
561
562 ! Fill in the clusters and their prices
563 CALL section_vals_val_get(distribution_section, "COST_MODEL", i_val=cost_model)
564 IF (.NOT. molecular_distribution) THEN
565 DO iatom = 1, natom
566 IF (iatom .GT. nclusters) &
567 cpabort("Bounds error")
568 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
569 cluster_list(iatom) = iatom
570 SELECT CASE (cost_model)
571 CASE (model_block_count)
572 CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf)
573 cluster_price = nsgf
574 CASE (model_block_lmax)
575 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
576 CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
577 cluster_price = maxval(lmax_basis)
578 CASE default
579 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
580 CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
581 cluster_price = 8 + (maxval(lmax_basis)**2)
582 END SELECT
583 cluster_prices(iatom) = cluster_price
584 END DO
585 ELSE
586 imol = 0
587 DO imolecule_kind = 1, nmolecule_kind
588 molecule_kind => molecule_kind_set(imolecule_kind)
589 CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, natom=natom_mol)
590 DO imolecule = 1, SIZE(molecule_list)
591 imol = imol + 1
592 cluster_list(imol) = imol
593 cluster_price = 0
594 DO iatom_mol = 1, natom_mol
595 iatom = molecule_set(molecule_list(imolecule))%first_atom + iatom_mol - 1
596 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
597 SELECT CASE (cost_model)
598 CASE (model_block_count)
599 CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf)
600 cluster_price = cluster_price + nsgf
601 CASE (model_block_lmax)
602 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
603 CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
604 cluster_price = cluster_price + maxval(lmax_basis)
605 CASE default
606 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
607 CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
608 cluster_price = cluster_price + 8 + (maxval(lmax_basis)**2)
609 END SELECT
610 END DO
611 cluster_prices(imol) = cluster_price
612 END DO
613 END DO
614 END IF
615
616 ! And distribute
617 IF (basic_optimization) THEN
618 CALL make_basic_distribution(cluster_list, cluster_prices, &
619 nprow, cluster_row_distribution(:, 1), npcol, cluster_col_distribution(:, 1))
620 ELSE
621 IF (basic_cluster_optimization) THEN
622 IF (molecular_distribution) &
623 cpabort("clustering and molecular blocking NYI")
624 ALLOCATE (pbc_scaled_coords(3, natom), coords(3, natom))
625 DO iatom = 1, natom
626 CALL real_to_scaled(pbc_scaled_coords(:, iatom), pbc(particle_set(iatom)%r(:), cell), cell)
627 coords(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
628 END DO
629 CALL make_cluster_distribution(coords, pbc_scaled_coords, cell, cluster_prices, &
630 nprow, cluster_row_distribution, npcol, cluster_col_distribution)
631 ELSE ! basic_spatial_optimization
632 ALLOCATE (pbc_scaled_coords(3, nclusters))
633 IF (.NOT. molecular_distribution) THEN
634 ! just scaled coords
635 DO iatom = 1, natom
636 CALL real_to_scaled(pbc_scaled_coords(:, iatom), pbc(particle_set(iatom)%r(:), cell), cell)
637 END DO
638 ELSE
639 ! use scaled coords of geometric center, folding when appropriate
640 imol = 0
641 DO imolecule_kind = 1, nmolecule_kind
642 molecule_kind => molecule_kind_set(imolecule_kind)
643 CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, natom=natom_mol)
644 DO imolecule = 1, SIZE(molecule_list)
645 imol = imol + 1
646 iatom_one = molecule_set(molecule_list(imolecule))%first_atom
647 center = 0.0_dp
648 DO iatom_mol = 1, natom_mol
649 iatom = molecule_set(molecule_list(imolecule))%first_atom + iatom_mol - 1
650 center = center + &
651 pbc(particle_set(iatom)%r(:) - particle_set(iatom_one)%r(:), cell) + particle_set(iatom_one)%r(:)
652 END DO
653 center = center/natom_mol
654 CALL real_to_scaled(pbc_scaled_coords(:, imol), pbc(center, cell), cell)
655 END DO
656 END DO
657 END IF
658
659 CALL make_basic_spatial_distribution(pbc_scaled_coords, cluster_prices, &
660 nprow, cluster_row_distribution(:, 1), npcol, cluster_col_distribution(:, 1))
661
662 DEALLOCATE (pbc_scaled_coords)
663 END IF
664 END IF
665
666 ! And assign back
667 IF (.NOT. molecular_distribution) THEN
668 row_distribution = cluster_row_distribution
669 col_distribution = cluster_col_distribution
670 ELSE
671 imol = 0
672 DO imolecule_kind = 1, nmolecule_kind
673 molecule_kind => molecule_kind_set(imolecule_kind)
674 CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, natom=natom_mol)
675 DO imolecule = 1, SIZE(molecule_list)
676 imol = imol + 1
677 DO iatom_mol = 1, natom_mol
678 iatom = molecule_set(molecule_list(imolecule))%first_atom + iatom_mol - 1
679 row_distribution(iatom, :) = cluster_row_distribution(imol, :)
680 col_distribution(iatom, :) = cluster_col_distribution(imol, :)
681 END DO
682 END DO
683 END DO
684 END IF
685
686 ! cleanup
687 DEALLOCATE (cluster_list)
688 DEALLOCATE (cluster_prices)
689 DEALLOCATE (cluster_row_distribution)
690 DEALLOCATE (cluster_col_distribution)
691
692 ELSE
693 ! expects nothing else
694 cpabort("")
695 END IF
696
697 ! prepare the lists of local particles
698
699 ! count local particles of a given kind
700 nparticle_local_col = 0
701 nparticle_local_row = 0
702 DO iatom = 1, natom
703 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, kind_number=kind_a)
704 IF (row_distribution(iatom, 1) == myprow) nparticle_local_row(kind_a) = nparticle_local_row(kind_a) + 1
705 IF (col_distribution(iatom, 1) == mypcol) nparticle_local_col(kind_a) = nparticle_local_col(kind_a) + 1
706 END DO
707
708 ! allocate space
709 DO iparticle_kind = 1, nparticle_kind
710 n = nparticle_local_row(iparticle_kind)
711 ALLOCATE (local_particle_row(iparticle_kind)%array(n))
712
713 n = nparticle_local_col(iparticle_kind)
714 ALLOCATE (local_particle_col(iparticle_kind)%array(n))
715 END DO
716
717 ! store
718 nparticle_local_col = 0
719 nparticle_local_row = 0
720 DO iatom = 1, natom
721 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, kind_number=kind_a)
722 IF (row_distribution(iatom, 1) == myprow) THEN
723 nparticle_local_row(kind_a) = nparticle_local_row(kind_a) + 1
724 local_particle_row(kind_a)%array(nparticle_local_row(kind_a)) = iatom
725 END IF
726 IF (col_distribution(iatom, 1) == mypcol) THEN
727 nparticle_local_col(kind_a) = nparticle_local_col(kind_a) + 1
728 local_particle_col(kind_a)%array(nparticle_local_col(kind_a)) = iatom
729 END IF
730 END DO
731
732! *** Generate the 2d distribution structure but take care of the zero offsets required
733 row_distribution(:, 1) = row_distribution(:, 1) - 1
734 col_distribution(:, 1) = col_distribution(:, 1) - 1
735 CALL distribution_2d_create(distribution_2d, &
736 row_distribution_ptr=row_distribution, &
737 col_distribution_ptr=col_distribution, &
738 local_rows_ptr=local_particle_row, &
739 local_cols_ptr=local_particle_col, &
740 blacs_env=blacs_env)
741
742 NULLIFY (local_particle_row)
743 NULLIFY (local_particle_col)
744 NULLIFY (row_distribution)
745 NULLIFY (col_distribution)
746
747! *** Print distribution, if requested ***
748 IF (btest(cp_print_key_should_output(logger%iter_info, &
749 force_env_section, "PRINT%DISTRIBUTION"), cp_p_file)) THEN
750
751 output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%DISTRIBUTION", &
752 extension=".Log")
753
754! *** Print row distribution ***
755
756 ALLOCATE (work(nprow))
757 work(:) = 0
758
759 IF (mypcol == 1) work(myprow) = sum(distribution_2d%n_local_rows)
760
761 CALL group%sum(work)
762
763 IF (output_unit > 0) THEN
764 WRITE (unit=output_unit, &
765 fmt="(/, T2, A, /, T15, A, /, (T16, I10, T41, I10, T71, I10))") &
766 "DISTRIBUTION OF THE PARTICLES (ROWS)", &
767 "Process row Number of particles Number of matrix rows", &
768 (iprow - 1, work(iprow), -1, iprow=1, nprow)
769 WRITE (unit=output_unit, fmt="(T23, A3, T41, I10, T71, I10)") &
770 "Sum", sum(work), -1
771 CALL m_flush(output_unit)
772 END IF
773
774 DEALLOCATE (work)
775
776! *** Print column distribution ***
777
778 ALLOCATE (work(npcol))
779 work(:) = 0
780
781 IF (myprow == 1) work(mypcol) = sum(distribution_2d%n_local_cols)
782
783 CALL group%sum(work)
784
785 IF (output_unit > 0) THEN
786 WRITE (unit=output_unit, &
787 fmt="(/, T2, A, /, T15, A, /, (T16, I10, T41, I10, T71, I10))") &
788 "DISTRIBUTION OF THE PARTICLES (COLUMNS)", &
789 "Process col Number of particles Number of matrix columns", &
790 (ipcol - 1, work(ipcol), -1, ipcol=1, npcol)
791 WRITE (unit=output_unit, fmt="(T23, A3, T41, I10, T71, I10)") &
792 "Sum", sum(work), -1
793 CALL m_flush(output_unit)
794 END IF
795
796 DEALLOCATE (work)
797
798 CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
799 "PRINT%DISTRIBUTION")
800 END IF
801 END associate
802
803 IF (btest(cp_print_key_should_output(logger%iter_info, &
804 force_env_section, "PRINT%DISTRIBUTION2D"), cp_p_file)) THEN
805
806 iw = cp_logger_get_default_unit_nr(logger, local=.true.)
807 CALL distribution_2d_write(distribution_2d, &
808 unit_nr=iw, &
809 local=.true., &
810 long_description=.true.)
811
812 END IF
813
814! *** Release work storage ***
815
816 DEALLOCATE (nparticle_local_row)
817
818 DEALLOCATE (nparticle_local_col)
819
820 CALL timestop(handle)
821
822 END SUBROUTINE distribute_molecules_2d
823
824! **************************************************************************************************
825!> \brief Creates a basic distribution
826!> \param cluster_list ...
827!> \param cluster_prices ...
828!> \param nprows ...
829!> \param row_distribution ...
830!> \param npcols ...
831!> \param col_distribution ...
832!> \par History
833!> - Created 2010-08-06 UB
834! **************************************************************************************************
835 SUBROUTINE make_basic_distribution(cluster_list, cluster_prices, &
836 nprows, row_distribution, npcols, col_distribution)
837 INTEGER, DIMENSION(:), INTENT(INOUT) :: cluster_list, cluster_prices
838 INTEGER, INTENT(IN) :: nprows
839 INTEGER, DIMENSION(:), INTENT(OUT) :: row_distribution
840 INTEGER, INTENT(IN) :: npcols
841 INTEGER, DIMENSION(:), INTENT(OUT) :: col_distribution
842
843 CHARACTER(len=*), PARAMETER :: routinen = 'make_basic_distribution'
844
845 INTEGER :: bin, cluster, cluster_index, &
846 cluster_price, nbins, nclusters, pcol, &
847 pgrid_gcd, prow, timing_handle
848 INTEGER(int_8) :: bin_price
849 LOGICAL :: found
850 TYPE(cp_heap_type) :: bin_heap
851
852! ---------------------------------------------------------------------------
853
854 CALL timeset(routinen, timing_handle)
855 nbins = lcm(nprows, npcols)
856 pgrid_gcd = gcd(nprows, npcols)
857 CALL sort(cluster_prices, SIZE(cluster_list), cluster_list)
858 CALL cp_heap_new(bin_heap, nbins)
859 CALL cp_heap_fill(bin_heap, (/(0_int_8, bin=1, nbins)/))
860 !
861 nclusters = SIZE(cluster_list)
862 ! Put the most expensive cluster in the bin with the smallest
863 ! price and repeat.
864 DO cluster_index = nclusters, 1, -1
865 cluster = cluster_list(cluster_index)
866 CALL cp_heap_get_first(bin_heap, bin, bin_price, found)
867 IF (.NOT. found) &
868 cpabort("No topmost heap element found.")
869 !
870 prow = int((bin - 1)*pgrid_gcd/npcols)
871 IF (prow .GE. nprows) &
872 cpabort("Invalid process row.")
873 pcol = int((bin - 1)*pgrid_gcd/nprows)
874 IF (pcol .GE. npcols) &
875 cpabort("Invalid process column.")
876 row_distribution(cluster) = prow + 1
877 col_distribution(cluster) = pcol + 1
878 !
879 cluster_price = cluster_prices(cluster_index)
880 bin_price = bin_price + cluster_price
881 CALL cp_heap_reset_first(bin_heap, bin_price)
882 END DO
883 CALL cp_heap_release(bin_heap)
884 CALL timestop(timing_handle)
885 END SUBROUTINE make_basic_distribution
886
887! **************************************************************************************************
888!> \brief Creates a basic spatial distribution
889!> that tries to make the corresponding blocks as homogeneous as possible
890!> \param pbc_scaled_coords ...
891!> \param costs ...
892!> \param nprows ...
893!> \param row_distribution ...
894!> \param npcols ...
895!> \param col_distribution ...
896!> \par History
897!> - Created 2010-11-11 Joost VandeVondele
898! **************************************************************************************************
899 SUBROUTINE make_basic_spatial_distribution(pbc_scaled_coords, costs, &
900 nprows, row_distribution, npcols, col_distribution)
901 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pbc_scaled_coords
902 INTEGER, DIMENSION(:), INTENT(IN) :: costs
903 INTEGER, INTENT(IN) :: nprows
904 INTEGER, DIMENSION(:), INTENT(OUT) :: row_distribution
905 INTEGER, INTENT(IN) :: npcols
906 INTEGER, DIMENSION(:), INTENT(OUT) :: col_distribution
907
908 CHARACTER(len=*), PARAMETER :: routinen = 'make_basic_spatial_distribution'
909
910 INTEGER :: handle, iatom, natoms, nbins, pgrid_gcd
911 INTEGER, ALLOCATABLE, DIMENSION(:) :: bin_costs, distribution
912
913 CALL timeset(routinen, handle)
914
915 natoms = SIZE(costs)
916 nbins = lcm(nprows, npcols)
917 pgrid_gcd = gcd(nprows, npcols)
918 ALLOCATE (bin_costs(nbins), distribution(natoms))
919 bin_costs = 0
920
921 CALL spatial_recurse(pbc_scaled_coords, costs, (/(iatom, iatom=1, natoms)/), bin_costs, distribution, 0)
922
923 ! WRITE(*, *) "Final bin costs: ", bin_costs
924
925 ! final row_distribution / col_distribution
926 DO iatom = 1, natoms
927 row_distribution(iatom) = (distribution(iatom) - 1)*pgrid_gcd/npcols + 1
928 col_distribution(iatom) = (distribution(iatom) - 1)*pgrid_gcd/nprows + 1
929 END DO
930
931 DEALLOCATE (bin_costs, distribution)
932
933 CALL timestop(handle)
934
935 END SUBROUTINE make_basic_spatial_distribution
936
937! **************************************************************************************************
938!> \brief ...
939!> \param pbc_scaled_coords ...
940!> \param costs ...
941!> \param indices ...
942!> \param bin_costs ...
943!> \param distribution ...
944!> \param level ...
945! **************************************************************************************************
946 RECURSIVE SUBROUTINE spatial_recurse(pbc_scaled_coords, costs, indices, bin_costs, distribution, level)
947 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pbc_scaled_coords
948 INTEGER, DIMENSION(:), INTENT(IN) :: costs, indices
949 INTEGER, DIMENSION(:), INTENT(INOUT) :: bin_costs, distribution
950 INTEGER, INTENT(IN) :: level
951
952 INTEGER :: iatom, ibin, natoms, nbins, nhalf
953 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_costs_sorted, atom_permutation, &
954 bin_costs_sorted, permutation
955 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coord
956
957 natoms = SIZE(costs)
958 nbins = SIZE(bin_costs)
959 nhalf = (natoms + 1)/2
960
961 IF (natoms <= nbins) THEN
962 ! assign the most expensive atom to the least costly bin
963 ALLOCATE (bin_costs_sorted(nbins), permutation(nbins))
964 bin_costs_sorted(:) = bin_costs
965 CALL sort(bin_costs_sorted, nbins, permutation)
966 ALLOCATE (atom_costs_sorted(natoms), atom_permutation(natoms))
967 atom_costs_sorted(:) = costs
968 CALL sort(atom_costs_sorted, natoms, atom_permutation)
969 ibin = 0
970 ! WRITE(*, *) "Dealing with a new bunch of atoms "
971 DO iatom = natoms, 1, -1
972 ibin = ibin + 1
973 ! WRITE(*, *) "atom", indices(atom_permutation(iatom)), "cost", atom_costs_sorted(iatom), &
974 ! "bin", permutation(ibin), "its cost", bin_costs(permutation(ibin))
975 ! WRITE(100, '(A, I0, 3F12.6)') "A", permutation(ibin), pbc_scaled_coords(:, atom_permutation(iatom))
976 bin_costs(permutation(ibin)) = bin_costs(permutation(ibin)) + atom_costs_sorted(iatom)
977 distribution(indices(atom_permutation(iatom))) = permutation(ibin)
978 END DO
979 DEALLOCATE (bin_costs_sorted, permutation, atom_costs_sorted, atom_permutation)
980 ELSE
981 ! divide atoms in two subsets, sorting according to their coordinates, alternatively x, y, z
982 ! recursively do this for both subsets
983 ALLOCATE (coord(natoms), permutation(natoms))
984 coord(:) = pbc_scaled_coords(mod(level, 3) + 1, :)
985 CALL sort(coord, natoms, permutation)
986 CALL spatial_recurse(pbc_scaled_coords(:, permutation(1:nhalf)), costs(permutation(1:nhalf)), &
987 indices(permutation(1:nhalf)), bin_costs, distribution, level + 1)
988 CALL spatial_recurse(pbc_scaled_coords(:, permutation(nhalf + 1:)), costs(permutation(nhalf + 1:)), &
989 indices(permutation(nhalf + 1:)), bin_costs, distribution, level + 1)
990 DEALLOCATE (coord, permutation)
991 END IF
992
993 END SUBROUTINE spatial_recurse
994
995! **************************************************************************************************
996!> \brief creates a distribution placing close by atoms into clusters and
997!> putting them on the same processors. Load balancing is
998!> performed by balancing sum of the cluster costs per processor
999!> \param coords coordinates of the system
1000!> \param scaled_coords scaled coordinates
1001!> \param cell the cell_type
1002!> \param costs costs per atomic block
1003!> \param nprows number of precessors per row on the 2d grid
1004!> \param row_distribution the resulting distribution over proc_rows of atomic blocks
1005!> \param npcols number of precessors per col on the 2d grid
1006!> \param col_distribution the resulting distribution over proc_cols of atomic blocks
1007! **************************************************************************************************
1008 SUBROUTINE make_cluster_distribution(coords, scaled_coords, cell, costs, &
1009 nprows, row_distribution, npcols, col_distribution)
1010 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: coords, scaled_coords
1011 TYPE(cell_type), POINTER :: cell
1012 INTEGER, DIMENSION(:), INTENT(IN) :: costs
1013 INTEGER, INTENT(IN) :: nprows
1014 INTEGER, DIMENSION(:, :), INTENT(OUT) :: row_distribution
1015 INTEGER, INTENT(IN) :: npcols
1016 INTEGER, DIMENSION(:, :), INTENT(OUT) :: col_distribution
1017
1018 CHARACTER(len=*), PARAMETER :: routinen = 'make_cluster_distribution'
1019
1020 INTEGER :: handle, i, icluster, level, natom, &
1021 output_unit
1022 INTEGER(KIND=int_8) :: ncluster
1023 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_to_cluster, cluster_cost, &
1024 cluster_count, cluster_to_col, &
1025 cluster_to_row, piv_cost, proc_cost, &
1026 sorted_cost
1027 REAL(kind=dp) :: fold(3)
1028 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cluster_center, cluster_high, cluster_low
1029
1030 CALL timeset(routinen, handle)
1031
1032 output_unit = cp_logger_get_default_io_unit()
1033
1034 natom = SIZE(costs)
1035 ncluster = dbcsr_distribution_get_num_images(sum(costs), natom, nprows, npcols)
1036 ALLOCATE (atom_to_cluster(natom))
1037 ALLOCATE (cluster_cost(ncluster))
1038 ALLOCATE (cluster_to_row(ncluster))
1039 ALLOCATE (cluster_to_col(ncluster))
1040 ALLOCATE (sorted_cost(ncluster))
1041 ALLOCATE (piv_cost(ncluster))
1042 cluster_cost(:) = 0
1043
1044 icluster = 0
1045 CALL cluster_recurse(coords, scaled_coords, cell, costs, atom_to_cluster, ncluster, icluster, cluster_cost)
1046
1047 sorted_cost(:) = cluster_cost(:)
1048 CALL sort(sorted_cost, int(ncluster), piv_cost)
1049
1050 ALLOCATE (proc_cost(nprows))
1051 proc_cost = 0; level = 1
1052 CALL assign_clusters(cluster_cost, piv_cost, proc_cost, cluster_to_row, nprows)
1053
1054 DEALLOCATE (proc_cost); ALLOCATE (proc_cost(npcols))
1055 proc_cost = 0; level = 1
1056 CALL assign_clusters(cluster_cost, piv_cost, proc_cost, cluster_to_col, npcols)
1057
1058 DO i = 1, natom
1059 row_distribution(i, 1) = cluster_to_row(atom_to_cluster(i))
1060 row_distribution(i, 2) = atom_to_cluster(i)
1061 col_distribution(i, 1) = cluster_to_col(atom_to_cluster(i))
1062 col_distribution(i, 2) = atom_to_cluster(i)
1063 END DO
1064
1065 ! generate some statistics on clusters
1066 ALLOCATE (cluster_center(3, ncluster))
1067 ALLOCATE (cluster_low(3, ncluster))
1068 ALLOCATE (cluster_high(3, ncluster))
1069 ALLOCATE (cluster_count(ncluster))
1070 cluster_count = 0
1071 DO i = 1, natom
1072 cluster_count(atom_to_cluster(i)) = cluster_count(atom_to_cluster(i)) + 1
1073 cluster_center(:, atom_to_cluster(i)) = coords(:, i)
1074 END DO
1075 cluster_low = huge(0.0_dp)/2
1076 cluster_high = -huge(0.0_dp)/2
1077 DO i = 1, natom
1078 fold = pbc(coords(:, i) - cluster_center(:, atom_to_cluster(i)), cell) + cluster_center(:, atom_to_cluster(i))
1079 cluster_low(:, atom_to_cluster(i)) = min(cluster_low(:, atom_to_cluster(i)), fold(:))
1080 cluster_high(:, atom_to_cluster(i)) = max(cluster_high(:, atom_to_cluster(i)), fold(:))
1081 END DO
1082 IF (output_unit > 0) THEN
1083 WRITE (output_unit, *)
1084 WRITE (output_unit, '(T2,A)') "Cluster distribution information"
1085 WRITE (output_unit, '(T2,A,T48,I8)') "Number of atoms", natom
1086 WRITE (output_unit, '(T2,A,T48,I8)') "Number of clusters", ncluster
1087 WRITE (output_unit, '(T2,A,T48,I8)') "Largest cluster in atoms", maxval(cluster_count)
1088 WRITE (output_unit, '(T2,A,T48,I8)') "Smallest cluster in atoms", minval(cluster_count)
1089 WRITE (output_unit, '(T2,A,T48,F8.3,I8)') "Largest cartesian extend [a.u.]/cluster x=", &
1090 maxval(cluster_high(1, :) - cluster_low(1, :), mask=(cluster_count > 0)), &
1091 maxloc(cluster_high(1, :) - cluster_low(1, :), mask=(cluster_count > 0))
1092 WRITE (output_unit, '(T2,A,T48,F8.3,I8)') "Largest cartesian extend [a.u.]/cluster y=", &
1093 maxval(cluster_high(2, :) - cluster_low(2, :), mask=(cluster_count > 0)), &
1094 maxloc(cluster_high(2, :) - cluster_low(2, :), mask=(cluster_count > 0))
1095 WRITE (output_unit, '(T2,A,T48,F8.3,I8)') "Largest cartesian extend [a.u.]/cluster z=", &
1096 maxval(cluster_high(3, :) - cluster_low(3, :), mask=(cluster_count > 0)), &
1097 maxloc(cluster_high(3, :) - cluster_low(3, :), mask=(cluster_count > 0))
1098 END IF
1099
1100 DEALLOCATE (atom_to_cluster, cluster_cost, cluster_to_row, cluster_to_col, sorted_cost, piv_cost, proc_cost)
1101 CALL timestop(handle)
1102
1103 END SUBROUTINE make_cluster_distribution
1104
1105! **************************************************************************************************
1106!> \brief assigns the clusters to processors, tryimg to balance the cost on the nodes
1107!> \param cluster_cost vector with the cost of each cluster
1108!> \param piv_cost pivoting vector sorting the cluster_cost
1109!> \param proc_cost cost per processor, on input 0 everywhere
1110!> \param cluster_assign assgnment of clusters on proc
1111!> \param nproc number of processor over which clusters are distributed
1112! **************************************************************************************************
1113 SUBROUTINE assign_clusters(cluster_cost, piv_cost, proc_cost, cluster_assign, nproc)
1114 INTEGER, ALLOCATABLE, DIMENSION(:) :: cluster_cost, piv_cost, proc_cost, &
1115 cluster_assign
1116 INTEGER :: nproc
1117
1118 CHARACTER(len=*), PARAMETER :: routinen = 'assign_clusters'
1119
1120 INTEGER :: handle, i, ilevel, offset, &
1121 piv_pcost(nproc), sort_proc_cost(nproc)
1122
1123 CALL timeset(routinen, handle)
1124
1125 DO ilevel = 1, SIZE(cluster_cost)/nproc
1126 sort_proc_cost(:) = proc_cost(:)
1127 CALL sort(sort_proc_cost, nproc, piv_pcost)
1128
1129 offset = (SIZE(cluster_cost)/nproc - ilevel + 1)*nproc + 1
1130 DO i = 1, nproc
1131 cluster_assign(piv_cost(offset - i)) = piv_pcost(i)
1132 proc_cost(piv_pcost(i)) = proc_cost(piv_pcost(i)) + cluster_cost(piv_cost(offset - i))
1133 END DO
1134 END DO
1135
1136 CALL timestop(handle)
1137
1138 END SUBROUTINE assign_clusters
1139
1140! **************************************************************************************************
1141!> \brief recursive routine to cluster atoms.
1142!> Low level uses a modified KMEANS algorithm
1143!> recursion is used to reduce cost.
1144!> each level will subdivide a cluster into smaller clusters
1145!> If only a single split is necessary atoms are assigned to the current cluster
1146!> \param coord coordinates of the system
1147!> \param scaled_coord scaled coordinates
1148!> \param cell the cell_type
1149!> \param costs costs per atomic block
1150!> \param cluster_inds the atom_to cluster mapping
1151!> \param ncluster number of clusters still to be created on a given recursion level
1152!> \param icluster the index of the current cluster to be created
1153!> \param fin_cluster_cost total cost of the final clusters
1154! **************************************************************************************************
1155 RECURSIVE SUBROUTINE cluster_recurse(coord, scaled_coord, cell, costs, cluster_inds, ncluster, icluster, fin_cluster_cost)
1156 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: coord, scaled_coord
1157 TYPE(cell_type), POINTER :: cell
1158 INTEGER, DIMENSION(:), INTENT(IN) :: costs
1159 INTEGER, DIMENSION(:), INTENT(INOUT) :: cluster_inds
1160 INTEGER(KIND=int_8), INTENT(INOUT) :: ncluster
1161 INTEGER, INTENT(INOUT) :: icluster
1162 INTEGER, DIMENSION(:), INTENT(INOUT) :: fin_cluster_cost
1163
1164 INTEGER :: i, ibeg, iend, maxv(1), min_seed, &
1165 natoms, nleft, nsplits, seed, tot_cost
1166 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: ncluster_new
1167 INTEGER, ALLOCATABLE, DIMENSION(:) :: cluster_cost, inds_tmp, nat_cluster, piv
1168 LOGICAL :: found
1169 REAL(kind=dp) :: balance, balance_new, conv
1170
1171 natoms = SIZE(coord, 2)
1172 ! This is a bit of an arbitrary choice, simply a try to avoid too many clusters on large systems and too few for balancing on
1173 ! small systems or subclusters
1174 IF (natoms .LE. 1) THEN
1175 nsplits = 1
1176 ELSE
1177 nsplits = min(int(min(int(max(6, int(60.00/log(real(natoms, kind=dp)))), kind=int_8), ncluster)), natoms)
1178 END IF
1179 IF (nsplits == 1) THEN
1180 icluster = icluster + 1
1181 cluster_inds = icluster
1182 fin_cluster_cost(icluster) = sum(costs)
1183 ELSE
1184 ALLOCATE (cluster_cost(nsplits), ncluster_new(nsplits), inds_tmp(natoms), piv(natoms), nat_cluster(nsplits))
1185 ! initialise some values
1186 cluster_cost = 0; seed = 300; found = .true.; min_seed = seed
1187 CALL kmeans(nsplits, coord, scaled_coord, cell, cluster_inds, nat_cluster, seed, conv)
1188 balance = maxval(real(nat_cluster, kind=dp))/minval(real(nat_cluster, kind=dp))
1189
1190 ! If the system is small enough try to do better in terms of balancing number of atoms per cluster
1191 ! by changing the seed for the initial guess
1192 IF (natoms .LT. 1000 .AND. balance .GT. 1.1) THEN
1193 found = .false.
1194 DO i = 1, 5
1195 IF (balance .GT. 1.1) THEN
1196 CALL kmeans(nsplits, coord, scaled_coord, cell, cluster_inds, nat_cluster, seed + i*40, conv)
1197 balance_new = maxval(real(nat_cluster, kind=dp))/minval(real(nat_cluster, kind=dp))
1198 IF (balance_new .LT. balance) THEN
1199 balance = balance_new
1200 min_seed = seed + i*40;
1201 END IF
1202 ELSE
1203 found = .true.
1204 EXIT
1205 END IF
1206 END DO
1207 END IF
1208 !If we do not match the convergence than recompute at least the best assignment
1209 IF (.NOT. found) CALL kmeans(nsplits, coord, scaled_coord, cell, cluster_inds, nat_cluster, min_seed, conv)
1210
1211 ! compute the cost of each cluster to decide how many splits have to be performed on the next lower level
1212 DO i = 1, natoms
1213 cluster_cost(cluster_inds(i)) = cluster_cost(cluster_inds(i)) + costs(i)
1214 END DO
1215 tot_cost = sum(cluster_cost)
1216 ! compute new splitting, can be done more elegant
1217 ncluster_new(:) = ncluster*cluster_cost(:)/tot_cost
1218 nleft = int(ncluster - sum(ncluster_new))
1219 ! As we won't have empty clusters, we can not have 0 as new size, so we correct for this at first
1220 DO i = 1, nsplits
1221 IF (ncluster_new(i) == 0) THEN
1222 ncluster_new(i) = 1
1223 nleft = nleft - 1
1224 END IF
1225 END DO
1226 ! now comes the next part that the number of clusters will not match anymore, so try to correct in a meaningful way without
1227 ! introducing 0 sized blocks again
1228 IF (nleft .NE. 0) THEN
1229 DO i = 1, abs(nleft)
1230 IF (nleft < 0) THEN
1231 maxv = minloc(cluster_cost/ncluster_new)
1232 IF (ncluster_new(maxv(1)) .NE. 1) THEN
1233 ncluster_new(maxv) = ncluster_new(maxv) - 1
1234 ELSE
1235 maxv = maxloc(ncluster_new)
1236 ncluster_new(maxv) = ncluster_new(maxv) - 1
1237 END IF
1238 ELSE
1239 maxv = maxloc(cluster_cost/ncluster_new)
1240 ncluster_new(maxv) = ncluster_new(maxv) + 1
1241 END IF
1242 END DO
1243 END IF
1244
1245 !Now get the permutations to sort the atoms in the nsplits clusters for the next level of iteration
1246 inds_tmp(:) = cluster_inds(:)
1247 CALL sort(inds_tmp, natoms, piv)
1248
1249 ibeg = 1; iend = 0
1250 DO i = 1, nsplits
1251 IF (nat_cluster(i) == 0) cycle
1252 iend = iend + nat_cluster(i)
1253 CALL cluster_recurse(coord(:, piv(ibeg:iend)), scaled_coord(:, piv(ibeg:iend)), cell, costs(piv(ibeg:iend)), &
1254 inds_tmp(ibeg:iend), ncluster_new(i), icluster, fin_cluster_cost)
1255 ibeg = ibeg + nat_cluster(i)
1256 END DO
1257 ! copy the sorted cluster IDs on the old layout, inds_tmp gets set at the lowest level of recursion
1258 cluster_inds(piv(:)) = inds_tmp
1259 DEALLOCATE (cluster_cost, ncluster_new, inds_tmp, piv, nat_cluster)
1260
1261 END IF
1262
1263 END SUBROUTINE cluster_recurse
1264
1265! **************************************************************************************************
1266!> \brief A modified version of the kmeans algorithm.
1267!> The assignment has a penalty function in case clusters become
1268!> larger than average. Like this more even sized clusters are created
1269!> trading it for locality
1270!> \param ncent number of centers to be created
1271!> \param coord coordinates
1272!> \param scaled_coord scaled coord
1273!> \param cell the cell_type
1274!> \param cluster atom to cluster assignment
1275!> \param nat_cl atoms per cluster
1276!> \param seed seed for the RNG. Algorithm might need multiple tries to deliver best results
1277!> \param tot_var the total variance of the clusters around the centers
1278! **************************************************************************************************
1279 SUBROUTINE kmeans(ncent, coord, scaled_coord, cell, cluster, nat_cl, seed, tot_var)
1280 INTEGER :: ncent
1281 REAL(kind=dp), DIMENSION(:, :) :: coord, scaled_coord
1282 TYPE(cell_type), POINTER :: cell
1283 INTEGER, DIMENSION(:) :: cluster, nat_cl
1284 INTEGER :: seed
1285 REAL(kind=dp) :: tot_var
1286
1287 CHARACTER(len=*), PARAMETER :: routinen = 'kmeans'
1288
1289 INTEGER :: handle, i, ind, itn, j, nat, oldc
1290 LOGICAL :: changed
1291 REAL(kind=dp) :: average(3, ncent, 2), cent_coord(3, ncent), devi, deviat(ncent), dist, &
1292 dvec(3), old_var, rn, scaled_cent(3, ncent), var_cl(ncent)
1293 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dmat
1294 REAL(kind=dp), DIMENSION(3, 2) :: initial_seed
1295 TYPE(rng_stream_type) :: rng_stream
1296
1297 CALL timeset(routinen, handle)
1298
1299 initial_seed = real(seed, dp); nat = SIZE(coord, 2)
1300 ALLOCATE (dmat(ncent, nat))
1301
1302 rng_stream = rng_stream_type(name="kmeans uniform distribution [0,1]", &
1303 distribution_type=uniform, seed=initial_seed)
1304
1305! try to find a clever initial guess with centers being somewhat distributed
1306 rn = rng_stream%next()
1307 ind = ceiling(rn*nat)
1308 cent_coord(:, 1) = coord(:, ind)
1309 DO i = 2, ncent
1310 DO
1311 rn = rng_stream%next()
1312 ind = ceiling(rn*nat)
1313 cent_coord(:, i) = coord(:, ind)
1314 devi = huge(1.0_dp)
1315 DO j = 1, i - 1
1316 dvec = pbc(cent_coord(:, j), cent_coord(:, i), cell)
1317 dist = sqrt(dot_product(dvec, dvec))
1318 IF (dist .LT. devi) devi = dist
1319 END DO
1320 rn = rng_stream%next()
1321 IF (rn .LT. devi**2/169.0) EXIT
1322 END DO
1323 END DO
1324
1325! Now start the KMEANS but penalise it in case it starts packing too many atoms into a single set
1326! Unfoirtunatelz as this is dependent on what happened before it cant be parallel
1327 cluster = 0; old_var = huge(1.0_dp)
1328 DO itn = 1, 1000
1329 changed = .false.; var_cl = 0.0_dp; tot_var = 0.0_dp; nat_cl = 0; deviat = 0.0_dp
1330! !$OMP PARALLEL DO PRIVATE(i,j,dvec)
1331 DO i = 1, nat
1332 DO j = 1, ncent
1333 dvec = pbc(cent_coord(:, j), coord(:, i), cell)
1334 dmat(j, i) = dot_product(dvec, dvec)
1335 END DO
1336 END DO
1337 DO i = 1, nat
1338 devi = huge(1.0_dp); oldc = cluster(i)
1339 DO j = 1, ncent
1340 dist = dmat(j, i) + max(nat_cl(j)**2/nat*ncent, nat/ncent)
1341 IF (dist .LT. devi) THEN
1342 devi = dist; cluster(i) = j
1343 END IF
1344 END DO
1345 deviat(cluster(i)) = deviat(cluster(i)) + sqrt(devi)
1346 nat_cl(cluster(i)) = nat_cl(cluster(i)) + 1
1347 tot_var = tot_var + devi
1348 IF (oldc .NE. cluster(i)) changed = .true.
1349 END DO
1350 ! get the update of the centers done, add a new one in case one center lost all its atoms
1351 ! the algorithm would survive, but its nice to really create what you demand
1352 IF (tot_var .GE. old_var) EXIT
1353 IF (changed) THEN
1354 ! Here misery of computing the center of geometry of the clusters in PBC.
1355 ! The mapping on the unit circle allows to circumvent all problems
1356 average = 0.0_dp
1357 DO i = 1, SIZE(coord, 2)
1358 average(:, cluster(i), 1) = average(:, cluster(i), 1) + cos(scaled_coord(:, i)*2.0_dp*pi)
1359 average(:, cluster(i), 2) = average(:, cluster(i), 2) + sin(scaled_coord(:, i)*2.0_dp*pi)
1360 END DO
1361
1362 DO i = 1, ncent
1363 IF (nat_cl(i) == 0) THEN
1364 rn = rng_stream%next()
1365 scaled_cent(:, i) = scaled_coord(:, ceiling(rn*nat))
1366 ELSE
1367 average(:, i, 1) = average(:, i, 1)/real(nat_cl(i), dp)
1368 average(:, i, 2) = average(:, i, 2)/real(nat_cl(i), dp)
1369 scaled_cent(:, i) = (atan2(-average(:, i, 2), -average(:, i, 1)) + pi)/(2.0_dp*pi)
1370 CALL scaled_to_real(cent_coord(:, i), scaled_cent(:, i), cell)
1371 END IF
1372 END DO
1373 ELSE
1374 EXIT
1375 END IF
1376 END DO
1377
1378 CALL timestop(handle)
1379
1380 END SUBROUTINE kmeans
1381
1382END MODULE distribution_methods
static int gcd(const int a, const int b)
Private routine for computing greatest common divisor of two numbers.
static int lcm(const int a, const int b)
Private routine for computing least common multiple of two numbers.
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:516
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:486
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
methods related to the blacs parallel environment
various routines to log and control the output. The idea is that decisions about where to log should ...
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
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...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
subroutine, public cp_heap_fill(heap, values)
Fill heap with given values.
subroutine, public cp_heap_new(heap, n)
...
subroutine, public cp_heap_get_first(heap, key, value, found)
Returns the first heap element without removing it.
subroutine, public cp_heap_release(heap)
...
subroutine, public cp_heap_reset_first(heap, value)
Changes the value of the minimum heap element and rebalances the heap.
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
subroutine, public distribution_1d_create(distribution_1d, para_env, listbased_distribution, n_el, n_lists)
creates a local list
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
subroutine, public distribution_2d_create(distribution_2d, blacs_env, local_rows_ptr, n_local_rows, local_cols_ptr, row_distribution_ptr, col_distribution_ptr, n_local_cols, n_row_distribution, n_col_distribution)
initializes the distribution_2d
subroutine, public distribution_2d_write(distribution_2d, unit_nr, local, long_description)
writes out the given distribution
Distribution methods for atoms, particles, or molecules.
subroutine, public distribute_molecules_1d(atomic_kind_set, particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, force_env_section, prev_molecule_kind_set, prev_local_molecules)
Distribute molecules and particles.
subroutine, public distribute_molecules_2d(cell, atomic_kind_set, particle_set, qs_kind_set, molecule_kind_set, molecule_set, distribution_2d, blacs_env, force_env_section)
Distributes the particle pairs creating a 2d distribution optimally suited for quickstep.
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public model_block_count
integer, parameter, public model_block_lmax
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:106
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
elemental integer function, public lcm(a, b)
computes the least common multiplier of two numbers
Definition mathlib.F:1326
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1291
Define the molecule kind structure types and the corresponding functionality.
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
subroutine, public get_molecule_kind_set(molecule_kind_set, maxatom, natom, nbond, nbend, nub, ntorsion, nimpr, nopbend, nconstraint, nconstraint_fixd, nmolecule, nrestraints)
Get informations about a molecule kind set.
Define the data structure for the molecule information.
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
integer, parameter, public uniform
Define the data structure for the particle information.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
All kind of helpful little routines.
Definition util.F:14
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
represent a pointer to a 1d array
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
type of a logger, at the moment it contains just a print level starting at which level it should be l...
structure to store local (to a processor) ordered lists of integers.
distributes pairs on a 2d grid of processors
Provides all information about a quickstep kind.