Line data Source code
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 Routines for a Kim-Gordon-like partitioning into molecular subunits
10 : !> \par History
11 : !> 2012.07 created [Martin Haeufel]
12 : !> \author Martin Haeufel
13 : ! **************************************************************************************************
14 : MODULE kg_environment
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE bibliography, ONLY: Andermatt2016,&
20 : cite_reference
21 : USE cell_types, ONLY: cell_type
22 : USE cp_files, ONLY: close_file,&
23 : open_file
24 : USE cp_log_handling, ONLY: cp_get_default_logger,&
25 : cp_logger_get_default_io_unit,&
26 : cp_logger_type
27 : USE distribution_1d_types, ONLY: distribution_1d_type
28 : USE distribution_2d_types, ONLY: distribution_2d_type
29 : USE external_potential_types, ONLY: get_potential,&
30 : local_potential_type
31 : USE input_constants, ONLY: kg_tnadd_atomic,&
32 : kg_tnadd_embed,&
33 : kg_tnadd_embed_ri,&
34 : kg_tnadd_none
35 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
36 : section_vals_type,&
37 : section_vals_val_get
38 : USE integration_grid_types, ONLY: allocate_intgrid,&
39 : integration_grid_type
40 : USE kg_environment_types, ONLY: kg_environment_type
41 : USE kg_vertex_coloring_methods, ONLY: kg_vertex_coloring
42 : USE kinds, ONLY: dp,&
43 : int_4,&
44 : int_4_size,&
45 : int_8
46 : USE lri_environment_init, ONLY: lri_env_basis,&
47 : lri_env_init
48 : USE message_passing, ONLY: mp_para_env_type
49 : USE molecule_types, ONLY: molecule_type
50 : USE particle_types, ONLY: particle_type
51 : USE qs_environment_types, ONLY: get_qs_env,&
52 : qs_environment_type
53 : USE qs_grid_atom, ONLY: initialize_atomic_grid
54 : USE qs_kind_types, ONLY: get_qs_kind,&
55 : qs_kind_type
56 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
57 : neighbor_list_iterate,&
58 : neighbor_list_iterator_create,&
59 : neighbor_list_iterator_p_type,&
60 : neighbor_list_iterator_release,&
61 : neighbor_list_set_p_type,&
62 : release_neighbor_list_sets
63 : USE qs_neighbor_lists, ONLY: atom2d_build,&
64 : atom2d_cleanup,&
65 : build_neighbor_lists,&
66 : local_atoms_type,&
67 : pair_radius_setup,&
68 : write_neighbor_lists
69 : USE string_utilities, ONLY: uppercase
70 : USE task_list_types, ONLY: deallocate_task_list
71 : USE util, ONLY: sort
72 : #include "./base/base_uses.f90"
73 :
74 : IMPLICIT NONE
75 :
76 : PRIVATE
77 :
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kg_environment'
79 :
80 : PUBLIC :: kg_env_create, kg_build_neighborlist, kg_build_subsets
81 :
82 : CONTAINS
83 :
84 : ! **************************************************************************************************
85 : !> \brief Allocates and intitializes kg_env
86 : !> \param qs_env ...
87 : !> \param kg_env the object to create
88 : !> \param qs_kind_set ...
89 : !> \param input ...
90 : !> \par History
91 : !> 2012.07 created [Martin Haeufel]
92 : !> \author Martin Haeufel
93 : ! **************************************************************************************************
94 66 : SUBROUTINE kg_env_create(qs_env, kg_env, qs_kind_set, input)
95 : TYPE(qs_environment_type), POINTER :: qs_env
96 : TYPE(kg_environment_type), POINTER :: kg_env
97 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
98 : TYPE(section_vals_type), OPTIONAL, POINTER :: input
99 :
100 66 : ALLOCATE (kg_env)
101 66 : CALL init_kg_env(qs_env, kg_env, qs_kind_set, input)
102 66 : END SUBROUTINE kg_env_create
103 :
104 : ! **************************************************************************************************
105 : !> \brief Initializes kg_env
106 : !> \param qs_env ...
107 : !> \param kg_env ...
108 : !> \param qs_kind_set ...
109 : !> \param input ...
110 : !> \par History
111 : !> 2012.07 created [Martin Haeufel]
112 : !> 2018.01 TNADD correction {JGH]
113 : !> \author Martin Haeufel
114 : ! **************************************************************************************************
115 66 : SUBROUTINE init_kg_env(qs_env, kg_env, qs_kind_set, input)
116 : TYPE(qs_environment_type), POINTER :: qs_env
117 : TYPE(kg_environment_type), POINTER :: kg_env
118 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
119 : TYPE(section_vals_type), OPTIONAL, POINTER :: input
120 :
121 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_kg_env'
122 :
123 : CHARACTER(LEN=10) :: intgrid
124 : INTEGER :: handle, i, iatom, ib, ikind, iunit, n, &
125 : na, natom, nbatch, nkind, np, nr
126 66 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bid
127 : REAL(KIND=dp) :: load, radb, rmax
128 : TYPE(cp_logger_type), POINTER :: logger
129 : TYPE(gto_basis_set_type), POINTER :: lri_aux_basis
130 : TYPE(integration_grid_type), POINTER :: ig_full, ig_mol
131 : TYPE(mp_para_env_type), POINTER :: para_env
132 : TYPE(qs_kind_type), POINTER :: qs_kind
133 : TYPE(section_vals_type), POINTER :: lri_section
134 :
135 66 : CALL timeset(routineN, handle)
136 :
137 66 : CALL cite_reference(Andermatt2016)
138 :
139 66 : NULLIFY (para_env)
140 66 : NULLIFY (kg_env%sab_orb_full)
141 66 : NULLIFY (kg_env%sac_kin)
142 66 : NULLIFY (kg_env%subset_of_mol)
143 66 : NULLIFY (kg_env%subset)
144 66 : NULLIFY (kg_env%tnadd_mat)
145 66 : NULLIFY (kg_env%lri_env)
146 66 : NULLIFY (kg_env%lri_env1)
147 66 : NULLIFY (kg_env%int_grid_atom)
148 66 : NULLIFY (kg_env%int_grid_molecules)
149 66 : NULLIFY (kg_env%int_grid_full)
150 66 : NULLIFY (kg_env%lri_density)
151 66 : NULLIFY (kg_env%lri_rho1)
152 :
153 66 : kg_env%nsubsets = 0
154 :
155 : ! get coloring method settings
156 66 : CALL section_vals_val_get(input, "DFT%KG_METHOD%COLORING_METHOD", i_val=kg_env%coloring_method)
157 : ! get method for nonadditive kinetic energy embedding potential
158 66 : CALL section_vals_val_get(input, "DFT%KG_METHOD%TNADD_METHOD", i_val=kg_env%tnadd_method)
159 : !
160 118 : SELECT CASE (kg_env%tnadd_method)
161 : CASE (kg_tnadd_embed, kg_tnadd_embed_ri)
162 : ! kinetic energy functional
163 52 : kg_env%xc_section_kg => section_vals_get_subs_vals(input, "DFT%KG_METHOD%XC")
164 52 : IF (.NOT. ASSOCIATED(kg_env%xc_section_kg)) THEN
165 : CALL cp_abort(__LOCATION__, &
166 0 : "KG runs require a kinetic energy functional set in &KG_METHOD")
167 : END IF
168 : CASE (kg_tnadd_atomic, kg_tnadd_none)
169 14 : NULLIFY (kg_env%xc_section_kg)
170 : CASE DEFAULT
171 66 : CPABORT("KG:TNADD METHOD")
172 : END SELECT
173 :
174 66 : IF (kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
175 : ! initialize the LRI environment
176 : ! Check if LRI_AUX basis is available
177 2 : rmax = 0.0_dp
178 2 : nkind = SIZE(qs_kind_set)
179 6 : DO ikind = 1, nkind
180 4 : qs_kind => qs_kind_set(ikind)
181 4 : NULLIFY (lri_aux_basis)
182 4 : CALL get_qs_kind(qs_kind, basis_set=lri_aux_basis, basis_type="LRI_AUX")
183 4 : CPASSERT(ASSOCIATED(lri_aux_basis))
184 4 : CALL get_gto_basis_set(gto_basis_set=lri_aux_basis, kind_radius=radb)
185 6 : rmax = MAX(rmax, radb)
186 : END DO
187 2 : rmax = 1.25_dp*rmax
188 2 : lri_section => section_vals_get_subs_vals(input, "DFT%KG_METHOD%LRIGPW")
189 2 : CALL lri_env_init(kg_env%lri_env, lri_section)
190 2 : CALL lri_env_basis("LRI", qs_env, kg_env%lri_env, qs_kind_set)
191 : !
192 : ! if energy correction is performed with force calculation,
193 : ! then response calculation requires
194 : ! perturbation density for kernel calculation
195 2 : CALL lri_env_init(kg_env%lri_env1, lri_section)
196 2 : CALL lri_env_basis("LRI", qs_env, kg_env%lri_env1, qs_kind_set)
197 : !
198 : ! integration grid
199 : !
200 2 : CALL section_vals_val_get(input, "DFT%KG_METHOD%INTEGRATION_GRID", c_val=intgrid)
201 2 : CALL uppercase(intgrid)
202 2 : SELECT CASE (intgrid)
203 : CASE ("SMALL")
204 2 : nr = 50
205 2 : na = 38
206 : CASE ("MEDIUM")
207 0 : nr = 100
208 0 : na = 110
209 : CASE ("LARGE")
210 0 : nr = 200
211 0 : na = 302
212 : CASE ("HUGE")
213 0 : nr = 400
214 0 : na = 590
215 : CASE DEFAULT
216 2 : CPABORT("KG:INTEGRATION_GRID")
217 : END SELECT
218 2 : NULLIFY (logger)
219 2 : logger => cp_get_default_logger()
220 2 : iunit = cp_logger_get_default_io_unit(logger)
221 2 : CALL initialize_atomic_grid(kg_env%int_grid_atom, nr, na, rmax, iunit=iunit)
222 : ! load balancing
223 2 : CALL get_qs_env(qs_env=qs_env, natom=natom, para_env=para_env)
224 2 : np = para_env%num_pe
225 2 : load = REAL(natom, KIND=dp)*kg_env%int_grid_atom%ntot/REAL(np, KIND=dp)
226 : !
227 2 : CALL allocate_intgrid(kg_env%int_grid_full)
228 2 : ig_full => kg_env%int_grid_full
229 2 : CALL allocate_intgrid(kg_env%int_grid_molecules)
230 2 : ig_mol => kg_env%int_grid_molecules
231 2 : nbatch = (natom*kg_env%int_grid_atom%nbatch)/np
232 2 : nbatch = NINT((nbatch + 1)*1.2_dp)
233 6 : ALLOCATE (bid(2, nbatch))
234 12 : nbatch = 0
235 12 : DO iatom = 1, natom
236 892 : DO ib = 1, kg_env%int_grid_atom%nbatch
237 890 : IF (para_env%mepos == MOD(iatom + ib, np)) THEN
238 440 : nbatch = nbatch + 1
239 440 : CPASSERT(nbatch <= SIZE(bid, 2))
240 440 : bid(1, nbatch) = iatom
241 440 : bid(2, nbatch) = ib
242 : END IF
243 : END DO
244 : END DO
245 : !
246 2 : ig_full%nbatch = nbatch
247 446 : ALLOCATE (ig_full%grid_batch(nbatch))
248 : !
249 2 : ig_mol%nbatch = nbatch
250 446 : ALLOCATE (ig_mol%grid_batch(nbatch))
251 : !
252 442 : DO i = 1, nbatch
253 440 : iatom = bid(1, i)
254 440 : ib = bid(2, i)
255 : !
256 440 : ig_full%grid_batch(i)%ref_atom = iatom
257 440 : ig_full%grid_batch(i)%ibatch = ib
258 440 : ig_full%grid_batch(i)%np = kg_env%int_grid_atom%batch(ib)%np
259 440 : ig_full%grid_batch(i)%radius = kg_env%int_grid_atom%batch(ib)%rad
260 3080 : ig_full%grid_batch(i)%rcenter(1:3) = kg_env%int_grid_atom%batch(ib)%rcenter(1:3)
261 440 : n = ig_full%grid_batch(i)%np
262 1320 : ALLOCATE (ig_full%grid_batch(i)%rco(3, n))
263 1320 : ALLOCATE (ig_full%grid_batch(i)%weight(n))
264 1320 : ALLOCATE (ig_full%grid_batch(i)%wref(n))
265 1320 : ALLOCATE (ig_full%grid_batch(i)%wsum(n))
266 19440 : ig_full%grid_batch(i)%weight(:) = kg_env%int_grid_atom%batch(ib)%weight(:)
267 : !
268 440 : ig_mol%grid_batch(i)%ref_atom = iatom
269 440 : ig_mol%grid_batch(i)%ibatch = ib
270 440 : ig_mol%grid_batch(i)%np = kg_env%int_grid_atom%batch(ib)%np
271 440 : ig_mol%grid_batch(i)%radius = kg_env%int_grid_atom%batch(ib)%rad
272 3080 : ig_mol%grid_batch(i)%rcenter(1:3) = kg_env%int_grid_atom%batch(ib)%rcenter(1:3)
273 440 : n = ig_mol%grid_batch(i)%np
274 1320 : ALLOCATE (ig_mol%grid_batch(i)%rco(3, n))
275 1320 : ALLOCATE (ig_mol%grid_batch(i)%weight(n))
276 1320 : ALLOCATE (ig_mol%grid_batch(i)%wref(n))
277 1320 : ALLOCATE (ig_mol%grid_batch(i)%wsum(n))
278 19442 : ig_mol%grid_batch(i)%weight(:) = kg_env%int_grid_atom%batch(ib)%weight(:)
279 : END DO
280 : !
281 2 : DEALLOCATE (bid)
282 : END IF
283 :
284 66 : CALL timestop(handle)
285 :
286 66 : END SUBROUTINE init_kg_env
287 :
288 : ! **************************************************************************************************
289 : !> \brief builds either the full neighborlist or neighborlists of molecular
290 : !> \brief subsets, depending on parameter values
291 : !> \param qs_env ...
292 : !> \param sab_orb the return type, a neighborlist
293 : !> \param sac_kin ...
294 : !> \param molecular if false, the full neighborlist is build
295 : !> \param subset_of_mol the molecular subsets
296 : !> \param current_subset the subset of which the neighborlist is to be build
297 : !> \par History
298 : !> 2012.07 created [Martin Haeufel]
299 : !> \author Martin Haeufel
300 : ! **************************************************************************************************
301 306 : SUBROUTINE kg_build_neighborlist(qs_env, sab_orb, sac_kin, &
302 : molecular, subset_of_mol, current_subset)
303 : TYPE(qs_environment_type), POINTER :: qs_env
304 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
305 : OPTIONAL, POINTER :: sab_orb, sac_kin
306 : LOGICAL, OPTIONAL :: molecular
307 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: subset_of_mol
308 : INTEGER, OPTIONAL :: current_subset
309 :
310 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_build_neighborlist'
311 :
312 : INTEGER :: handle, ikind, nkind
313 : LOGICAL :: mic, molecule_only
314 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: orb_present, tpot_present
315 : REAL(dp) :: subcells
316 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: orb_radius, tpot_radius
317 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
318 306 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
319 : TYPE(cell_type), POINTER :: cell
320 : TYPE(distribution_1d_type), POINTER :: distribution_1d
321 : TYPE(distribution_2d_type), POINTER :: distribution_2d
322 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
323 306 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
324 : TYPE(local_potential_type), POINTER :: tnadd_potential
325 306 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
326 : TYPE(mp_para_env_type), POINTER :: para_env
327 306 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
328 306 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
329 : TYPE(section_vals_type), POINTER :: neighbor_list_section
330 :
331 306 : CALL timeset(routineN, handle)
332 306 : NULLIFY (para_env)
333 :
334 : ! restrict lists to molecular subgroups
335 306 : molecule_only = .FALSE.
336 306 : IF (PRESENT(molecular)) molecule_only = molecular
337 : ! enforce minimum image convention if we use molecules
338 306 : mic = molecule_only
339 :
340 : CALL get_qs_env(qs_env=qs_env, &
341 : atomic_kind_set=atomic_kind_set, &
342 : qs_kind_set=qs_kind_set, &
343 : cell=cell, &
344 : distribution_2d=distribution_2d, &
345 : molecule_set=molecule_set, &
346 : local_particles=distribution_1d, &
347 : particle_set=particle_set, &
348 306 : para_env=para_env)
349 :
350 306 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
351 :
352 : ! Allocate work storage
353 306 : nkind = SIZE(atomic_kind_set)
354 1530 : ALLOCATE (orb_radius(nkind), tpot_radius(nkind))
355 770 : orb_radius(:) = 0.0_dp
356 770 : tpot_radius(:) = 0.0_dp
357 1530 : ALLOCATE (orb_present(nkind), tpot_present(nkind))
358 1224 : ALLOCATE (pair_radius(nkind, nkind))
359 918 : ALLOCATE (atom2d(nkind))
360 :
361 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
362 306 : molecule_set, molecule_only, particle_set=particle_set)
363 :
364 770 : DO ikind = 1, nkind
365 464 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
366 464 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
367 770 : IF (ASSOCIATED(orb_basis_set)) THEN
368 464 : orb_present(ikind) = .TRUE.
369 464 : IF (PRESENT(subset_of_mol)) THEN
370 258 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, kind_radius=orb_radius(ikind))
371 : ELSE
372 206 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, short_kind_radius=orb_radius(ikind))
373 : END IF
374 : ELSE
375 0 : orb_present(ikind) = .FALSE.
376 0 : orb_radius(ikind) = 0.0_dp
377 : END IF
378 : END DO
379 :
380 306 : IF (PRESENT(sab_orb)) THEN
381 : ! Build the orbital-orbital overlap neighbor list
382 284 : CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
383 284 : IF (PRESENT(subset_of_mol)) THEN
384 : CALL build_neighbor_lists(sab_orb, particle_set, atom2d, cell, pair_radius, &
385 : mic=mic, subcells=subcells, molecular=molecule_only, subset_of_mol=subset_of_mol, &
386 170 : current_subset=current_subset, nlname="sab_orb")
387 : ELSE
388 : CALL build_neighbor_lists(sab_orb, particle_set, atom2d, cell, pair_radius, &
389 114 : mic=mic, subcells=subcells, molecular=molecule_only, nlname="sab_orb")
390 : END IF
391 : ! Print out the neighborlist
392 284 : neighbor_list_section => section_vals_get_subs_vals(qs_env%input, "DFT%KG_METHOD%PRINT%NEIGHBOR_LISTS")
393 284 : IF (molecule_only) THEN
394 : CALL write_neighbor_lists(sab_orb, particle_set, cell, para_env, neighbor_list_section, &
395 172 : "/SAB_ORB_MOLECULAR", "sab_orb", "MOLECULAR SUBSET NEIGHBORLIST")
396 : ELSE
397 : CALL write_neighbor_lists(sab_orb, particle_set, cell, para_env, neighbor_list_section, &
398 112 : "/SAB_ORB_FULL", "sab_orb", "FULL NEIGHBORLIST")
399 : END IF
400 : END IF
401 :
402 306 : IF (PRESENT(sac_kin)) THEN
403 56 : DO ikind = 1, nkind
404 34 : tpot_present(ikind) = .FALSE.
405 34 : CALL get_qs_kind(qs_kind_set(ikind), tnadd_potential=tnadd_potential)
406 56 : IF (ASSOCIATED(tnadd_potential)) THEN
407 34 : CALL get_potential(potential=tnadd_potential, radius=tpot_radius(ikind))
408 34 : tpot_present(ikind) = .TRUE.
409 : END IF
410 : END DO
411 22 : CALL pair_radius_setup(orb_present, tpot_present, orb_radius, tpot_radius, pair_radius)
412 : CALL build_neighbor_lists(sac_kin, particle_set, atom2d, cell, pair_radius, &
413 22 : subcells=subcells, operator_type="ABC", nlname="sac_kin")
414 : neighbor_list_section => section_vals_get_subs_vals(qs_env%input, &
415 22 : "DFT%KG_METHOD%PRINT%NEIGHBOR_LISTS")
416 : CALL write_neighbor_lists(sac_kin, particle_set, cell, para_env, neighbor_list_section, &
417 22 : "/SAC_KIN", "sac_kin", "ORBITAL kin energy potential")
418 : END IF
419 :
420 : ! Release work storage
421 306 : CALL atom2d_cleanup(atom2d)
422 306 : DEALLOCATE (atom2d)
423 306 : DEALLOCATE (orb_present, tpot_present)
424 306 : DEALLOCATE (orb_radius, tpot_radius)
425 306 : DEALLOCATE (pair_radius)
426 :
427 306 : CALL timestop(handle)
428 :
429 918 : END SUBROUTINE kg_build_neighborlist
430 :
431 : ! **************************************************************************************************
432 : !> \brief Removes all replicated pairs from a 2d integer buffer array
433 : !> \param pairs_buffer the array, assumed to have the shape (2,:)
434 : !> \param n number of pairs (in), number of disjunct pairs (out)
435 : !> \par History
436 : !> 2012.07 created [Martin Haeufel]
437 : !> 2014.11 simplified [Ole Schuett]
438 : !> \author Martin Haeufel
439 : ! **************************************************************************************************
440 123 : SUBROUTINE kg_remove_duplicates(pairs_buffer, n)
441 : INTEGER(KIND=int_4), DIMENSION(:, :), &
442 : INTENT(INOUT) :: pairs_buffer
443 : INTEGER, INTENT(INOUT) :: n
444 :
445 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_remove_duplicates'
446 :
447 : INTEGER :: handle, i, npairs
448 246 : INTEGER, DIMENSION(n) :: ind
449 246 : INTEGER(KIND=int_8), DIMENSION(n) :: sort_keys
450 246 : INTEGER(KIND=int_4), DIMENSION(2, n) :: work
451 :
452 123 : CALL timeset(routineN, handle)
453 :
454 123 : IF (n > 0) THEN
455 : ! represent a pair of int_4 as a single int_8 number, simplifies sorting.
456 4489 : sort_keys(1:n) = ISHFT(INT(pairs_buffer(1, 1:n), KIND=int_8), 8*int_4_size)
457 4489 : sort_keys(1:n) = sort_keys(1:n) + pairs_buffer(2, 1:n) !upper + lower bytes
458 111 : CALL sort(sort_keys, n, ind)
459 :
460 : ! add first pair, the case npairs==0 was excluded above
461 111 : npairs = 1
462 333 : work(:, 1) = pairs_buffer(:, ind(1))
463 :
464 : ! remove duplicates from the sorted list
465 4378 : DO i = 2, n
466 4378 : IF (sort_keys(i) /= sort_keys(i - 1)) THEN
467 195 : npairs = npairs + 1
468 585 : work(:, npairs) = pairs_buffer(:, ind(i))
469 : END IF
470 : END DO
471 :
472 111 : n = npairs
473 1029 : pairs_buffer(:, :n) = work(:, :n)
474 : END IF
475 :
476 123 : CALL timestop(handle)
477 :
478 123 : END SUBROUTINE kg_remove_duplicates
479 :
480 : ! **************************************************************************************************
481 : !> \brief writes the graph to file using the DIMACS standard format
482 : !> for a definition of the file format see
483 : !> mat.gsia.cmu.edu?COLOR/general/ccformat.ps
484 : !> c comment line
485 : !> p edge NODES EDGES
486 : !> with NODES - number of nodes
487 : !> EDGES - numer of edges
488 : !> e W V
489 : !> ...
490 : !> there is one edge descriptor line for each edge in the graph
491 : !> for an edge (w,v) the fields W and V specify its endpoints
492 : !> \param pairs ...
493 : !> \param nnodes the total number of nodes
494 : ! **************************************************************************************************
495 0 : SUBROUTINE write_to_file(pairs, nnodes)
496 : INTEGER(KIND=int_4), ALLOCATABLE, &
497 : DIMENSION(:, :), INTENT(IN) :: pairs
498 : INTEGER, INTENT(IN) :: nnodes
499 :
500 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_to_file'
501 :
502 : INTEGER :: handle, i, imol, jmol, npairs, unit_nr
503 : INTEGER(KIND=int_4), ALLOCATABLE, DIMENSION(:, :) :: sorted_pairs
504 :
505 0 : CALL timeset(routineN, handle)
506 :
507 : ! get the number of disjunct pairs
508 0 : npairs = SIZE(pairs, 2)
509 :
510 0 : ALLOCATE (sorted_pairs(2, npairs))
511 :
512 : ! reorder pairs such that pairs(1,*) < pairs(2,*)
513 0 : DO i = 1, npairs
514 : ! get molecular ids
515 0 : imol = pairs(1, i)
516 0 : jmol = pairs(2, i)
517 0 : IF (imol > jmol) THEN
518 : ! switch pair and store
519 0 : sorted_pairs(1, i) = jmol
520 0 : sorted_pairs(2, i) = imol
521 : ELSE
522 : ! keep ordering just copy
523 0 : sorted_pairs(1, i) = imol
524 0 : sorted_pairs(2, i) = jmol
525 : END IF
526 : END DO
527 :
528 : ! remove duplicates and get the number of disjunct pairs (number of edges)
529 0 : CALL kg_remove_duplicates(sorted_pairs, npairs)
530 :
531 : ! should now be half as much pairs as before
532 0 : CPASSERT(npairs == SIZE(pairs, 2)/2)
533 :
534 0 : CALL open_file(unit_number=unit_nr, file_name="graph.col")
535 :
536 0 : WRITE (unit_nr, '(A6,1X,I8,1X,I8)') "p edge", nnodes, npairs
537 :
538 : ! only write out the first npairs entries
539 0 : DO i = 1, npairs
540 0 : WRITE (unit_nr, '(A1,1X,I8,1X,I8)') "e", sorted_pairs(1, i), sorted_pairs(2, i)
541 : END DO
542 :
543 0 : CALL close_file(unit_nr)
544 :
545 0 : DEALLOCATE (sorted_pairs)
546 :
547 0 : CALL timestop(handle)
548 :
549 0 : END SUBROUTINE write_to_file
550 :
551 : ! **************************************************************************************************
552 : !> \brief ...
553 : !> \param kg_env ...
554 : !> \param para_env ...
555 : ! **************************************************************************************************
556 82 : SUBROUTINE kg_build_subsets(kg_env, para_env)
557 : TYPE(kg_environment_type), POINTER :: kg_env
558 : TYPE(mp_para_env_type), POINTER :: para_env
559 :
560 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_build_subsets'
561 :
562 : INTEGER :: color, handle, i, iatom, imol, isub, &
563 : jatom, jmol, nmol, npairs, npairs_local
564 : INTEGER(KIND=int_4) :: ncolors
565 82 : INTEGER(KIND=int_4), ALLOCATABLE, DIMENSION(:) :: color_of_node
566 82 : INTEGER(KIND=int_4), ALLOCATABLE, DIMENSION(:, :) :: msg_gather, pairs, pairs_buffer
567 82 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nnodes_of_color
568 : TYPE(neighbor_list_iterator_p_type), &
569 82 : DIMENSION(:), POINTER :: nl_iterator
570 :
571 82 : CALL timeset(routineN, handle)
572 :
573 : ! first: get a (local) list of pairs from the (local) neighbor list data
574 82 : nmol = SIZE(kg_env%molecule_set)
575 :
576 82 : npairs = 0
577 82 : CALL neighbor_list_iterator_create(nl_iterator, kg_env%sab_orb_full)
578 4976 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
579 4894 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom)
580 :
581 4894 : imol = kg_env%atom_to_molecule(iatom)
582 4894 : jmol = kg_env%atom_to_molecule(jatom)
583 :
584 : !IF (imol<jmol) THEN
585 4976 : IF (imol .NE. jmol) THEN
586 :
587 2087 : npairs = npairs + 2
588 :
589 : END IF
590 :
591 : END DO
592 82 : CALL neighbor_list_iterator_release(nl_iterator)
593 :
594 238 : ALLOCATE (pairs_buffer(2, npairs))
595 :
596 82 : npairs = 0
597 82 : CALL neighbor_list_iterator_create(nl_iterator, kg_env%sab_orb_full)
598 4976 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
599 4894 : CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom)
600 :
601 4894 : imol = kg_env%atom_to_molecule(iatom)
602 4894 : jmol = kg_env%atom_to_molecule(jatom)
603 :
604 4976 : IF (imol .NE. jmol) THEN
605 :
606 : ! add pair to the local list
607 :
608 : ! add both orderings - makes it easier to build the neighborlist
609 2087 : npairs = npairs + 1
610 :
611 2087 : pairs_buffer(1, npairs) = imol
612 2087 : pairs_buffer(2, npairs) = jmol
613 :
614 2087 : npairs = npairs + 1
615 :
616 2087 : pairs_buffer(2, npairs) = imol
617 2087 : pairs_buffer(1, npairs) = jmol
618 :
619 : END IF
620 :
621 : END DO
622 82 : CALL neighbor_list_iterator_release(nl_iterator)
623 :
624 : ! remove duplicates
625 82 : CALL kg_remove_duplicates(pairs_buffer, npairs)
626 :
627 : ! get the maximum number of local pairs on all nodes (size of the mssg)
628 : ! remember how many pairs we have local
629 82 : npairs_local = npairs
630 82 : CALL para_env%max(npairs)
631 :
632 : ! allocate message
633 238 : ALLOCATE (pairs(2, npairs))
634 :
635 694 : pairs(:, 1:npairs_local) = pairs_buffer(:, 1:npairs_local)
636 82 : pairs(:, npairs_local + 1:) = 0
637 :
638 82 : DEALLOCATE (pairs_buffer)
639 :
640 : ! second: gather all data on the master node
641 : ! better would be a scheme where duplicates are removed in a tree-like reduction scheme.
642 : ! this step can be needlessly memory intensive in the current implementation.
643 :
644 82 : IF (para_env%is_source()) THEN
645 119 : ALLOCATE (msg_gather(2, npairs*para_env%num_pe))
646 : ELSE
647 41 : ALLOCATE (msg_gather(2, 1))
648 : END IF
649 :
650 817 : msg_gather = 0
651 :
652 82 : CALL para_env%gather(pairs, msg_gather)
653 :
654 82 : DEALLOCATE (pairs)
655 :
656 82 : IF (para_env%is_source()) THEN
657 :
658 : ! shift all non-zero entries to the beginning of the array and count the number of actual pairs
659 41 : npairs = 0
660 :
661 245 : DO i = 1, SIZE(msg_gather, 2)
662 245 : IF (msg_gather(1, i) .NE. 0) THEN
663 204 : npairs = npairs + 1
664 612 : msg_gather(:, npairs) = msg_gather(:, i)
665 : END IF
666 : END DO
667 :
668 : ! remove duplicates
669 41 : CALL kg_remove_duplicates(msg_gather, npairs)
670 :
671 119 : ALLOCATE (pairs(2, npairs))
672 :
673 347 : pairs(:, 1:npairs) = msg_gather(:, 1:npairs)
674 :
675 41 : DEALLOCATE (msg_gather)
676 :
677 : !WRITE(*,'(A48,5X,I10,4X,A2,1X,I10)') " KG| Total number of overlapping molecular pairs",npairs/2,"of",nmol*(nmol-1)/2
678 :
679 : ! write to file, nnodes = number of molecules
680 : IF (.FALSE.) THEN
681 : CALL write_to_file(pairs, SIZE(kg_env%molecule_set))
682 : END IF
683 :
684 : ! vertex coloring algorithm
685 41 : CALL kg_vertex_coloring(kg_env, pairs, ncolors, color_of_node)
686 :
687 41 : DEALLOCATE (pairs)
688 :
689 : ELSE
690 :
691 41 : DEALLOCATE (msg_gather)
692 :
693 : END IF
694 :
695 : !WRITE(*,'(A27,40X,I6,1X,A6)') ' KG| Vertex coloring result', ncolors, 'colors'
696 :
697 : ! broadcast the number of colors to all nodes
698 82 : CALL para_env%bcast(ncolors)
699 :
700 164 : IF (.NOT. ALLOCATED(color_of_node)) ALLOCATE (color_of_node(nmol))
701 :
702 : ! broadcast the resulting coloring to all nodes.....
703 82 : CALL para_env%bcast(color_of_node)
704 :
705 82 : IF ((kg_env%nsubsets .NE. 0) .AND. (ncolors .NE. kg_env%nsubsets)) THEN
706 : ! number of subsets has changed
707 :
708 : ! deallocate stuff if necessary
709 0 : IF (ASSOCIATED(kg_env%subset)) THEN
710 0 : DO isub = 1, kg_env%nsubsets
711 0 : CALL release_neighbor_list_sets(kg_env%subset(isub)%sab_orb)
712 0 : CALL deallocate_task_list(kg_env%subset(isub)%task_list)
713 : END DO
714 0 : DEALLOCATE (kg_env%subset)
715 0 : NULLIFY (kg_env%subset)
716 : END IF
717 :
718 : END IF
719 :
720 : ! allocate and nullify some stuff
721 82 : IF (.NOT. ASSOCIATED(kg_env%subset)) THEN
722 :
723 150 : ALLOCATE (kg_env%subset(ncolors))
724 :
725 156 : DO i = 1, ncolors
726 106 : NULLIFY (kg_env%subset(i)%sab_orb)
727 156 : NULLIFY (kg_env%subset(i)%task_list)
728 : END DO
729 : END IF
730 :
731 : ! set the number of subsets
732 82 : kg_env%nsubsets = ncolors
733 :
734 : ! counting loop
735 246 : ALLOCATE (nnodes_of_color(ncolors))
736 252 : nnodes_of_color = 0
737 268 : DO i = 1, nmol ! nmol=nnodes
738 186 : color = color_of_node(i)
739 186 : kg_env%subset_of_mol(i) = color
740 268 : nnodes_of_color(color) = nnodes_of_color(color) + 1
741 : END DO
742 :
743 82 : DEALLOCATE (nnodes_of_color)
744 82 : DEALLOCATE (color_of_node)
745 :
746 82 : CALL timestop(handle)
747 :
748 164 : END SUBROUTINE kg_build_subsets
749 :
750 : END MODULE kg_environment
|
