1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Build up the plane wave density by collocating the primitive Gaussian 8!> functions (pgf). 9!> \par History 10!> Joost VandeVondele (02.2002) 11!> 1) rewrote collocate_pgf for increased accuracy and speed 12!> 2) collocate_core hack for PGI compiler 13!> 3) added multiple grid feature 14!> 4) new way to go over the grid 15!> Joost VandeVondele (05.2002) 16!> 1) prelim. introduction of the real space grid type 17!> JGH [30.08.02] multigrid arrays independent from potential 18!> JGH [17.07.03] distributed real space code 19!> JGH [23.11.03] refactoring and new loop ordering 20!> JGH [04.12.03] OpneMP parallelization of main loops 21!> Joost VandeVondele (12.2003) 22!> 1) modified to compute tau 23!> Joost removed incremental build feature 24!> Joost introduced map consistent 25!> Rewrote grid integration/collocation routines, [Joost VandeVondele,03.2007] 26!> JGH [26.06.15] modification to allow for k-points 27!> \author Matthias Krack (03.04.2001) 28! ************************************************************************************************** 29MODULE qs_integrate_potential_product 30 USE admm_types, ONLY: admm_type 31 USE atomic_kind_types, ONLY: atomic_kind_type,& 32 get_atomic_kind_set 33 USE basis_set_types, ONLY: get_gto_basis_set,& 34 gto_basis_set_type 35 USE cell_types, ONLY: cell_type,& 36 pbc 37 USE cp_control_types, ONLY: dft_control_type 38 USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set 39 USE cube_utils, ONLY: cube_info_type 40 USE dbcsr_api, ONLY: & 41 dbcsr_add_block_node, dbcsr_copy, dbcsr_create, dbcsr_distribution_get, & 42 dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, dbcsr_p_type, & 43 dbcsr_type, dbcsr_work_create 44 USE gaussian_gridlevels, ONLY: gridlevel_info_type 45 USE input_constants, ONLY: do_admm_exch_scaling_merlot 46 USE kinds, ONLY: default_string_length,& 47 dp,& 48 int_8 49 USE memory_utilities, ONLY: reallocate 50 USE orbital_pointers, ONLY: ncoset 51 USE particle_types, ONLY: particle_type 52 USE pw_env_types, ONLY: pw_env_get,& 53 pw_env_type 54 USE pw_types, ONLY: pw_p_type 55 USE qs_environment_types, ONLY: get_qs_env,& 56 qs_environment_type 57 USE qs_force_types, ONLY: qs_force_type 58 USE qs_integrate_potential_low, ONLY: integrate_pgf_product_rspace 59 USE qs_kind_types, ONLY: get_qs_kind,& 60 get_qs_kind_set,& 61 qs_kind_type 62 USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,& 63 realspace_grid_p_type,& 64 rs_grid_release,& 65 rs_grid_retain 66 USE rs_pw_interface, ONLY: potential_pw2rs 67 USE task_list_methods, ONLY: int2pair,& 68 rs_distribute_matrix 69 USE task_list_types, ONLY: task_list_type 70 USE virial_types, ONLY: virial_type 71 72!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads 73 74#include "./base/base_uses.f90" 75 76 IMPLICIT NONE 77 78 PRIVATE 79 80 INTEGER :: debug_count = 0 81 82 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE. 83 84 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_integrate_potential_product' 85 86! *** Public subroutines *** 87! *** Don't include this routines directly, use the interface to 88! *** qs_integrate_potential 89 90 PUBLIC :: integrate_v_rspace 91 92CONTAINS 93 94! ************************************************************************************************** 95!> \brief computes matrix elements corresponding to a given potential 96!> \param v_rspace ... 97!> \param hmat ... 98!> \param hmat_kp ... 99!> \param pmat ... 100!> \param pmat_kp ... 101!> \param qs_env ... 102!> \param calculate_forces ... 103!> \param force_adm whether force of in aux. dens. matrix is calculated 104!> \param ispin ... 105!> \param compute_tau ... 106!> \param gapw ... 107!> \param basis_type ... 108!> \param pw_env_external ... 109!> \param task_list_external ... 110!> \par History 111!> IAB (29-Apr-2010): Added OpenMP parallelisation to task loop 112!> (c) The Numerical Algorithms Group (NAG) Ltd, 2010 on behalf of the HECToR project 113!> Some refactoring, get priorities for options correct (JGH, 04.2014) 114!> Added options to allow for k-points 115!> For a smooth transition we allow for old and new (vector) matrices (hmat, pmat) (JGH, 06.2015) 116!> \note 117!> integrates a given potential (or other object on a real 118!> space grid) = v_rspace using a multi grid technique (mgrid_*) 119!> over the basis set producing a number for every element of h 120!> (should have the same sparsity structure of S) 121!> additional screening is available using the magnitude of the 122!> elements in p (? I'm not sure this is a very good idea) 123!> this argument is optional 124!> derivatives of these matrix elements with respect to the ionic 125!> coordinates can be computed as well 126! ************************************************************************************************** 127 SUBROUTINE integrate_v_rspace(v_rspace, hmat, hmat_kp, pmat, pmat_kp, & 128 qs_env, calculate_forces, force_adm, ispin, & 129 compute_tau, gapw, basis_type, pw_env_external, task_list_external) 130 131 TYPE(pw_p_type) :: v_rspace 132 TYPE(dbcsr_p_type), INTENT(INOUT), OPTIONAL :: hmat 133 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & 134 POINTER :: hmat_kp 135 TYPE(dbcsr_p_type), INTENT(IN), OPTIONAL :: pmat 136 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & 137 POINTER :: pmat_kp 138 TYPE(qs_environment_type), POINTER :: qs_env 139 LOGICAL, INTENT(IN) :: calculate_forces 140 LOGICAL, INTENT(IN), OPTIONAL :: force_adm 141 INTEGER, INTENT(IN), OPTIONAL :: ispin 142 LOGICAL, INTENT(IN), OPTIONAL :: compute_tau, gapw 143 CHARACTER(len=*), INTENT(IN), OPTIONAL :: basis_type 144 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external 145 TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external 146 147 CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace', & 148 routineP = moduleN//':'//routineN 149 150 CHARACTER(len=default_string_length) :: my_basis_type 151 INTEGER :: atom_a, atom_b, bcol, brow, handle, i, iatom, igrid_level, ikind, ikind_old, & 152 ilevel, img, ipair, ipgf, ipgf_new, iset, iset_new, iset_old, itask, ithread, jatom, & 153 jkind, jkind_old, jpgf, jpgf_new, jset, jset_new, jset_old, maxco, maxpgf, maxset, & 154 maxsgf_set, na1, na2, natom, nb1, nb2, ncoa, ncob, nimages, nkind, nseta, nsetb, nthread, & 155 offs_dv, sgfa, sgfb 156 INTEGER(KIND=int_8), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send 157 INTEGER(kind=int_8), DIMENSION(:, :), POINTER :: tasks 158 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 159 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: block_touched 160 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 161 npgfb, nsgfa, nsgfb 162 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 163 LOGICAL :: atom_pair_changed, atom_pair_done, distributed_grids, do_kp, found, h_duplicated, & 164 has_threads, map_consistent, my_compute_tau, my_force_adm, my_gapw, new_set_pair_coming, & 165 p_duplicated, pab_required, scatter, use_subpatch, use_virial 166 REAL(KIND=dp) :: admm_scal_fac, dab, eps_gvg_rspace, & 167 rab2, scalef, zetp 168 REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra, rab, rab_inv, rb 169 REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b 170 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 171 REAL(KIND=dp), DIMENSION(:, :), POINTER :: dist_ab, h_block, hab, p_block, pab, & 172 rpgfa, rpgfb, sphi_a, sphi_b, work, & 173 zeta, zetb 174 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: habt, hadb, hdab, pabt, workt 175 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: hadbt, hdabt 176 TYPE(admm_type), POINTER :: admm_env 177 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 178 TYPE(cell_type), POINTER :: cell 179 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info 180 TYPE(dbcsr_distribution_type) :: dist 181 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: deltap, dhmat, htemp 182 TYPE(dbcsr_type), POINTER :: href 183 TYPE(dft_control_type), POINTER :: dft_control 184 TYPE(gridlevel_info_type), POINTER :: gridlevel_info 185 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 186 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 187 TYPE(pw_env_type), POINTER :: pw_env 188 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 189 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 190 TYPE(realspace_grid_desc_p_type), DIMENSION(:), & 191 POINTER :: rs_descs 192 TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: rs_v 193 TYPE(task_list_type), POINTER :: task_list, task_list_soft 194 TYPE(virial_type), POINTER :: virial 195 196 CALL timeset(routineN, handle) 197 198 ! we test here if the provided operator matrices are consistent 199 CPASSERT(PRESENT(hmat) .OR. PRESENT(hmat_kp)) 200 do_kp = .FALSE. 201 IF (PRESENT(hmat_kp)) do_kp = .TRUE. 202 IF (PRESENT(pmat)) THEN 203 CPASSERT(PRESENT(hmat)) 204 ELSE IF (PRESENT(pmat_kp)) THEN 205 CPASSERT(PRESENT(hmat_kp)) 206 END IF 207 208 NULLIFY (pw_env, rs_descs, tasks, dist_ab, admm_env) 209 210 debug_count = debug_count + 1 211 212 offs_dv = 0 213 214 ! this routine works in two modes: 215 ! normal mode : <a| V | b> 216 ! tau mode : < nabla a| V | nabla b> 217 my_compute_tau = .FALSE. 218 IF (PRESENT(compute_tau)) my_compute_tau = compute_tau 219 220 my_force_adm = .FALSE. 221 IF (PRESENT(force_adm)) my_force_adm = force_adm 222 223 ! this sets the basis set to be used. GAPW(==soft basis) overwrites basis_type 224 ! default is "ORB" 225 my_gapw = .FALSE. 226 IF (PRESENT(gapw)) my_gapw = gapw 227 IF (PRESENT(basis_type)) THEN 228 my_basis_type = basis_type 229 ELSE 230 my_basis_type = "ORB" 231 END IF 232 233 ! get the task lists 234 ! task lists have to be in sync with basis sets 235 ! there is an option to provide the task list from outside (not through qs_env) 236 ! outside option has highest priority 237 SELECT CASE (my_basis_type) 238 CASE ("ORB") 239 CALL get_qs_env(qs_env=qs_env, & 240 task_list=task_list, & 241 task_list_soft=task_list_soft) 242 CASE ("AUX_FIT") 243 CALL get_qs_env(qs_env=qs_env, & 244 task_list_aux_fit=task_list, & 245 task_list_soft=task_list_soft) 246 END SELECT 247 IF (my_gapw) task_list => task_list_soft 248 IF (PRESENT(task_list_external)) task_list => task_list_external 249 CPASSERT(ASSOCIATED(task_list)) 250 251 ! the information on the grids is provided through pw_env 252 ! pw_env has to be the parent env for the potential grid (input) 253 ! there is an option to provide an external grid 254 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 255 IF (PRESENT(pw_env_external)) pw_env => pw_env_external 256 257 ! get all the general information on the system we are working on 258 CALL get_qs_env(qs_env=qs_env, & 259 atomic_kind_set=atomic_kind_set, & 260 qs_kind_set=qs_kind_set, & 261 cell=cell, & 262 dft_control=dft_control, & 263 particle_set=particle_set, & 264 force=force, & 265 virial=virial) 266 267 admm_scal_fac = 1.0_dp 268 IF (my_force_adm) THEN 269 CALL get_qs_env(qs_env=qs_env, admm_env=admm_env) 270 ! Calculate bare scaling of force according to Merlot, 1. IF: ADMMP, 2. IF: ADMMS, 271 IF ((.NOT. admm_env%charge_constrain) .AND. & 272 (admm_env%scaling_model == do_admm_exch_scaling_merlot)) THEN 273 admm_scal_fac = admm_env%gsi(ispin)**2 274 ELSE IF (admm_env%charge_constrain .AND. & 275 (admm_env%scaling_model == do_admm_exch_scaling_merlot)) THEN 276 admm_scal_fac = (admm_env%gsi(ispin))**(2.0_dp/3.0_dp) 277 END IF 278 END IF 279 280 ! short cuts to task list variables 281 tasks => task_list%tasks 282 dist_ab => task_list%dist_ab 283 atom_pair_send => task_list%atom_pair_send 284 atom_pair_recv => task_list%atom_pair_recv 285 286 CPASSERT(ASSOCIATED(pw_env)) 287 CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_v) 288 DO i = 1, SIZE(rs_v) 289 CALL rs_grid_retain(rs_v(i)%rs_grid) 290 END DO 291 292 ! assign from pw_env 293 gridlevel_info => pw_env%gridlevel_info 294 cube_info => pw_env%cube_info 295 296 ! transform the potential on the rs_multigrids 297 CALL potential_pw2rs(rs_v, v_rspace, pw_env) 298 299 nimages = dft_control%nimages 300 IF (nimages > 1) THEN 301 CPASSERT(do_kp) 302 END IF 303 nkind = SIZE(qs_kind_set) 304 natom = SIZE(particle_set) 305 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 306 307 IF (calculate_forces) THEN 308 ALLOCATE (atom_of_kind(natom)) 309 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind) 310 END IF 311 312 map_consistent = dft_control%qs_control%map_consistent 313 IF (map_consistent) THEN 314 ! needs to be consistent with rho_rspace 315 eps_gvg_rspace = dft_control%qs_control%eps_rho_rspace 316 ELSE 317 eps_gvg_rspace = dft_control%qs_control%eps_gvg_rspace 318 ENDIF 319 320 pab_required = (PRESENT(pmat) .OR. PRESENT(pmat_kp)) & 321 .AND. (calculate_forces .OR. .NOT. map_consistent) 322 323 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, & 324 maxco=maxco, & 325 maxsgf_set=maxsgf_set, & 326 basis_type=my_basis_type) 327 328 distributed_grids = .FALSE. 329 DO igrid_level = 1, gridlevel_info%ngrid_levels 330 IF (rs_v(igrid_level)%rs_grid%desc%distributed) THEN 331 distributed_grids = .TRUE. 332 ENDIF 333 ENDDO 334 335 ! initialize the working hmat structures 336 h_duplicated = .FALSE. 337 ALLOCATE (dhmat(nimages)) 338 IF (do_kp) THEN 339 DO img = 1, nimages 340 dhmat(img)%matrix => hmat_kp(img)%matrix 341 END DO 342 ELSE 343 dhmat(1)%matrix => hmat%matrix 344 END IF 345 IF (distributed_grids) THEN 346 h_duplicated = .TRUE. 347 href => dhmat(1)%matrix 348 DO img = 1, nimages 349 NULLIFY (dhmat(img)%matrix) 350 ALLOCATE (dhmat(img)%matrix) 351 CALL dbcsr_create(dhmat(img)%matrix, template=href, name='LocalH') 352 END DO 353 END IF 354 355 p_duplicated = .FALSE. 356 IF (pab_required) THEN 357 ! initialize the working pmat structures 358 ALLOCATE (deltap(nimages)) 359 IF (do_kp) THEN 360 DO img = 1, nimages 361 deltap(img)%matrix => pmat_kp(img)%matrix 362 END DO 363 ELSE 364 deltap(1)%matrix => pmat%matrix 365 END IF 366 IF (distributed_grids) THEN 367 p_duplicated = .TRUE. 368 DO img = 1, nimages 369 NULLIFY (deltap(img)%matrix) 370 ALLOCATE (deltap(img)%matrix) 371 END DO 372 IF (do_kp) THEN 373 DO img = 1, nimages 374 CALL dbcsr_copy(deltap(img)%matrix, pmat_kp(img)%matrix, name="LocalP") 375 END DO 376 ELSE 377 CALL dbcsr_copy(deltap(1)%matrix, pmat%matrix, name="LocalP") 378 END IF 379 END IF 380 END IF 381 382 nthread = 1 383!$ nthread = omp_get_max_threads() 384 385 ! *** Allocate work storage *** 386 NULLIFY (pabt, habt, workt) 387 CALL reallocate(habt, 1, maxco, 1, maxco, 0, nthread) 388 CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread) 389 IF (pab_required) THEN 390 CALL reallocate(pabt, 1, maxco, 1, maxco, 0, nthread) 391 ELSE 392 CPASSERT(.NOT. calculate_forces) 393 END IF 394 395 NULLIFY (hdabt, hadbt, hdab, hadb) 396 397 ! get maximum numbers 398 natom = SIZE(particle_set) 399 maxset = 0 400 maxpgf = 0 401 DO ikind = 1, nkind 402 CALL get_qs_kind(qs_kind_set(ikind), & 403 softb=my_gapw, & 404 basis_set=orb_basis_set, basis_type=my_basis_type) 405 406 IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE 407 408 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 409 npgf=npgfa, nset=nseta) 410 411 maxset = MAX(nseta, maxset) 412 maxpgf = MAX(MAXVAL(npgfa), maxpgf) 413 END DO 414 415 IF (distributed_grids .AND. pab_required) THEN 416 CALL rs_distribute_matrix(rs_descs, deltap, atom_pair_send, atom_pair_recv, & 417 natom, nimages, scatter=.TRUE.) 418 ENDIF 419 420 IF (debug_this_module) THEN 421 ALLOCATE (block_touched(natom, natom)) 422 END IF 423 424!$OMP PARALLEL DEFAULT(NONE), & 425!$OMP SHARED(workt,habt,hdabt,hadbt,pabt,tasks,particle_set,natom,maxset), & 426!$OMP SHARED(maxpgf,my_basis_type,my_gapw,dhmat,deltap,use_virial,admm_scal_fac), & 427!$OMP SHARED(pab_required,calculate_forces,ncoset,rs_v,cube_info,my_compute_tau), & 428!$OMP SHARED(map_consistent,eps_gvg_rspace,force,virial,cell,atom_of_kind,dist_ab), & 429!$OMP SHARED(gridlevel_info,task_list,block_touched,nthread,qs_kind_set), & 430!$OMP SHARED(nimages,do_kp), & 431!$OMP PRIVATE(ithread,work,hab,hdab,hadb,pab,iset_old,jset_old), & 432!$OMP PRIVATE(ikind_old,jkind_old,iatom,jatom,iset,jset,ikind,jkind,ilevel,ipgf,jpgf), & 433!$OMP PRIVATE(img,brow,bcol,orb_basis_set,first_sgfa,la_max,la_min,npgfa,nseta,nsgfa), & 434!$OMP PRIVATE(rpgfa,set_radius_a,sphi_a,zeta,first_sgfb,lb_max,lb_min,npgfb), & 435!$OMP PRIVATE(nsetb,nsgfb,rpgfb,set_radius_b,sphi_b,zetb,found,atom_a,atom_b), & 436!$OMP PRIVATE(force_a,force_b,my_virial_a,my_virial_b,atom_pair_changed,h_block), & 437!$OMP PRIVATE(p_block,ncoa,sgfa,ncob,sgfb,rab,rab2,ra,rb,zetp,dab,igrid_level), & 438!$OMP PRIVATE(na1,na2,nb1,nb2,use_subpatch,rab_inv,new_set_pair_coming,atom_pair_done), & 439!$OMP PRIVATE(iset_new,jset_new,ipgf_new,jpgf_new,scalef), & 440!$OMP PRIVATE(itask,dist,has_threads) 441 442 ithread = 0 443!$ ithread = omp_get_thread_num() 444 work => workt(:, :, ithread) 445 hab => habt(:, :, ithread) 446 IF (pab_required) THEN 447 pab => pabt(:, :, ithread) 448 END IF 449 450 iset_old = -1; jset_old = -1 451 ikind_old = -1; jkind_old = -1 452 453 ! Here we loop over gridlevels first, finalising the matrix after each grid level is 454 ! completed. On each grid level, we loop over atom pairs, which will only access 455 ! a single block of each matrix, so with OpenMP, each matrix block is only touched 456 ! by a single thread for each grid level 457 loop_gridlevels: DO igrid_level = 1, gridlevel_info%ngrid_levels 458 459 DO img = 1, nimages 460 CALL dbcsr_work_create(dhmat(img)%matrix, work_mutable=.TRUE., n=nthread) 461 CALL dbcsr_get_info(dhmat(img)%matrix, distribution=dist) 462 CALL dbcsr_distribution_get(dist, has_threads=has_threads) 463!$ IF (.NOT. has_threads) & 464!$ CPABORT("No thread distribution defined.") 465 END DO 466!$OMP BARRIER 467 468 IF (debug_this_module) THEN 469!$OMP SINGLE 470 block_touched = -1 471!$OMP END SINGLE 472!$OMP FLUSH 473 END IF 474 475!$OMP DO schedule (dynamic, MAX(1,task_list%npairs(igrid_level)/(nthread*50))) 476 loop_pairs: DO ipair = 1, task_list%npairs(igrid_level) 477 loop_tasks: DO itask = task_list%taskstart(ipair, igrid_level), task_list%taskstop(ipair, igrid_level) 478 479 CALL int2pair(tasks(3, itask), ilevel, img, iatom, jatom, iset, jset, ipgf, jpgf, & 480 nimages, natom, maxset, maxpgf) 481 CPASSERT(img == 1 .OR. do_kp) 482 483 ! At the start of a block of tasks, get atom data (and kind data, if needed) 484 IF (itask .EQ. task_list%taskstart(ipair, igrid_level)) THEN 485 486 ikind = particle_set(iatom)%atomic_kind%kind_number 487 jkind = particle_set(jatom)%atomic_kind%kind_number 488 489 ra(:) = pbc(particle_set(iatom)%r, cell) 490 491 IF (iatom <= jatom) THEN 492 brow = iatom 493 bcol = jatom 494 ELSE 495 brow = jatom 496 bcol = iatom 497 END IF 498 499 IF (ikind .NE. ikind_old) THEN 500 CALL get_qs_kind(qs_kind_set(ikind), & 501 softb=my_gapw, & 502 basis_set=orb_basis_set, basis_type=my_basis_type) 503 504 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 505 first_sgf=first_sgfa, & 506 lmax=la_max, & 507 lmin=la_min, & 508 npgf=npgfa, & 509 nset=nseta, & 510 nsgf_set=nsgfa, & 511 pgf_radius=rpgfa, & 512 set_radius=set_radius_a, & 513 sphi=sphi_a, & 514 zet=zeta) 515 ENDIF 516 517 IF (jkind .NE. jkind_old) THEN 518 CALL get_qs_kind(qs_kind_set(jkind), & 519 softb=my_gapw, & 520 basis_set=orb_basis_set, basis_type=my_basis_type) 521 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 522 first_sgf=first_sgfb, & 523 lmax=lb_max, & 524 lmin=lb_min, & 525 npgf=npgfb, & 526 nset=nsetb, & 527 nsgf_set=nsgfb, & 528 pgf_radius=rpgfb, & 529 set_radius=set_radius_b, & 530 sphi=sphi_b, & 531 zet=zetb) 532 533 ENDIF 534 535 IF (debug_this_module) THEN 536!$OMP CRITICAL (block_touched_critical) 537 IF ((block_touched(brow, bcol) .NE. ithread) .AND. (block_touched(brow, bcol) .NE. -1)) THEN 538 CPABORT("Block has been modified by another thread") 539 END IF 540 block_touched(brow, bcol) = ithread 541!$OMP END CRITICAL (block_touched_critical) 542 END IF 543 544 NULLIFY (h_block) 545 CALL dbcsr_get_block_p(dhmat(img)%matrix, brow, bcol, h_block, found) 546 IF (.NOT. ASSOCIATED(h_block)) THEN 547 CALL dbcsr_add_block_node(dhmat(img)%matrix, brow, bcol, h_block) 548 END IF 549 550 IF (pab_required) THEN 551 CALL dbcsr_get_block_p(matrix=deltap(img)%matrix, & 552 row=brow, col=bcol, BLOCK=p_block, found=found) 553 CPASSERT(found) 554 END IF 555 556 IF (calculate_forces) THEN 557 atom_a = atom_of_kind(iatom) 558 atom_b = atom_of_kind(jatom) 559 force_a(:) = 0.0_dp 560 force_b(:) = 0.0_dp 561 ENDIF 562 IF (use_virial) THEN 563 my_virial_a = 0.0_dp 564 my_virial_b = 0.0_dp 565 ENDIF 566 567 ikind_old = ikind 568 jkind_old = jkind 569 570 atom_pair_changed = .TRUE. 571 572 ELSE 573 574 atom_pair_changed = .FALSE. 575 576 ENDIF 577 578 IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN 579 580 ncoa = npgfa(iset)*ncoset(la_max(iset)) 581 sgfa = first_sgfa(1, iset) 582 ncob = npgfb(jset)*ncoset(lb_max(jset)) 583 sgfb = first_sgfb(1, jset) 584 IF (pab_required) THEN 585 IF (iatom <= jatom) THEN 586 CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), & 587 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 588 p_block(sgfa, sgfb), SIZE(p_block, 1), & 589 0.0_dp, work(1, 1), SIZE(work, 1)) 590 CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), & 591 1.0_dp, work(1, 1), SIZE(work, 1), & 592 sphi_b(1, sgfb), SIZE(sphi_b, 1), & 593 0.0_dp, pab(1, 1), SIZE(pab, 1)) 594 ELSE 595 CALL dgemm("N", "N", ncob, nsgfa(iset), nsgfb(jset), & 596 1.0_dp, sphi_b(1, sgfb), SIZE(sphi_b, 1), & 597 p_block(sgfb, sgfa), SIZE(p_block, 1), & 598 0.0_dp, work(1, 1), SIZE(work, 1)) 599 CALL dgemm("N", "T", ncob, ncoa, nsgfa(iset), & 600 1.0_dp, work(1, 1), SIZE(work, 1), & 601 sphi_a(1, sgfa), SIZE(sphi_a, 1), & 602 0.0_dp, pab(1, 1), SIZE(pab, 1)) 603 END IF 604 END IF 605 606 IF (iatom <= jatom) THEN 607 hab(1:ncoa, 1:ncob) = 0._dp 608 ELSE 609 hab(1:ncob, 1:ncoa) = 0._dp 610 ENDIF 611 612 iset_old = iset 613 jset_old = jset 614 615 ENDIF 616 617 rab(1) = dist_ab(1, itask) 618 rab(2) = dist_ab(2, itask) 619 rab(3) = dist_ab(3, itask) 620 rab2 = DOT_PRODUCT(rab, rab) 621 rb(1) = ra(1) + rab(1) 622 rb(2) = ra(2) + rab(2) 623 rb(3) = ra(3) + rab(3) 624 zetp = zeta(ipgf, iset) + zetb(jpgf, jset) 625 dab = SQRT(rab2) 626 627 na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1 628 na2 = ipgf*ncoset(la_max(iset)) 629 nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1 630 nb2 = jpgf*ncoset(lb_max(jset)) 631 632 ! check whether we need to use fawzi's generalised collocation scheme 633 IF (rs_v(igrid_level)%rs_grid%desc%distributed) THEN 634 !tasks(4,:) is 0 for replicated, 1 for distributed 2 for exceptional distributed tasks 635 IF (tasks(4, itask) .EQ. 2) THEN 636 use_subpatch = .TRUE. 637 ELSE 638 use_subpatch = .FALSE. 639 ENDIF 640 ELSE 641 use_subpatch = .FALSE. 642 ENDIF 643 644 IF (pab_required) THEN 645 IF (iatom <= jatom) THEN 646 CALL integrate_pgf_product_rspace( & 647 la_max(iset), zeta(ipgf, iset), la_min(iset), & 648 lb_max(jset), zetb(jpgf, jset), lb_min(jset), & 649 ra, rab, rab2, rs_v(igrid_level)%rs_grid, cell, & 650 cube_info(igrid_level), & 651 hab, pab=pab, o1=na1 - 1, o2=nb1 - 1, & 652 eps_gvg_rspace=eps_gvg_rspace, & 653 calculate_forces=calculate_forces, & 654 force_a=force_a, force_b=force_b, & 655 compute_tau=my_compute_tau, map_consistent=map_consistent, & 656 use_virial=use_virial, my_virial_a=my_virial_a, & 657 my_virial_b=my_virial_b, use_subpatch=use_subpatch, subpatch_pattern=tasks(6, itask)) 658 ELSE 659 rab_inv = -rab 660 CALL integrate_pgf_product_rspace( & 661 lb_max(jset), zetb(jpgf, jset), lb_min(jset), & 662 la_max(iset), zeta(ipgf, iset), la_min(iset), & 663 rb, rab_inv, rab2, rs_v(igrid_level)%rs_grid, cell, & 664 cube_info(igrid_level), & 665 hab, pab=pab, o1=nb1 - 1, o2=na1 - 1, & 666 eps_gvg_rspace=eps_gvg_rspace, & 667 calculate_forces=calculate_forces, & 668 force_a=force_b, force_b=force_a, & 669 compute_tau=my_compute_tau, map_consistent=map_consistent, & 670 use_virial=use_virial, my_virial_a=my_virial_b, & 671 my_virial_b=my_virial_a, use_subpatch=use_subpatch, subpatch_pattern=tasks(6, itask)) 672 END IF 673 ELSE 674 IF (iatom <= jatom) THEN 675 CALL integrate_pgf_product_rspace( & 676 la_max(iset), zeta(ipgf, iset), la_min(iset), & 677 lb_max(jset), zetb(jpgf, jset), lb_min(jset), & 678 ra, rab, rab2, rs_v(igrid_level)%rs_grid, cell, & 679 cube_info(igrid_level), & 680 hab, o1=na1 - 1, o2=nb1 - 1, & 681 eps_gvg_rspace=eps_gvg_rspace, & 682 calculate_forces=calculate_forces, & 683 force_a=force_a, force_b=force_b, & 684 compute_tau=my_compute_tau, & 685 map_consistent=map_consistent, use_subpatch=use_subpatch, subpatch_pattern=tasks(6, itask)) 686 ELSE 687 rab_inv = -rab 688 CALL integrate_pgf_product_rspace( & 689 lb_max(jset), zetb(jpgf, jset), lb_min(jset), & 690 la_max(iset), zeta(ipgf, iset), la_min(iset), & 691 rb, rab_inv, rab2, rs_v(igrid_level)%rs_grid, cell, & 692 cube_info(igrid_level), & 693 hab, o1=nb1 - 1, o2=na1 - 1, & 694 eps_gvg_rspace=eps_gvg_rspace, & 695 calculate_forces=calculate_forces, & 696 force_a=force_b, force_b=force_a, & 697 compute_tau=my_compute_tau, & 698 map_consistent=map_consistent, use_subpatch=use_subpatch, subpatch_pattern=tasks(6, itask)) 699 END IF 700 END IF 701 702 new_set_pair_coming = .FALSE. 703 atom_pair_done = .FALSE. 704 IF (itask < task_list%taskstop(ipair, igrid_level)) THEN 705 CALL int2pair(tasks(3, itask + 1), ilevel, img, iatom, jatom, iset_new, jset_new, ipgf_new, jpgf_new, & 706 nimages, natom, maxset, maxpgf) 707 IF (iset_new .NE. iset .OR. jset_new .NE. jset) THEN 708 new_set_pair_coming = .TRUE. 709 ENDIF 710 ELSE 711 ! do not forget the last block 712 new_set_pair_coming = .TRUE. 713 atom_pair_done = .TRUE. 714 ENDIF 715 716 ! contract the block into h if we're done with the current set pair 717 IF (new_set_pair_coming) THEN 718 IF (iatom <= jatom) THEN 719 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, & 720 1.0_dp, hab(1, 1), SIZE(hab, 1), & 721 sphi_b(1, sgfb), SIZE(sphi_b, 1), & 722 0.0_dp, work(1, 1), SIZE(work, 1)) 723 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, & 724 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 725 work(1, 1), SIZE(work, 1), & 726 1.0_dp, h_block(sgfa, sgfb), SIZE(h_block, 1)) 727 ELSE 728 CALL dgemm("N", "N", ncob, nsgfa(iset), ncoa, & 729 1.0_dp, hab(1, 1), SIZE(hab, 1), & 730 sphi_a(1, sgfa), SIZE(sphi_a, 1), & 731 0.0_dp, work(1, 1), SIZE(work, 1)) 732 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncob, & 733 1.0_dp, sphi_b(1, sgfb), SIZE(sphi_b, 1), & 734 work(1, 1), SIZE(work, 1), & 735 1.0_dp, h_block(sgfb, sgfa), SIZE(h_block, 1)) 736 END IF 737 END IF 738 739 IF (atom_pair_done) THEN 740!$OMP CRITICAL(force_critical) 741 IF (iatom == jatom) THEN 742 scalef = 1.0_dp 743 ELSE 744 scalef = 2.0_dp 745 END IF 746 IF (calculate_forces) THEN 747 force(ikind)%rho_elec(:, atom_a) = & 748 force(ikind)%rho_elec(:, atom_a) + scalef*admm_scal_fac*force_a(:) 749 force(jkind)%rho_elec(:, atom_b) = & 750 force(jkind)%rho_elec(:, atom_b) + scalef*admm_scal_fac*force_b(:) 751 ENDIF 752 IF (use_virial) THEN 753 IF (use_virial .AND. calculate_forces) THEN 754 virial%pv_virial = virial%pv_virial + scalef*admm_scal_fac*my_virial_a 755 virial%pv_virial = virial%pv_virial + scalef*admm_scal_fac*my_virial_b 756 END IF 757 END IF 758!$OMP END CRITICAL (force_critical) 759 ENDIF 760 END DO loop_tasks 761 END DO loop_pairs 762!$OMP END DO 763 764 DO img = 1, nimages 765 CALL dbcsr_finalize(dhmat(img)%matrix) 766 END DO 767 768 END DO loop_gridlevels 769 770!$OMP END PARALLEL 771 772 IF (debug_this_module) THEN 773 DEALLOCATE (block_touched) 774 END IF 775 776 IF (h_duplicated) THEN 777 ! Reconstruct H matrix if using distributed RS grids 778 ! note send and recv direction reversed WRT collocate 779 scatter = .FALSE. 780 IF (do_kp) THEN 781 CALL rs_distribute_matrix(rs_descs, dhmat, atom_pair_recv, atom_pair_send, & 782 natom, nimages, scatter, hmats=hmat_kp) 783 ELSE 784 ALLOCATE (htemp(1)) 785 htemp(1)%matrix => hmat%matrix 786 787 CALL rs_distribute_matrix(rs_descs, dhmat, atom_pair_recv, atom_pair_send, & 788 natom, nimages, scatter, hmats=htemp) 789 790 NULLIFY (htemp(1)%matrix) 791 DEALLOCATE (htemp) 792 END IF 793 CALL dbcsr_deallocate_matrix_set(dhmat) 794 ELSE 795 DO img = 1, nimages 796 NULLIFY (dhmat(img)%matrix) 797 END DO 798 DEALLOCATE (dhmat) 799 END IF 800 801 IF (pab_required) THEN 802 IF (p_duplicated) THEN 803 CALL dbcsr_deallocate_matrix_set(deltap) 804 ELSE 805 DO img = 1, nimages 806 NULLIFY (deltap(img)%matrix) 807 END DO 808 DEALLOCATE (deltap) 809 END IF 810 END IF 811 812 ! *** Release work storage *** 813 814 DEALLOCATE (habt, workt) 815 816 IF (pab_required) THEN 817 DEALLOCATE (pabt) 818 END IF 819 820 IF (ASSOCIATED(rs_v)) THEN 821 DO i = 1, SIZE(rs_v) 822 CALL rs_grid_release(rs_v(i)%rs_grid) 823 END DO 824 END IF 825 826 IF (calculate_forces) THEN 827 DEALLOCATE (atom_of_kind) 828 END IF 829 830 CALL timestop(handle) 831 832 END SUBROUTINE integrate_v_rspace 833 834END MODULE qs_integrate_potential_product 835