(git:b2ae3e3)
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-2026 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,&
27 USE cp_dbcsr_api, ONLY: dbcsr_distribution_get_num_images
32 USE cp_min_heap, ONLY: cp_heap_fill,&
38 USE cp_output_handling, ONLY: cp_p_file,&
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) THEN
211 cpabort("No topmost heap element found.")
212 END IF
213
214 ipe = bin
215 IF (bin_price /= workload_count(ipe)) THEN
216 cpabort("inconsistent heap")
217 END IF
218
219 workload_count(ipe) = workload_count(ipe) + nload
220 IF (ipe == mype) THEN
221 nmolecule_local(imolecule_kind) = nmolecule_local(imolecule_kind) + 1
222 END IF
223
224 bin_price = workload_count(ipe)
225 CALL cp_heap_reset_first(bin_heap_count, bin_price)
226 END IF
227 END DO
228
229! *** Distribute the molecules ***
230 n = nmolecule_local(imolecule_kind)
231
232 IF (n > 0) THEN
233 ALLOCATE (local_molecule(imolecule_kind)%array(n))
234 ELSE
235 NULLIFY (local_molecule(imolecule_kind)%array)
236 END IF
237
238 imolecule_local = 0
239 DO imolecule = 1, nmolecule
240 is_local = .false.
241 IF (has_prev_subsys_info) THEN
242 DO imolecule_prev_kind = 1, SIZE(prev_molecule_kind_set)
243 IF (any(prev_local_molecules%list(imolecule_prev_kind)%array( &
244 1:prev_local_molecules%n_el(imolecule_prev_kind)) == molecule_list(imolecule))) THEN
245 is_local = .true.
246 END IF
247 END DO
248 ELSE
249 CALL cp_heap_get_first(bin_heap_fill, bin, bin_price, found)
250 IF (.NOT. found) THEN
251 cpabort("No topmost heap element found.")
252 END IF
253
254 ipe = bin
255 IF (bin_price /= workload_fill(ipe)) THEN
256 cpabort("inconsistent heap")
257 END IF
258
259 workload_fill(ipe) = workload_fill(ipe) + nload
260 is_local = (ipe == mype)
261 END IF
262 IF (is_local) THEN
263 imolecule_local = imolecule_local + 1
264 molecule_a = molecule_list(imolecule)
265 local_molecule(imolecule_kind)%array(imolecule_local) = molecule_a
266 DO iatom = 1, natom
267 atom_a = molecule_set(molecule_a)%first_atom + iatom - 1
268
269 CALL get_atomic_kind(atomic_kind=particle_set(atom_a)%atomic_kind, &
270 kind_number=kind_a)
271 nparticle_local(kind_a) = nparticle_local(kind_a) + 1
272 END DO
273 END IF
274 IF (.NOT. has_prev_subsys_info) THEN
275 bin_price = workload_fill(ipe)
276 CALL cp_heap_reset_first(bin_heap_fill, bin_price)
277 END IF
278 END DO
279
280 END DO
281
282 IF (any(workload_fill /= workload_count)) THEN
283 cpabort("Inconsistent heaps encountered")
284 END IF
285
286 CALL cp_heap_release(bin_heap_count)
287 CALL cp_heap_release(bin_heap_fill)
288
289! *** Create the local molecule structure ***
290
291 CALL distribution_1d_create(local_molecules, &
292 n_el=nmolecule_local, &
293 para_env=logger%para_env)
294
295! *** Create the local particle structure ***
296
297 CALL distribution_1d_create(local_particles, &
298 n_el=nparticle_local, &
299 para_env=logger%para_env)
300
301! *** Store the generated local molecule and particle distributions ***
302
303 nparticle_local(:) = 0
304
305 DO imolecule_kind = 1, nmolecule_kind
306
307 IF (nmolecule_local(imolecule_kind) == 0) cycle
308
309 local_molecules%list(imolecule_kind)%array(:) = &
310 local_molecule(imolecule_kind)%array(:)
311
312 molecule_kind => molecule_kind_set(imolecule_kind)
313
314 CALL get_molecule_kind(molecule_kind=molecule_kind, &
315 natom=natom)
316
317 DO imolecule = 1, nmolecule_local(imolecule_kind)
318 molecule_a = local_molecule(imolecule_kind)%array(imolecule)
319 DO iatom = 1, natom
320 atom_a = molecule_set(molecule_a)%first_atom + iatom - 1
321 CALL get_atomic_kind(atomic_kind=particle_set(atom_a)%atomic_kind, &
322 kind_number=kind_a)
323 nparticle_local(kind_a) = nparticle_local(kind_a) + 1
324 local_particles%list(kind_a)%array(nparticle_local(kind_a)) = atom_a
325 END DO
326 END DO
327
328 END DO
329
330! *** Print distribution, if requested ***
331
332 IF (btest(cp_print_key_should_output(logger%iter_info, &
333 force_env_section, "PRINT%DISTRIBUTION1D"), cp_p_file)) THEN
334
335 output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%DISTRIBUTION1D", &
336 extension=".Log")
337
338 iw = output_unit
339 IF (output_unit < 0) iw = cp_logger_get_default_unit_nr(logger, local=.true.)
340
341! *** Print molecule distribution ***
342
343 ALLOCATE (work(npe))
344 work(:) = 0
345
346 work(mype) = sum(nmolecule_local)
347 CALL group%sum(work)
348
349 IF (output_unit > 0) THEN
350 WRITE (unit=output_unit, &
351 fmt="(/, T2, A, T51, A, /, (T52, I6, T73, I8))") &
352 "DISTRIBUTION OF THE MOLECULES", &
353 "Process Number of molecules", &
354 (ipe - 1, work(ipe), ipe=1, npe)
355 WRITE (unit=output_unit, fmt="(T55, A3, T73, I8)") &
356 "Sum", sum(work)
357 CALL m_flush(output_unit)
358 END IF
359
360 CALL group%sync()
361
362 DO ipe = 1, npe
363 IF (ipe == mype) THEN
364 WRITE (unit=iw, fmt="(/, T3, A)") &
365 "Process Kind Local molecules (global indices)"
366 DO imolecule_kind = 1, nmolecule_kind
367 IF (imolecule_kind == 1) THEN
368 WRITE (unit=iw, fmt="(T4, I6, 2X, I5, (T21, 10I6))") &
369 ipe - 1, imolecule_kind, &
370 (local_molecules%list(imolecule_kind)%array(imolecule), &
371 imolecule=1, nmolecule_local(imolecule_kind))
372 ELSE
373 WRITE (unit=iw, fmt="(T12, I5, (T21, 10I6))") &
374 imolecule_kind, &
375 (local_molecules%list(imolecule_kind)%array(imolecule), &
376 imolecule=1, nmolecule_local(imolecule_kind))
377 END IF
378 END DO
379 END IF
380 CALL m_flush(iw)
381 CALL group%sync()
382 END DO
383
384! *** Print particle distribution ***
385
386 work(:) = 0
387
388 work(mype) = sum(nparticle_local)
389 CALL group%sum(work)
390
391 IF (output_unit > 0) THEN
392 WRITE (unit=output_unit, &
393 fmt="(/, T2, A, T51, A, /, (T52, I6, T73, I8))") &
394 "DISTRIBUTION OF THE PARTICLES", &
395 "Process Number of particles", &
396 (ipe - 1, work(ipe), ipe=1, npe)
397 WRITE (unit=output_unit, fmt="(T55, A3, T73, I8)") &
398 "Sum", sum(work)
399 CALL m_flush(output_unit)
400 END IF
401
402 CALL group%sync()
403
404 DO ipe = 1, npe
405 IF (ipe == mype) THEN
406 WRITE (unit=iw, fmt="(/, T3, A)") &
407 "Process Kind Local particles (global indices)"
408 DO iparticle_kind = 1, nparticle_kind
409 IF (iparticle_kind == 1) THEN
410 WRITE (unit=iw, fmt="(T4, I6, 2X, I5, (T20, 10I6))") &
411 ipe - 1, iparticle_kind, &
412 (local_particles%list(iparticle_kind)%array(iatom), &
413 iatom=1, nparticle_local(iparticle_kind))
414 ELSE
415 WRITE (unit=iw, fmt="(T12, I5, (T20, 10I6))") &
416 iparticle_kind, &
417 (local_particles%list(iparticle_kind)%array(iatom), &
418 iatom=1, nparticle_local(iparticle_kind))
419 END IF
420 END DO
421 END IF
422 CALL m_flush(iw)
423 CALL group%sync()
424 END DO
425 DEALLOCATE (work)
426
427 CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
428 "PRINT%DISTRIBUTION1D")
429 END IF
430 END associate
431! *** Release work storage ***
432
433 DEALLOCATE (workload_count)
434
435 DEALLOCATE (workload_fill)
436
437 DEALLOCATE (nmolecule_local)
438
439 DEALLOCATE (nparticle_local)
440
441 DO imolecule_kind = 1, nmolecule_kind
442 IF (ASSOCIATED(local_molecule(imolecule_kind)%array)) THEN
443 DEALLOCATE (local_molecule(imolecule_kind)%array)
444 END IF
445 END DO
446 DEALLOCATE (local_molecule)
447
448 CALL timestop(handle)
449
450 END SUBROUTINE distribute_molecules_1d
451
452! **************************************************************************************************
453!> \brief Distributes the particle pairs creating a 2d distribution optimally
454!> suited for quickstep
455!> \param cell ...
456!> \param atomic_kind_set ...
457!> \param particle_set ...
458!> \param qs_kind_set ...
459!> \param molecule_kind_set ...
460!> \param molecule_set ...
461!> \param distribution_2d the distribution that will be created by this
462!> method
463!> \param blacs_env the parallel environment at the basis of the
464!> distribution
465!> \param force_env_section ...
466!> \par History
467!> - local_rows & cols blocksize optimizations (Aug. 2003, MK)
468!> - cleanup of distribution_2d (Sep. 2003, fawzi)
469!> - update for molecules (Oct. 2003, MK)
470!> \author fawzi (Feb. 2003)
471!> \note
472!> Intermediate generation of a 2d distribution of the molecules, but
473!> only the corresponding particle (atomic) distribution is currently
474!> used. The 2d distribution of the molecules is deleted, but may easily
475!> be recovered (MK).
476! **************************************************************************************************
477 SUBROUTINE distribute_molecules_2d(cell, atomic_kind_set, particle_set, &
478 qs_kind_set, molecule_kind_set, molecule_set, &
479 distribution_2d, blacs_env, force_env_section)
480 TYPE(cell_type), POINTER :: cell
481 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
482 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
483 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
484 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
485 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
486 TYPE(distribution_2d_type), POINTER :: distribution_2d
487 TYPE(cp_blacs_env_type), POINTER :: blacs_env
488 TYPE(section_vals_type), POINTER :: force_env_section
489
490 CHARACTER(len=*), PARAMETER :: routinen = 'distribute_molecules_2d'
491
492 INTEGER :: cluster_price, cost_model, handle, iatom, iatom_mol, iatom_one, ikind, imol, &
493 imolecule, imolecule_kind, iparticle_kind, ipcol, iprow, iw, kind_a, n, natom, natom_mol, &
494 nclusters, nmolecule, nmolecule_kind, nparticle_kind, nsgf, output_unit
495 INTEGER, ALLOCATABLE, DIMENSION(:) :: cluster_list, cluster_prices, &
496 nparticle_local_col, &
497 nparticle_local_row, work
498 INTEGER, DIMENSION(:), POINTER :: lmax_basis, molecule_list
499 INTEGER, DIMENSION(:, :), POINTER :: cluster_col_distribution, &
500 cluster_row_distribution, &
501 col_distribution, row_distribution
502 LOGICAL :: basic_cluster_optimization, basic_optimization, basic_spatial_optimization, &
503 molecular_distribution, skip_optimization
504 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coords, pbc_scaled_coords
505 REAL(kind=dp), DIMENSION(3) :: center
506 TYPE(cp_1d_i_p_type), DIMENSION(:), POINTER :: local_particle_col, local_particle_row
507 TYPE(cp_logger_type), POINTER :: logger
508 TYPE(gto_basis_set_type), POINTER :: orb_basis_set
509 TYPE(molecule_kind_type), POINTER :: molecule_kind
510 TYPE(section_vals_type), POINTER :: distribution_section
511
512!...
513
514 CALL timeset(routinen, handle)
515
516 logger => cp_get_default_logger()
517
518 distribution_section => section_vals_get_subs_vals(force_env_section, "DFT%QS%DISTRIBUTION")
519
520 CALL section_vals_val_get(distribution_section, "2D_MOLECULAR_DISTRIBUTION", l_val=molecular_distribution)
521 CALL section_vals_val_get(distribution_section, "SKIP_OPTIMIZATION", l_val=skip_optimization)
522 CALL section_vals_val_get(distribution_section, "BASIC_OPTIMIZATION", l_val=basic_optimization)
523 CALL section_vals_val_get(distribution_section, "BASIC_SPATIAL_OPTIMIZATION", l_val=basic_spatial_optimization)
524 CALL section_vals_val_get(distribution_section, "BASIC_CLUSTER_OPTIMIZATION", l_val=basic_cluster_optimization)
525
526 CALL section_vals_val_get(distribution_section, "COST_MODEL", i_val=cost_model)
527 !
528
529 associate(group => blacs_env%para_env, myprow => blacs_env%mepos(1) + 1, mypcol => blacs_env%mepos(2) + 1, &
530 nprow => blacs_env%num_pe(1), npcol => blacs_env%num_pe(2))
531
532 nmolecule_kind = SIZE(molecule_kind_set)
533 CALL get_molecule_kind_set(molecule_kind_set, nmolecule=nmolecule)
534
535 nparticle_kind = SIZE(atomic_kind_set)
536 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
537
538 !
539 ! we need to generate two representations of the distribution, one as a straight array with global particles
540 ! one ordered wrt to kinds and only listing the local particles
541 !
542 ALLOCATE (row_distribution(natom, 2))
543 ALLOCATE (col_distribution(natom, 2))
544 ! Initialize the distributions to -1, as the second dimension only gets set with cluster optimization
545 ! but the information is needed by dbcsr
546 row_distribution = -1; col_distribution = -1
547
548 ALLOCATE (local_particle_col(nparticle_kind))
549 ALLOCATE (local_particle_row(nparticle_kind))
550 ALLOCATE (nparticle_local_row(nparticle_kind))
551 ALLOCATE (nparticle_local_col(nparticle_kind))
552
553 IF (basic_optimization .OR. basic_spatial_optimization .OR. basic_cluster_optimization) THEN
554
555 IF (molecular_distribution) THEN
556 nclusters = nmolecule
557 ELSE
558 nclusters = natom
559 END IF
560
561 ALLOCATE (cluster_list(nclusters))
562 ALLOCATE (cluster_prices(nclusters))
563 ALLOCATE (cluster_row_distribution(nclusters, 2))
564 ALLOCATE (cluster_col_distribution(nclusters, 2))
565 cluster_row_distribution = -1; cluster_col_distribution = -1
566
567 ! Fill in the clusters and their prices
568 CALL section_vals_val_get(distribution_section, "COST_MODEL", i_val=cost_model)
569 IF (.NOT. molecular_distribution) THEN
570 DO iatom = 1, natom
571 IF (iatom > nclusters) THEN
572 cpabort("Bounds error")
573 END IF
574 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
575 cluster_list(iatom) = iatom
576 SELECT CASE (cost_model)
577 CASE (model_block_count)
578 CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf)
579 cluster_price = nsgf
580 CASE (model_block_lmax)
581 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
582 CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
583 cluster_price = maxval(lmax_basis)
584 CASE default
585 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
586 CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
587 cluster_price = 8 + (maxval(lmax_basis)**2)
588 END SELECT
589 cluster_prices(iatom) = cluster_price
590 END DO
591 ELSE
592 imol = 0
593 DO imolecule_kind = 1, nmolecule_kind
594 molecule_kind => molecule_kind_set(imolecule_kind)
595 CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, natom=natom_mol)
596 DO imolecule = 1, SIZE(molecule_list)
597 imol = imol + 1
598 cluster_list(imol) = imol
599 cluster_price = 0
600 DO iatom_mol = 1, natom_mol
601 iatom = molecule_set(molecule_list(imolecule))%first_atom + iatom_mol - 1
602 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
603 SELECT CASE (cost_model)
604 CASE (model_block_count)
605 CALL get_qs_kind(qs_kind_set(ikind), nsgf=nsgf)
606 cluster_price = cluster_price + nsgf
607 CASE (model_block_lmax)
608 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
609 CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
610 cluster_price = cluster_price + maxval(lmax_basis)
611 CASE default
612 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
613 CALL get_gto_basis_set(orb_basis_set, lmax=lmax_basis)
614 cluster_price = cluster_price + 8 + (maxval(lmax_basis)**2)
615 END SELECT
616 END DO
617 cluster_prices(imol) = cluster_price
618 END DO
619 END DO
620 END IF
621
622 ! And distribute
623 IF (basic_optimization) THEN
624 CALL make_basic_distribution(cluster_list, cluster_prices, &
625 nprow, cluster_row_distribution(:, 1), npcol, cluster_col_distribution(:, 1))
626 ELSE
627 IF (basic_cluster_optimization) THEN
628 IF (molecular_distribution) THEN
629 cpabort("clustering and molecular blocking NYI")
630 END IF
631 ALLOCATE (pbc_scaled_coords(3, natom), coords(3, natom))
632 DO iatom = 1, natom
633 CALL real_to_scaled(pbc_scaled_coords(:, iatom), pbc(particle_set(iatom)%r(:), cell), cell)
634 coords(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
635 END DO
636 CALL make_cluster_distribution(coords, pbc_scaled_coords, cell, cluster_prices, &
637 nprow, cluster_row_distribution, npcol, cluster_col_distribution)
638 ELSE ! basic_spatial_optimization
639 ALLOCATE (pbc_scaled_coords(3, nclusters))
640 IF (.NOT. molecular_distribution) THEN
641 ! just scaled coords
642 DO iatom = 1, natom
643 CALL real_to_scaled(pbc_scaled_coords(:, iatom), pbc(particle_set(iatom)%r(:), cell), cell)
644 END DO
645 ELSE
646 ! use scaled coords of geometric center, folding when appropriate
647 imol = 0
648 DO imolecule_kind = 1, nmolecule_kind
649 molecule_kind => molecule_kind_set(imolecule_kind)
650 CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, natom=natom_mol)
651 DO imolecule = 1, SIZE(molecule_list)
652 imol = imol + 1
653 iatom_one = molecule_set(molecule_list(imolecule))%first_atom
654 center = 0.0_dp
655 DO iatom_mol = 1, natom_mol
656 iatom = molecule_set(molecule_list(imolecule))%first_atom + iatom_mol - 1
657 center = center + &
658 pbc(particle_set(iatom)%r(:) - particle_set(iatom_one)%r(:), cell) + particle_set(iatom_one)%r(:)
659 END DO
660 center = center/natom_mol
661 CALL real_to_scaled(pbc_scaled_coords(:, imol), pbc(center, cell), cell)
662 END DO
663 END DO
664 END IF
665
666 CALL make_basic_spatial_distribution(pbc_scaled_coords, cluster_prices, &
667 nprow, cluster_row_distribution(:, 1), npcol, cluster_col_distribution(:, 1))
668
669 DEALLOCATE (pbc_scaled_coords)
670 END IF
671 END IF
672
673 ! And assign back
674 IF (.NOT. molecular_distribution) THEN
675 row_distribution = cluster_row_distribution
676 col_distribution = cluster_col_distribution
677 ELSE
678 imol = 0
679 DO imolecule_kind = 1, nmolecule_kind
680 molecule_kind => molecule_kind_set(imolecule_kind)
681 CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, natom=natom_mol)
682 DO imolecule = 1, SIZE(molecule_list)
683 imol = imol + 1
684 DO iatom_mol = 1, natom_mol
685 iatom = molecule_set(molecule_list(imolecule))%first_atom + iatom_mol - 1
686 row_distribution(iatom, :) = cluster_row_distribution(imol, :)
687 col_distribution(iatom, :) = cluster_col_distribution(imol, :)
688 END DO
689 END DO
690 END DO
691 END IF
692
693 ! cleanup
694 DEALLOCATE (cluster_list)
695 DEALLOCATE (cluster_prices)
696 DEALLOCATE (cluster_row_distribution)
697 DEALLOCATE (cluster_col_distribution)
698
699 ELSE
700 ! expects nothing else
701 cpabort("Invalid optimization for DFT%QS%DISTRIBUTION")
702 END IF
703
704 ! prepare the lists of local particles
705
706 ! count local particles of a given kind
707 nparticle_local_col = 0
708 nparticle_local_row = 0
709 DO iatom = 1, natom
710 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, kind_number=kind_a)
711 IF (row_distribution(iatom, 1) == myprow) nparticle_local_row(kind_a) = nparticle_local_row(kind_a) + 1
712 IF (col_distribution(iatom, 1) == mypcol) nparticle_local_col(kind_a) = nparticle_local_col(kind_a) + 1
713 END DO
714
715 ! allocate space
716 DO iparticle_kind = 1, nparticle_kind
717 n = nparticle_local_row(iparticle_kind)
718 ALLOCATE (local_particle_row(iparticle_kind)%array(n))
719
720 n = nparticle_local_col(iparticle_kind)
721 ALLOCATE (local_particle_col(iparticle_kind)%array(n))
722 END DO
723
724 ! store
725 nparticle_local_col = 0
726 nparticle_local_row = 0
727 DO iatom = 1, natom
728 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, kind_number=kind_a)
729 IF (row_distribution(iatom, 1) == myprow) THEN
730 nparticle_local_row(kind_a) = nparticle_local_row(kind_a) + 1
731 local_particle_row(kind_a)%array(nparticle_local_row(kind_a)) = iatom
732 END IF
733 IF (col_distribution(iatom, 1) == mypcol) THEN
734 nparticle_local_col(kind_a) = nparticle_local_col(kind_a) + 1
735 local_particle_col(kind_a)%array(nparticle_local_col(kind_a)) = iatom
736 END IF
737 END DO
738
739! *** Generate the 2d distribution structure but take care of the zero offsets required
740 row_distribution(:, 1) = row_distribution(:, 1) - 1
741 col_distribution(:, 1) = col_distribution(:, 1) - 1
742 CALL distribution_2d_create(distribution_2d, &
743 row_distribution_ptr=row_distribution, &
744 col_distribution_ptr=col_distribution, &
745 local_rows_ptr=local_particle_row, &
746 local_cols_ptr=local_particle_col, &
747 blacs_env=blacs_env)
748
749 NULLIFY (local_particle_row)
750 NULLIFY (local_particle_col)
751 NULLIFY (row_distribution)
752 NULLIFY (col_distribution)
753
754! *** Print distribution, if requested ***
755 IF (btest(cp_print_key_should_output(logger%iter_info, &
756 force_env_section, "PRINT%DISTRIBUTION"), cp_p_file)) THEN
757
758 output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%DISTRIBUTION", &
759 extension=".Log")
760
761! *** Print row distribution ***
762
763 ALLOCATE (work(nprow))
764 work(:) = 0
765
766 IF (mypcol == 1) work(myprow) = sum(distribution_2d%n_local_rows)
767
768 CALL group%sum(work)
769
770 IF (output_unit > 0) THEN
771 WRITE (unit=output_unit, &
772 fmt="(/, T2, A, /, T15, A, /, (T16, I10, T41, I10, T71, I10))") &
773 "DISTRIBUTION OF THE PARTICLES (ROWS)", &
774 "Process row Number of particles Number of matrix rows", &
775 (iprow - 1, work(iprow), -1, iprow=1, nprow)
776 WRITE (unit=output_unit, fmt="(T23, A3, T41, I10, T71, I10)") &
777 "Sum", sum(work), -1
778 CALL m_flush(output_unit)
779 END IF
780
781 DEALLOCATE (work)
782
783! *** Print column distribution ***
784
785 ALLOCATE (work(npcol))
786 work(:) = 0
787
788 IF (myprow == 1) work(mypcol) = sum(distribution_2d%n_local_cols)
789
790 CALL group%sum(work)
791
792 IF (output_unit > 0) THEN
793 WRITE (unit=output_unit, &
794 fmt="(/, T2, A, /, T15, A, /, (T16, I10, T41, I10, T71, I10))") &
795 "DISTRIBUTION OF THE PARTICLES (COLUMNS)", &
796 "Process col Number of particles Number of matrix columns", &
797 (ipcol - 1, work(ipcol), -1, ipcol=1, npcol)
798 WRITE (unit=output_unit, fmt="(T23, A3, T41, I10, T71, I10)") &
799 "Sum", sum(work), -1
800 CALL m_flush(output_unit)
801 END IF
802
803 DEALLOCATE (work)
804
805 CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
806 "PRINT%DISTRIBUTION")
807 END IF
808 END associate
809
810 IF (btest(cp_print_key_should_output(logger%iter_info, &
811 force_env_section, "PRINT%DISTRIBUTION2D"), cp_p_file)) THEN
812
813 iw = cp_logger_get_default_unit_nr(logger, local=.true.)
814 CALL distribution_2d_write(distribution_2d, &
815 unit_nr=iw, &
816 local=.true., &
817 long_description=.true.)
818
819 END IF
820
821! *** Release work storage ***
822
823 DEALLOCATE (nparticle_local_row)
824
825 DEALLOCATE (nparticle_local_col)
826
827 CALL timestop(handle)
828
829 END SUBROUTINE distribute_molecules_2d
830
831! **************************************************************************************************
832!> \brief Creates a basic distribution
833!> \param cluster_list ...
834!> \param cluster_prices ...
835!> \param nprows ...
836!> \param row_distribution ...
837!> \param npcols ...
838!> \param col_distribution ...
839!> \par History
840!> - Created 2010-08-06 UB
841! **************************************************************************************************
842 SUBROUTINE make_basic_distribution(cluster_list, cluster_prices, &
843 nprows, row_distribution, npcols, col_distribution)
844 INTEGER, DIMENSION(:), INTENT(INOUT) :: cluster_list, cluster_prices
845 INTEGER, INTENT(IN) :: nprows
846 INTEGER, DIMENSION(:), INTENT(OUT) :: row_distribution
847 INTEGER, INTENT(IN) :: npcols
848 INTEGER, DIMENSION(:), INTENT(OUT) :: col_distribution
849
850 CHARACTER(len=*), PARAMETER :: routinen = 'make_basic_distribution'
851
852 INTEGER :: bin, cluster, cluster_index, &
853 cluster_price, nbins, nclusters, pcol, &
854 pgrid_gcd, prow, timing_handle
855 INTEGER(int_8) :: bin_price
856 LOGICAL :: found
857 TYPE(cp_heap_type) :: bin_heap
858
859! ---------------------------------------------------------------------------
860
861 CALL timeset(routinen, timing_handle)
862 nbins = lcm(nprows, npcols)
863 pgrid_gcd = gcd(nprows, npcols)
864 CALL sort(cluster_prices, SIZE(cluster_list), cluster_list)
865 CALL cp_heap_new(bin_heap, nbins)
866 CALL cp_heap_fill(bin_heap, [(0_int_8, bin=1, nbins)])
867 !
868 nclusters = SIZE(cluster_list)
869 ! Put the most expensive cluster in the bin with the smallest
870 ! price and repeat.
871 DO cluster_index = nclusters, 1, -1
872 cluster = cluster_list(cluster_index)
873 CALL cp_heap_get_first(bin_heap, bin, bin_price, found)
874 IF (.NOT. found) THEN
875 cpabort("No topmost heap element found.")
876 END IF
877 !
878 prow = int((bin - 1)*pgrid_gcd/npcols)
879 IF (prow >= nprows) THEN
880 cpabort("Invalid process row.")
881 END IF
882 pcol = int((bin - 1)*pgrid_gcd/nprows)
883 IF (pcol >= npcols) THEN
884 cpabort("Invalid process column.")
885 END IF
886 row_distribution(cluster) = prow + 1
887 col_distribution(cluster) = pcol + 1
888 !
889 cluster_price = cluster_prices(cluster_index)
890 bin_price = bin_price + cluster_price
891 CALL cp_heap_reset_first(bin_heap, bin_price)
892 END DO
893 CALL cp_heap_release(bin_heap)
894 CALL timestop(timing_handle)
895 END SUBROUTINE make_basic_distribution
896
897! **************************************************************************************************
898!> \brief Creates a basic spatial distribution
899!> that tries to make the corresponding blocks as homogeneous as possible
900!> \param pbc_scaled_coords ...
901!> \param costs ...
902!> \param nprows ...
903!> \param row_distribution ...
904!> \param npcols ...
905!> \param col_distribution ...
906!> \par History
907!> - Created 2010-11-11 Joost VandeVondele
908! **************************************************************************************************
909 SUBROUTINE make_basic_spatial_distribution(pbc_scaled_coords, costs, &
910 nprows, row_distribution, npcols, col_distribution)
911 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pbc_scaled_coords
912 INTEGER, DIMENSION(:), INTENT(IN) :: costs
913 INTEGER, INTENT(IN) :: nprows
914 INTEGER, DIMENSION(:), INTENT(OUT) :: row_distribution
915 INTEGER, INTENT(IN) :: npcols
916 INTEGER, DIMENSION(:), INTENT(OUT) :: col_distribution
917
918 CHARACTER(len=*), PARAMETER :: routinen = 'make_basic_spatial_distribution'
919
920 INTEGER :: handle, iatom, natoms, nbins, pgrid_gcd
921 INTEGER, ALLOCATABLE, DIMENSION(:) :: bin_costs, distribution
922
923 CALL timeset(routinen, handle)
924
925 natoms = SIZE(costs)
926 nbins = lcm(nprows, npcols)
927 pgrid_gcd = gcd(nprows, npcols)
928 ALLOCATE (bin_costs(nbins), distribution(natoms))
929 bin_costs = 0
930
931 CALL spatial_recurse(pbc_scaled_coords, costs, [(iatom, iatom=1, natoms)], bin_costs, distribution, 0)
932
933 ! WRITE(*, *) "Final bin costs: ", bin_costs
934
935 ! final row_distribution / col_distribution
936 DO iatom = 1, natoms
937 row_distribution(iatom) = (distribution(iatom) - 1)*pgrid_gcd/npcols + 1
938 col_distribution(iatom) = (distribution(iatom) - 1)*pgrid_gcd/nprows + 1
939 END DO
940
941 DEALLOCATE (bin_costs, distribution)
942
943 CALL timestop(handle)
944
945 END SUBROUTINE make_basic_spatial_distribution
946
947! **************************************************************************************************
948!> \brief ...
949!> \param pbc_scaled_coords ...
950!> \param costs ...
951!> \param indices ...
952!> \param bin_costs ...
953!> \param distribution ...
954!> \param level ...
955! **************************************************************************************************
956 RECURSIVE SUBROUTINE spatial_recurse(pbc_scaled_coords, costs, indices, bin_costs, distribution, level)
957 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pbc_scaled_coords
958 INTEGER, DIMENSION(:), INTENT(IN) :: costs, indices
959 INTEGER, DIMENSION(:), INTENT(INOUT) :: bin_costs, distribution
960 INTEGER, INTENT(IN) :: level
961
962 INTEGER :: iatom, ibin, natoms, nbins, nhalf
963 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_costs_sorted, atom_permutation, &
964 bin_costs_sorted, permutation
965 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coord
966
967 natoms = SIZE(costs)
968 nbins = SIZE(bin_costs)
969 nhalf = (natoms + 1)/2
970
971 IF (natoms <= nbins) THEN
972 ! assign the most expensive atom to the least costly bin
973 ALLOCATE (bin_costs_sorted(nbins), permutation(nbins))
974 bin_costs_sorted(:) = bin_costs
975 CALL sort(bin_costs_sorted, nbins, permutation)
976 ALLOCATE (atom_costs_sorted(natoms), atom_permutation(natoms))
977 atom_costs_sorted(:) = costs
978 CALL sort(atom_costs_sorted, natoms, atom_permutation)
979 ibin = 0
980 ! WRITE(*, *) "Dealing with a new bunch of atoms "
981 DO iatom = natoms, 1, -1
982 ibin = ibin + 1
983 ! WRITE(*, *) "atom", indices(atom_permutation(iatom)), "cost", atom_costs_sorted(iatom), &
984 ! "bin", permutation(ibin), "its cost", bin_costs(permutation(ibin))
985 ! WRITE(100, '(A, I0, 3F12.6)') "A", permutation(ibin), pbc_scaled_coords(:, atom_permutation(iatom))
986 bin_costs(permutation(ibin)) = bin_costs(permutation(ibin)) + atom_costs_sorted(iatom)
987 distribution(indices(atom_permutation(iatom))) = permutation(ibin)
988 END DO
989 DEALLOCATE (bin_costs_sorted, permutation, atom_costs_sorted, atom_permutation)
990 ELSE
991 ! divide atoms in two subsets, sorting according to their coordinates, alternatively x, y, z
992 ! recursively do this for both subsets
993 ALLOCATE (coord(natoms), permutation(natoms))
994 coord(:) = pbc_scaled_coords(mod(level, 3) + 1, :)
995 CALL sort(coord, natoms, permutation)
996 CALL spatial_recurse(pbc_scaled_coords(:, permutation(1:nhalf)), costs(permutation(1:nhalf)), &
997 indices(permutation(1:nhalf)), bin_costs, distribution, level + 1)
998 CALL spatial_recurse(pbc_scaled_coords(:, permutation(nhalf + 1:)), costs(permutation(nhalf + 1:)), &
999 indices(permutation(nhalf + 1:)), bin_costs, distribution, level + 1)
1000 DEALLOCATE (coord, permutation)
1001 END IF
1002
1003 END SUBROUTINE spatial_recurse
1004
1005! **************************************************************************************************
1006!> \brief creates a distribution placing close by atoms into clusters and
1007!> putting them on the same processors. Load balancing is
1008!> performed by balancing sum of the cluster costs per processor
1009!> \param coords coordinates of the system
1010!> \param scaled_coords scaled coordinates
1011!> \param cell the cell_type
1012!> \param costs costs per atomic block
1013!> \param nprows number of precessors per row on the 2d grid
1014!> \param row_distribution the resulting distribution over proc_rows of atomic blocks
1015!> \param npcols number of precessors per col on the 2d grid
1016!> \param col_distribution the resulting distribution over proc_cols of atomic blocks
1017! **************************************************************************************************
1018 SUBROUTINE make_cluster_distribution(coords, scaled_coords, cell, costs, &
1019 nprows, row_distribution, npcols, col_distribution)
1020 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: coords, scaled_coords
1021 TYPE(cell_type), POINTER :: cell
1022 INTEGER, DIMENSION(:), INTENT(IN) :: costs
1023 INTEGER, INTENT(IN) :: nprows
1024 INTEGER, DIMENSION(:, :), INTENT(OUT) :: row_distribution
1025 INTEGER, INTENT(IN) :: npcols
1026 INTEGER, DIMENSION(:, :), INTENT(OUT) :: col_distribution
1027
1028 CHARACTER(len=*), PARAMETER :: routinen = 'make_cluster_distribution'
1029
1030 INTEGER :: handle, i, icluster, level, natom, &
1031 output_unit
1032 INTEGER(KIND=int_8) :: ncluster
1033 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_to_cluster, cluster_cost, &
1034 cluster_count, cluster_to_col, &
1035 cluster_to_row, piv_cost, proc_cost, &
1036 sorted_cost
1037 REAL(kind=dp) :: fold(3)
1038 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cluster_center, cluster_high, cluster_low
1039
1040 CALL timeset(routinen, handle)
1041
1042 output_unit = cp_logger_get_default_io_unit()
1043
1044 natom = SIZE(costs)
1045 ncluster = dbcsr_distribution_get_num_images(sum(costs), natom, nprows, npcols)
1046 ALLOCATE (atom_to_cluster(natom))
1047 ALLOCATE (cluster_cost(ncluster))
1048 ALLOCATE (cluster_to_row(ncluster))
1049 ALLOCATE (cluster_to_col(ncluster))
1050 ALLOCATE (sorted_cost(ncluster))
1051 ALLOCATE (piv_cost(ncluster))
1052 cluster_cost(:) = 0
1053
1054 icluster = 0
1055 CALL cluster_recurse(coords, scaled_coords, cell, costs, atom_to_cluster, ncluster, icluster, cluster_cost)
1056
1057 sorted_cost(:) = cluster_cost(:)
1058 CALL sort(sorted_cost, int(ncluster), piv_cost)
1059
1060 ALLOCATE (proc_cost(nprows))
1061 proc_cost = 0; level = 1
1062 CALL assign_clusters(cluster_cost, piv_cost, proc_cost, cluster_to_row, nprows)
1063
1064 DEALLOCATE (proc_cost); ALLOCATE (proc_cost(npcols))
1065 proc_cost = 0; level = 1
1066 CALL assign_clusters(cluster_cost, piv_cost, proc_cost, cluster_to_col, npcols)
1067
1068 DO i = 1, natom
1069 row_distribution(i, 1) = cluster_to_row(atom_to_cluster(i))
1070 row_distribution(i, 2) = atom_to_cluster(i)
1071 col_distribution(i, 1) = cluster_to_col(atom_to_cluster(i))
1072 col_distribution(i, 2) = atom_to_cluster(i)
1073 END DO
1074
1075 ! generate some statistics on clusters
1076 ALLOCATE (cluster_center(3, ncluster))
1077 ALLOCATE (cluster_low(3, ncluster))
1078 ALLOCATE (cluster_high(3, ncluster))
1079 ALLOCATE (cluster_count(ncluster))
1080 cluster_count = 0
1081 DO i = 1, natom
1082 cluster_count(atom_to_cluster(i)) = cluster_count(atom_to_cluster(i)) + 1
1083 cluster_center(:, atom_to_cluster(i)) = coords(:, i)
1084 END DO
1085 cluster_low = huge(0.0_dp)/2
1086 cluster_high = -huge(0.0_dp)/2
1087 DO i = 1, natom
1088 fold = pbc(coords(:, i) - cluster_center(:, atom_to_cluster(i)), cell) + cluster_center(:, atom_to_cluster(i))
1089 cluster_low(:, atom_to_cluster(i)) = min(cluster_low(:, atom_to_cluster(i)), fold(:))
1090 cluster_high(:, atom_to_cluster(i)) = max(cluster_high(:, atom_to_cluster(i)), fold(:))
1091 END DO
1092 IF (output_unit > 0) THEN
1093 WRITE (output_unit, *)
1094 WRITE (output_unit, '(T2,A)') "Cluster distribution information"
1095 WRITE (output_unit, '(T2,A,T48,I8)') "Number of atoms", natom
1096 WRITE (output_unit, '(T2,A,T48,I8)') "Number of clusters", ncluster
1097 WRITE (output_unit, '(T2,A,T48,I8)') "Largest cluster in atoms", maxval(cluster_count)
1098 WRITE (output_unit, '(T2,A,T48,I8)') "Smallest cluster in atoms", minval(cluster_count)
1099 WRITE (output_unit, '(T2,A,T48,F8.3,I8)') "Largest cartesian extend [a.u.]/cluster x=", &
1100 maxval(cluster_high(1, :) - cluster_low(1, :), mask=(cluster_count > 0)), &
1101 maxloc(cluster_high(1, :) - cluster_low(1, :), mask=(cluster_count > 0))
1102 WRITE (output_unit, '(T2,A,T48,F8.3,I8)') "Largest cartesian extend [a.u.]/cluster y=", &
1103 maxval(cluster_high(2, :) - cluster_low(2, :), mask=(cluster_count > 0)), &
1104 maxloc(cluster_high(2, :) - cluster_low(2, :), mask=(cluster_count > 0))
1105 WRITE (output_unit, '(T2,A,T48,F8.3,I8)') "Largest cartesian extend [a.u.]/cluster z=", &
1106 maxval(cluster_high(3, :) - cluster_low(3, :), mask=(cluster_count > 0)), &
1107 maxloc(cluster_high(3, :) - cluster_low(3, :), mask=(cluster_count > 0))
1108 END IF
1109
1110 DEALLOCATE (atom_to_cluster, cluster_cost, cluster_to_row, cluster_to_col, sorted_cost, piv_cost, proc_cost)
1111 CALL timestop(handle)
1112
1113 END SUBROUTINE make_cluster_distribution
1114
1115! **************************************************************************************************
1116!> \brief assigns the clusters to processors, tryimg to balance the cost on the nodes
1117!> \param cluster_cost vector with the cost of each cluster
1118!> \param piv_cost pivoting vector sorting the cluster_cost
1119!> \param proc_cost cost per processor, on input 0 everywhere
1120!> \param cluster_assign assgnment of clusters on proc
1121!> \param nproc number of processor over which clusters are distributed
1122! **************************************************************************************************
1123 SUBROUTINE assign_clusters(cluster_cost, piv_cost, proc_cost, cluster_assign, nproc)
1124 INTEGER, ALLOCATABLE, DIMENSION(:) :: cluster_cost, piv_cost, proc_cost, &
1125 cluster_assign
1126 INTEGER :: nproc
1127
1128 CHARACTER(len=*), PARAMETER :: routinen = 'assign_clusters'
1129
1130 INTEGER :: handle, i, ilevel, offset, &
1131 piv_pcost(nproc), sort_proc_cost(nproc)
1132
1133 CALL timeset(routinen, handle)
1134
1135 DO ilevel = 1, SIZE(cluster_cost)/nproc
1136 sort_proc_cost(:) = proc_cost(:)
1137 CALL sort(sort_proc_cost, nproc, piv_pcost)
1138
1139 offset = (SIZE(cluster_cost)/nproc - ilevel + 1)*nproc + 1
1140 DO i = 1, nproc
1141 cluster_assign(piv_cost(offset - i)) = piv_pcost(i)
1142 proc_cost(piv_pcost(i)) = proc_cost(piv_pcost(i)) + cluster_cost(piv_cost(offset - i))
1143 END DO
1144 END DO
1145
1146 CALL timestop(handle)
1147
1148 END SUBROUTINE assign_clusters
1149
1150! **************************************************************************************************
1151!> \brief recursive routine to cluster atoms.
1152!> Low level uses a modified KMEANS algorithm
1153!> recursion is used to reduce cost.
1154!> each level will subdivide a cluster into smaller clusters
1155!> If only a single split is necessary atoms are assigned to the current cluster
1156!> \param coord coordinates of the system
1157!> \param scaled_coord scaled coordinates
1158!> \param cell the cell_type
1159!> \param costs costs per atomic block
1160!> \param cluster_inds the atom_to cluster mapping
1161!> \param ncluster number of clusters still to be created on a given recursion level
1162!> \param icluster the index of the current cluster to be created
1163!> \param fin_cluster_cost total cost of the final clusters
1164! **************************************************************************************************
1165 RECURSIVE SUBROUTINE cluster_recurse(coord, scaled_coord, cell, costs, cluster_inds, ncluster, icluster, fin_cluster_cost)
1166 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: coord, scaled_coord
1167 TYPE(cell_type), POINTER :: cell
1168 INTEGER, DIMENSION(:), INTENT(IN) :: costs
1169 INTEGER, DIMENSION(:), INTENT(INOUT) :: cluster_inds
1170 INTEGER(KIND=int_8), INTENT(INOUT) :: ncluster
1171 INTEGER, INTENT(INOUT) :: icluster
1172 INTEGER, DIMENSION(:), INTENT(INOUT) :: fin_cluster_cost
1173
1174 INTEGER :: i, ibeg, iend, maxv(1), min_seed, &
1175 natoms, nleft, nsplits, seed, tot_cost
1176 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: ncluster_new
1177 INTEGER, ALLOCATABLE, DIMENSION(:) :: cluster_cost, inds_tmp, nat_cluster, piv
1178 LOGICAL :: found
1179 REAL(kind=dp) :: balance, balance_new, conv
1180
1181 natoms = SIZE(coord, 2)
1182 ! 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
1183 ! small systems or subclusters
1184 IF (natoms <= 1) THEN
1185 nsplits = 1
1186 ELSE
1187 nsplits = min(int(min(int(max(6, int(60.00/log(real(natoms, kind=dp)))), kind=int_8), ncluster)), natoms)
1188 END IF
1189 IF (nsplits == 1) THEN
1190 icluster = icluster + 1
1191 cluster_inds = icluster
1192 fin_cluster_cost(icluster) = sum(costs)
1193 ELSE
1194 ALLOCATE (cluster_cost(nsplits), ncluster_new(nsplits), inds_tmp(natoms), piv(natoms), nat_cluster(nsplits))
1195 ! initialise some values
1196 cluster_cost = 0; seed = 300; found = .true.; min_seed = seed
1197 CALL kmeans(nsplits, coord, scaled_coord, cell, cluster_inds, nat_cluster, seed, conv)
1198 balance = maxval(real(nat_cluster, kind=dp))/minval(real(nat_cluster, kind=dp))
1199
1200 ! If the system is small enough try to do better in terms of balancing number of atoms per cluster
1201 ! by changing the seed for the initial guess
1202 IF (natoms < 1000 .AND. balance > 1.1) THEN
1203 found = .false.
1204 DO i = 1, 5
1205 IF (balance > 1.1) THEN
1206 CALL kmeans(nsplits, coord, scaled_coord, cell, cluster_inds, nat_cluster, seed + i*40, conv)
1207 balance_new = maxval(real(nat_cluster, kind=dp))/minval(real(nat_cluster, kind=dp))
1208 IF (balance_new < balance) THEN
1209 balance = balance_new
1210 min_seed = seed + i*40
1211 END IF
1212 ELSE
1213 found = .true.
1214 EXIT
1215 END IF
1216 END DO
1217 END IF
1218 !If we do not match the convergence than recompute at least the best assignment
1219 IF (.NOT. found) CALL kmeans(nsplits, coord, scaled_coord, cell, cluster_inds, nat_cluster, min_seed, conv)
1220
1221 ! compute the cost of each cluster to decide how many splits have to be performed on the next lower level
1222 DO i = 1, natoms
1223 cluster_cost(cluster_inds(i)) = cluster_cost(cluster_inds(i)) + costs(i)
1224 END DO
1225 tot_cost = sum(cluster_cost)
1226 ! compute new splitting, can be done more elegant
1227 ncluster_new(:) = ncluster*cluster_cost(:)/tot_cost
1228 nleft = int(ncluster - sum(ncluster_new))
1229 ! As we won't have empty clusters, we can not have 0 as new size, so we correct for this at first
1230 DO i = 1, nsplits
1231 IF (ncluster_new(i) == 0) THEN
1232 ncluster_new(i) = 1
1233 nleft = nleft - 1
1234 END IF
1235 END DO
1236 ! now comes the next part that the number of clusters will not match anymore, so try to correct in a meaningful way without
1237 ! introducing 0 sized blocks again
1238 IF (nleft /= 0) THEN
1239 DO i = 1, abs(nleft)
1240 IF (nleft < 0) THEN
1241 maxv = minloc(cluster_cost/ncluster_new)
1242 IF (ncluster_new(maxv(1)) /= 1) THEN
1243 ncluster_new(maxv) = ncluster_new(maxv) - 1
1244 ELSE
1245 maxv = maxloc(ncluster_new)
1246 ncluster_new(maxv) = ncluster_new(maxv) - 1
1247 END IF
1248 ELSE
1249 maxv = maxloc(cluster_cost/ncluster_new)
1250 ncluster_new(maxv) = ncluster_new(maxv) + 1
1251 END IF
1252 END DO
1253 END IF
1254
1255 !Now get the permutations to sort the atoms in the nsplits clusters for the next level of iteration
1256 inds_tmp(:) = cluster_inds(:)
1257 CALL sort(inds_tmp, natoms, piv)
1258
1259 ibeg = 1; iend = 0
1260 DO i = 1, nsplits
1261 IF (nat_cluster(i) == 0) cycle
1262 iend = iend + nat_cluster(i)
1263 CALL cluster_recurse(coord(:, piv(ibeg:iend)), scaled_coord(:, piv(ibeg:iend)), cell, costs(piv(ibeg:iend)), &
1264 inds_tmp(ibeg:iend), ncluster_new(i), icluster, fin_cluster_cost)
1265 ibeg = ibeg + nat_cluster(i)
1266 END DO
1267 ! copy the sorted cluster IDs on the old layout, inds_tmp gets set at the lowest level of recursion
1268 cluster_inds(piv(:)) = inds_tmp
1269 DEALLOCATE (cluster_cost, ncluster_new, inds_tmp, piv, nat_cluster)
1270
1271 END IF
1272
1273 END SUBROUTINE cluster_recurse
1274
1275! **************************************************************************************************
1276!> \brief A modified version of the kmeans algorithm.
1277!> The assignment has a penalty function in case clusters become
1278!> larger than average. Like this more even sized clusters are created
1279!> trading it for locality
1280!> \param ncent number of centers to be created
1281!> \param coord coordinates
1282!> \param scaled_coord scaled coord
1283!> \param cell the cell_type
1284!> \param cluster atom to cluster assignment
1285!> \param nat_cl atoms per cluster
1286!> \param seed seed for the RNG. Algorithm might need multiple tries to deliver best results
1287!> \param tot_var the total variance of the clusters around the centers
1288! **************************************************************************************************
1289 SUBROUTINE kmeans(ncent, coord, scaled_coord, cell, cluster, nat_cl, seed, tot_var)
1290 INTEGER :: ncent
1291 REAL(kind=dp), DIMENSION(:, :) :: coord, scaled_coord
1292 TYPE(cell_type), POINTER :: cell
1293 INTEGER, DIMENSION(:) :: cluster, nat_cl
1294 INTEGER :: seed
1295 REAL(kind=dp) :: tot_var
1296
1297 CHARACTER(len=*), PARAMETER :: routinen = 'kmeans'
1298
1299 INTEGER :: handle, i, ind, itn, j, nat, oldc
1300 LOGICAL :: changed
1301 REAL(kind=dp) :: average(3, ncent, 2), cent_coord(3, ncent), devi, deviat(ncent), dist, &
1302 dvec(3), old_var, rn, scaled_cent(3, ncent), var_cl(ncent)
1303 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dmat
1304 REAL(kind=dp), DIMENSION(3, 2) :: initial_seed
1305 TYPE(rng_stream_type) :: rng_stream
1306
1307 CALL timeset(routinen, handle)
1308
1309 initial_seed = real(seed, dp); nat = SIZE(coord, 2)
1310 ALLOCATE (dmat(ncent, nat))
1311
1312 rng_stream = rng_stream_type(name="kmeans uniform distribution [0,1]", &
1313 distribution_type=uniform, seed=initial_seed)
1314
1315! try to find a clever initial guess with centers being somewhat distributed
1316 rn = rng_stream%next()
1317 ind = ceiling(rn*nat)
1318 cent_coord(:, 1) = coord(:, ind)
1319 DO i = 2, ncent
1320 DO
1321 rn = rng_stream%next()
1322 ind = ceiling(rn*nat)
1323 cent_coord(:, i) = coord(:, ind)
1324 devi = huge(1.0_dp)
1325 DO j = 1, i - 1
1326 dvec = pbc(cent_coord(:, j), cent_coord(:, i), cell)
1327 dist = norm2(dvec)
1328 IF (dist < devi) devi = dist
1329 END DO
1330 rn = rng_stream%next()
1331 IF (rn < devi**2/169.0) EXIT
1332 END DO
1333 END DO
1334
1335! Now start the KMEANS but penalise it in case it starts packing too many atoms into a single set
1336! Unfoirtunatelz as this is dependent on what happened before it cant be parallel
1337 cluster = 0; old_var = huge(1.0_dp)
1338 DO itn = 1, 1000
1339 changed = .false.; var_cl = 0.0_dp; tot_var = 0.0_dp; nat_cl = 0; deviat = 0.0_dp
1340! !$OMP PARALLEL DO PRIVATE(i,j,dvec)
1341 DO i = 1, nat
1342 DO j = 1, ncent
1343 dvec = pbc(cent_coord(:, j), coord(:, i), cell)
1344 dmat(j, i) = dot_product(dvec, dvec)
1345 END DO
1346 END DO
1347 DO i = 1, nat
1348 devi = huge(1.0_dp); oldc = cluster(i)
1349 DO j = 1, ncent
1350 dist = dmat(j, i) + max(nat_cl(j)**2/nat*ncent, nat/ncent)
1351 IF (dist < devi) THEN
1352 devi = dist; cluster(i) = j
1353 END IF
1354 END DO
1355 deviat(cluster(i)) = deviat(cluster(i)) + sqrt(devi)
1356 nat_cl(cluster(i)) = nat_cl(cluster(i)) + 1
1357 tot_var = tot_var + devi
1358 IF (oldc /= cluster(i)) changed = .true.
1359 END DO
1360 ! get the update of the centers done, add a new one in case one center lost all its atoms
1361 ! the algorithm would survive, but its nice to really create what you demand
1362 IF (tot_var >= old_var) EXIT
1363 IF (changed) THEN
1364 ! Here misery of computing the center of geometry of the clusters in PBC.
1365 ! The mapping on the unit circle allows to circumvent all problems
1366 average = 0.0_dp
1367 DO i = 1, SIZE(coord, 2)
1368 average(:, cluster(i), 1) = average(:, cluster(i), 1) + cos(scaled_coord(:, i)*2.0_dp*pi)
1369 average(:, cluster(i), 2) = average(:, cluster(i), 2) + sin(scaled_coord(:, i)*2.0_dp*pi)
1370 END DO
1371
1372 DO i = 1, ncent
1373 IF (nat_cl(i) == 0) THEN
1374 rn = rng_stream%next()
1375 scaled_cent(:, i) = scaled_coord(:, ceiling(rn*nat))
1376 ELSE
1377 average(:, i, 1) = average(:, i, 1)/real(nat_cl(i), dp)
1378 average(:, i, 2) = average(:, i, 2)/real(nat_cl(i), dp)
1379 scaled_cent(:, i) = (atan2(-average(:, i, 2), -average(:, i, 1)) + pi)/(2.0_dp*pi)
1380 CALL scaled_to_real(cent_coord(:, i), scaled_cent(:, i), cell)
1381 END IF
1382 END DO
1383 ELSE
1384 EXIT
1385 END IF
1386 END DO
1387
1388 CALL timestop(handle)
1389
1390 END SUBROUTINE kmeans
1391
1392END 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, npgf_seg_sum, ccon)
...
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:565
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:535
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:124
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:1324
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1289
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, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, 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_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, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, 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:60
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.