Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_operations, ONLY: dbcsr_allocate_matrix_set
28 : USE dbcsr_api, ONLY: dbcsr_filter,&
29 : dbcsr_finalize,&
30 : dbcsr_get_block_p,&
31 : dbcsr_p_type,&
32 : dbcsr_type
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 17989 : 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 17989 : 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 47890 : nderivative)
125 :
126 17989 : 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 17989 : 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 17989 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
170 : INTEGER, DIMENSION(3) :: cell
171 17989 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
172 17989 : npgfb, nsgfa, nsgfb
173 17989 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
174 17989 : 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 17989 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, qab
179 17989 : 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 35978 : REAL(KIND=dp), DIMENSION(3, natom) :: force_thread
183 17989 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
184 17989 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
185 17989 : zeta, zetb
186 17989 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dab
187 17989 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
188 17989 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: k_block
189 : TYPE(dft_control_type), POINTER :: dft_control
190 17989 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
191 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
192 : TYPE(kpoint_type), POINTER :: kpoints
193 17989 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
194 17989 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
195 : TYPE(virial_type), POINTER :: virial
196 :
197 : !$ INTEGER(kind=omp_lock_kind), &
198 17989 : !$ ALLOCATABLE, DIMENSION(:) :: locks
199 : !$ INTEGER :: lock_num, hash, hash1, hash2
200 : !$ INTEGER(KIND=int_8) :: iatom8
201 : !$ INTEGER, PARAMETER :: nlock = 501
202 :
203 : MARK_USED(int_8)
204 :
205 17989 : CALL timeset(routineN, handle)
206 :
207 : ! test for matrices (kpoints or standard gamma point)
208 17989 : IF (PRESENT(matrix_t)) THEN
209 : dokp = .FALSE.
210 : use_cell_mapping = .FALSE.
211 16673 : ELSEIF (PRESENT(matrixkp_t)) THEN
212 16673 : dokp = .TRUE.
213 16673 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
214 16673 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
215 16673 : use_cell_mapping = (SIZE(cell_to_index) > 1)
216 : ELSE
217 0 : CPABORT("")
218 : END IF
219 :
220 17989 : NULLIFY (atomic_kind_set, qs_kind_set, p_block, dft_control)
221 : CALL get_ks_env(ks_env, &
222 : dft_control=dft_control, &
223 : atomic_kind_set=atomic_kind_set, &
224 17989 : qs_kind_set=qs_kind_set)
225 :
226 17989 : nimg = dft_control%nimages
227 17989 : nkind = SIZE(atomic_kind_set)
228 :
229 17989 : do_forces = .FALSE.
230 17989 : IF (PRESENT(calculate_forces)) do_forces = calculate_forces
231 :
232 17989 : nder = 0
233 17989 : IF (PRESENT(nderivative)) nder = nderivative
234 17989 : maxder = ncoset(nder)
235 :
236 : ! check for symmetry
237 17989 : CPASSERT(SIZE(sab_nl) > 0)
238 17989 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
239 :
240 : ! prepare basis set
241 85862 : ALLOCATE (basis_set_list(nkind))
242 17989 : CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
243 :
244 17989 : IF (dokp) THEN
245 16673 : CALL dbcsr_allocate_matrix_set(matrixkp_t, 1, nimg)
246 : CALL create_sab_matrix(ks_env, matrixkp_t, matrix_name, basis_set_list, basis_set_list, &
247 16673 : sab_nl, do_symmetric)
248 : ELSE
249 1316 : CALL dbcsr_allocate_matrix_set(matrix_t, maxder)
250 : CALL create_sab_matrix(ks_env, matrix_t, matrix_name, basis_set_list, basis_set_list, &
251 1316 : sab_nl, do_symmetric)
252 : END IF
253 :
254 17989 : use_virial = .FALSE.
255 17989 : IF (do_forces) THEN
256 : ! if forces -> maybe virial too
257 7369 : CALL get_ks_env(ks_env=ks_env, force=force, virial=virial)
258 7369 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
259 : ! we need density matrix for forces
260 7369 : IF (dokp) THEN
261 6077 : CPASSERT(PRESENT(matrixkp_p))
262 : ELSE
263 1292 : CPASSERT(PRESENT(matrix_p))
264 : END IF
265 : END IF
266 :
267 284365 : force_thread = 0.0_dp
268 17989 : pv_thread = 0.0_dp
269 :
270 : ! *** Allocate work storage ***
271 17989 : ldsab = get_memory_usage(qs_kind_set, basis_type)
272 :
273 : !$OMP PARALLEL DEFAULT(NONE) &
274 : !$OMP SHARED (do_forces, ldsab, use_cell_mapping, do_symmetric, dokp,&
275 : !$OMP sab_nl, ncoset, maxder, nder, ndod, use_virial, matrix_t, matrixkp_t,&
276 : !$OMP matrix_p, basis_set_list, atom_of_kind, cell_to_index, matrixkp_p, locks, natom) &
277 : !$OMP PRIVATE (k_block, kab, qab, pab, ikind, jkind, iatom, jatom, rab, cell, basis_set_a, basis_set_b, f, &
278 : !$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
279 : !$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, tab, &
280 : !$OMP slot, zetb, scon_a, scon_b, i, ic, irow, icol, f0, ff, found, trans, rab2, sgfa, sgfb, iset, jset, &
281 : !$OMP hash, hash1, hash2, iatom8) &
282 17989 : !$OMP REDUCTION (+ : pv_thread, force_thread )
283 :
284 : !$OMP SINGLE
285 : !$ ALLOCATE (locks(nlock))
286 : !$OMP END SINGLE
287 :
288 : !$OMP DO
289 : !$ DO lock_num = 1, nlock
290 : !$ call omp_init_lock(locks(lock_num))
291 : !$ END DO
292 : !$OMP END DO
293 :
294 : ALLOCATE (kab(ldsab, ldsab, maxder), qab(ldsab, ldsab))
295 : IF (do_forces) THEN
296 : ALLOCATE (dab(ldsab, ldsab, 3), pab(ldsab, ldsab))
297 : END IF
298 :
299 : ALLOCATE (k_block(maxder))
300 : DO i = 1, maxder
301 : NULLIFY (k_block(i)%block)
302 : END DO
303 :
304 : !$OMP DO SCHEDULE(GUIDED)
305 : DO slot = 1, sab_nl(1)%nl_size
306 :
307 : ikind = sab_nl(1)%nlist_task(slot)%ikind
308 : jkind = sab_nl(1)%nlist_task(slot)%jkind
309 : iatom = sab_nl(1)%nlist_task(slot)%iatom
310 : jatom = sab_nl(1)%nlist_task(slot)%jatom
311 : cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
312 : rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
313 :
314 : basis_set_a => basis_set_list(ikind)%gto_basis_set
315 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
316 : basis_set_b => basis_set_list(jkind)%gto_basis_set
317 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
318 :
319 : !$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
320 : !$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
321 :
322 : ! basis ikind
323 : first_sgfa => basis_set_a%first_sgf
324 : la_max => basis_set_a%lmax
325 : la_min => basis_set_a%lmin
326 : npgfa => basis_set_a%npgf
327 : nseta = basis_set_a%nset
328 : nsgfa => basis_set_a%nsgf_set
329 : rpgfa => basis_set_a%pgf_radius
330 : set_radius_a => basis_set_a%set_radius
331 : scon_a => basis_set_a%scon
332 : zeta => basis_set_a%zet
333 : ! basis jkind
334 : first_sgfb => basis_set_b%first_sgf
335 : lb_max => basis_set_b%lmax
336 : lb_min => basis_set_b%lmin
337 : npgfb => basis_set_b%npgf
338 : nsetb = basis_set_b%nset
339 : nsgfb => basis_set_b%nsgf_set
340 : rpgfb => basis_set_b%pgf_radius
341 : set_radius_b => basis_set_b%set_radius
342 : scon_b => basis_set_b%scon
343 : zetb => basis_set_b%zet
344 :
345 : IF (use_cell_mapping) THEN
346 : ic = cell_to_index(cell(1), cell(2), cell(3))
347 : CPASSERT(ic > 0)
348 : ELSE
349 : ic = 1
350 : END IF
351 :
352 : IF (do_symmetric) THEN
353 : IF (iatom <= jatom) THEN
354 : irow = iatom
355 : icol = jatom
356 : ELSE
357 : irow = jatom
358 : icol = iatom
359 : END IF
360 : f0 = 2.0_dp
361 : IF (iatom == jatom) f0 = 1.0_dp
362 : ff = 2.0_dp
363 : ELSE
364 : irow = iatom
365 : icol = jatom
366 : f0 = 1.0_dp
367 : ff = 1.0_dp
368 : END IF
369 : IF (dokp) THEN
370 : CALL dbcsr_get_block_p(matrix=matrixkp_t(1, ic)%matrix, &
371 : row=irow, col=icol, BLOCK=k_block(1)%block, found=found)
372 : CPASSERT(found)
373 : ELSE
374 : DO i = 1, maxder
375 : NULLIFY (k_block(i)%block)
376 : CALL dbcsr_get_block_p(matrix=matrix_t(i)%matrix, &
377 : row=irow, col=icol, BLOCK=k_block(i)%block, found=found)
378 : CPASSERT(found)
379 : END DO
380 : END IF
381 :
382 : IF (do_forces) THEN
383 : NULLIFY (p_block)
384 : IF (dokp) THEN
385 : CALL dbcsr_get_block_p(matrix=matrixkp_p(1, ic)%matrix, &
386 : row=irow, col=icol, block=p_block, found=found)
387 : CPASSERT(found)
388 : ELSE
389 : CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
390 : block=p_block, found=found)
391 : CPASSERT(found)
392 : END IF
393 : END IF
394 :
395 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
396 : tab = SQRT(rab2)
397 : trans = do_symmetric .AND. (iatom > jatom)
398 :
399 : DO iset = 1, nseta
400 :
401 : ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
402 : sgfa = first_sgfa(1, iset)
403 :
404 : DO jset = 1, nsetb
405 :
406 : IF (set_radius_a(iset) + set_radius_b(jset) < tab) CYCLE
407 :
408 : !$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
409 : !$ hash = MOD(hash1 + hash2, nlock) + 1
410 :
411 : ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
412 : sgfb = first_sgfb(1, jset)
413 :
414 : IF (do_forces .AND. ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
415 : ! Decontract P matrix block
416 : kab = 0.0_dp
417 : CALL block_add("OUT", kab(:, :, 1), nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
418 : CALL decontraction(kab(:, :, 1), pab, scon_a(:, sgfa:), ncoa, nsgfa(iset), scon_b(:, sgfb:), ncob, nsgfb(jset), &
419 : trans=trans)
420 : ! calculate integrals and derivatives
421 : CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
422 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
423 : rab, kab(:, :, 1), dab)
424 : CALL force_trace(force_a, dab, pab, ncoa, ncob, 3)
425 : force_thread(:, iatom) = force_thread(:, iatom) + ff*force_a(:)
426 : force_thread(:, jatom) = force_thread(:, jatom) - ff*force_a(:)
427 : IF (use_virial) THEN
428 : CALL virial_pair_force(pv_thread, f0, force_a, rab)
429 : END IF
430 : ELSE
431 : ! calclulate integrals
432 : IF (nder == 0) THEN
433 : CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
434 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
435 : rab, kab=kab(:, :, 1))
436 : ELSE IF (nder == 1) THEN
437 : CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
438 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
439 : rab, kab=kab(:, :, 1), dab=kab(:, :, 2:4))
440 : END IF
441 : END IF
442 : DO i = 1, maxder
443 : f = 1.0_dp
444 : IF (ndod(i) == 1 .AND. trans) f = -1.0_dp
445 : ! Contraction step
446 : CALL contraction(kab(:, :, i), qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), &
447 : cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), fscale=f, &
448 : trans=trans)
449 : !$ CALL omp_set_lock(locks(hash))
450 : CALL block_add("IN", qab, nsgfa(iset), nsgfb(jset), k_block(i)%block, sgfa, sgfb, trans=trans)
451 : !$ CALL omp_unset_lock(locks(hash))
452 : END DO
453 : END DO
454 : END DO
455 :
456 : END DO
457 : DEALLOCATE (kab, qab)
458 : IF (do_forces) DEALLOCATE (pab, dab)
459 : !$OMP DO
460 : !$ DO lock_num = 1, nlock
461 : !$ call omp_destroy_lock(locks(lock_num))
462 : !$ END DO
463 : !$OMP END DO
464 :
465 : !$OMP SINGLE
466 : !$ DEALLOCATE (locks)
467 : !$OMP END SINGLE NOWAIT
468 :
469 : !$OMP END PARALLEL
470 :
471 17989 : IF (do_forces) THEN
472 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, &
473 7369 : kind_of=kind_of)
474 :
475 : !$OMP DO
476 : DO iatom = 1, natom
477 25915 : atom_a = atom_of_kind(iatom)
478 25915 : ikind = kind_of(iatom)
479 103660 : force(ikind)%kinetic(:, atom_a) = force(ikind)%kinetic(:, atom_a) + force_thread(:, iatom)
480 : END DO
481 : !$OMP END DO
482 : END IF
483 17989 : IF (do_forces .AND. use_virial) THEN
484 10530 : virial%pv_ekinetic = virial%pv_ekinetic + pv_thread
485 10530 : virial%pv_virial = virial%pv_virial + pv_thread
486 : END IF
487 :
488 17989 : IF (dokp) THEN
489 64630 : DO ic = 1, nimg
490 47957 : CALL dbcsr_finalize(matrixkp_t(1, ic)%matrix)
491 64630 : IF (PRESENT(eps_filter)) THEN
492 47179 : CALL dbcsr_filter(matrixkp_t(1, ic)%matrix, eps_filter)
493 : END IF
494 : END DO
495 : ELSE
496 1316 : CALL dbcsr_finalize(matrix_t(1)%matrix)
497 1316 : IF (PRESENT(eps_filter)) THEN
498 12 : CALL dbcsr_filter(matrix_t(1)%matrix, eps_filter)
499 : END IF
500 : END IF
501 :
502 17989 : DEALLOCATE (basis_set_list)
503 :
504 17989 : CALL timestop(handle)
505 :
506 53967 : END SUBROUTINE build_kinetic_matrix_low
507 :
508 : END MODULE qs_kinetic
509 :
|
