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