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 kinetic energy matrix and forces 8!> \par History 9!> JGH: from core_hamiltonian 10!> simplify further [7.2014] 11!> \author Juerg Hutter 12! ************************************************************************************************** 13MODULE qs_kinetic 14 USE ai_contraction, ONLY: block_add,& 15 contraction,& 16 decontraction,& 17 force_trace 18 USE ai_kinetic, ONLY: kinetic 19 USE atomic_kind_types, ONLY: atomic_kind_type,& 20 get_atomic_kind_set 21 USE basis_set_types, ONLY: gto_basis_set_p_type,& 22 gto_basis_set_type 23 USE cp_control_types, ONLY: dft_control_type 24 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set 25 USE dbcsr_api, ONLY: dbcsr_filter,& 26 dbcsr_finalize,& 27 dbcsr_get_block_p,& 28 dbcsr_p_type,& 29 dbcsr_type 30 USE kinds, ONLY: dp 31 USE kpoint_types, ONLY: get_kpoint_info,& 32 kpoint_type 33 USE orbital_pointers, ONLY: ncoset 34 USE qs_force_types, ONLY: qs_force_type 35 USE qs_integral_utils, ONLY: basis_set_list_setup,& 36 get_memory_usage 37 USE qs_kind_types, ONLY: qs_kind_type 38 USE qs_ks_types, ONLY: get_ks_env,& 39 qs_ks_env_type 40 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 41 get_neighbor_list_set_p,& 42 neighbor_list_iterate,& 43 neighbor_list_iterator_create,& 44 neighbor_list_iterator_p_type,& 45 neighbor_list_iterator_release,& 46 neighbor_list_set_p_type 47 USE qs_overlap, ONLY: create_sab_matrix 48 USE virial_methods, ONLY: virial_pair_force 49 USE virial_types, ONLY: virial_type 50 51!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num 52#include "./base/base_uses.f90" 53 54 IMPLICIT NONE 55 56 PRIVATE 57 58! *** Global parameters *** 59 60 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kinetic' 61 62! *** Public subroutines *** 63 64 PUBLIC :: build_kinetic_matrix 65 66CONTAINS 67 68! ************************************************************************************************** 69!> \brief Calculation of the kinetic energy matrix over Cartesian Gaussian functions. 70!> \param ks_env the QS environment 71!> \param matrix_t The kinetic energy matrix to be calculated (optional) 72!> \param matrixkp_t The kinetic energy matrices to be calculated (kpoints,optional) 73!> \param matrix_name The name of the matrix (i.e. for output) 74!> \param basis_type basis set to be used 75!> \param sab_nl pair list (must be consistent with basis sets!) 76!> \param calculate_forces (optional) ... 77!> \param matrix_p density matrix for force calculation (optional) 78!> \param matrixkp_p density matrix for force calculation with kpoints (optional) 79!> \param eps_filter Filter final matrix (optional) 80!> \date 11.10.2010 81!> \par History 82!> Ported from qs_overlap, replaces code in build_core_hamiltonian 83!> Refactoring [07.2014] JGH 84!> Simplify options and use new kinetic energy integral routine 85!> kpoints [08.2014] JGH 86!> \author JGH 87!> \version 1.0 88! ************************************************************************************************** 89 SUBROUTINE build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, & 90 basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, & 91 eps_filter) 92 93 TYPE(qs_ks_env_type), POINTER :: ks_env 94 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & 95 POINTER :: matrix_t 96 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & 97 POINTER :: matrixkp_t 98 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name 99 CHARACTER(LEN=*), INTENT(IN) :: basis_type 100 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 101 POINTER :: sab_nl 102 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces 103 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p 104 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, & 105 POINTER :: matrixkp_p 106 REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter 107 108 CHARACTER(len=*), PARAMETER :: routineN = 'build_kinetic_matrix', & 109 routineP = moduleN//':'//routineN 110 111 INTEGER :: atom_a, atom_b, handle, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, & 112 ldsab, mepos, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, nthread, sgfa, sgfb 113 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind 114 INTEGER, DIMENSION(3) :: cell 115 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 116 npgfb, nsgfa, nsgfb 117 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 118 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 119 LOGICAL :: do_forces, do_symmetric, dokp, found, & 120 trans, use_cell_mapping, use_virial 121 REAL(KIND=dp) :: f0, ff, rab2, tab 122 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: kab, pab, qab 123 REAL(KIND=dp), DIMENSION(3) :: force_a, rab 124 REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc 125 REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 126 REAL(KIND=dp), DIMENSION(:, :), POINTER :: k_block, p_block, rpgfa, rpgfb, scon_a, & 127 scon_b, zeta, zetb 128 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dab 129 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 130 TYPE(dft_control_type), POINTER :: dft_control 131 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list 132 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b 133 TYPE(kpoint_type), POINTER :: kpoints 134 TYPE(neighbor_list_iterator_p_type), & 135 DIMENSION(:), POINTER :: nl_iterator 136 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 137 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 138 TYPE(virial_type), POINTER :: virial 139 140 CALL timeset(routineN, handle) 141 142 ! test for matrices (kpoints or standard gamma point) 143 IF (PRESENT(matrix_t)) THEN 144 dokp = .FALSE. 145 use_cell_mapping = .FALSE. 146 ELSEIF (PRESENT(matrixkp_t)) THEN 147 dokp = .TRUE. 148 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) 149 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 150 use_cell_mapping = (SIZE(cell_to_index) > 1) 151 ELSE 152 CPABORT("") 153 END IF 154 155 NULLIFY (atomic_kind_set, qs_kind_set, p_block, dft_control) 156 CALL get_ks_env(ks_env, & 157 dft_control=dft_control, & 158 atomic_kind_set=atomic_kind_set, & 159 natom=natom, & 160 qs_kind_set=qs_kind_set) 161 162 nimg = dft_control%nimages 163 nkind = SIZE(atomic_kind_set) 164 165 ALLOCATE (atom_of_kind(natom)) 166 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind) 167 168 do_forces = .FALSE. 169 IF (PRESENT(calculate_forces)) do_forces = calculate_forces 170 171 ! check for symmetry 172 CPASSERT(SIZE(sab_nl) > 0) 173 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric) 174 175 ! prepare basis set 176 ALLOCATE (basis_set_list(nkind)) 177 CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set) 178 179 IF (dokp) THEN 180 CALL dbcsr_allocate_matrix_set(matrixkp_t, 1, nimg) 181 CALL create_sab_matrix(ks_env, matrixkp_t, matrix_name, basis_set_list, basis_set_list, & 182 sab_nl, do_symmetric) 183 ELSE 184 CALL dbcsr_allocate_matrix_set(matrix_t, 1) 185 CALL create_sab_matrix(ks_env, matrix_t, matrix_name, basis_set_list, basis_set_list, & 186 sab_nl, do_symmetric) 187 END IF 188 189 IF (do_forces) THEN 190 ! if forces -> maybe virial too 191 CALL get_ks_env(ks_env=ks_env, force=force, virial=virial) 192 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 193 pv_loc = virial%pv_virial 194 ! we need density matrix for forces 195 IF (dokp) THEN 196 CPASSERT(PRESENT(matrixkp_p)) 197 ELSE 198 CPASSERT(PRESENT(matrix_p)) 199 END IF 200 END IF 201 ! *** Allocate work storage *** 202 ldsab = get_memory_usage(qs_kind_set, basis_type) 203 204 nthread = 1 205!$ nthread = omp_get_max_threads() 206 ! Iterate of neighbor list 207 CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread) 208 209!$OMP PARALLEL DEFAULT(NONE) & 210!$OMP SHARED (nthread,do_forces,ldsab,nl_iterator, use_cell_mapping, do_symmetric,dokp,& 211!$OMP ncoset,use_virial,force,virial,matrix_t,matrixkp_t,& 212!$OMP matrix_p,basis_set_list, atom_of_kind, cell_to_index, matrixkp_p) & 213!$OMP PRIVATE (k_block,mepos,kab,qab,pab,ikind,jkind,iatom,jatom,rab,cell,basis_set_a,basis_set_b,atom_a,atom_b,& 214!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, & 215!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, tab, & 216!$OMP zetb, scon_a, scon_b, ic, irow, icol, f0, ff, found, trans, rab2, sgfa, sgfb, iset, jset ) 217 218 mepos = 0 219!$ mepos = omp_get_thread_num() 220 221 ALLOCATE (kab(ldsab, ldsab), qab(ldsab, ldsab)) 222 IF (do_forces) THEN 223 ALLOCATE (dab(ldsab, ldsab, 3), pab(ldsab, ldsab)) 224 END IF 225 226 DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0) 227 CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, & 228 iatom=iatom, jatom=jatom, r=rab, cell=cell) 229 atom_a = atom_of_kind(iatom) 230 atom_b = atom_of_kind(jatom) 231 basis_set_a => basis_set_list(ikind)%gto_basis_set 232 IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE 233 basis_set_b => basis_set_list(jkind)%gto_basis_set 234 IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE 235 ! basis ikind 236 first_sgfa => basis_set_a%first_sgf 237 la_max => basis_set_a%lmax 238 la_min => basis_set_a%lmin 239 npgfa => basis_set_a%npgf 240 nseta = basis_set_a%nset 241 nsgfa => basis_set_a%nsgf_set 242 rpgfa => basis_set_a%pgf_radius 243 set_radius_a => basis_set_a%set_radius 244 scon_a => basis_set_a%scon 245 zeta => basis_set_a%zet 246 ! basis jkind 247 first_sgfb => basis_set_b%first_sgf 248 lb_max => basis_set_b%lmax 249 lb_min => basis_set_b%lmin 250 npgfb => basis_set_b%npgf 251 nsetb = basis_set_b%nset 252 nsgfb => basis_set_b%nsgf_set 253 rpgfb => basis_set_b%pgf_radius 254 set_radius_b => basis_set_b%set_radius 255 scon_b => basis_set_b%scon 256 zetb => basis_set_b%zet 257 258 IF (use_cell_mapping) THEN 259 ic = cell_to_index(cell(1), cell(2), cell(3)) 260 CPASSERT(ic > 0) 261 ELSE 262 ic = 1 263 END IF 264 265 IF (do_symmetric) THEN 266 IF (iatom <= jatom) THEN 267 irow = iatom 268 icol = jatom 269 ELSE 270 irow = jatom 271 icol = iatom 272 END IF 273 f0 = 2.0_dp 274 IF (iatom == jatom) f0 = 1.0_dp 275 ff = 2.0_dp 276 ELSE 277 irow = iatom 278 icol = jatom 279 f0 = 1.0_dp 280 ff = 1.0_dp 281 END IF 282 NULLIFY (k_block) 283 IF (dokp) THEN 284 CALL dbcsr_get_block_p(matrix=matrixkp_t(1, ic)%matrix, & 285 row=irow, col=icol, BLOCK=k_block, found=found) 286 CPASSERT(found) 287 ELSE 288 CALL dbcsr_get_block_p(matrix=matrix_t(1)%matrix, & 289 row=irow, col=icol, BLOCK=k_block, found=found) 290 CPASSERT(found) 291 END IF 292 293 IF (do_forces) THEN 294 NULLIFY (p_block) 295 IF (dokp) THEN 296 CALL dbcsr_get_block_p(matrix=matrixkp_p(1, ic)%matrix, & 297 row=irow, col=icol, block=p_block, found=found) 298 CPASSERT(found) 299 ELSE 300 CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, & 301 block=p_block, found=found) 302 CPASSERT(found) 303 END IF 304 END IF 305 306 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 307 tab = SQRT(rab2) 308 trans = do_symmetric .AND. (iatom > jatom) 309 310 DO iset = 1, nseta 311 312 ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1)) 313 sgfa = first_sgfa(1, iset) 314 315 DO jset = 1, nsetb 316 317 IF (set_radius_a(iset) + set_radius_b(jset) < tab) CYCLE 318 319 ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)) 320 sgfb = first_sgfb(1, jset) 321 322 IF (do_forces .AND. ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN 323 ! Decontract P matrix block 324 kab = 0.0_dp 325 CALL block_add("OUT", kab, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans) 326 CALL decontraction(kab, pab, scon_a(:, sgfa:), ncoa, nsgfa(iset), scon_b(:, sgfb:), ncob, nsgfb(jset), & 327 trans=trans) 328 ! calculate integrals and derivatives 329 CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 330 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 331 rab, kab, dab) 332 CALL force_trace(force_a, dab, pab, ncoa, ncob, 3) 333!$OMP CRITICAL(forceupate) 334 force(ikind)%kinetic(:, atom_a) = force(ikind)%kinetic(:, atom_a) + ff*force_a(:) 335 force(jkind)%kinetic(:, atom_b) = force(jkind)%kinetic(:, atom_b) - ff*force_a(:) 336 IF (use_virial) THEN 337 CALL virial_pair_force(virial%pv_virial, f0, force_a, rab) 338 END IF 339!$OMP END CRITICAL(forceupate) 340 ELSE 341 ! calclulate integrals 342 CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), & 343 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), & 344 rab, kab) 345 END IF 346 ! Contraction step 347 CALL contraction(kab, qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), & 348 cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), & 349 trans=trans) 350!$OMP CRITICAL(blockadd) 351 CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), k_block, sgfa, sgfb, trans=trans) 352!$OMP END CRITICAL(blockadd) 353 354 END DO 355 END DO 356 357 END DO 358 DEALLOCATE (kab, qab) 359 IF (do_forces) DEALLOCATE (pab, dab) 360!$OMP END PARALLEL 361 CALL neighbor_list_iterator_release(nl_iterator) 362 363 IF (do_forces .AND. use_virial) THEN 364 virial%pv_ekin = virial%pv_virial - pv_loc 365 END IF 366 367 IF (dokp) THEN 368 DO ic = 1, nimg 369 CALL dbcsr_finalize(matrixkp_t(1, ic)%matrix) 370 IF (PRESENT(eps_filter)) THEN 371 CALL dbcsr_filter(matrixkp_t(1, ic)%matrix, eps_filter) 372 END IF 373 END DO 374 ELSE 375 CALL dbcsr_finalize(matrix_t(1)%matrix) 376 IF (PRESENT(eps_filter)) THEN 377 CALL dbcsr_filter(matrix_t(1)%matrix, eps_filter) 378 END IF 379 END IF 380 381 ! Release work storage 382 DEALLOCATE (atom_of_kind) 383 DEALLOCATE (basis_set_list) 384 385 CALL timestop(handle) 386 387 END SUBROUTINE build_kinetic_matrix 388 389END MODULE qs_kinetic 390 391