1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of overlap matrix, its derivatives and forces 8!> \par History 9!> JGH: removed printing routines 10!> JGH: upgraded to unique routine for overlaps 11!> JGH: Add specific routine for 'forces only' 12!> Major refactoring for new overlap routines 13!> JGH: Kpoints 14!> \author Matthias Krack (03.09.2001,25.06.2003) 15! ************************************************************************************************** 16MODULE qs_overlap 17 USE ai_contraction, ONLY: block_add,& 18 contraction,& 19 decontraction,& 20 force_trace 21 USE ai_overlap, ONLY: overlap_ab 22 USE atomic_kind_types, ONLY: atomic_kind_type,& 23 get_atomic_kind_set 24 USE basis_set_types, ONLY: get_gto_basis_set,& 25 gto_basis_set_p_type,& 26 gto_basis_set_type 27 USE block_p_types, ONLY: block_p_type 28 USE cp_control_types, ONLY: dft_control_type 29 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 30 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set 31 USE dbcsr_api, ONLY: & 32 dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_get_block_p, & 33 dbcsr_p_type, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, & 34 dbcsr_type_symmetric 35 USE kinds, ONLY: default_string_length,& 36 dp 37 USE kpoint_types, ONLY: get_kpoint_info,& 38 kpoint_type 39 USE orbital_pointers, ONLY: indco,& 40 ncoset 41 USE orbital_symbols, ONLY: cgf_symbol 42 USE particle_methods, ONLY: get_particle_set 43 USE particle_types, ONLY: particle_type 44 USE qs_force_types, ONLY: qs_force_type 45 USE qs_integral_utils, ONLY: basis_set_list_setup,& 46 get_memory_usage 47 USE qs_kind_types, ONLY: qs_kind_type 48 USE qs_ks_types, ONLY: get_ks_env,& 49 qs_ks_env_type 50 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 51 get_neighbor_list_set_p,& 52 neighbor_list_iterate,& 53 neighbor_list_iterator_create,& 54 neighbor_list_iterator_p_type,& 55 neighbor_list_iterator_release,& 56 neighbor_list_set_p_type 57 USE string_utilities, ONLY: compress,& 58 uppercase 59 USE virial_methods, ONLY: virial_pair_force 60 USE virial_types, ONLY: virial_type 61 62!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num 63#include "./base/base_uses.f90" 64 65 IMPLICIT NONE 66 67 PRIVATE 68 69! *** Global parameters *** 70 71 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_overlap' 72 73 ! should be a parameter, but this triggers a bug in OMPed code with gfortran 4.9 74 INTEGER, DIMENSION(1:56), SAVE :: ndod = (/0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, & 75 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, & 76 1, 1, 1, 1, 1, 1, 1/) 77 78 INTERFACE create_sab_matrix 79 MODULE PROCEDURE create_sab_matrix_1d, create_sab_matrix_2d 80 END INTERFACE 81 82! *** Public subroutines *** 83 84 PUBLIC :: build_overlap_matrix, build_overlap_matrix_simple, & 85 build_overlap_force, create_sab_matrix 86 87CONTAINS 88 89! ************************************************************************************************** 90 91!> \brief Calculation of the overlap matrix over Cartesian Gaussian functions. 92!> \param ks_env the QS env 93!> \param matrix_s The overlap matrix to be calculated (optional) 94!> \param matrixkp_s The overlap matrices to be calculated (kpoints, optional) 95!> \param matrix_name The name of the overlap matrix (i.e. for output) 96!> \param nderivative Derivative with respect to basis origin 97!> \param basis_type_a basis set to be used for bra in <a|b> 98!> \param basis_type_b basis set to be used for ket in <a|b> 99!> \param sab_nl pair list (must be consistent with basis sets!) 100!> \param calculate_forces (optional) ... 101!> \param matrix_p density matrix for force calculation (optional) 102!> \param matrixkp_p density matrix for force calculation with k_points (optional) 103!> \date 11.03.2002 104!> \par History 105!> Enlarged functionality of this routine. Now overlap matrices based 106!> on different basis sets can be calculated, taking into account also 107!> mixed overlaps NOTE: the pointer to the overlap matrix must now be 108!> put into its corresponding env outside of this routine 109!> [Manuel Guidon] 110!> Generalized for derivatives and force calculations [JHU] 111!> Kpoints, returns overlap matrices in real space index form 112!> \author MK 113!> \version 1.0 114! ************************************************************************************************** 115 SUBROUTINE build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, & 116 nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, & 117 matrix_p, matrixkp_p) 118 119 TYPE(qs_ks_env_type), POINTER :: ks_env 120 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & 121 POINTER :: matrix_s 122 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & 123 POINTER :: matrixkp_s 124 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name 125 INTEGER, INTENT(IN), OPTIONAL :: nderivative 126 CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b 127 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 128 POINTER :: sab_nl 129 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces 130 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p 131 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & 132 POINTER :: matrixkp_p 133 134 CHARACTER(len=*), PARAMETER :: routineN = 'build_overlap_matrix', & 135 routineP = moduleN//':'//routineN 136 137 INTEGER :: atom_a, atom_b, handle, i, iatom, ic, icol, ikind, irow, iset, jatom, jkind, & 138 jset, ldsab, maxder, maxs, mepos, n1, n2, natom, ncoa, ncob, nder, nimg, nkind, nseta, & 139 nsetb, nthread, sgfa, sgfb 140 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 141 INTEGER, DIMENSION(3) :: cell 142 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 143 npgfb, nsgfa, nsgfb 144 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 145 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 146 LOGICAL :: do_forces, do_symmetric, dokp, found, & 147 trans, use_cell_mapping, use_virial 148 REAL(KIND=dp) :: dab, f, f0, ff, rab2 149 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork, pmat 150 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint 151 REAL(KIND=dp), DIMENSION(3) :: force_a, rab 152 REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc 153 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 154 REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, & 155 zeta, zetb 156 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 157 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint 158 TYPE(dft_control_type), POINTER :: dft_control 159 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b 160 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 161 TYPE(kpoint_type), POINTER :: kpoints 162 TYPE(neighbor_list_iterator_p_type), & 163 DIMENSION(:), POINTER :: nl_iterator 164 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 165 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 166 TYPE(virial_type), POINTER :: virial 167 168 NULLIFY (dft_control) 169 170 CALL timeset(routineN, handle) 171 172 ! test for matrices (kpoints or standard gamma point) 173 IF (PRESENT(matrix_s)) THEN 174 dokp = .FALSE. 175 use_cell_mapping = .FALSE. 176 ELSEIF (PRESENT(matrixkp_s)) THEN 177 dokp = .TRUE. 178 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) 179 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 180 use_cell_mapping = (SIZE(cell_to_index) > 1) 181 ELSE 182 CPABORT("") 183 END IF 184 185 NULLIFY (atomic_kind_set) 186 CALL get_ks_env(ks_env, & 187 atomic_kind_set=atomic_kind_set, & 188 natom=natom, & 189 qs_kind_set=qs_kind_set, & 190 dft_control=dft_control) 191 192 nimg = dft_control%nimages 193 nkind = SIZE(qs_kind_set) 194 195 ALLOCATE (atom_of_kind(natom)) 196 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind) 197 198 IF (PRESENT(calculate_forces)) THEN 199 do_forces = calculate_forces 200 ELSE 201 do_forces = .FALSE. 202 END IF 203 204 IF (PRESENT(nderivative)) THEN 205 nder = nderivative 206 ELSE 207 nder = 0 208 END IF 209 maxder = ncoset(nder) 210 211 ! check for symmetry 212 CPASSERT(SIZE(sab_nl) > 0) 213 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric) 214 IF (do_symmetric) THEN 215 CPASSERT(basis_type_a == basis_type_b) 216 END IF 217 218 ! set up basis set lists 219 ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind)) 220 CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set) 221 CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set) 222 223 IF (dokp) THEN 224 CALL dbcsr_allocate_matrix_set(matrixkp_s, maxder, nimg) 225 CALL create_sab_matrix(ks_env, matrixkp_s, matrix_name, basis_set_list_a, basis_set_list_b, & 226 sab_nl, do_symmetric) 227 ELSE 228 CALL dbcsr_allocate_matrix_set(matrix_s, maxder) 229 CALL create_sab_matrix(ks_env, matrix_s, matrix_name, basis_set_list_a, basis_set_list_b, & 230 sab_nl, do_symmetric) 231 END IF 232 maxs = maxder 233 234 IF (do_forces) THEN 235 CALL get_ks_env(ks_env=ks_env, force=force, virial=virial) 236 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 237 pv_loc = virial%pv_virial 238 END IF 239 240 ldsab = get_memory_usage(qs_kind_set, basis_type_a, basis_type_b) 241 IF (do_forces) THEN 242 ! we need density matrix for forces 243 IF (dokp) THEN 244 CPASSERT(PRESENT(matrixkp_p)) 245 ELSE 246 CPASSERT(PRESENT(matrix_p)) 247 END IF 248 nder = MAX(nder, 1) 249 END IF 250 maxder = ncoset(nder) 251 252 nthread = 1 253!$ nthread = omp_get_max_threads() 254 CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread) 255 256!$OMP PARALLEL DEFAULT(NONE) & 257!$OMP SHARED (nthread,do_forces,ldsab,maxder,nl_iterator, use_cell_mapping, do_symmetric,maxs,dokp,& 258!$OMP ncoset,nder,use_virial,force,virial,ndod,& 259!$OMP matrix_s, matrix_p,basis_set_list_a, basis_set_list_b, atom_of_kind, cell_to_index, matrixkp_s, matrixkp_p) & 260!$OMP PRIVATE (mepos,oint,owork,pmat,sint,ikind,jkind,iatom,jatom,rab,cell,basis_set_a,basis_set_b,atom_a,atom_b,& 261!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, & 262!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, f, & 263!$OMP zetb, scon_a, scon_b, ic, irow, icol, f0, ff, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset ) 264 265 mepos = 0 266!$ mepos = omp_get_thread_num() 267 268 NULLIFY (p_block) 269 270 ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab)) 271 IF (do_forces) ALLOCATE (pmat(ldsab, ldsab)) 272 ALLOCATE (sint(maxs)) 273 DO i = 1, maxs 274 NULLIFY (sint(i)%block) 275 END DO 276 277 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 278 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, & 279 iatom=iatom, jatom=jatom, r=rab, cell=cell) 280 basis_set_a => basis_set_list_a(ikind)%gto_basis_set 281 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 282 basis_set_b => basis_set_list_b(jkind)%gto_basis_set 283 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 284 atom_a = atom_of_kind(iatom) 285 atom_b = atom_of_kind(jatom) 286 ! basis ikind 287 first_sgfa => basis_set_a%first_sgf 288 la_max => basis_set_a%lmax 289 la_min => basis_set_a%lmin 290 npgfa => basis_set_a%npgf 291 nseta = basis_set_a%nset 292 nsgfa => basis_set_a%nsgf_set 293 rpgfa => basis_set_a%pgf_radius 294 set_radius_a => basis_set_a%set_radius 295 scon_a => basis_set_a%scon 296 zeta => basis_set_a%zet 297 ! basis jkind 298 first_sgfb => basis_set_b%first_sgf 299 lb_max => basis_set_b%lmax 300 lb_min => basis_set_b%lmin 301 npgfb => basis_set_b%npgf 302 nsetb = basis_set_b%nset 303 nsgfb => basis_set_b%nsgf_set 304 rpgfb => basis_set_b%pgf_radius 305 set_radius_b => basis_set_b%set_radius 306 scon_b => basis_set_b%scon 307 zetb => basis_set_b%zet 308 309 IF (use_cell_mapping) THEN 310 ic = cell_to_index(cell(1), cell(2), cell(3)) 311 CPASSERT(ic > 0) 312 ELSE 313 ic = 1 314 END IF 315 316 IF (do_symmetric) THEN 317 IF (iatom <= jatom) THEN 318 irow = iatom 319 icol = jatom 320 ELSE 321 irow = jatom 322 icol = iatom 323 END IF 324 f0 = 2.0_dp 325 ff = 2.0_dp 326 IF (iatom == jatom) f0 = 1.0_dp 327 ELSE 328 irow = iatom 329 icol = jatom 330 f0 = 1.0_dp 331 ff = 1.0_dp 332 END IF 333 DO i = 1, maxs 334 NULLIFY (sint(i)%block) 335 IF (dokp) THEN 336 CALL dbcsr_get_block_p(matrix=matrixkp_s(i, ic)%matrix, & 337 row=irow, col=icol, BLOCK=sint(i)%block, found=found) 338 CPASSERT(found) 339 ELSE 340 CALL dbcsr_get_block_p(matrix=matrix_s(i)%matrix, & 341 row=irow, col=icol, BLOCK=sint(i)%block, found=found) 342 CPASSERT(found) 343 END IF 344 END DO 345 IF (do_forces) THEN 346 NULLIFY (p_block) 347 IF (dokp) THEN 348 CALL dbcsr_get_block_p(matrix=matrixkp_p(1, ic)%matrix, & 349 row=irow, col=icol, block=p_block, found=found) 350 CPASSERT(found) 351 ELSE 352 CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, & 353 block=p_block, found=found) 354 CPASSERT(found) 355 END IF 356 END IF 357 trans = do_symmetric .AND. (iatom > jatom) 358 359 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 360 dab = SQRT(rab2) 361 362 DO iset = 1, nseta 363 364 ncoa = npgfa(iset)*ncoset(la_max(iset)) 365 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1)) 366 sgfa = first_sgfa(1, iset) 367 368 DO jset = 1, nsetb 369 370 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 371 372 ncob = npgfb(jset)*ncoset(lb_max(jset)) 373 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)) 374 sgfb = first_sgfb(1, jset) 375 376 ! calculate integrals and derivatives 377 SELECT CASE (nder) 378 CASE (0) 379 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 380 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 381 rab, sab=oint(:, :, 1)) 382 CASE (1) 383 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 384 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 385 rab, sab=oint(:, :, 1), dab=oint(:, :, 2:4)) 386 CASE (2) 387 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 388 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 389 rab, sab=oint(:, :, 1), dab=oint(:, :, 2:4), ddab=oint(:, :, 5:10)) 390 CASE DEFAULT 391 CPABORT("") 392 END SELECT 393 IF (do_forces .AND. ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN 394 ! Decontract P matrix block 395 owork = 0.0_dp 396 CALL block_add("OUT", owork, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans) 397 CALL decontraction(owork, pmat, scon_a(:, sgfa:), n1, nsgfa(iset), scon_b(:, sgfb:), n2, nsgfb(jset), & 398 trans=trans) 399 CALL force_trace(force_a, oint(:, :, 2:4), pmat, n1, n2, 3) 400!$OMP CRITICAL(forceupdate) 401 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - ff*force_a(:) 402 force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + ff*force_a(:) 403 IF (use_virial) THEN 404 CALL virial_pair_force(virial%pv_virial, -f0, force_a, rab) 405 END IF 406!$OMP END CRITICAL(forceupdate) 407 END IF 408 ! Contraction 409 DO i = 1, maxs 410 f = 1.0_dp 411 IF (ndod(i) == 1 .AND. trans) f = -1.0_dp 412 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), & 413 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=f, trans=trans) 414!$OMP CRITICAL(blockadd) 415 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(i)%block, & 416 sgfa, sgfb, trans=trans) 417!$OMP END CRITICAL(blockadd) 418 END DO 419 420 END DO 421 END DO 422 423 END DO 424 IF (do_forces) DEALLOCATE (pmat) 425 DEALLOCATE (oint, owork) 426 DEALLOCATE (sint) 427!$OMP END PARALLEL 428 CALL neighbor_list_iterator_release(nl_iterator) 429 430 IF (do_forces .AND. use_virial) THEN 431 virial%pv_overlap = virial%pv_virial - pv_loc 432 ENDIF 433 434 IF (dokp) THEN 435 DO i = 1, maxs 436 DO ic = 1, nimg 437 CALL dbcsr_finalize(matrixkp_s(i, ic)%matrix) 438 CALL dbcsr_filter(matrixkp_s(i, ic)%matrix, & 439 dft_control%qs_control%eps_filter_matrix) 440 ENDDO 441 ENDDO 442 ELSE 443 DO i = 1, maxs 444 CALL dbcsr_finalize(matrix_s(i)%matrix) 445 CALL dbcsr_filter(matrix_s(i)%matrix, & 446 dft_control%qs_control%eps_filter_matrix) 447 ENDDO 448 END IF 449 450 ! *** Release work storage *** 451 DEALLOCATE (atom_of_kind) 452 DEALLOCATE (basis_set_list_a, basis_set_list_b) 453 454 CALL timestop(handle) 455 456 END SUBROUTINE build_overlap_matrix 457 458! ************************************************************************************************** 459!> \brief Calculation of the overlap matrix over Cartesian Gaussian functions. 460!> \param ks_env the QS env 461!> \param matrix_s The overlap matrix to be calculated 462!> \param basis_set_list_a basis set list to be used for bra in <a|b> 463!> \param basis_set_list_b basis set list to be used for ket in <a|b> 464!> \param sab_nl pair list (must be consistent with basis sets!) 465!> \date 11.03.2016 466!> \par History 467!> Simplified version of build_overlap_matrix 468!> \author MK 469!> \version 1.0 470! ************************************************************************************************** 471 SUBROUTINE build_overlap_matrix_simple(ks_env, matrix_s, & 472 basis_set_list_a, basis_set_list_b, sab_nl) 473 474 TYPE(qs_ks_env_type), POINTER :: ks_env 475 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 476 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b 477 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 478 POINTER :: sab_nl 479 480 CHARACTER(len=*), PARAMETER :: routineN = 'build_overlap_matrix_simple', & 481 routineP = moduleN//':'//routineN 482 483 INTEGER :: atom_a, atom_b, handle, iatom, icol, ikind, irow, iset, jatom, jkind, jset, & 484 ldsab, m1, m2, mepos, n1, n2, natom, ncoa, ncob, nkind, nseta, nsetb, nthread, sgfa, sgfb 485 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 486 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 487 npgfb, nsgfa, nsgfb 488 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 489 LOGICAL :: do_symmetric, found, trans 490 REAL(KIND=dp) :: dab, rab2 491 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork 492 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint 493 REAL(KIND=dp), DIMENSION(3) :: rab 494 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 495 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb 496 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 497 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint 498 TYPE(dft_control_type), POINTER :: dft_control 499 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 500 TYPE(neighbor_list_iterator_p_type), & 501 DIMENSION(:), POINTER :: nl_iterator 502 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 503 504 NULLIFY (dft_control) 505 506 CALL timeset(routineN, handle) 507 508 NULLIFY (atomic_kind_set) 509 CALL get_ks_env(ks_env, & 510 atomic_kind_set=atomic_kind_set, & 511 natom=natom, & 512 qs_kind_set=qs_kind_set, & 513 dft_control=dft_control) 514 515 ! check for symmetry 516 CPASSERT(SIZE(sab_nl) > 0) 517 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric) 518 519 nkind = SIZE(qs_kind_set) 520 521 ALLOCATE (atom_of_kind(natom)) 522 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind) 523 524 CALL dbcsr_allocate_matrix_set(matrix_s, 1) 525 CALL create_sab_matrix(ks_env, matrix_s, "Matrix", basis_set_list_a, basis_set_list_b, & 526 sab_nl, do_symmetric) 527 528 ldsab = 0 529 DO ikind = 1, nkind 530 basis_set_a => basis_set_list_a(ikind)%gto_basis_set 531 CALL get_gto_basis_set(gto_basis_set=basis_set_a, maxco=m1, nsgf=m2) 532 ldsab = MAX(m1, m2, ldsab) 533 basis_set_b => basis_set_list_b(ikind)%gto_basis_set 534 CALL get_gto_basis_set(gto_basis_set=basis_set_b, maxco=m1, nsgf=m2) 535 ldsab = MAX(m1, m2, ldsab) 536 END DO 537 538 nthread = 1 539!$ nthread = omp_get_max_threads() 540 CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread) 541 542!$OMP PARALLEL DEFAULT(NONE) & 543!$OMP SHARED (nthread,ldsab,nl_iterator,do_symmetric,ncoset,& 544!$OMP matrix_s,basis_set_list_a,basis_set_list_b,atom_of_kind) & 545!$OMP PRIVATE (mepos,oint,owork,sint,ikind,jkind,iatom,jatom,rab,basis_set_a,basis_set_b,atom_a,atom_b,& 546!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, & 547!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, dab, & 548!$OMP zetb, scon_a, scon_b, irow, icol, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset ) 549 550 mepos = 0 551!$ mepos = omp_get_thread_num() 552 553 ALLOCATE (oint(ldsab, ldsab, 1), owork(ldsab, ldsab)) 554 ALLOCATE (sint(1)) 555 NULLIFY (sint(1)%block) 556 557 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 558 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, & 559 iatom=iatom, jatom=jatom, r=rab) 560 basis_set_a => basis_set_list_a(ikind)%gto_basis_set 561 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 562 basis_set_b => basis_set_list_b(jkind)%gto_basis_set 563 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 564 atom_a = atom_of_kind(iatom) 565 atom_b = atom_of_kind(jatom) 566 ! basis ikind 567 first_sgfa => basis_set_a%first_sgf 568 la_max => basis_set_a%lmax 569 la_min => basis_set_a%lmin 570 npgfa => basis_set_a%npgf 571 nseta = basis_set_a%nset 572 nsgfa => basis_set_a%nsgf_set 573 rpgfa => basis_set_a%pgf_radius 574 set_radius_a => basis_set_a%set_radius 575 scon_a => basis_set_a%scon 576 zeta => basis_set_a%zet 577 ! basis jkind 578 first_sgfb => basis_set_b%first_sgf 579 lb_max => basis_set_b%lmax 580 lb_min => basis_set_b%lmin 581 npgfb => basis_set_b%npgf 582 nsetb = basis_set_b%nset 583 nsgfb => basis_set_b%nsgf_set 584 rpgfb => basis_set_b%pgf_radius 585 set_radius_b => basis_set_b%set_radius 586 scon_b => basis_set_b%scon 587 zetb => basis_set_b%zet 588 589 IF (do_symmetric) THEN 590 IF (iatom <= jatom) THEN 591 irow = iatom 592 icol = jatom 593 ELSE 594 irow = jatom 595 icol = iatom 596 END IF 597 ELSE 598 irow = iatom 599 icol = jatom 600 END IF 601 602 NULLIFY (sint(1)%block) 603 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, & 604 row=irow, col=icol, BLOCK=sint(1)%block, found=found) 605 CPASSERT(found) 606 trans = do_symmetric .AND. (iatom > jatom) 607 608 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 609 dab = SQRT(rab2) 610 611 DO iset = 1, nseta 612 613 ncoa = npgfa(iset)*ncoset(la_max(iset)) 614 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1)) 615 sgfa = first_sgfa(1, iset) 616 617 DO jset = 1, nsetb 618 619 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 620 621 ncob = npgfb(jset)*ncoset(lb_max(jset)) 622 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)) 623 sgfb = first_sgfb(1, jset) 624 625 ! calculate integrals and derivatives 626 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 627 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 628 rab, sab=oint(:, :, 1)) 629 ! Contraction 630 CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), & 631 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=trans) 632!$OMP CRITICAL(blockadd) 633 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(1)%block, & 634 sgfa, sgfb, trans=trans) 635!$OMP END CRITICAL(blockadd) 636 637 END DO 638 END DO 639 640 END DO 641 DEALLOCATE (oint, owork) 642 DEALLOCATE (sint) 643!$OMP END PARALLEL 644 CALL neighbor_list_iterator_release(nl_iterator) 645 646 CALL dbcsr_finalize(matrix_s(1)%matrix) 647 CALL dbcsr_filter(matrix_s(1)%matrix, dft_control%qs_control%eps_filter_matrix) 648 649 ! *** Release work storage *** 650 DEALLOCATE (atom_of_kind) 651 652 CALL timestop(handle) 653 654 END SUBROUTINE build_overlap_matrix_simple 655 656! ************************************************************************************************** 657 658!> \brief Calculation of the force contribution from an overlap matrix 659!> over Cartesian Gaussian functions. 660!> \param ks_env the QS environment 661!> \param force holds the calcuated force Tr(P dS/dR) 662!> \param basis_type_a basis set to be used for bra in <a|b> 663!> \param basis_type_b basis set to be used for ket in <a|b> 664!> \param sab_nl pair list (must be consistent with basis sets!) 665!> \param matrix_p density matrix for force calculation 666!> \date 11.03.2002 667!> \par History 668!> Enlarged functionality of this routine. Now overlap matrices based 669!> on different basis sets can be calculated, taking into account also 670!> mixed overlaps NOTE: the pointer to the overlap matrix must now be 671!> put into its corresponding env outside of this routine 672!> [Manuel Guidon] 673!> Generalized for derivatives and force calculations [JHU] 674!> This special version is only for forces [07.07.2014, JGH] 675!> \author MK/JGH 676!> \version 1.0 677! ************************************************************************************************** 678 SUBROUTINE build_overlap_force(ks_env, force, basis_type_a, basis_type_b, & 679 sab_nl, matrix_p) 680 681 TYPE(qs_ks_env_type), POINTER :: ks_env 682 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force 683 CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b 684 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 685 POINTER :: sab_nl 686 TYPE(dbcsr_type), POINTER :: matrix_p 687 688 CHARACTER(len=*), PARAMETER :: routineN = 'build_overlap_force', & 689 routineP = moduleN//':'//routineN 690 691 INTEGER :: handle, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, & 692 ldsab, n1, n2, ncoa, ncob, nder, nkind, nseta, nsetb, sgfa, sgfb 693 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 694 npgfb, nsgfa, nsgfb 695 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 696 LOGICAL :: do_symmetric, found, new_atom_b, trans, & 697 use_virial 698 REAL(KIND=dp) :: dab, f0, ff, rab2 699 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, sab 700 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: drab 701 REAL(KIND=dp), DIMENSION(3) :: force_a, rab 702 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 703 REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, & 704 zeta, zetb 705 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b 706 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 707 TYPE(neighbor_list_iterator_p_type), & 708 DIMENSION(:), POINTER :: nl_iterator 709 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 710 TYPE(virial_type), POINTER :: virial 711 712 CALL timeset(routineN, handle) 713 714 NULLIFY (qs_kind_set) 715 CALL get_ks_env(ks_env=ks_env, qs_kind_set=qs_kind_set) 716 717 nkind = SIZE(qs_kind_set) 718 nder = 1 719 720 ! check for symmetry 721 CPASSERT(SIZE(sab_nl) > 0) 722 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric) 723 724 CALL get_ks_env(ks_env=ks_env, virial=virial) 725 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 726 727 ! *** Allocate work storage *** 728 ldsab = get_memory_usage(qs_kind_set, basis_type_a, basis_type_b) 729 ALLOCATE (sab(ldsab, ldsab), pab(ldsab, ldsab)) 730 ALLOCATE (drab(ldsab, ldsab, 3)) 731 732 ! set up basis sets 733 ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind)) 734 CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set) 735 CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set) 736 737 ! Loop over neighbor list 738 CALL neighbor_list_iterator_create(nl_iterator, sab_nl) 739 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 740 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, & 741 iatom=iatom, jatom=jatom, r=rab) 742 basis_set_a => basis_set_list_a(ikind)%gto_basis_set 743 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 744 basis_set_b => basis_set_list_b(jkind)%gto_basis_set 745 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 746 ! basis ikind 747 first_sgfa => basis_set_a%first_sgf 748 la_max => basis_set_a%lmax 749 la_min => basis_set_a%lmin 750 npgfa => basis_set_a%npgf 751 nseta = basis_set_a%nset 752 nsgfa => basis_set_a%nsgf_set 753 rpgfa => basis_set_a%pgf_radius 754 set_radius_a => basis_set_a%set_radius 755 scon_a => basis_set_a%scon 756 zeta => basis_set_a%zet 757 ! basis jkind 758 first_sgfb => basis_set_b%first_sgf 759 lb_max => basis_set_b%lmax 760 lb_min => basis_set_b%lmin 761 npgfb => basis_set_b%npgf 762 nsetb = basis_set_b%nset 763 nsgfb => basis_set_b%nsgf_set 764 rpgfb => basis_set_b%pgf_radius 765 set_radius_b => basis_set_b%set_radius 766 scon_b => basis_set_b%scon 767 zetb => basis_set_b%zet 768 769 IF (inode == 1) last_jatom = 0 770 771 IF (jatom /= last_jatom) THEN 772 new_atom_b = .TRUE. 773 last_jatom = jatom 774 ELSE 775 new_atom_b = .FALSE. 776 END IF 777 778 IF (new_atom_b) THEN 779 IF (do_symmetric) THEN 780 IF (iatom <= jatom) THEN 781 irow = iatom 782 icol = jatom 783 ELSE 784 irow = jatom 785 icol = iatom 786 END IF 787 f0 = 2.0_dp 788 IF (iatom == jatom) f0 = 1.0_dp 789 ff = 2.0_dp 790 ELSE 791 irow = iatom 792 icol = jatom 793 f0 = 1.0_dp 794 ff = 1.0_dp 795 END IF 796 NULLIFY (p_block) 797 CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, & 798 block=p_block, found=found) 799 END IF 800 trans = do_symmetric .AND. (iatom > jatom) 801 802 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 803 dab = SQRT(rab2) 804 805 DO iset = 1, nseta 806 807 ncoa = npgfa(iset)*ncoset(la_max(iset)) 808 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1)) 809 sgfa = first_sgfa(1, iset) 810 811 DO jset = 1, nsetb 812 813 IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE 814 815 ncob = npgfb(jset)*ncoset(lb_max(jset)) 816 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)) 817 sgfb = first_sgfb(1, jset) 818 819 IF (ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN 820 ! Decontract P matrix block 821 sab = 0.0_dp 822 CALL block_add("OUT", sab, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans) 823 CALL decontraction(sab, pab, scon_a(:, sgfa:), n1, nsgfa(iset), scon_b(:, sgfb:), n2, nsgfb(jset), & 824 trans=trans) 825 ! calculate integrals and derivatives 826 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 827 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 828 rab, dab=drab) 829 CALL force_trace(force_a, drab, pab, n1, n2, 3) 830 force(1:3, iatom) = force(1:3, iatom) - ff*force_a(1:3) 831 force(1:3, jatom) = force(1:3, jatom) + ff*force_a(1:3) 832 IF (use_virial) THEN 833 CALL virial_pair_force(virial%pv_virial, -f0, force_a, rab) 834 END IF 835 END IF 836 837 END DO 838 END DO 839 840 END DO 841 CALL neighbor_list_iterator_release(nl_iterator) 842 843 ! *** Release work storage *** 844 DEALLOCATE (sab, pab, drab) 845 DEALLOCATE (basis_set_list_a, basis_set_list_b) 846 847 CALL timestop(handle) 848 849 END SUBROUTINE build_overlap_force 850 851! ************************************************************************************************** 852!> \brief Setup the structure of a sparse matrix based on the overlap 853!> neighbor list 854!> \param ks_env The QS environment 855!> \param matrix_s Matrices to be constructed 856!> \param matrix_name Matrix base name 857!> \param basis_set_list_a Basis set used for <a| 858!> \param basis_set_list_b Basis set used for |b> 859!> \param sab_nl Overlap neighbor list 860!> \param symmetric Is symmetry used in the neighbor list? 861! ************************************************************************************************** 862 SUBROUTINE create_sab_matrix_1d(ks_env, matrix_s, matrix_name, & 863 basis_set_list_a, basis_set_list_b, sab_nl, symmetric) 864 865 TYPE(qs_ks_env_type), POINTER :: ks_env 866 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 867 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name 868 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b 869 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 870 POINTER :: sab_nl 871 LOGICAL, INTENT(IN) :: symmetric 872 873 CHARACTER(len=*), PARAMETER :: routineN = 'create_sab_matrix_1d', & 874 routineP = moduleN//':'//routineN 875 876 CHARACTER(LEN=12) :: cgfsym 877 CHARACTER(LEN=32) :: symmetry_string 878 CHARACTER(LEN=default_string_length) :: mname, name 879 INTEGER :: i, maxs, natom 880 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes 881 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 882 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 883 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 884 885 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, & 886 qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist) 887 888 natom = SIZE(particle_set) 889 890 IF (PRESENT(matrix_name)) THEN 891 mname = matrix_name 892 ELSE 893 mname = "DUMMY" 894 END IF 895 896 maxs = SIZE(matrix_s) 897 898 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom)) 899 900 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, & 901 basis=basis_set_list_a) 902 CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, & 903 basis=basis_set_list_b) 904 905 ! prepare for allocation 906 IF (symmetric) THEN 907 symmetry_string = dbcsr_type_symmetric 908 ELSE 909 symmetry_string = dbcsr_type_no_symmetry 910 END IF 911 912 DO i = 1, maxs 913 IF (symmetric) THEN 914 IF (ndod(i) == 1) THEN 915 ! odd derivatives are anti-symmetric 916 symmetry_string = dbcsr_type_antisymmetric 917 ELSE 918 symmetry_string = dbcsr_type_symmetric 919 END IF 920 ELSE 921 symmetry_string = dbcsr_type_no_symmetry 922 END IF 923 cgfsym = cgf_symbol(1, indco(1:3, i)) 924 IF (i == 1) THEN 925 name = mname 926 ELSE 927 name = TRIM(cgfsym(4:))//" DERIVATIVE OF THE "//TRIM(mname)// & 928 " W.R.T. THE NUCLEAR COORDINATES" 929 END IF 930 CALL compress(name) 931 CALL uppercase(name) 932 ALLOCATE (matrix_s(i)%matrix) 933 CALL dbcsr_create(matrix=matrix_s(i)%matrix, & 934 name=TRIM(name), & 935 dist=dbcsr_dist, matrix_type=symmetry_string, & 936 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, & 937 nze=0) 938 CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i)%matrix, sab_nl) 939 END DO 940 941 DEALLOCATE (row_blk_sizes, col_blk_sizes) 942 943 END SUBROUTINE create_sab_matrix_1d 944 945! ************************************************************************************************** 946!> \brief Setup the structure of a sparse matrix based on the overlap 947!> neighbor list, 2d version 948!> \param ks_env The QS environment 949!> \param matrix_s Matrices to be constructed 950!> \param matrix_name Matrix base name 951!> \param basis_set_list_a Basis set used for <a| 952!> \param basis_set_list_b Basis set used for |b> 953!> \param sab_nl Overlap neighbor list 954!> \param symmetric Is symmetry used in the neighbor list? 955! ************************************************************************************************** 956 SUBROUTINE create_sab_matrix_2d(ks_env, matrix_s, matrix_name, & 957 basis_set_list_a, basis_set_list_b, sab_nl, symmetric) 958 959 TYPE(qs_ks_env_type), POINTER :: ks_env 960 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s 961 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name 962 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b 963 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 964 POINTER :: sab_nl 965 LOGICAL, INTENT(IN) :: symmetric 966 967 CHARACTER(len=*), PARAMETER :: routineN = 'create_sab_matrix_2d', & 968 routineP = moduleN//':'//routineN 969 970 CHARACTER(LEN=12) :: cgfsym 971 CHARACTER(LEN=32) :: symmetry_string 972 CHARACTER(LEN=default_string_length) :: mname, name 973 INTEGER :: i1, i2, natom 974 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes 975 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 976 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 977 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 978 979 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, & 980 qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist) 981 982 natom = SIZE(particle_set) 983 984 IF (PRESENT(matrix_name)) THEN 985 mname = matrix_name 986 ELSE 987 mname = "DUMMY" 988 END IF 989 990 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom)) 991 992 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, & 993 basis=basis_set_list_a) 994 CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, & 995 basis=basis_set_list_b) 996 997 ! prepare for allocation 998 IF (symmetric) THEN 999 symmetry_string = dbcsr_type_symmetric 1000 ELSE 1001 symmetry_string = dbcsr_type_no_symmetry 1002 END IF 1003 1004 DO i2 = 1, SIZE(matrix_s, 2) 1005 DO i1 = 1, SIZE(matrix_s, 1) 1006 IF (symmetric) THEN 1007 IF (ndod(i1) == 1) THEN 1008 ! odd derivatives are anti-symmetric 1009 symmetry_string = dbcsr_type_antisymmetric 1010 ELSE 1011 symmetry_string = dbcsr_type_symmetric 1012 END IF 1013 ELSE 1014 symmetry_string = dbcsr_type_no_symmetry 1015 END IF 1016 cgfsym = cgf_symbol(1, indco(1:3, i1)) 1017 IF (i1 == 1) THEN 1018 name = mname 1019 ELSE 1020 name = TRIM(cgfsym(4:))//" DERIVATIVE OF THE "//TRIM(mname)// & 1021 " W.R.T. THE NUCLEAR COORDINATES" 1022 END IF 1023 CALL compress(name) 1024 CALL uppercase(name) 1025 ALLOCATE (matrix_s(i1, i2)%matrix) 1026 CALL dbcsr_create(matrix=matrix_s(i1, i2)%matrix, & 1027 name=TRIM(name), & 1028 dist=dbcsr_dist, matrix_type=symmetry_string, & 1029 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, & 1030 nze=0) 1031 CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i1, i2)%matrix, sab_nl) 1032 END DO 1033 END DO 1034 1035 DEALLOCATE (row_blk_sizes, col_blk_sizes) 1036 1037 END SUBROUTINE create_sab_matrix_2d 1038 1039! ************************************************************************************************** 1040 1041END MODULE qs_overlap 1042 1043