Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 eps_filter Filter final matrix (optional)
87 : !> \param nderivative The number of calculated derivatives
88 : !> \date 11.10.2010
89 : !> \par History
90 : !> Ported from qs_overlap, replaces code in build_core_hamiltonian
91 : !> Refactoring [07.2014] JGH
92 : !> Simplify options and use new kinetic energy integral routine
93 : !> kpoints [08.2014] JGH
94 : !> Include the derivatives [2021] SL, ED
95 : !> \author JGH
96 : !> \version 1.0
97 : ! **************************************************************************************************
98 18829 : SUBROUTINE build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, &
99 : basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, &
100 : eps_filter, nderivative)
101 :
102 : TYPE(qs_ks_env_type), POINTER :: ks_env
103 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
104 : POINTER :: matrix_t
105 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
106 : POINTER :: matrixkp_t
107 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
108 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
109 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
110 : POINTER :: sab_nl
111 : LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
112 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p
113 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
114 : POINTER :: matrixkp_p
115 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
116 : INTEGER, INTENT(IN), OPTIONAL :: nderivative
117 :
118 : INTEGER :: natom
119 :
120 18829 : CALL get_ks_env(ks_env, natom=natom)
121 :
122 : CALL build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, &
123 : sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, natom, &
124 19569 : nderivative)
125 :
126 18829 : END SUBROUTINE build_kinetic_matrix
127 :
128 : ! **************************************************************************************************
129 : !> \brief Implementation of build_kinetic_matrix with the additional natom argument passed in to
130 : !> allow for explicitly shaped force_thread array which is needed for OMP REDUCTION.
131 : !> \param ks_env ...
132 : !> \param matrix_t ...
133 : !> \param matrixkp_t ...
134 : !> \param matrix_name ...
135 : !> \param basis_type ...
136 : !> \param sab_nl ...
137 : !> \param calculate_forces ...
138 : !> \param matrix_p ...
139 : !> \param matrixkp_p ...
140 : !> \param eps_filter ...
141 : !> \param natom ...
142 : !> \param nderivative ...
143 : ! **************************************************************************************************
144 18829 : SUBROUTINE build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, &
145 : sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, natom, &
146 : nderivative)
147 :
148 : TYPE(qs_ks_env_type), POINTER :: ks_env
149 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
150 : POINTER :: matrix_t
151 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
152 : POINTER :: matrixkp_t
153 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
154 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
155 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
156 : POINTER :: sab_nl
157 : LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
158 : TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p
159 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
160 : POINTER :: matrixkp_p
161 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_filter
162 : INTEGER, INTENT(IN) :: natom
163 : INTEGER, INTENT(IN), OPTIONAL :: nderivative
164 :
165 : CHARACTER(len=*), PARAMETER :: routineN = 'build_kinetic_matrix_low'
166 :
167 : INTEGER :: atom_a, handle, i, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, &
168 : maxder, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
169 18829 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
170 : INTEGER, DIMENSION(3) :: cell
171 18829 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
172 18829 : npgfb, nsgfa, nsgfb
173 18829 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
174 18829 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
175 : LOGICAL :: do_forces, do_symmetric, dokp, found, &
176 : trans, use_cell_mapping, use_virial
177 : REAL(KIND=dp) :: f, f0, ff, rab2, tab
178 18829 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, qab
179 18829 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kab
180 : REAL(KIND=dp), DIMENSION(3) :: force_a, rab
181 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_thread
182 37658 : REAL(KIND=dp), DIMENSION(3, natom) :: force_thread
183 18829 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
184 18829 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
185 18829 : zeta, zetb
186 18829 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dab
187 18829 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
188 18829 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: k_block
189 : TYPE(dbcsr_type), POINTER :: matp
190 : TYPE(dft_control_type), POINTER :: dft_control
191 18829 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
192 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
193 : TYPE(kpoint_type), POINTER :: kpoints
194 18829 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
195 18829 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
196 : TYPE(virial_type), POINTER :: virial
197 :
198 : !$ INTEGER(kind=omp_lock_kind), &
199 18829 : !$ ALLOCATABLE, DIMENSION(:) :: locks
200 : !$ INTEGER :: lock_num, hash, hash1, hash2
201 : !$ INTEGER(KIND=int_8) :: iatom8
202 : !$ INTEGER, PARAMETER :: nlock = 501
203 :
204 : MARK_USED(int_8)
205 :
206 18829 : CALL timeset(routineN, handle)
207 :
208 : ! test for matrices (kpoints or standard gamma point)
209 18829 : IF (PRESENT(matrix_t)) THEN
210 : dokp = .FALSE.
211 : use_cell_mapping = .FALSE.
212 18743 : ELSEIF (PRESENT(matrixkp_t)) THEN
213 18743 : dokp = .TRUE.
214 18743 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
215 18743 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
216 74972 : use_cell_mapping = (SIZE(cell_to_index) > 1)
217 : ELSE
218 0 : CPABORT("")
219 : END IF
220 :
221 18829 : NULLIFY (atomic_kind_set, qs_kind_set, p_block, dft_control)
222 : CALL get_ks_env(ks_env, &
223 : dft_control=dft_control, &
224 : atomic_kind_set=atomic_kind_set, &
225 18829 : qs_kind_set=qs_kind_set)
226 :
227 18829 : nimg = dft_control%nimages
228 18829 : nkind = SIZE(atomic_kind_set)
229 :
230 18829 : do_forces = .FALSE.
231 18829 : IF (PRESENT(calculate_forces)) do_forces = calculate_forces
232 :
233 18829 : nder = 0
234 18829 : IF (PRESENT(nderivative)) nder = nderivative
235 18829 : maxder = ncoset(nder)
236 :
237 : ! check for symmetry
238 18829 : CPASSERT(SIZE(sab_nl) > 0)
239 18829 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
240 :
241 : ! prepare basis set
242 89908 : ALLOCATE (basis_set_list(nkind))
243 18829 : CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
244 :
245 18829 : IF (dokp) THEN
246 18743 : CALL dbcsr_allocate_matrix_set(matrixkp_t, 1, nimg)
247 : CALL create_sab_matrix(ks_env, matrixkp_t, matrix_name, basis_set_list, basis_set_list, &
248 19483 : sab_nl, do_symmetric)
249 : ELSE
250 86 : CALL dbcsr_allocate_matrix_set(matrix_t, maxder)
251 : CALL create_sab_matrix(ks_env, matrix_t, matrix_name, basis_set_list, basis_set_list, &
252 86 : sab_nl, do_symmetric)
253 : END IF
254 :
255 18829 : use_virial = .FALSE.
256 18829 : IF (do_forces) THEN
257 : ! if forces -> maybe virial too
258 7697 : CALL get_ks_env(ks_env=ks_env, force=force, virial=virial)
259 7697 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
260 : ! we need density matrix for forces
261 7697 : IF (dokp) THEN
262 7647 : CPASSERT(PRESENT(matrixkp_p))
263 : ELSE
264 50 : IF (PRESENT(matrix_p)) THEN
265 0 : matp => matrix_p
266 50 : ELSEIF (PRESENT(matrixkp_p)) THEN
267 50 : matp => matrixkp_p(1, 1)%matrix
268 : ELSE
269 0 : CPABORT("Missing density matrix")
270 : END IF
271 : END IF
272 : END IF
273 :
274 296949 : force_thread = 0.0_dp
275 18829 : pv_thread = 0.0_dp
276 :
277 : ! *** Allocate work storage ***
278 18829 : ldsab = get_memory_usage(qs_kind_set, basis_type)
279 :
280 : !$OMP PARALLEL DEFAULT(NONE) &
281 : !$OMP SHARED (do_forces, ldsab, use_cell_mapping, do_symmetric, dokp,&
282 : !$OMP sab_nl, ncoset, maxder, nder, ndod, use_virial, matrix_t, matrixkp_t,&
283 : !$OMP matp, matrix_p, basis_set_list, atom_of_kind, cell_to_index, matrixkp_p, locks, natom) &
284 : !$OMP PRIVATE (k_block, kab, qab, pab, ikind, jkind, iatom, jatom, rab, cell, basis_set_a, basis_set_b, f, &
285 : !$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
286 : !$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, tab, &
287 : !$OMP slot, zetb, scon_a, scon_b, i, ic, irow, icol, f0, ff, found, trans, rab2, sgfa, sgfb, iset, jset, &
288 : !$OMP hash, hash1, hash2, iatom8, lock_num) &
289 18829 : !$OMP REDUCTION (+ : pv_thread, force_thread )
290 :
291 : !$OMP SINGLE
292 : !$ ALLOCATE (locks(nlock))
293 : !$OMP END SINGLE
294 :
295 : !$OMP DO
296 : !$ DO lock_num = 1, nlock
297 : !$ call omp_init_lock(locks(lock_num))
298 : !$ END DO
299 : !$OMP END DO
300 :
301 : ALLOCATE (kab(ldsab, ldsab, maxder), qab(ldsab, ldsab))
302 : IF (do_forces) THEN
303 : ALLOCATE (dab(ldsab, ldsab, 3), pab(ldsab, ldsab))
304 : END IF
305 :
306 : ALLOCATE (k_block(maxder))
307 : DO i = 1, maxder
308 : NULLIFY (k_block(i)%block)
309 : END DO
310 :
311 : !$OMP DO SCHEDULE(GUIDED)
312 : DO slot = 1, sab_nl(1)%nl_size
313 :
314 : ikind = sab_nl(1)%nlist_task(slot)%ikind
315 : jkind = sab_nl(1)%nlist_task(slot)%jkind
316 : iatom = sab_nl(1)%nlist_task(slot)%iatom
317 : jatom = sab_nl(1)%nlist_task(slot)%jatom
318 : cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
319 : rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
320 :
321 : basis_set_a => basis_set_list(ikind)%gto_basis_set
322 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
323 : basis_set_b => basis_set_list(jkind)%gto_basis_set
324 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
325 :
326 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
327 : !$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
328 :
329 : ! basis ikind
330 : first_sgfa => basis_set_a%first_sgf
331 : la_max => basis_set_a%lmax
332 : la_min => basis_set_a%lmin
333 : npgfa => basis_set_a%npgf
334 : nseta = basis_set_a%nset
335 : nsgfa => basis_set_a%nsgf_set
336 : rpgfa => basis_set_a%pgf_radius
337 : set_radius_a => basis_set_a%set_radius
338 : scon_a => basis_set_a%scon
339 : zeta => basis_set_a%zet
340 : ! basis jkind
341 : first_sgfb => basis_set_b%first_sgf
342 : lb_max => basis_set_b%lmax
343 : lb_min => basis_set_b%lmin
344 : npgfb => basis_set_b%npgf
345 : nsetb = basis_set_b%nset
346 : nsgfb => basis_set_b%nsgf_set
347 : rpgfb => basis_set_b%pgf_radius
348 : set_radius_b => basis_set_b%set_radius
349 : scon_b => basis_set_b%scon
350 : zetb => basis_set_b%zet
351 :
352 : IF (use_cell_mapping) THEN
353 : ic = cell_to_index(cell(1), cell(2), cell(3))
354 : CPASSERT(ic > 0)
355 : ELSE
356 : ic = 1
357 : END IF
358 :
359 : IF (do_symmetric) THEN
360 : IF (iatom <= jatom) THEN
361 : irow = iatom
362 : icol = jatom
363 : ELSE
364 : irow = jatom
365 : icol = iatom
366 : END IF
367 : f0 = 2.0_dp
368 : IF (iatom == jatom) f0 = 1.0_dp
369 : ff = 2.0_dp
370 : ELSE
371 : irow = iatom
372 : icol = jatom
373 : f0 = 1.0_dp
374 : ff = 1.0_dp
375 : END IF
376 : IF (dokp) THEN
377 : CALL dbcsr_get_block_p(matrix=matrixkp_t(1, ic)%matrix, &
378 : row=irow, col=icol, BLOCK=k_block(1)%block, found=found)
379 : CPASSERT(found)
380 : ELSE
381 : DO i = 1, maxder
382 : NULLIFY (k_block(i)%block)
383 : CALL dbcsr_get_block_p(matrix=matrix_t(i)%matrix, &
384 : row=irow, col=icol, BLOCK=k_block(i)%block, found=found)
385 : CPASSERT(found)
386 : END DO
387 : END IF
388 :
389 : IF (do_forces) THEN
390 : NULLIFY (p_block)
391 : IF (dokp) THEN
392 : CALL dbcsr_get_block_p(matrix=matrixkp_p(1, ic)%matrix, &
393 : row=irow, col=icol, block=p_block, found=found)
394 : CPASSERT(found)
395 : ELSE
396 : !deb CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
397 : !deb block=p_block, found=found)
398 : CALL dbcsr_get_block_p(matrix=matp, row=irow, col=icol, &
399 : block=p_block, found=found)
400 : CPASSERT(found)
401 : END IF
402 : END IF
403 :
404 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
405 : tab = SQRT(rab2)
406 : trans = do_symmetric .AND. (iatom > jatom)
407 :
408 : DO iset = 1, nseta
409 :
410 : ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
411 : sgfa = first_sgfa(1, iset)
412 :
413 : DO jset = 1, nsetb
414 :
415 : IF (set_radius_a(iset) + set_radius_b(jset) < tab) CYCLE
416 :
417 : !$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
418 : !$ hash = MOD(hash1 + hash2, nlock) + 1
419 :
420 : ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
421 : sgfb = first_sgfb(1, jset)
422 :
423 : IF (do_forces .AND. ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
424 : ! Decontract P matrix block
425 : kab = 0.0_dp
426 : CALL block_add("OUT", kab(:, :, 1), nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
427 : CALL decontraction(kab(:, :, 1), pab, scon_a(:, sgfa:), ncoa, nsgfa(iset), scon_b(:, sgfb:), ncob, nsgfb(jset), &
428 : trans=trans)
429 : ! calculate integrals and derivatives
430 : CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
431 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
432 : rab, kab(:, :, 1), dab)
433 : CALL force_trace(force_a, dab, pab, ncoa, ncob, 3)
434 : force_thread(:, iatom) = force_thread(:, iatom) + ff*force_a(:)
435 : force_thread(:, jatom) = force_thread(:, jatom) - ff*force_a(:)
436 : IF (use_virial) THEN
437 : CALL virial_pair_force(pv_thread, f0, force_a, rab)
438 : END IF
439 : ELSE
440 : ! calclulate integrals
441 : IF (nder == 0) THEN
442 : CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
443 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
444 : rab, kab=kab(:, :, 1))
445 : ELSE IF (nder == 1) THEN
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=kab(:, :, 1), dab=kab(:, :, 2:4))
449 : END IF
450 : END IF
451 : DO i = 1, maxder
452 : f = 1.0_dp
453 : IF (ndod(i) == 1 .AND. trans) f = -1.0_dp
454 : ! Contraction step
455 : CALL contraction(kab(:, :, i), qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), &
456 : cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), fscale=f, &
457 : trans=trans)
458 : !$ CALL omp_set_lock(locks(hash))
459 : CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), k_block(i)%block, sgfa, sgfb, trans=trans)
460 : !$ CALL omp_unset_lock(locks(hash))
461 : END DO
462 : END DO
463 : END DO
464 :
465 : END DO
466 : DEALLOCATE (kab, qab)
467 : IF (do_forces) DEALLOCATE (pab, dab)
468 : !$OMP DO
469 : !$ DO lock_num = 1, nlock
470 : !$ call omp_destroy_lock(locks(lock_num))
471 : !$ END DO
472 : !$OMP END DO
473 :
474 : !$OMP SINGLE
475 : !$ DEALLOCATE (locks)
476 : !$OMP END SINGLE NOWAIT
477 :
478 : !$OMP END PARALLEL
479 :
480 18829 : IF (do_forces) THEN
481 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, &
482 7697 : kind_of=kind_of)
483 :
484 : !$OMP DO
485 : DO iatom = 1, natom
486 26999 : atom_a = atom_of_kind(iatom)
487 26999 : ikind = kind_of(iatom)
488 107996 : force(ikind)%kinetic(:, atom_a) = force(ikind)%kinetic(:, atom_a) + force_thread(:, iatom)
489 : END DO
490 : !$OMP END DO
491 : END IF
492 7697 : IF (do_forces .AND. use_virial) THEN
493 10972 : virial%pv_ekinetic = virial%pv_ekinetic + pv_thread
494 10972 : virial%pv_virial = virial%pv_virial + pv_thread
495 : END IF
496 :
497 18829 : IF (dokp) THEN
498 70490 : DO ic = 1, nimg
499 51747 : CALL dbcsr_finalize(matrixkp_t(1, ic)%matrix)
500 70490 : IF (PRESENT(eps_filter)) THEN
501 49629 : CALL dbcsr_filter(matrixkp_t(1, ic)%matrix, eps_filter)
502 : END IF
503 : END DO
504 : ELSE
505 86 : CALL dbcsr_finalize(matrix_t(1)%matrix)
506 86 : IF (PRESENT(eps_filter)) THEN
507 24 : CALL dbcsr_filter(matrix_t(1)%matrix, eps_filter)
508 : END IF
509 : END IF
510 :
511 18829 : DEALLOCATE (basis_set_list)
512 :
513 18829 : CALL timestop(handle)
514 :
515 56487 : END SUBROUTINE build_kinetic_matrix_low
516 :
517 : END MODULE qs_kinetic
518 :
|
