(git:5557c94)
Loading...
Searching...
No Matches
qs_overlap.F
Go to the documentation of this file.
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 overlap matrix, its derivatives and forces
10!> \par History
11!> JGH: removed printing routines
12!> JGH: upgraded to unique routine for overlaps
13!> JGH: Add specific routine for 'forces only'
14!> Major refactoring for new overlap routines
15!> JGH: Kpoints
16!> JGH: openMP refactoring
17!> \author Matthias Krack (03.09.2001,25.06.2003)
18! **************************************************************************************************
19MODULE qs_overlap
20 USE ai_contraction, ONLY: block_add,&
21 contraction,&
22 decontraction,&
23 force_trace
24 USE ai_overlap, ONLY: overlap_ab
25 USE atomic_kind_types, ONLY: atomic_kind_type,&
26 get_atomic_kind_set
27 USE basis_set_types, ONLY: get_gto_basis_set,&
28 gto_basis_set_p_type,&
29 gto_basis_set_type
30 USE block_p_types, ONLY: block_p_type
31 USE cp_control_types, ONLY: dft_control_type
32 USE cp_dbcsr_api, ONLY: &
33 dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_get_block_p, &
34 dbcsr_p_type, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
35 dbcsr_type_symmetric
36 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
37 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
38 USE kinds, ONLY: default_string_length,&
39 dp,&
40 int_8
41 USE kpoint_types, ONLY: get_kpoint_info,&
42 kpoint_type
43 USE orbital_pointers, ONLY: indco,&
44 ncoset
45 USE orbital_symbols, ONLY: cgf_symbol
46 USE particle_methods, ONLY: get_particle_set
47 USE particle_types, ONLY: particle_type
48 USE qs_force_types, ONLY: qs_force_type
49 USE qs_integral_utils, ONLY: basis_set_list_setup,&
50 get_memory_usage
51 USE qs_kind_types, ONLY: qs_kind_type
52 USE qs_ks_types, ONLY: get_ks_env,&
53 qs_ks_env_type
54 USE qs_neighbor_list_types, ONLY: get_neighbor_list_set_p,&
55 neighbor_list_set_p_type
56 USE string_utilities, ONLY: compress,&
57 uppercase
58 USE virial_methods, ONLY: virial_pair_force
59 USE virial_types, ONLY: virial_type
60
61!$ USE OMP_LIB, ONLY: omp_lock_kind, &
62!$ omp_init_lock, omp_set_lock, &
63!$ omp_unset_lock, omp_destroy_lock
64
65#include "./base/base_uses.f90"
66
67 IMPLICIT NONE
68
69 PRIVATE
70
71! *** Global parameters ***
72
73 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_overlap'
74
75 ! should be a parameter, but this triggers a bug in OMPed code with gfortran 4.9
76 INTEGER, DIMENSION(1:56), SAVE :: ndod = [0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, &
77 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, &
78 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, &
79 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
80
81 INTERFACE create_sab_matrix
82 MODULE PROCEDURE create_sab_matrix_1d, create_sab_matrix_2d
83 END INTERFACE
84
85! *** Public subroutines ***
86
87 PUBLIC :: build_overlap_matrix, build_overlap_matrix_simple, &
88 build_overlap_force, create_sab_matrix
89
90CONTAINS
91
92! **************************************************************************************************
93!> \brief Calculation of the overlap matrix over Cartesian Gaussian functions.
94!> \param ks_env the QS env
95!> \param matrix_s The overlap matrix to be calculated (optional)
96!> \param matrixkp_s The overlap matrices to be calculated (kpoints, optional)
97!> \param matrix_name The name of the overlap matrix (i.e. for output)
98!> \param nderivative Derivative with respect to basis origin
99!> \param basis_type_a basis set to be used for bra in <a|b>
100!> \param basis_type_b basis set to be used for ket in <a|b>
101!> \param sab_nl pair list (must be consistent with basis sets!)
102!> \param calculate_forces (optional) ...
103!> \param matrix_p density matrix for force calculation (optional)
104!> \param matrixkp_p density matrix for force calculation with k_points (optional)
105!> \param ext_kpoints ...
106!> \date 11.03.2002
107!> \par History
108!> Enlarged functionality of this routine. Now overlap matrices based
109!> on different basis sets can be calculated, taking into account also
110!> mixed overlaps NOTE: the pointer to the overlap matrix must now be
111!> put into its corresponding env outside of this routine
112!> [Manuel Guidon]
113!> Generalized for derivatives and force calculations [JHU]
114!> Kpoints, returns overlap matrices in real space index form
115!> \author MK
116!> \version 1.0
117! **************************************************************************************************
118 SUBROUTINE build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, &
119 nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, &
120 matrix_p, matrixkp_p, ext_kpoints)
121
122 TYPE(qs_ks_env_type), POINTER :: ks_env
123 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
124 POINTER :: matrix_s
125 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
126 POINTER :: matrixkp_s
127 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
128 INTEGER, INTENT(IN), OPTIONAL :: nderivative
129 CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b
130 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
131 POINTER :: sab_nl
132 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
133 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p
134 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
135 POINTER :: matrixkp_p
136 TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
137
138 INTEGER :: natom
139
140 mark_used(int_8)
141
142 CALL get_ks_env(ks_env, natom=natom)
143
144 CALL build_overlap_matrix_low(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, &
145 basis_type_a, basis_type_b, sab_nl, calculate_forces, &
146 matrix_p, matrixkp_p, ext_kpoints, natom)
147
148 END SUBROUTINE build_overlap_matrix
149
150! **************************************************************************************************
151!> \brief Implementation of build_overlap_matrix with the additional natom argument passed in to
152!> allow for explicitly shaped force_thread array which is needed for OMP REDUCTION.
153!> \param ks_env ...
154!> \param matrix_s ...
155!> \param matrixkp_s ...
156!> \param matrix_name ...
157!> \param nderivative ...
158!> \param basis_type_a ...
159!> \param basis_type_b ...
160!> \param sab_nl ...
161!> \param calculate_forces ...
162!> \param matrix_p ...
163!> \param matrixkp_p ...
164!> \param ext_kpoints ...
165!> \param natom ...
166! **************************************************************************************************
167 SUBROUTINE build_overlap_matrix_low(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, &
168 basis_type_a, basis_type_b, sab_nl, calculate_forces, &
169 matrix_p, matrixkp_p, ext_kpoints, natom)
170
171 TYPE(qs_ks_env_type), POINTER :: ks_env
172 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
173 POINTER :: matrix_s
174 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
175 POINTER :: matrixkp_s
176 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
177 INTEGER, INTENT(IN), OPTIONAL :: nderivative
178 CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b
179 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
180 POINTER :: sab_nl
181 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
182 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p
183 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
184 POINTER :: matrixkp_p
185 TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
186 INTEGER, INTENT(IN) :: natom
187
188 CHARACTER(len=*), PARAMETER :: routinen = 'build_overlap_matrix_low'
189
190 INTEGER :: atom_a, handle, i, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, &
191 maxder, maxs, n1, n2, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
192 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
193 INTEGER, DIMENSION(3) :: cell
194 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
195 npgfb, nsgfa, nsgfb
196 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
197 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
198 LOGICAL :: do_forces, do_symmetric, dokp, found, &
199 trans, use_cell_mapping, use_virial
200 REAL(kind=dp) :: dab, f, f0, ff, rab2
201 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: owork, pmat
202 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint
203 REAL(kind=dp), DIMENSION(3) :: force_a, rab
204 REAL(kind=dp), DIMENSION(3, 3) :: pv_thread
205 REAL(kind=dp), DIMENSION(3, natom) :: force_thread
206 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
207 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
208 zeta, zetb
209 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
210 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint
211 TYPE(dft_control_type), POINTER :: dft_control
212 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
213 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
214 TYPE(kpoint_type), POINTER :: kpoints
215 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
216 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
217 TYPE(virial_type), POINTER :: virial
218
219!$ INTEGER(kind=omp_lock_kind), &
220!$ ALLOCATABLE, DIMENSION(:) :: locks
221!$ INTEGER :: lock_num, hash, hash1, hash2
222!$ INTEGER(KIND=int_8) :: iatom8
223!$ INTEGER, PARAMETER :: nlock = 501
224
225 CALL timeset(routinen, handle)
226
227 NULLIFY (atomic_kind_set)
228 CALL get_ks_env(ks_env, &
229 atomic_kind_set=atomic_kind_set, &
230 qs_kind_set=qs_kind_set, &
231 dft_control=dft_control)
232
233 ! test for matrices (kpoints or standard gamma point)
234 IF (PRESENT(matrix_s)) THEN
235 dokp = .false.
236 use_cell_mapping = .false.
237 nimg = dft_control%nimages
238 ELSEIF (PRESENT(matrixkp_s)) THEN
239 dokp = .true.
240 IF (PRESENT(ext_kpoints)) THEN
241 IF (ASSOCIATED(ext_kpoints)) THEN
242 CALL get_kpoint_info(kpoint=ext_kpoints, cell_to_index=cell_to_index)
243 use_cell_mapping = (SIZE(cell_to_index) > 1)
244 nimg = SIZE(ext_kpoints%index_to_cell, 2)
245 ELSE
246 use_cell_mapping = .false.
247 nimg = 1
248 END IF
249 ELSE
250 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
251 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
252 use_cell_mapping = (SIZE(cell_to_index) > 1)
253 nimg = dft_control%nimages
254 END IF
255 ELSE
256 cpabort("")
257 END IF
258
259 nkind = SIZE(qs_kind_set)
260
261 IF (PRESENT(calculate_forces)) THEN
262 do_forces = calculate_forces
263 ELSE
264 do_forces = .false.
265 END IF
266
267 IF (PRESENT(nderivative)) THEN
268 nder = nderivative
269 ELSE
270 nder = 0
271 END IF
272 maxder = ncoset(nder)
273
274 ! check for symmetry
275 cpassert(SIZE(sab_nl) > 0)
276 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
277 IF (do_symmetric) THEN
278 cpassert(basis_type_a == basis_type_b)
279 END IF
280
281 ! set up basis set lists
282 ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
283 CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set)
284 CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set)
285
286 IF (dokp) THEN
287 CALL dbcsr_allocate_matrix_set(matrixkp_s, maxder, nimg)
288 CALL create_sab_matrix(ks_env, matrixkp_s, matrix_name, basis_set_list_a, basis_set_list_b, &
289 sab_nl, do_symmetric)
290 ELSE
291 CALL dbcsr_allocate_matrix_set(matrix_s, maxder)
292 CALL create_sab_matrix(ks_env, matrix_s, matrix_name, basis_set_list_a, basis_set_list_b, &
293 sab_nl, do_symmetric)
294 END IF
295 maxs = maxder
296
297 use_virial = .false.
298 IF (do_forces) THEN
299 CALL get_ks_env(ks_env=ks_env, force=force, virial=virial)
300 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
301 END IF
302
303 force_thread = 0.0_dp
304 pv_thread = 0.0_dp
305
306 ldsab = get_memory_usage(qs_kind_set, basis_type_a, basis_type_b)
307 IF (do_forces) THEN
308 ! we need density matrix for forces
309 IF (dokp) THEN
310 cpassert(PRESENT(matrixkp_p))
311 ELSE
312 cpassert(PRESENT(matrix_p))
313 END IF
314 nder = max(nder, 1)
315 END IF
316 maxder = ncoset(nder)
317
318!$OMP PARALLEL DEFAULT(NONE) &
319!$OMP SHARED (do_forces, ldsab, maxder, use_cell_mapping, do_symmetric, maxs, dokp, &
320!$OMP ncoset, nder, use_virial, ndod, sab_nl, nimg,&
321!$OMP matrix_s, matrix_p,basis_set_list_a, basis_set_list_b, cell_to_index, &
322!$OMP matrixkp_s, matrixkp_p, locks, natom) &
323!$OMP PRIVATE (oint, owork, pmat, sint, ikind, jkind, iatom, jatom, rab, cell, &
324!$OMP basis_set_a, basis_set_b, lock_num, &
325!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
326!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, f, &
327!$OMP zetb, scon_a, scon_b, ic, irow, icol, f0, ff, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
328!$OMP hash, hash1, hash2, iatom8, slot ) &
329!$OMP REDUCTION (+ : pv_thread, force_thread )
330
331!$OMP SINGLE
332!$ ALLOCATE (locks(nlock))
333!$OMP END SINGLE
334
335!$OMP DO
336!$ DO lock_num = 1, nlock
337!$ call omp_init_lock(locks(lock_num))
338!$ END DO
339!$OMP END DO
340
341 NULLIFY (p_block)
342 ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
343 IF (do_forces) ALLOCATE (pmat(ldsab, ldsab))
344 ALLOCATE (sint(maxs))
345 DO i = 1, maxs
346 NULLIFY (sint(i)%block)
347 END DO
348
349!$OMP DO SCHEDULE(GUIDED)
350 DO slot = 1, sab_nl(1)%nl_size
351
352 ikind = sab_nl(1)%nlist_task(slot)%ikind
353 jkind = sab_nl(1)%nlist_task(slot)%jkind
354 iatom = sab_nl(1)%nlist_task(slot)%iatom
355 jatom = sab_nl(1)%nlist_task(slot)%jatom
356 cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
357 rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
358
359 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
360 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
361 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
362 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
363
364!$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
365!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
366
367 ! basis ikind
368 first_sgfa => basis_set_a%first_sgf
369 la_max => basis_set_a%lmax
370 la_min => basis_set_a%lmin
371 npgfa => basis_set_a%npgf
372 nseta = basis_set_a%nset
373 nsgfa => basis_set_a%nsgf_set
374 rpgfa => basis_set_a%pgf_radius
375 set_radius_a => basis_set_a%set_radius
376 scon_a => basis_set_a%scon
377 zeta => basis_set_a%zet
378 ! basis jkind
379 first_sgfb => basis_set_b%first_sgf
380 lb_max => basis_set_b%lmax
381 lb_min => basis_set_b%lmin
382 npgfb => basis_set_b%npgf
383 nsetb = basis_set_b%nset
384 nsgfb => basis_set_b%nsgf_set
385 rpgfb => basis_set_b%pgf_radius
386 set_radius_b => basis_set_b%set_radius
387 scon_b => basis_set_b%scon
388 zetb => basis_set_b%zet
389
390 IF (use_cell_mapping) THEN
391 ic = cell_to_index(cell(1), cell(2), cell(3))
392 IF (ic < 1 .OR. ic > nimg) cycle
393 ELSE
394 ic = 1
395 END IF
396
397 IF (do_symmetric) THEN
398 IF (iatom <= jatom) THEN
399 irow = iatom
400 icol = jatom
401 ELSE
402 irow = jatom
403 icol = iatom
404 END IF
405 f0 = 2.0_dp
406 ff = 2.0_dp
407 IF (iatom == jatom) f0 = 1.0_dp
408 ELSE
409 irow = iatom
410 icol = jatom
411 f0 = 1.0_dp
412 ff = 1.0_dp
413 END IF
414 DO i = 1, maxs
415 NULLIFY (sint(i)%block)
416 IF (dokp) THEN
417 CALL dbcsr_get_block_p(matrix=matrixkp_s(i, ic)%matrix, &
418 row=irow, col=icol, block=sint(i)%block, found=found)
419 cpassert(found)
420 ELSE
421 CALL dbcsr_get_block_p(matrix=matrix_s(i)%matrix, &
422 row=irow, col=icol, block=sint(i)%block, found=found)
423 cpassert(found)
424 END IF
425 END DO
426 IF (do_forces) THEN
427 NULLIFY (p_block)
428 IF (dokp) THEN
429 CALL dbcsr_get_block_p(matrix=matrixkp_p(1, ic)%matrix, &
430 row=irow, col=icol, block=p_block, found=found)
431 cpassert(found)
432 ELSE
433 CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
434 block=p_block, found=found)
435 cpassert(found)
436 END IF
437 END IF
438 trans = do_symmetric .AND. (iatom > jatom)
439
440 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
441 dab = sqrt(rab2)
442
443 DO iset = 1, nseta
444
445 ncoa = npgfa(iset)*ncoset(la_max(iset))
446 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
447 sgfa = first_sgfa(1, iset)
448
449 DO jset = 1, nsetb
450
451 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
452
453!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
454!$ hash = MOD(hash1 + hash2, nlock) + 1
455
456 ncob = npgfb(jset)*ncoset(lb_max(jset))
457 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
458 sgfb = first_sgfb(1, jset)
459
460 ! calculate integrals and derivatives
461 SELECT CASE (nder)
462 CASE (0)
463 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
464 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
465 rab, sab=oint(:, :, 1))
466 CASE (1)
467 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
468 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
469 rab, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
470 CASE (2)
471 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
472 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
473 rab, sab=oint(:, :, 1), dab=oint(:, :, 2:4), ddab=oint(:, :, 5:10))
474 CASE DEFAULT
475 cpabort("")
476 END SELECT
477 IF (do_forces .AND. ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
478 ! Decontract P matrix block
479 owork = 0.0_dp
480 CALL block_add("OUT", owork, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
481 CALL decontraction(owork, pmat, scon_a(:, sgfa:), n1, nsgfa(iset), scon_b(:, sgfb:), n2, &
482 nsgfb(jset), trans=trans)
483 CALL force_trace(force_a, oint(:, :, 2:4), pmat, n1, n2, 3)
484 force_thread(:, iatom) = force_thread(:, iatom) - ff*force_a(:)
485 force_thread(:, jatom) = force_thread(:, jatom) + ff*force_a(:)
486 IF (use_virial) THEN
487 CALL virial_pair_force(pv_thread, -f0, force_a, rab)
488 END IF
489 END IF
490 ! Contraction
491 DO i = 1, maxs
492 f = 1.0_dp
493 IF (ndod(i) == 1 .AND. trans) f = -1.0_dp
494 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
495 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=f, trans=trans)
496!$ CALL omp_set_lock(locks(hash))
497 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(i)%block, &
498 sgfa, sgfb, trans=trans)
499!$ CALL omp_unset_lock(locks(hash))
500 END DO
501
502 END DO
503 END DO
504
505 END DO
506 IF (do_forces) DEALLOCATE (pmat)
507 DEALLOCATE (oint, owork)
508 DEALLOCATE (sint)
509!$OMP DO
510!$ DO lock_num = 1, nlock
511!$ call omp_destroy_lock(locks(lock_num))
512!$ END DO
513!$OMP END DO
514
515!$OMP SINGLE
516!$ DEALLOCATE (locks)
517!$OMP END SINGLE NOWAIT
518
519!$OMP END PARALLEL
520
521 IF (do_forces) THEN
522 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
523!$OMP DO
524 DO iatom = 1, natom
525 atom_a = atom_of_kind(iatom)
526 ikind = kind_of(iatom)
527 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) + force_thread(:, iatom)
528 END DO
529!$OMP END DO
530 END IF
531 IF (do_forces .AND. use_virial) THEN
532 virial%pv_overlap = virial%pv_overlap + pv_thread
533 virial%pv_virial = virial%pv_virial + pv_thread
534 END IF
535
536 IF (dokp) THEN
537 DO i = 1, maxs
538 DO ic = 1, nimg
539 CALL dbcsr_finalize(matrixkp_s(i, ic)%matrix)
540 CALL dbcsr_filter(matrixkp_s(i, ic)%matrix, &
541 dft_control%qs_control%eps_filter_matrix)
542 END DO
543 END DO
544 ELSE
545 DO i = 1, maxs
546 CALL dbcsr_finalize(matrix_s(i)%matrix)
547 CALL dbcsr_filter(matrix_s(i)%matrix, &
548 dft_control%qs_control%eps_filter_matrix)
549 END DO
550 END IF
551
552 ! *** Release work storage ***
553 DEALLOCATE (basis_set_list_a, basis_set_list_b)
554
555 CALL timestop(handle)
556
557 END SUBROUTINE build_overlap_matrix_low
558
559! **************************************************************************************************
560!> \brief Calculation of the overlap matrix over Cartesian Gaussian functions.
561!> \param ks_env the QS env
562!> \param matrix_s The overlap matrix to be calculated
563!> \param basis_set_list_a basis set list to be used for bra in <a|b>
564!> \param basis_set_list_b basis set list to be used for ket in <a|b>
565!> \param sab_nl pair list (must be consistent with basis sets!)
566!> \date 11.03.2016
567!> \par History
568!> Simplified version of build_overlap_matrix
569!> \author MK
570!> \version 1.0
571! **************************************************************************************************
572 SUBROUTINE build_overlap_matrix_simple(ks_env, matrix_s, &
573 basis_set_list_a, basis_set_list_b, sab_nl)
574
575 TYPE(qs_ks_env_type), POINTER :: ks_env
576 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
577 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
578 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
579 POINTER :: sab_nl
580
581 CHARACTER(len=*), PARAMETER :: routinen = 'build_overlap_matrix_simple'
582
583 INTEGER :: handle, iatom, icol, ikind, irow, iset, &
584 jatom, jkind, jset, ldsab, m1, m2, n1, &
585 n2, natom, ncoa, ncob, nkind, nseta, &
586 nsetb, sgfa, sgfb, slot
587 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
588 npgfb, nsgfa, nsgfb
589 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
590 LOGICAL :: do_symmetric, found, trans
591 REAL(kind=dp) :: dab, rab2
592 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
593 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint
594 REAL(kind=dp), DIMENSION(3) :: rab
595 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
596 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
597 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
598 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint
599 TYPE(dft_control_type), POINTER :: dft_control
600 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
601 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
602
603!$ INTEGER(kind=omp_lock_kind), &
604!$ ALLOCATABLE, DIMENSION(:) :: locks
605!$ INTEGER :: lock_num, hash, hash1, hash2
606!$ INTEGER(KIND=int_8) :: iatom8
607!$ INTEGER, PARAMETER :: nlock = 501
608
609 NULLIFY (dft_control)
610
611 CALL timeset(routinen, handle)
612
613 NULLIFY (atomic_kind_set)
614 CALL get_ks_env(ks_env, &
615 atomic_kind_set=atomic_kind_set, &
616 natom=natom, &
617 qs_kind_set=qs_kind_set, &
618 dft_control=dft_control)
619
620 ! check for symmetry
621 cpassert(SIZE(sab_nl) > 0)
622 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
623
624 nkind = SIZE(qs_kind_set)
625
626 CALL dbcsr_allocate_matrix_set(matrix_s, 1)
627 CALL create_sab_matrix(ks_env, matrix_s, "Matrix", basis_set_list_a, basis_set_list_b, &
628 sab_nl, do_symmetric)
629
630 ldsab = 0
631 DO ikind = 1, nkind
632 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
633 CALL get_gto_basis_set(gto_basis_set=basis_set_a, maxco=m1, nsgf=m2)
634 ldsab = max(m1, m2, ldsab)
635 basis_set_b => basis_set_list_b(ikind)%gto_basis_set
636 CALL get_gto_basis_set(gto_basis_set=basis_set_b, maxco=m1, nsgf=m2)
637 ldsab = max(m1, m2, ldsab)
638 END DO
639
640!$OMP PARALLEL DEFAULT(NONE) &
641!$OMP SHARED (ldsab,sab_nl,do_symmetric,ncoset,natom,&
642!$OMP matrix_s,basis_set_list_a,basis_set_list_b,locks) &
643!$OMP PRIVATE (oint,owork,sint,ikind,jkind,iatom,jatom,rab,basis_set_a,basis_set_b,&
644!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, &
645!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, dab, &
646!$OMP zetb, scon_a, scon_b, irow, icol, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
647!$OMP slot, lock_num, hash, hash1, hash2, iatom8 )
648
649!$OMP SINGLE
650!$ ALLOCATE (locks(nlock))
651!$OMP END SINGLE
652
653!$OMP DO
654!$ DO lock_num = 1, nlock
655!$ call omp_init_lock(locks(lock_num))
656!$ END DO
657!$OMP END DO
658
659 ALLOCATE (oint(ldsab, ldsab, 1), owork(ldsab, ldsab))
660 ALLOCATE (sint(1))
661 NULLIFY (sint(1)%block)
662
663!$OMP DO SCHEDULE(GUIDED)
664 DO slot = 1, sab_nl(1)%nl_size
665 ikind = sab_nl(1)%nlist_task(slot)%ikind
666 jkind = sab_nl(1)%nlist_task(slot)%jkind
667 iatom = sab_nl(1)%nlist_task(slot)%iatom
668 jatom = sab_nl(1)%nlist_task(slot)%jatom
669 rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
670!$ iatom8 = (iatom - 1)*natom + jatom
671!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
672
673 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
674 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
675 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
676 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
677 ! basis ikind
678 first_sgfa => basis_set_a%first_sgf
679 la_max => basis_set_a%lmax
680 la_min => basis_set_a%lmin
681 npgfa => basis_set_a%npgf
682 nseta = basis_set_a%nset
683 nsgfa => basis_set_a%nsgf_set
684 rpgfa => basis_set_a%pgf_radius
685 set_radius_a => basis_set_a%set_radius
686 scon_a => basis_set_a%scon
687 zeta => basis_set_a%zet
688 ! basis jkind
689 first_sgfb => basis_set_b%first_sgf
690 lb_max => basis_set_b%lmax
691 lb_min => basis_set_b%lmin
692 npgfb => basis_set_b%npgf
693 nsetb = basis_set_b%nset
694 nsgfb => basis_set_b%nsgf_set
695 rpgfb => basis_set_b%pgf_radius
696 set_radius_b => basis_set_b%set_radius
697 scon_b => basis_set_b%scon
698 zetb => basis_set_b%zet
699
700 IF (do_symmetric) THEN
701 IF (iatom <= jatom) THEN
702 irow = iatom
703 icol = jatom
704 ELSE
705 irow = jatom
706 icol = iatom
707 END IF
708 ELSE
709 irow = iatom
710 icol = jatom
711 END IF
712
713 NULLIFY (sint(1)%block)
714 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
715 row=irow, col=icol, block=sint(1)%block, found=found)
716 cpassert(found)
717 trans = do_symmetric .AND. (iatom > jatom)
718
719 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
720 dab = sqrt(rab2)
721
722 DO iset = 1, nseta
723
724 ncoa = npgfa(iset)*ncoset(la_max(iset))
725 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
726 sgfa = first_sgfa(1, iset)
727
728 DO jset = 1, nsetb
729
730 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
731
732!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
733!$ hash = MOD(hash1 + hash2, nlock) + 1
734
735 ncob = npgfb(jset)*ncoset(lb_max(jset))
736 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
737 sgfb = first_sgfb(1, jset)
738
739 ! calculate integrals and derivatives
740 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
741 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
742 rab, sab=oint(:, :, 1))
743 ! Contraction
744 CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
745 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=trans)
746!$OMP CRITICAL(blockadd)
747 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(1)%block, &
748 sgfa, sgfb, trans=trans)
749!$OMP END CRITICAL(blockadd)
750
751 END DO
752 END DO
753
754 END DO
755 DEALLOCATE (oint, owork)
756 DEALLOCATE (sint)
757!$OMP DO
758!$ DO lock_num = 1, nlock
759!$ call omp_destroy_lock(locks(lock_num))
760!$ END DO
761!$OMP END DO
762
763!$OMP SINGLE
764!$ DEALLOCATE (locks)
765!$OMP END SINGLE NOWAIT
766
767!$OMP END PARALLEL
768
769 CALL dbcsr_finalize(matrix_s(1)%matrix)
770 CALL dbcsr_filter(matrix_s(1)%matrix, dft_control%qs_control%eps_filter_matrix)
771
772 CALL timestop(handle)
773
774 END SUBROUTINE build_overlap_matrix_simple
775
776! **************************************************************************************************
777!> \brief Calculation of the force contribution from an overlap matrix
778!> over Cartesian Gaussian functions.
779!> \param ks_env the QS environment
780!> \param force holds the calcuated force Tr(P dS/dR)
781!> \param basis_type_a basis set to be used for bra in <a|b>
782!> \param basis_type_b basis set to be used for ket in <a|b>
783!> \param sab_nl pair list (must be consistent with basis sets!)
784!> \param matrix_p density matrix for force calculation
785!> \param matrixkp_p ...
786!> \date 11.03.2002
787!> \par History
788!> Enlarged functionality of this routine. Now overlap matrices based
789!> on different basis sets can be calculated, taking into account also
790!> mixed overlaps NOTE: the pointer to the overlap matrix must now be
791!> put into its corresponding env outside of this routine
792!> [Manuel Guidon]
793!> Generalized for derivatives and force calculations [JHU]
794!> This special version is only for forces [07.07.2014, JGH]
795!> \author MK/JGH
796!> \version 1.0
797! **************************************************************************************************
798 SUBROUTINE build_overlap_force(ks_env, force, basis_type_a, basis_type_b, &
799 sab_nl, matrix_p, matrixkp_p)
800
801 TYPE(qs_ks_env_type), POINTER :: ks_env
802 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
803 CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b
804 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
805 POINTER :: sab_nl
806 TYPE(dbcsr_type), OPTIONAL :: matrix_p
807 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL :: matrixkp_p
808
809 CHARACTER(len=*), PARAMETER :: routinen = 'build_overlap_force'
810
811 INTEGER :: handle, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, n1, n2, &
812 natom, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
813 INTEGER, DIMENSION(3) :: cell
814 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
815 npgfb, nsgfa, nsgfb
816 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
817 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
818 LOGICAL :: do_symmetric, dokp, found, trans, &
819 use_cell_mapping, use_virial
820 REAL(kind=dp) :: dab, f0, ff, rab2
821 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, sab
822 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: drab
823 REAL(kind=dp), DIMENSION(3) :: force_a, rab
824 REAL(kind=dp), DIMENSION(3, 3) :: virial_thread
825 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
826 REAL(kind=dp), DIMENSION(3, SIZE(force, 2)) :: force_thread
827 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
828 zeta, zetb
829 TYPE(dft_control_type), POINTER :: dft_control
830 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
831 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
832 TYPE(kpoint_type), POINTER :: kpoints
833 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
834 TYPE(virial_type), POINTER :: virial
835
836 CALL timeset(routinen, handle)
837
838 NULLIFY (qs_kind_set)
839 CALL get_ks_env(ks_env=ks_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
840 nimg = dft_control%nimages
841
842 ! test for matrices (kpoints or standard gamma point)
843 IF (PRESENT(matrix_p)) THEN
844 dokp = .false.
845 use_cell_mapping = .false.
846 ELSEIF (PRESENT(matrixkp_p)) THEN
847 dokp = .true.
848 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
849 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
850 use_cell_mapping = (SIZE(cell_to_index) > 1)
851 ELSE
852 cpabort("")
853 END IF
854
855 nkind = SIZE(qs_kind_set)
856 nder = 1
857
858 ! check for symmetry
859 cpassert(SIZE(sab_nl) > 0)
860 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
861
862 CALL get_ks_env(ks_env=ks_env, virial=virial)
863 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
864 virial_thread = 0.0_dp
865
866 ! set up basis sets
867 ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
868 CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set)
869 CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set)
870 ldsab = get_memory_usage(qs_kind_set, basis_type_a, basis_type_b)
871
872 natom = SIZE(force, 2)
873 force_thread = 0.0_dp
874
875!$OMP PARALLEL DEFAULT(NONE) &
876!$OMP SHARED (ldsab, do_symmetric, ncoset, nder, use_virial, force, virial, ndod, sab_nl, cell_to_index, &
877!$OMP matrix_p, basis_set_list_a, basis_set_list_b, dokp, use_cell_mapping, nimg, matrixkp_p) &
878!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, ic, &
879!$OMP basis_set_a, basis_set_b, sab, pab, drab, &
880!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
881!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, &
882!$OMP zetb, scon_a, scon_b, irow, icol, f0, ff, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
883!$OMP slot, cell) &
884!$OMP REDUCTION (+ : virial_thread, force_thread )
885
886 ! *** Allocate work storage ***
887 ALLOCATE (sab(ldsab, ldsab), pab(ldsab, ldsab))
888 ALLOCATE (drab(ldsab, ldsab, 3))
889
890 ! Loop over neighbor list
891!$OMP DO SCHEDULE(GUIDED)
892 DO slot = 1, sab_nl(1)%nl_size
893 ikind = sab_nl(1)%nlist_task(slot)%ikind
894 jkind = sab_nl(1)%nlist_task(slot)%jkind
895 iatom = sab_nl(1)%nlist_task(slot)%iatom
896 jatom = sab_nl(1)%nlist_task(slot)%jatom
897 cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
898 rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
899
900 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
901 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
902 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
903 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
904 ! basis ikind
905 first_sgfa => basis_set_a%first_sgf
906 la_max => basis_set_a%lmax
907 la_min => basis_set_a%lmin
908 npgfa => basis_set_a%npgf
909 nseta = basis_set_a%nset
910 nsgfa => basis_set_a%nsgf_set
911 rpgfa => basis_set_a%pgf_radius
912 set_radius_a => basis_set_a%set_radius
913 scon_a => basis_set_a%scon
914 zeta => basis_set_a%zet
915 ! basis jkind
916 first_sgfb => basis_set_b%first_sgf
917 lb_max => basis_set_b%lmax
918 lb_min => basis_set_b%lmin
919 npgfb => basis_set_b%npgf
920 nsetb = basis_set_b%nset
921 nsgfb => basis_set_b%nsgf_set
922 rpgfb => basis_set_b%pgf_radius
923 set_radius_b => basis_set_b%set_radius
924 scon_b => basis_set_b%scon
925 zetb => basis_set_b%zet
926
927 IF (use_cell_mapping) THEN
928 ic = cell_to_index(cell(1), cell(2), cell(3))
929 IF (ic < 1 .OR. ic > nimg) cycle
930 ELSE
931 ic = 1
932 END IF
933
934 IF (do_symmetric) THEN
935 IF (iatom <= jatom) THEN
936 irow = iatom
937 icol = jatom
938 ELSE
939 irow = jatom
940 icol = iatom
941 END IF
942 f0 = 2.0_dp
943 IF (iatom == jatom) f0 = 1.0_dp
944 ff = 2.0_dp
945 ELSE
946 irow = iatom
947 icol = jatom
948 f0 = 1.0_dp
949 ff = 1.0_dp
950 END IF
951 NULLIFY (p_block)
952 IF (dokp) THEN
953 CALL dbcsr_get_block_p(matrix=matrixkp_p(ic)%matrix, row=irow, col=icol, &
954 block=p_block, found=found)
955 ELSE
956 CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
957 block=p_block, found=found)
958 END IF
959
960 trans = do_symmetric .AND. (iatom > jatom)
961
962 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
963 dab = sqrt(rab2)
964
965 DO iset = 1, nseta
966
967 ncoa = npgfa(iset)*ncoset(la_max(iset))
968 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
969 sgfa = first_sgfa(1, iset)
970
971 DO jset = 1, nsetb
972
973 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
974
975 ncob = npgfb(jset)*ncoset(lb_max(jset))
976 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
977 sgfb = first_sgfb(1, jset)
978
979 IF (ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
980 ! Decontract P matrix block
981 sab = 0.0_dp
982 CALL block_add("OUT", sab, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
983 CALL decontraction(sab, pab, scon_a(:, sgfa:), n1, nsgfa(iset), scon_b(:, sgfb:), n2, &
984 nsgfb(jset), trans=trans)
985 ! calculate integrals and derivatives
986 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
987 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
988 rab, dab=drab)
989 CALL force_trace(force_a, drab, pab, n1, n2, 3)
990 force_thread(1:3, iatom) = force_thread(1:3, iatom) - ff*force_a(1:3)
991 force_thread(1:3, jatom) = force_thread(1:3, jatom) + ff*force_a(1:3)
992 IF (use_virial) THEN
993 CALL virial_pair_force(virial_thread, -f0, force_a, rab)
994 END IF
995 END IF
996
997 END DO
998 END DO
999
1000 END DO
1001 DEALLOCATE (sab, pab, drab)
1002!$OMP END PARALLEL
1003
1004!$OMP PARALLEL DEFAULT(NONE) &
1005!$OMP SHARED(force,force_thread,natom)
1006!$OMP WORKSHARE
1007 force(1:3, 1:natom) = force(1:3, 1:natom) + force_thread(1:3, 1:natom)
1008!$OMP END WORKSHARE
1009!$OMP END PARALLEL
1010 IF (use_virial) THEN
1011 virial%pv_virial = virial%pv_virial + virial_thread
1012 virial%pv_overlap = virial%pv_overlap + virial_thread
1013 END IF
1014
1015 DEALLOCATE (basis_set_list_a, basis_set_list_b)
1016
1017 CALL timestop(handle)
1018
1019 END SUBROUTINE build_overlap_force
1020
1021! **************************************************************************************************
1022!> \brief Setup the structure of a sparse matrix based on the overlap
1023!> neighbor list
1024!> \param ks_env The QS environment
1025!> \param matrix_s Matrices to be constructed
1026!> \param matrix_name Matrix base name
1027!> \param basis_set_list_a Basis set used for <a|
1028!> \param basis_set_list_b Basis set used for |b>
1029!> \param sab_nl Overlap neighbor list
1030!> \param symmetric Is symmetry used in the neighbor list?
1031! **************************************************************************************************
1032 SUBROUTINE create_sab_matrix_1d(ks_env, matrix_s, matrix_name, &
1033 basis_set_list_a, basis_set_list_b, sab_nl, symmetric)
1034
1035 TYPE(qs_ks_env_type), POINTER :: ks_env
1036 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1037 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
1038 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
1039 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1040 POINTER :: sab_nl
1041 LOGICAL, INTENT(IN) :: symmetric
1042
1043 CHARACTER(LEN=12) :: cgfsym
1044 CHARACTER(LEN=32) :: symmetry_string
1045 CHARACTER(LEN=default_string_length) :: mname, name
1046 INTEGER :: i, maxs, natom
1047 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
1048 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1049 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1050 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1051
1052 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, &
1053 qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist)
1054
1055 natom = SIZE(particle_set)
1056
1057 IF (PRESENT(matrix_name)) THEN
1058 mname = matrix_name
1059 ELSE
1060 mname = "DUMMY"
1061 END IF
1062
1063 maxs = SIZE(matrix_s)
1064
1065 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
1066
1067 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1068 basis=basis_set_list_a)
1069 CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, &
1070 basis=basis_set_list_b)
1071
1072 ! prepare for allocation
1073 IF (symmetric) THEN
1074 symmetry_string = dbcsr_type_symmetric
1075 ELSE
1076 symmetry_string = dbcsr_type_no_symmetry
1077 END IF
1078
1079 DO i = 1, maxs
1080 IF (symmetric) THEN
1081 IF (ndod(i) == 1) THEN
1082 ! odd derivatives are anti-symmetric
1083 symmetry_string = dbcsr_type_antisymmetric
1084 ELSE
1085 symmetry_string = dbcsr_type_symmetric
1086 END IF
1087 ELSE
1088 symmetry_string = dbcsr_type_no_symmetry
1089 END IF
1090 cgfsym = cgf_symbol(1, indco(1:3, i))
1091 IF (i == 1) THEN
1092 name = mname
1093 ELSE
1094 name = trim(cgfsym(4:))//" DERIVATIVE OF THE "//trim(mname)// &
1095 " W.R.T. THE NUCLEAR COORDINATES"
1096 END IF
1097 CALL compress(name)
1098 CALL uppercase(name)
1099 ALLOCATE (matrix_s(i)%matrix)
1100 CALL dbcsr_create(matrix=matrix_s(i)%matrix, &
1101 name=trim(name), &
1102 dist=dbcsr_dist, matrix_type=symmetry_string, &
1103 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
1104 CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i)%matrix, sab_nl)
1105 END DO
1106
1107 DEALLOCATE (row_blk_sizes, col_blk_sizes)
1108
1109 END SUBROUTINE create_sab_matrix_1d
1110
1111! **************************************************************************************************
1112!> \brief Setup the structure of a sparse matrix based on the overlap
1113!> neighbor list, 2d version
1114!> \param ks_env The QS environment
1115!> \param matrix_s Matrices to be constructed
1116!> \param matrix_name Matrix base name
1117!> \param basis_set_list_a Basis set used for <a|
1118!> \param basis_set_list_b Basis set used for |b>
1119!> \param sab_nl Overlap neighbor list
1120!> \param symmetric Is symmetry used in the neighbor list?
1121! **************************************************************************************************
1122 SUBROUTINE create_sab_matrix_2d(ks_env, matrix_s, matrix_name, &
1123 basis_set_list_a, basis_set_list_b, sab_nl, symmetric)
1124
1125 TYPE(qs_ks_env_type), POINTER :: ks_env
1126 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1127 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
1128 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
1129 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1130 POINTER :: sab_nl
1131 LOGICAL, INTENT(IN) :: symmetric
1132
1133 CHARACTER(LEN=12) :: cgfsym
1134 CHARACTER(LEN=32) :: symmetry_string
1135 CHARACTER(LEN=default_string_length) :: mname, name
1136 INTEGER :: i1, i2, natom
1137 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
1138 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1139 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1140 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1141
1142 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, &
1143 qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist)
1144
1145 natom = SIZE(particle_set)
1146
1147 IF (PRESENT(matrix_name)) THEN
1148 mname = matrix_name
1149 ELSE
1150 mname = "DUMMY"
1151 END IF
1152
1153 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
1154
1155 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1156 basis=basis_set_list_a)
1157 CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, &
1158 basis=basis_set_list_b)
1159
1160 ! prepare for allocation
1161 IF (symmetric) THEN
1162 symmetry_string = dbcsr_type_symmetric
1163 ELSE
1164 symmetry_string = dbcsr_type_no_symmetry
1165 END IF
1166
1167 DO i2 = 1, SIZE(matrix_s, 2)
1168 DO i1 = 1, SIZE(matrix_s, 1)
1169 IF (symmetric) THEN
1170 IF (ndod(i1) == 1) THEN
1171 ! odd derivatives are anti-symmetric
1172 symmetry_string = dbcsr_type_antisymmetric
1173 ELSE
1174 symmetry_string = dbcsr_type_symmetric
1175 END IF
1176 ELSE
1177 symmetry_string = dbcsr_type_no_symmetry
1178 END IF
1179 cgfsym = cgf_symbol(1, indco(1:3, i1))
1180 IF (i1 == 1) THEN
1181 name = mname
1182 ELSE
1183 name = trim(cgfsym(4:))//" DERIVATIVE OF THE "//trim(mname)// &
1184 " W.R.T. THE NUCLEAR COORDINATES"
1185 END IF
1186 CALL compress(name)
1187 CALL uppercase(name)
1188 ALLOCATE (matrix_s(i1, i2)%matrix)
1189 CALL dbcsr_create(matrix=matrix_s(i1, i2)%matrix, &
1190 name=trim(name), &
1191 dist=dbcsr_dist, matrix_type=symmetry_string, &
1192 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
1193 CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i1, i2)%matrix, sab_nl)
1194 END DO
1195 END DO
1196
1197 DEALLOCATE (row_blk_sizes, col_blk_sizes)
1198
1199 END SUBROUTINE create_sab_matrix_2d
1200
1201END MODULE qs_overlap
ai_contraction::block_add
Definition ai_contraction.F:43
ai_contraction::contraction
Definition ai_contraction.F:31
ai_contraction::decontraction
Definition ai_contraction.F:35
ai_contraction::force_trace
Definition ai_contraction.F:39
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
qs_integral_utils::get_memory_usage
Definition qs_integral_utils.F:34
qs_overlap::create_sab_matrix
Definition qs_overlap.F:81
ai_contraction
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Definition ai_contraction.F:18
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:681
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:631
block_p_types
collect pointers to a block of reals
Definition block_p_types.F:14
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:692
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:657
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:672
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:372
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43
orbital_symbols
orbital_symbols
Definition orbital_symbols.F:14
orbital_symbols::cgf_symbol
character(len=12) function, public cgf_symbol(n, lxyz)
Build a Cartesian orbital symbol (orbital labels for printing).
Definition orbital_symbols.F:44
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition particle_methods.F:99
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_force_types
Definition qs_force_types.F:13
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:341
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::get_neighbor_list_set_p
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
Definition qs_neighbor_list_types.F:995
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
qs_overlap::build_overlap_force
subroutine, public build_overlap_force(ks_env, force, basis_type_a, basis_type_b, sab_nl, matrix_p, matrixkp_p)
Calculation of the force contribution from an overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:800
qs_overlap::build_overlap_matrix_simple
subroutine, public build_overlap_matrix_simple(ks_env, matrix_s, basis_set_list_a, basis_set_list_b, sab_nl)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:574
qs_overlap::build_overlap_matrix
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p, ext_kpoints)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:121
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::compress
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Definition string_utilities.F:3191
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:142
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
block_p_types::block_p_type
Definition block_p_types.F:23
cp_control_types::dft_control_type
Definition cp_control_types.F:631
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:157
particle_types::particle_type
Definition particle_types.F:35
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
virial_types::virial_type
Definition virial_types.F:26