Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of kinetic energy matrix and forces
10 : !> \par History
11 : !> JGH: from core_hamiltonian
12 : !> simplify further [7.2014]
13 : !> \author Juerg Hutter
14 : ! **************************************************************************************************
15 : MODULE qs_kinetic
16 : USE ai_contraction, ONLY: block_add,&
17 : contraction,&
18 : decontraction,&
19 : force_trace
20 : USE ai_kinetic, ONLY: kinetic
21 : USE atomic_kind_types, ONLY: atomic_kind_type,&
22 : get_atomic_kind_set
23 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
24 : gto_basis_set_type
25 : USE block_p_types, ONLY: block_p_type
26 : USE cp_control_types, ONLY: dft_control_type
27 : USE cp_dbcsr_api, ONLY: dbcsr_filter,&
28 : dbcsr_finalize,&
29 : dbcsr_get_block_p,&
30 : dbcsr_p_type,&
31 : dbcsr_type
32 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
33 : USE kinds, ONLY: dp,&
34 : int_8
35 : USE kpoint_types, ONLY: get_kpoint_info,&
36 : kpoint_type
37 : USE orbital_pointers, ONLY: ncoset
38 : USE qs_force_types, ONLY: qs_force_type
39 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
40 : get_memory_usage
41 : USE qs_kind_types, ONLY: qs_kind_type
42 : USE qs_ks_types, ONLY: get_ks_env,&
43 : qs_ks_env_type
44 : USE qs_neighbor_list_types, ONLY: get_neighbor_list_set_p,&
45 : neighbor_list_set_p_type
46 : USE qs_overlap, ONLY: create_sab_matrix
47 : USE virial_methods, ONLY: virial_pair_force
48 : USE virial_types, ONLY: virial_type
49 :
50 : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
51 : !$ omp_init_lock, omp_set_lock, &
52 : !$ omp_unset_lock, omp_destroy_lock
53 :
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 :
60 : ! *** Global parameters ***
61 :
62 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kinetic'
63 :
64 : ! *** Public subroutines ***
65 :
66 : PUBLIC :: build_kinetic_matrix
67 :
68 : INTEGER, DIMENSION(1:56), SAVE :: ndod = [0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, &
69 : 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, &
70 : 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, &
71 : 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
77 : !> \param ks_env the QS environment
78 : !> \param matrix_t The kinetic energy matrix to be calculated (optional)
79 : !> \param matrixkp_t The kinetic energy matrices to be calculated (kpoints,optional)
80 : !> \param matrix_name The name of the matrix (i.e. for output)
81 : !> \param basis_type basis set to be used
82 : !> \param sab_nl pair list (must be consistent with basis sets!)
83 : !> \param calculate_forces (optional) ...
84 : !> \param matrix_p density matrix for force calculation (optional)
85 : !> \param matrixkp_p density matrix for force calculation with kpoints (optional)
86 : !> \param ext_kpoints ...
87 : !> \param eps_filter Filter final matrix (optional)
88 : !> \param nderivative The number of calculated derivatives
89 : !> \date 11.10.2010
90 : !> \par History
91 : !> Ported from qs_overlap, replaces code in build_core_hamiltonian
92 : !> Refactoring [07.2014] JGH
93 : !> Simplify options and use new kinetic energy integral routine
94 : !> kpoints [08.2014] JGH
95 : !> Include the derivatives [2021] SL, ED
96 : !> \author JGH
97 : !> \version 1.0
98 : ! **************************************************************************************************
99 19583 : SUBROUTINE build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, &
100 : basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, &
101 : ext_kpoints, eps_filter, nderivative)
102 :
103 : TYPE(qs_ks_env_type), POINTER :: ks_env
104 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
105 : POINTER :: matrix_t
106 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
107 : POINTER :: matrixkp_t
108 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
109 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
110 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
111 : POINTER :: sab_nl
112 : LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
113 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p
114 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
115 : POINTER :: matrixkp_p
116 : TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
117 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
118 : INTEGER, INTENT(IN), OPTIONAL :: nderivative
119 :
120 : INTEGER :: natom
121 :
122 19583 : CALL get_ks_env(ks_env, natom=natom)
123 :
124 : CALL build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, &
125 : sab_nl, calculate_forces, matrix_p, matrixkp_p, &
126 20341 : ext_kpoints, eps_filter, natom, nderivative)
127 :
128 19583 : END SUBROUTINE build_kinetic_matrix
129 :
130 : ! **************************************************************************************************
131 : !> \brief Implementation of build_kinetic_matrix with the additional natom argument passed in to
132 : !> allow for explicitly shaped force_thread array which is needed for OMP REDUCTION.
133 : !> \param ks_env ...
134 : !> \param matrix_t ...
135 : !> \param matrixkp_t ...
136 : !> \param matrix_name ...
137 : !> \param basis_type ...
138 : !> \param sab_nl ...
139 : !> \param calculate_forces ...
140 : !> \param matrix_p ...
141 : !> \param matrixkp_p ...
142 : !> \param ext_kpoints ...
143 : !> \param eps_filter ...
144 : !> \param natom ...
145 : !> \param nderivative ...
146 : ! **************************************************************************************************
147 19583 : SUBROUTINE build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, &
148 : sab_nl, calculate_forces, matrix_p, matrixkp_p, &
149 : ext_kpoints, eps_filter, natom, nderivative)
150 :
151 : TYPE(qs_ks_env_type), POINTER :: ks_env
152 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
153 : POINTER :: matrix_t
154 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
155 : POINTER :: matrixkp_t
156 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
157 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
158 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
159 : POINTER :: sab_nl
160 : LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
161 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p
162 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
163 : POINTER :: matrixkp_p
164 : TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
165 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
166 : INTEGER, INTENT(IN) :: natom
167 : INTEGER, INTENT(IN), OPTIONAL :: nderivative
168 :
169 : CHARACTER(len=*), PARAMETER :: routineN = 'build_kinetic_matrix_low'
170 :
171 : INTEGER :: atom_a, handle, i, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, &
172 : maxder, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
173 19583 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
174 : INTEGER, DIMENSION(3) :: cell
175 19583 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
176 19583 : npgfb, nsgfa, nsgfb
177 19583 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
178 19583 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
179 : LOGICAL :: do_forces, do_symmetric, dokp, found, &
180 : trans, use_cell_mapping, use_virial
181 : REAL(KIND=dp) :: f, f0, ff, rab2, tab
182 19583 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, qab
183 19583 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kab
184 : REAL(KIND=dp), DIMENSION(3) :: force_a, rab
185 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
186 39166 : REAL(KIND=dp), DIMENSION(3, natom) :: force_thread
187 19583 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
188 19583 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
189 19583 : zeta, zetb
190 19583 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dab
191 19583 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
192 19583 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: k_block
193 : TYPE(dbcsr_type), POINTER :: matp
194 : TYPE(dft_control_type), POINTER :: dft_control
195 19583 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
196 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
197 : TYPE(kpoint_type), POINTER :: kpoints
198 19583 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
199 19583 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
200 : TYPE(virial_type), POINTER :: virial
201 :
202 : !$ INTEGER(kind=omp_lock_kind), &
203 19583 : !$ ALLOCATABLE, DIMENSION(:) :: locks
204 : !$ INTEGER :: lock_num, hash, hash1, hash2
205 : !$ INTEGER(KIND=int_8) :: iatom8
206 : !$ INTEGER, PARAMETER :: nlock = 501
207 :
208 : MARK_USED(int_8)
209 :
210 19583 : CALL timeset(routineN, handle)
211 :
212 19583 : NULLIFY (atomic_kind_set, qs_kind_set, p_block, dft_control)
213 : CALL get_ks_env(ks_env, &
214 : dft_control=dft_control, &
215 : atomic_kind_set=atomic_kind_set, &
216 19583 : qs_kind_set=qs_kind_set)
217 :
218 : ! test for matrices (kpoints or standard gamma point)
219 19583 : IF (PRESENT(matrix_t)) THEN
220 86 : dokp = .FALSE.
221 86 : use_cell_mapping = .FALSE.
222 86 : nimg = dft_control%nimages
223 19497 : ELSEIF (PRESENT(matrixkp_t)) THEN
224 19497 : dokp = .TRUE.
225 19497 : IF (PRESENT(ext_kpoints)) THEN
226 0 : IF (ASSOCIATED(ext_kpoints)) THEN
227 0 : CALL get_kpoint_info(kpoint=ext_kpoints, cell_to_index=cell_to_index)
228 0 : use_cell_mapping = (SIZE(cell_to_index) > 1)
229 0 : nimg = SIZE(ext_kpoints%index_to_cell, 2)
230 : ELSE
231 0 : use_cell_mapping = .FALSE.
232 0 : nimg = 1
233 : END IF
234 : ELSE
235 19497 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
236 19497 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
237 77988 : use_cell_mapping = (SIZE(cell_to_index) > 1)
238 19497 : nimg = dft_control%nimages
239 : END IF
240 : ELSE
241 0 : CPABORT("")
242 : END IF
243 :
244 19583 : nkind = SIZE(atomic_kind_set)
245 :
246 19583 : do_forces = .FALSE.
247 19583 : IF (PRESENT(calculate_forces)) do_forces = calculate_forces
248 :
249 19583 : nder = 0
250 19583 : IF (PRESENT(nderivative)) nder = nderivative
251 19583 : maxder = ncoset(nder)
252 :
253 : ! check for symmetry
254 19583 : CPASSERT(SIZE(sab_nl) > 0)
255 19583 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
256 :
257 : ! prepare basis set
258 93628 : ALLOCATE (basis_set_list(nkind))
259 19583 : CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
260 :
261 19583 : IF (dokp) THEN
262 19497 : CALL dbcsr_allocate_matrix_set(matrixkp_t, 1, nimg)
263 : CALL create_sab_matrix(ks_env, matrixkp_t, matrix_name, basis_set_list, basis_set_list, &
264 20255 : sab_nl, do_symmetric)
265 : ELSE
266 86 : CALL dbcsr_allocate_matrix_set(matrix_t, maxder)
267 : CALL create_sab_matrix(ks_env, matrix_t, matrix_name, basis_set_list, basis_set_list, &
268 86 : sab_nl, do_symmetric)
269 : END IF
270 :
271 19583 : use_virial = .FALSE.
272 19583 : IF (do_forces) THEN
273 : ! if forces -> maybe virial too
274 7973 : CALL get_ks_env(ks_env=ks_env, force=force, virial=virial)
275 7973 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
276 : ! we need density matrix for forces
277 7973 : IF (dokp) THEN
278 7923 : CPASSERT(PRESENT(matrixkp_p))
279 : ELSE
280 50 : IF (PRESENT(matrix_p)) THEN
281 0 : matp => matrix_p
282 50 : ELSEIF (PRESENT(matrixkp_p)) THEN
283 50 : matp => matrixkp_p(1, 1)%matrix
284 : ELSE
285 0 : CPABORT("Missing density matrix")
286 : END IF
287 : END IF
288 : END IF
289 :
290 307655 : force_thread = 0.0_dp
291 19583 : pv_thread = 0.0_dp
292 :
293 : ! *** Allocate work storage ***
294 19583 : ldsab = get_memory_usage(qs_kind_set, basis_type)
295 :
296 : !$OMP PARALLEL DEFAULT(NONE) &
297 : !$OMP SHARED (do_forces, ldsab, use_cell_mapping, do_symmetric, dokp,&
298 : !$OMP sab_nl, ncoset, maxder, nder, ndod, use_virial, matrix_t, matrixkp_t,&
299 : !$OMP matp, matrix_p, basis_set_list, atom_of_kind, cell_to_index, matrixkp_p, locks, natom) &
300 : !$OMP PRIVATE (k_block, kab, qab, pab, ikind, jkind, iatom, jatom, rab, cell, basis_set_a, basis_set_b, f, &
301 : !$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
302 : !$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, tab, &
303 : !$OMP slot, zetb, scon_a, scon_b, i, ic, irow, icol, f0, ff, found, trans, rab2, sgfa, sgfb, iset, jset, &
304 : !$OMP hash, hash1, hash2, iatom8, lock_num) &
305 19583 : !$OMP REDUCTION (+ : pv_thread, force_thread )
306 :
307 : !$OMP SINGLE
308 : !$ ALLOCATE (locks(nlock))
309 : !$OMP END SINGLE
310 :
311 : !$OMP DO
312 : !$ DO lock_num = 1, nlock
313 : !$ call omp_init_lock(locks(lock_num))
314 : !$ END DO
315 : !$OMP END DO
316 :
317 : ALLOCATE (kab(ldsab, ldsab, maxder), qab(ldsab, ldsab))
318 : IF (do_forces) THEN
319 : ALLOCATE (dab(ldsab, ldsab, 3), pab(ldsab, ldsab))
320 : END IF
321 :
322 : ALLOCATE (k_block(maxder))
323 : DO i = 1, maxder
324 : NULLIFY (k_block(i)%block)
325 : END DO
326 :
327 : !$OMP DO SCHEDULE(GUIDED)
328 : DO slot = 1, sab_nl(1)%nl_size
329 :
330 : ikind = sab_nl(1)%nlist_task(slot)%ikind
331 : jkind = sab_nl(1)%nlist_task(slot)%jkind
332 : iatom = sab_nl(1)%nlist_task(slot)%iatom
333 : jatom = sab_nl(1)%nlist_task(slot)%jatom
334 : cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
335 : rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
336 :
337 : basis_set_a => basis_set_list(ikind)%gto_basis_set
338 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
339 : basis_set_b => basis_set_list(jkind)%gto_basis_set
340 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
341 :
342 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
343 : !$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
344 :
345 : ! basis ikind
346 : first_sgfa => basis_set_a%first_sgf
347 : la_max => basis_set_a%lmax
348 : la_min => basis_set_a%lmin
349 : npgfa => basis_set_a%npgf
350 : nseta = basis_set_a%nset
351 : nsgfa => basis_set_a%nsgf_set
352 : rpgfa => basis_set_a%pgf_radius
353 : set_radius_a => basis_set_a%set_radius
354 : scon_a => basis_set_a%scon
355 : zeta => basis_set_a%zet
356 : ! basis jkind
357 : first_sgfb => basis_set_b%first_sgf
358 : lb_max => basis_set_b%lmax
359 : lb_min => basis_set_b%lmin
360 : npgfb => basis_set_b%npgf
361 : nsetb = basis_set_b%nset
362 : nsgfb => basis_set_b%nsgf_set
363 : rpgfb => basis_set_b%pgf_radius
364 : set_radius_b => basis_set_b%set_radius
365 : scon_b => basis_set_b%scon
366 : zetb => basis_set_b%zet
367 :
368 : IF (use_cell_mapping) THEN
369 : ic = cell_to_index(cell(1), cell(2), cell(3))
370 : CPASSERT(ic > 0)
371 : ELSE
372 : ic = 1
373 : END IF
374 :
375 : IF (do_symmetric) THEN
376 : IF (iatom <= jatom) THEN
377 : irow = iatom
378 : icol = jatom
379 : ELSE
380 : irow = jatom
381 : icol = iatom
382 : END IF
383 : f0 = 2.0_dp
384 : IF (iatom == jatom) f0 = 1.0_dp
385 : ff = 2.0_dp
386 : ELSE
387 : irow = iatom
388 : icol = jatom
389 : f0 = 1.0_dp
390 : ff = 1.0_dp
391 : END IF
392 : IF (dokp) THEN
393 : CALL dbcsr_get_block_p(matrix=matrixkp_t(1, ic)%matrix, &
394 : row=irow, col=icol, BLOCK=k_block(1)%block, found=found)
395 : CPASSERT(found)
396 : ELSE
397 : DO i = 1, maxder
398 : NULLIFY (k_block(i)%block)
399 : CALL dbcsr_get_block_p(matrix=matrix_t(i)%matrix, &
400 : row=irow, col=icol, BLOCK=k_block(i)%block, found=found)
401 : CPASSERT(found)
402 : END DO
403 : END IF
404 :
405 : IF (do_forces) THEN
406 : NULLIFY (p_block)
407 : IF (dokp) THEN
408 : CALL dbcsr_get_block_p(matrix=matrixkp_p(1, ic)%matrix, &
409 : row=irow, col=icol, block=p_block, found=found)
410 : CPASSERT(found)
411 : ELSE
412 : !deb CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
413 : !deb block=p_block, found=found)
414 : CALL dbcsr_get_block_p(matrix=matp, row=irow, col=icol, &
415 : block=p_block, found=found)
416 : CPASSERT(found)
417 : END IF
418 : END IF
419 :
420 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
421 : tab = SQRT(rab2)
422 : trans = do_symmetric .AND. (iatom > jatom)
423 :
424 : DO iset = 1, nseta
425 :
426 : ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
427 : sgfa = first_sgfa(1, iset)
428 :
429 : DO jset = 1, nsetb
430 :
431 : IF (set_radius_a(iset) + set_radius_b(jset) < tab) CYCLE
432 :
433 : !$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
434 : !$ hash = MOD(hash1 + hash2, nlock) + 1
435 :
436 : ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
437 : sgfb = first_sgfb(1, jset)
438 :
439 : IF (do_forces .AND. ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
440 : ! Decontract P matrix block
441 : kab = 0.0_dp
442 : CALL block_add("OUT", kab(:, :, 1), nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
443 : CALL decontraction(kab(:, :, 1), pab, scon_a(:, sgfa:), ncoa, nsgfa(iset), scon_b(:, sgfb:), ncob, nsgfb(jset), &
444 : trans=trans)
445 : ! calculate integrals and derivatives
446 : CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
447 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
448 : rab, kab(:, :, 1), dab)
449 : CALL force_trace(force_a, dab, pab, ncoa, ncob, 3)
450 : force_thread(:, iatom) = force_thread(:, iatom) + ff*force_a(:)
451 : force_thread(:, jatom) = force_thread(:, jatom) - ff*force_a(:)
452 : IF (use_virial) THEN
453 : CALL virial_pair_force(pv_thread, f0, force_a, rab)
454 : END IF
455 : ELSE
456 : ! calclulate integrals
457 : IF (nder == 0) THEN
458 : CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
459 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
460 : rab, kab=kab(:, :, 1))
461 : ELSE IF (nder == 1) THEN
462 : CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
463 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
464 : rab, kab=kab(:, :, 1), dab=kab(:, :, 2:4))
465 : END IF
466 : END IF
467 : DO i = 1, maxder
468 : f = 1.0_dp
469 : IF (ndod(i) == 1 .AND. trans) f = -1.0_dp
470 : ! Contraction step
471 : CALL contraction(kab(:, :, i), qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), &
472 : cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), fscale=f, &
473 : trans=trans)
474 : !$ CALL omp_set_lock(locks(hash))
475 : CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), k_block(i)%block, sgfa, sgfb, trans=trans)
476 : !$ CALL omp_unset_lock(locks(hash))
477 : END DO
478 : END DO
479 : END DO
480 :
481 : END DO
482 : DEALLOCATE (kab, qab)
483 : IF (do_forces) DEALLOCATE (pab, dab)
484 : !$OMP DO
485 : !$ DO lock_num = 1, nlock
486 : !$ call omp_destroy_lock(locks(lock_num))
487 : !$ END DO
488 : !$OMP END DO
489 :
490 : !$OMP SINGLE
491 : !$ DEALLOCATE (locks)
492 : !$OMP END SINGLE NOWAIT
493 :
494 : !$OMP END PARALLEL
495 :
496 19583 : IF (do_forces) THEN
497 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, &
498 7973 : kind_of=kind_of)
499 :
500 : !$OMP DO
501 : DO iatom = 1, natom
502 27939 : atom_a = atom_of_kind(iatom)
503 27939 : ikind = kind_of(iatom)
504 111756 : force(ikind)%kinetic(:, atom_a) = force(ikind)%kinetic(:, atom_a) + force_thread(:, iatom)
505 : END DO
506 : !$OMP END DO
507 : END IF
508 7973 : IF (do_forces .AND. use_virial) THEN
509 10972 : virial%pv_ekinetic = virial%pv_ekinetic + pv_thread
510 10972 : virial%pv_virial = virial%pv_virial + pv_thread
511 : END IF
512 :
513 19583 : IF (dokp) THEN
514 88892 : DO ic = 1, nimg
515 69395 : CALL dbcsr_finalize(matrixkp_t(1, ic)%matrix)
516 88892 : IF (PRESENT(eps_filter)) THEN
517 67159 : CALL dbcsr_filter(matrixkp_t(1, ic)%matrix, eps_filter)
518 : END IF
519 : END DO
520 : ELSE
521 86 : CALL dbcsr_finalize(matrix_t(1)%matrix)
522 86 : IF (PRESENT(eps_filter)) THEN
523 24 : CALL dbcsr_filter(matrix_t(1)%matrix, eps_filter)
524 : END IF
525 : END IF
526 :
527 19583 : DEALLOCATE (basis_set_list)
528 :
529 19583 : CALL timestop(handle)
530 :
531 58749 : END SUBROUTINE build_kinetic_matrix_low
532 :
533 : END MODULE qs_kinetic
534 :
|
