1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5! ************************************************************************************************** 6!> \brief Calculation of the non-local pseudopotential contribution to the core Hamiltonian 7!> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b> 8!> \par History 9!> - refactered from qs_core_hamiltian [Joost VandeVondele, 2008-11-01] 10!> - full rewrite [jhu, 2009-01-23] 11! ************************************************************************************************** 12MODULE core_ppnl 13 USE ai_overlap, ONLY: overlap 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind_set 16 USE basis_set_types, ONLY: gto_basis_set_p_type,& 17 gto_basis_set_type 18 USE dbcsr_api, ONLY: dbcsr_add,& 19 dbcsr_get_block_p,& 20 dbcsr_p_type 21 USE external_potential_types, ONLY: gth_potential_p_type,& 22 gth_potential_type,& 23 sgp_potential_p_type,& 24 sgp_potential_type 25 USE kinds, ONLY: dp 26 USE orbital_pointers, ONLY: init_orbital_pointers,& 27 nco,& 28 ncoset 29 USE particle_types, ONLY: particle_type 30 USE qs_force_types, ONLY: qs_force_type 31 USE qs_kind_types, ONLY: get_qs_kind,& 32 get_qs_kind_set,& 33 qs_kind_type 34 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 35 neighbor_list_iterate,& 36 neighbor_list_iterator_create,& 37 neighbor_list_iterator_p_type,& 38 neighbor_list_iterator_release,& 39 neighbor_list_set_p_type 40 USE sap_kind_types, ONLY: alist_type,& 41 clist_type,& 42 get_alist,& 43 release_sap_int,& 44 sap_int_type,& 45 sap_sort 46 USE virial_methods, ONLY: virial_pair_force 47 USE virial_types, ONLY: virial_type 48 49!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads 50 51#include "./base/base_uses.f90" 52 53 IMPLICIT NONE 54 55 PRIVATE 56 57 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'core_ppnl' 58 59 PUBLIC :: build_core_ppnl 60 61CONTAINS 62 63! ************************************************************************************************** 64!> \brief ... 65!> \param matrix_h ... 66!> \param matrix_p ... 67!> \param force ... 68!> \param virial ... 69!> \param calculate_forces ... 70!> \param use_virial ... 71!> \param nder ... 72!> \param qs_kind_set ... 73!> \param atomic_kind_set ... 74!> \param particle_set ... 75!> \param sab_orb ... 76!> \param sap_ppnl ... 77!> \param eps_ppnl ... 78!> \param nimages ... 79!> \param cell_to_index ... 80!> \param basis_type ... 81! ************************************************************************************************** 82 SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & 83 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, & 84 nimages, cell_to_index, basis_type) 85 86 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p 87 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 88 TYPE(virial_type), POINTER :: virial 89 LOGICAL, INTENT(IN) :: calculate_forces 90 LOGICAL :: use_virial 91 INTEGER :: nder 92 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 93 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 94 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 95 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 96 POINTER :: sab_orb, sap_ppnl 97 REAL(KIND=dp), INTENT(IN) :: eps_ppnl 98 INTEGER, INTENT(IN) :: nimages 99 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 100 CHARACTER(LEN=*), INTENT(IN) :: basis_type 101 102 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_ppnl', & 103 routineP = moduleN//':'//routineN 104 105 INTEGER :: atom_a, atom_b, atom_c, first_col, handle, i, iab, iac, iatom, ibc, icol, ikind, & 106 ilist, img, inode, irow, iset, j, jatom, jkind, jneighbor, kac, katom, kbc, kkind, l, & 107 lc_max, lc_min, ldai, ldsab, lppnl, maxco, maxder, maxl, maxlgto, maxlppnl, maxppnl, & 108 maxsgf, mepos, na, natom, nb, ncoa, ncoc, nkind, nlist, nneighbor, nnl, nnode, np, nppnl, & 109 nprjc, nseta, nsgfa, nthread, prjc, sgfa 110 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 111 INTEGER, DIMENSION(3) :: cell_b, cell_c 112 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nprj_ppnl, & 113 nsgf_seta 114 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa 115 LOGICAL :: dogth, dokp, found, ppnl_present 116 LOGICAL, DIMENSION(0:9) :: is_nonlocal 117 REAL(KIND=dp) :: dac, f0, ppnl_radius 118 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radp 119 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work 120 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work 121 REAL(KIND=dp), DIMENSION(1) :: rprjc, zetc 122 REAL(KIND=dp), DIMENSION(3) :: fa, fb, rab, rac, rbc 123 REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc 124 REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, alpha_ppnl, hprj, set_radius_a 125 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj, h_block, h_nl, p_block, rpgfa, & 126 sphi_a, vprj_ppnl, zeta 127 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: achint, acint, bchint, bcint, c_nl 128 TYPE(alist_type), POINTER :: alist_ac, alist_bc 129 TYPE(clist_type), POINTER :: clist 130 TYPE(gth_potential_p_type), DIMENSION(:), POINTER :: gpotential 131 TYPE(gth_potential_type), POINTER :: gth_potential 132 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set 133 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 134 TYPE(neighbor_list_iterator_p_type), & 135 DIMENSION(:), POINTER :: nl_iterator 136 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int 137 TYPE(sgp_potential_p_type), DIMENSION(:), POINTER :: spotential 138 TYPE(sgp_potential_type), POINTER :: sgp_potential 139 140 IF (calculate_forces) THEN 141 CALL timeset(routineN//"_forces", handle) 142 ELSE 143 CALL timeset(routineN, handle) 144 ENDIF 145 146 ppnl_present = ASSOCIATED(sap_ppnl) 147 148 IF (ppnl_present) THEN 149 150 nkind = SIZE(atomic_kind_set) 151 natom = SIZE(particle_set) 152 153 dokp = (nimages > 1) 154 155 ALLOCATE (atom_of_kind(natom)) 156 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind) 157 158 IF (calculate_forces) THEN 159 IF (SIZE(matrix_p, 1) == 2) THEN 160 DO img = 1, nimages 161 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & 162 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 163 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, & 164 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp) 165 END DO 166 END IF 167 END IF 168 169 IF (use_virial) THEN 170 pv_loc = virial%pv_virial 171 END IF 172 173 maxder = ncoset(nder) 174 175 CALL get_qs_kind_set(qs_kind_set, & 176 maxco=maxco, & 177 maxlgto=maxlgto, & 178 maxsgf=maxsgf, & 179 maxlppnl=maxlppnl, & 180 maxppnl=maxppnl, & 181 basis_type=basis_type) 182 183 maxl = MAX(maxlgto, maxlppnl) 184 CALL init_orbital_pointers(maxl + nder + 1) 185 186 ldsab = MAX(maxco, ncoset(maxlppnl), maxsgf, maxppnl) 187 ldai = ncoset(maxl + nder + 1) 188 189 !sap_int needs to be shared as multiple threads need to access this 190 ALLOCATE (sap_int(nkind*nkind)) 191 DO i = 1, nkind*nkind 192 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex) 193 sap_int(i)%nalist = 0 194 END DO 195 196 !set up direct access to basis and potential 197 ALLOCATE (basis_set(nkind), gpotential(nkind), spotential(nkind)) 198 DO ikind = 1, nkind 199 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type) 200 IF (ASSOCIATED(orb_basis_set)) THEN 201 basis_set(ikind)%gto_basis_set => orb_basis_set 202 ELSE 203 NULLIFY (basis_set(ikind)%gto_basis_set) 204 END IF 205 CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential) 206 NULLIFY (gpotential(ikind)%gth_potential) 207 NULLIFY (spotential(ikind)%sgp_potential) 208 IF (ASSOCIATED(gth_potential)) THEN 209 gpotential(ikind)%gth_potential => gth_potential 210 ELSE IF (ASSOCIATED(sgp_potential)) THEN 211 spotential(ikind)%sgp_potential => sgp_potential 212 END IF 213 END DO 214 215 nthread = 1 216!$ nthread = omp_get_max_threads() 217 218 !calculate the overlap integrals <a|p> 219 CALL neighbor_list_iterator_create(nl_iterator, sap_ppnl, nthread=nthread) 220!$OMP PARALLEL & 221!$OMP DEFAULT (NONE) & 222!$OMP SHARED (nl_iterator, basis_set, gpotential, spotential, maxder, ncoset, & 223!$OMP sap_int, nkind, ldsab, ldai, nder, nco ) & 224!$OMP PRIVATE (mepos, ikind, kkind, iatom, katom, nlist, ilist, nneighbor, jneighbor, & 225!$OMP cell_c, rac, iac, first_sgfa, la_max, la_min, npgfa, nseta, nsgfa, nsgf_seta, & 226!$OMP sphi_a, zeta, cprj, hprj, lppnl, nppnl, nprj_ppnl, & 227!$OMP clist, iset, ncoa, sgfa, prjc, work, sab, ai_work, nprjc, ppnl_radius, & 228!$OMP ncoc, rpgfa, first_col, vprj_ppnl, i, j, l, dogth, & 229!$OMP set_radius_a, rprjc, dac, lc_max, lc_min, zetc, alpha_ppnl, & 230!$OMP nnl, is_nonlocal, a_nl, h_nl, c_nl, radp) 231 mepos = 0 232!$ mepos = omp_get_thread_num() 233 234 ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder)) 235 sab = 0.0_dp 236 ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1))) 237 ai_work = 0.0_dp 238 239 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 240 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=kkind, iatom=iatom, & 241 jatom=katom, nlist=nlist, ilist=ilist, nnode=nneighbor, & 242 inode=jneighbor, cell=cell_c, r=rac) 243 iac = ikind + nkind*(kkind - 1) 244 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE 245 ! get definition of basis set 246 first_sgfa => basis_set(ikind)%gto_basis_set%first_sgf 247 la_max => basis_set(ikind)%gto_basis_set%lmax 248 la_min => basis_set(ikind)%gto_basis_set%lmin 249 npgfa => basis_set(ikind)%gto_basis_set%npgf 250 nseta = basis_set(ikind)%gto_basis_set%nset 251 nsgfa = basis_set(ikind)%gto_basis_set%nsgf 252 nsgf_seta => basis_set(ikind)%gto_basis_set%nsgf_set 253 rpgfa => basis_set(ikind)%gto_basis_set%pgf_radius 254 set_radius_a => basis_set(ikind)%gto_basis_set%set_radius 255 sphi_a => basis_set(ikind)%gto_basis_set%sphi 256 zeta => basis_set(ikind)%gto_basis_set%zet 257 ! get definition of PP projectors 258 IF (ASSOCIATED(gpotential(kkind)%gth_potential)) THEN 259 ! GTH potential 260 dogth = .TRUE. 261 alpha_ppnl => gpotential(kkind)%gth_potential%alpha_ppnl 262 cprj => gpotential(kkind)%gth_potential%cprj 263 lppnl = gpotential(kkind)%gth_potential%lppnl 264 nppnl = gpotential(kkind)%gth_potential%nppnl 265 nprj_ppnl => gpotential(kkind)%gth_potential%nprj_ppnl 266 ppnl_radius = gpotential(kkind)%gth_potential%ppnl_radius 267 vprj_ppnl => gpotential(kkind)%gth_potential%vprj_ppnl 268 ELSE IF (ASSOCIATED(spotential(kkind)%sgp_potential)) THEN 269 ! SGP potential 270 dogth = .FALSE. 271 nprjc = spotential(kkind)%sgp_potential%nppnl 272 IF (nprjc == 0) CYCLE 273 nnl = spotential(kkind)%sgp_potential%n_nonlocal 274 lppnl = spotential(kkind)%sgp_potential%lmax 275 is_nonlocal = .FALSE. 276 is_nonlocal(0:lppnl) = spotential(kkind)%sgp_potential%is_nonlocal(0:lppnl) 277 a_nl => spotential(kkind)%sgp_potential%a_nonlocal 278 h_nl => spotential(kkind)%sgp_potential%h_nonlocal 279 c_nl => spotential(kkind)%sgp_potential%c_nonlocal 280 ppnl_radius = spotential(kkind)%sgp_potential%ppnl_radius 281 ALLOCATE (radp(nnl)) 282 radp(:) = ppnl_radius 283 cprj => spotential(kkind)%sgp_potential%cprj_ppnl 284 hprj => spotential(kkind)%sgp_potential%vprj_ppnl 285 nppnl = SIZE(cprj, 2) 286 ELSE 287 CYCLE 288 END IF 289!$OMP CRITICAL(sap_int_critical) 290 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN 291 sap_int(iac)%a_kind = ikind 292 sap_int(iac)%p_kind = kkind 293 sap_int(iac)%nalist = nlist 294 ALLOCATE (sap_int(iac)%alist(nlist)) 295 DO i = 1, nlist 296 NULLIFY (sap_int(iac)%alist(i)%clist) 297 sap_int(iac)%alist(i)%aatom = 0 298 sap_int(iac)%alist(i)%nclist = 0 299 END DO 300 END IF 301 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN 302 sap_int(iac)%alist(ilist)%aatom = iatom 303 sap_int(iac)%alist(ilist)%nclist = nneighbor 304 ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor)) 305 DO i = 1, nneighbor 306 sap_int(iac)%alist(ilist)%clist(i)%catom = 0 307 END DO 308 END IF 309!$OMP END CRITICAL(sap_int_critical) 310 dac = SQRT(SUM(rac*rac)) 311 clist => sap_int(iac)%alist(ilist)%clist(jneighbor) 312 clist%catom = katom 313 clist%cell = cell_c 314 clist%rac = rac 315 ALLOCATE (clist%acint(nsgfa, nppnl, maxder), & 316 clist%achint(nsgfa, nppnl, maxder)) 317 clist%acint = 0._dp 318 clist%achint = 0._dp 319 clist%nsgf_cnt = 0 320 NULLIFY (clist%sgf_list) 321 DO iset = 1, nseta 322 ncoa = npgfa(iset)*ncoset(la_max(iset)) 323 sgfa = first_sgfa(1, iset) 324 IF (dogth) THEN 325 ! GTH potential 326 ! XXX fix, use correct bounds 327 prjc = 1 328 work = 0._dp 329 DO l = 0, lppnl 330 nprjc = nprj_ppnl(l)*nco(l) 331 IF (nprjc == 0) CYCLE 332 rprjc(1) = ppnl_radius 333 IF (set_radius_a(iset) + rprjc(1) < dac) CYCLE 334 lc_max = l + 2*(nprj_ppnl(l) - 1) 335 lc_min = l 336 zetc(1) = alpha_ppnl(l) 337 ncoc = ncoset(lc_max) 338 339 ! *** Calculate the primitive overlap integrals *** 340 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 341 lc_max, lc_min, 1, rprjc, zetc, rac, dac, sab, nder, .TRUE., ai_work, ldai) 342 ! *** Transformation step projector functions (cartesian->spherical) *** 343 DO i = 1, maxder 344 first_col = (i - 1)*ldsab 345 CALL dgemm("N", "N", ncoa, nprjc, ncoc, 1.0_dp, sab(1, first_col + 1), SIZE(sab, 1), & 346 cprj(1, prjc), SIZE(cprj, 1), 0.0_dp, work(1, first_col + prjc), ldsab) 347 END DO 348 prjc = prjc + nprjc 349 END DO 350 DO i = 1, maxder 351 first_col = (i - 1)*ldsab + 1 352 ! *** Contraction step (basis functions) *** 353 CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 354 work(1, first_col), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa) 355 ! *** Multiply with interaction matrix(h) *** 356 CALL dgemm("N", "N", nsgf_seta(iset), nppnl, nppnl, 1.0_dp, clist%acint(sgfa, 1, i), nsgfa, & 357 vprj_ppnl(1, 1), SIZE(vprj_ppnl, 1), 0.0_dp, clist%achint(sgfa, 1, i), nsgfa) 358 END DO 359 ELSE 360 ! SGP potential 361 ! *** Calculate the primitive overlap integrals *** 362 CALL overlap(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 363 lppnl, 0, nnl, radp, a_nl, rac, dac, sab, nder, .TRUE., ai_work, ldai) 364 DO i = 1, maxder 365 first_col = (i - 1)*ldsab + 1 366 ! *** Transformation step projector functions (cartesian->spherical) *** 367 CALL dgemm("N", "N", ncoa, nppnl, nprjc, 1.0_dp, sab(1, first_col), ldsab, & 368 cprj(1, 1), SIZE(cprj, 1), 0.0_dp, work(1, 1), ldsab) 369 ! *** Contraction step (basis functions) *** 370 CALL dgemm("T", "N", nsgf_seta(iset), nppnl, ncoa, 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), & 371 work(1, 1), ldsab, 0.0_dp, clist%acint(sgfa, 1, i), nsgfa) 372 ! *** Multiply with interaction matrix(h) *** 373 ncoc = sgfa + nsgf_seta(iset) - 1 374 DO j = 1, nppnl 375 clist%achint(sgfa:ncoc, j, i) = clist%acint(sgfa:ncoc, j, i)*hprj(j) 376 END DO 377 END DO 378 END IF 379 END DO 380 clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1))) 381 clist%maxach = MAXVAL(ABS(clist%achint(:, :, 1))) 382 IF (.NOT. dogth) DEALLOCATE (radp) 383 END DO 384 385 DEALLOCATE (sab, ai_work, work) 386!$OMP END PARALLEL 387 CALL neighbor_list_iterator_release(nl_iterator) 388 389 ! *** Set up a sorting index 390 CALL sap_sort(sap_int) 391 ! *** All integrals needed have been calculated and stored in sap_int 392 ! *** We now calculate the Hamiltonian matrix elements 393 CALL neighbor_list_iterator_create(nl_iterator, sab_orb, nthread=nthread) 394 395!$OMP PARALLEL & 396!$OMP DEFAULT (NONE) & 397!$OMP SHARED (nl_iterator, dokp, basis_set, atom_of_kind, matrix_h, cell_to_index,& 398!$OMP matrix_p, sap_int, nkind, eps_ppnl, force, virial, use_virial, calculate_forces) & 399!$OMP PRIVATE (mepos, ikind, jkind, iatom, jatom, nlist, ilist, nnode, inode, cell_b, rab, & 400!$OMP iab, atom_a, atom_b, f0, irow, icol, h_block, & 401!$OMP found,p_block, iac, ibc, alist_ac, alist_bc, acint, bcint, & 402!$OMP achint, bchint, na, np, nb, katom, atom_c, j, fa, fb, rbc, rac, & 403!$OMP kkind, kac, kbc, i, img) 404 405 mepos = 0 406!$ mepos = omp_get_thread_num() 407 408 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 409 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, & 410 jatom=jatom, nlist=nlist, ilist=ilist, nnode=nnode, inode=inode, cell=cell_b, r=rab) 411 IF (.NOT. ASSOCIATED(basis_set(ikind)%gto_basis_set)) CYCLE 412 IF (.NOT. ASSOCIATED(basis_set(jkind)%gto_basis_set)) CYCLE 413 iab = ikind + nkind*(jkind - 1) 414 atom_a = atom_of_kind(iatom) 415 atom_b = atom_of_kind(jatom) 416 417 ! *** Use the symmetry of the first derivatives *** 418 IF (iatom == jatom) THEN 419 f0 = 1.0_dp 420 ELSE 421 f0 = 2.0_dp 422 END IF 423 424 IF (dokp) THEN 425 img = cell_to_index(cell_b(1), cell_b(2), cell_b(3)) 426 ELSE 427 img = 1 428 END IF 429 430 ! *** Create matrix blocks for a new matrix block column *** 431 IF (iatom <= jatom) THEN 432 irow = iatom 433 icol = jatom 434 ELSE 435 irow = jatom 436 icol = iatom 437 END IF 438 NULLIFY (h_block) 439 CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found) 440 IF (calculate_forces) THEN 441 NULLIFY (p_block) 442 CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found) 443 END IF 444 445 ! loop over all kinds for projector atom 446 IF (ASSOCIATED(h_block)) THEN 447 DO kkind = 1, nkind 448 iac = ikind + nkind*(kkind - 1) 449 ibc = jkind + nkind*(kkind - 1) 450 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE 451 IF (.NOT. ASSOCIATED(sap_int(ibc)%alist)) CYCLE 452 CALL get_alist(sap_int(iac), alist_ac, iatom) 453 CALL get_alist(sap_int(ibc), alist_bc, jatom) 454 IF (.NOT. ASSOCIATED(alist_ac)) CYCLE 455 IF (.NOT. ASSOCIATED(alist_bc)) CYCLE 456 DO kac = 1, alist_ac%nclist 457 DO kbc = 1, alist_bc%nclist 458 IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE 459 IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN 460 IF (alist_ac%clist(kac)%maxac*alist_bc%clist(kbc)%maxach < eps_ppnl) CYCLE 461 acint => alist_ac%clist(kac)%acint 462 bcint => alist_bc%clist(kbc)%acint 463 achint => alist_ac%clist(kac)%achint 464 bchint => alist_bc%clist(kbc)%achint 465 na = SIZE(acint, 1) 466 np = SIZE(acint, 2) 467 nb = SIZE(bcint, 1) 468!$OMP CRITICAL(h_block_critical) 469 IF (iatom <= jatom) THEN 470 CALL dgemm("N", "T", na, nb, np, 1._dp, achint(1, 1, 1), na, & 471 bcint(1, 1, 1), nb, 1.0_dp, h_block, SIZE(h_block, 1)) 472 ELSE 473 CALL dgemm("N", "T", nb, na, np, 1.0_dp, bchint(1, 1, 1), nb, & 474 acint(1, 1, 1), na, 1.0_dp, h_block, SIZE(h_block, 1)) 475 END IF 476!$OMP END CRITICAL(h_block_critical) 477 IF (calculate_forces) THEN 478 IF (ASSOCIATED(p_block)) THEN 479 katom = alist_ac%clist(kac)%catom 480 atom_c = atom_of_kind(katom) 481 DO i = 1, 3 482 j = i + 1 483 IF (iatom <= jatom) THEN 484 fa(i) = SUM(p_block(1:na, 1:nb)* & 485 MATMUL(acint(1:na, 1:np, j), TRANSPOSE(bchint(1:nb, 1:np, 1)))) 486 fb(i) = SUM(p_block(1:na, 1:nb)* & 487 MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, j)))) 488 ELSE 489 fa(i) = SUM(p_block(1:nb, 1:na)* & 490 MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, j)))) 491 fb(i) = SUM(p_block(1:nb, 1:na)* & 492 MATMUL(bcint(1:nb, 1:np, j), TRANSPOSE(achint(1:na, 1:np, 1)))) 493 END IF 494!$OMP CRITICAL(force_critical) 495 force(ikind)%gth_ppnl(i, atom_a) = force(ikind)%gth_ppnl(i, atom_a) + f0*fa(i) 496 force(kkind)%gth_ppnl(i, atom_c) = force(kkind)%gth_ppnl(i, atom_c) - f0*fa(i) 497 force(jkind)%gth_ppnl(i, atom_b) = force(jkind)%gth_ppnl(i, atom_b) + f0*fb(i) 498 force(kkind)%gth_ppnl(i, atom_c) = force(kkind)%gth_ppnl(i, atom_c) - f0*fb(i) 499!$OMP END CRITICAL(force_critical) 500 END DO 501 502 IF (use_virial) THEN 503 rac = alist_ac%clist(kac)%rac 504 rbc = alist_bc%clist(kbc)%rac 505!$OMP CRITICAL(virial_critical) 506 CALL virial_pair_force(virial%pv_virial, f0, fa, rac) 507 CALL virial_pair_force(virial%pv_virial, f0, fb, rbc) 508!$OMP END CRITICAL(virial_critical) 509 END IF 510 ENDIF 511 END IF 512 EXIT ! We have found a match and there can be only one single match 513 END IF 514 END DO 515 END DO 516 END DO 517 ENDIF 518 END DO 519!$OMP END PARALLEL 520 CALL neighbor_list_iterator_release(nl_iterator) 521 522 CALL release_sap_int(sap_int) 523 524 DEALLOCATE (atom_of_kind) 525 DEALLOCATE (basis_set, gpotential, spotential) 526 527 IF (calculate_forces) THEN 528 ! *** If LSD, then recover alpha density and beta density *** 529 ! *** from the total density (1) and the spin density (2) *** 530 IF (SIZE(matrix_p, 1) == 2) THEN 531 DO img = 1, nimages 532 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & 533 alpha_scalar=0.5_dp, beta_scalar=0.5_dp) 534 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, & 535 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp) 536 END DO 537 END IF 538 END IF 539 540 IF (calculate_forces .AND. use_virial) THEN 541 virial%pv_ppnl = virial%pv_virial - pv_loc 542 END IF 543 544 END IF !ppnl_present 545 546 CALL timestop(handle) 547 548 END SUBROUTINE build_core_ppnl 549 550! ************************************************************************************************** 551 552END MODULE core_ppnl 553