(git:4cf809f)
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, lcart=.false.)
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, lcart=.false.)
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!> \param lcart ...
567!> \date 11.03.2016
568!> \par History
569!> Simplified version of build_overlap_matrix
570!> \author MK
571!> \version 1.0
572! **************************************************************************************************
573 SUBROUTINE build_overlap_matrix_simple(ks_env, matrix_s, &
574 basis_set_list_a, basis_set_list_b, sab_nl, lcart)
575
576 TYPE(qs_ks_env_type), POINTER :: ks_env
577 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
578 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
579 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
580 POINTER :: sab_nl
581 LOGICAL, OPTIONAL :: lcart
582
583 CHARACTER(len=*), PARAMETER :: routinen = 'build_overlap_matrix_simple'
584
585 INTEGER :: handle, iatom, icol, ikind, irow, iset, &
586 jatom, jkind, jset, ldsab, m1, m2, n1, &
587 n2, natom, ncoa, ncob, nkind, nseta, &
588 nsetb, sgfa, sgfb, slot
589 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
590 npgfb, nsgfa, nsgfb
591 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
592 LOGICAL :: do_symmetric, found, ldocart, trans
593 REAL(kind=dp) :: dab, rab2
594 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
595 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint
596 REAL(kind=dp), DIMENSION(3) :: rab
597 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
598 REAL(kind=dp), DIMENSION(:, :), POINTER :: ccona, cconb, rpgfa, rpgfb, scon_a, &
599 scon_b, zeta, zetb
600 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
601 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint
602 TYPE(dft_control_type), POINTER :: dft_control
603 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
604 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
605
606!$ INTEGER(kind=omp_lock_kind), &
607!$ ALLOCATABLE, DIMENSION(:) :: locks
608!$ INTEGER :: lock_num, hash, hash1, hash2
609!$ INTEGER(KIND=int_8) :: iatom8
610!$ INTEGER, PARAMETER :: nlock = 501
611
612 ldocart = .false.
613 IF (PRESENT(lcart)) ldocart = lcart
614 NULLIFY (dft_control)
615
616 CALL timeset(routinen, handle)
617
618 NULLIFY (atomic_kind_set)
619 CALL get_ks_env(ks_env, &
620 atomic_kind_set=atomic_kind_set, &
621 natom=natom, &
622 qs_kind_set=qs_kind_set, &
623 dft_control=dft_control)
624
625 ! check for symmetry
626 cpassert(SIZE(sab_nl) > 0)
627 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
628
629 nkind = SIZE(qs_kind_set)
630
631 CALL dbcsr_allocate_matrix_set(matrix_s, 1)
632 IF (.NOT. ldocart) THEN
633 CALL create_sab_matrix(ks_env, matrix_s, "Matrix", basis_set_list_a, basis_set_list_b, &
634 sab_nl, do_symmetric, lcart=ldocart)
635 ELSE
636 CALL create_sab_matrix(ks_env, matrix_s, "Cartesian Overlap Matrix", basis_set_list_a, basis_set_list_b, &
637 sab_nl, do_symmetric, lcart=ldocart)
638 END IF
639
640 ldsab = 0
641 DO ikind = 1, nkind
642 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
643 CALL get_gto_basis_set(gto_basis_set=basis_set_a, maxco=m1, nsgf=m2)
644 ldsab = max(m1, m2, ldsab)
645 basis_set_b => basis_set_list_b(ikind)%gto_basis_set
646 CALL get_gto_basis_set(gto_basis_set=basis_set_b, maxco=m1, nsgf=m2)
647 ldsab = max(m1, m2, ldsab)
648 END DO
649
650!$OMP PARALLEL DEFAULT(NONE) &
651!$OMP SHARED (ldsab,sab_nl,do_symmetric,ncoset,natom,&
652!$OMP matrix_s,basis_set_list_a,basis_set_list_b,locks,ldocart) &
653!$OMP PRIVATE (oint,owork,sint,ikind,jkind,iatom,jatom,rab,basis_set_a,basis_set_b,&
654!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, &
655!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, dab, &
656!$OMP zetb, scon_a, scon_b, irow, icol, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
657!$OMP slot, lock_num, hash, hash1, hash2, iatom8, ccona, cconb )
658
659!$OMP SINGLE
660!$ ALLOCATE (locks(nlock))
661!$OMP END SINGLE
662
663!$OMP DO
664!$ DO lock_num = 1, nlock
665!$ call omp_init_lock(locks(lock_num))
666!$ END DO
667!$OMP END DO
668
669 ALLOCATE (oint(ldsab, ldsab, 1), owork(ldsab, ldsab))
670 ALLOCATE (sint(1))
671 NULLIFY (sint(1)%block)
672
673!$OMP DO SCHEDULE(GUIDED)
674 DO slot = 1, sab_nl(1)%nl_size
675 ikind = sab_nl(1)%nlist_task(slot)%ikind
676 jkind = sab_nl(1)%nlist_task(slot)%jkind
677 iatom = sab_nl(1)%nlist_task(slot)%iatom
678 jatom = sab_nl(1)%nlist_task(slot)%jatom
679 rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
680!$ iatom8 = (iatom - 1)*natom + jatom
681!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
682
683 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
684 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
685 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
686 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
687 ! basis ikind
688 first_sgfa => basis_set_a%first_sgf
689 la_max => basis_set_a%lmax
690 la_min => basis_set_a%lmin
691 npgfa => basis_set_a%npgf
692 nseta = basis_set_a%nset
693 nsgfa => basis_set_a%nsgf_set
694 rpgfa => basis_set_a%pgf_radius
695 set_radius_a => basis_set_a%set_radius
696 scon_a => basis_set_a%scon
697 IF (ldocart) THEN
698 first_sgfa => basis_set_a%first_cgf
699 nsgfa => basis_set_a%ncgf_set
700 ccona => basis_set_a%ccon
701 END IF
702 zeta => basis_set_a%zet
703 ! basis jkind
704 first_sgfb => basis_set_b%first_sgf
705 lb_max => basis_set_b%lmax
706 lb_min => basis_set_b%lmin
707 npgfb => basis_set_b%npgf
708 nsetb = basis_set_b%nset
709 nsgfb => basis_set_b%nsgf_set
710 rpgfb => basis_set_b%pgf_radius
711 set_radius_b => basis_set_b%set_radius
712 scon_b => basis_set_b%scon
713 IF (ldocart) THEN
714 first_sgfb => basis_set_b%first_cgf
715 nsgfb => basis_set_b%ncgf_set
716 cconb => basis_set_b%ccon
717 END IF
718 zetb => basis_set_b%zet
719
720 IF (do_symmetric) THEN
721 IF (iatom <= jatom) THEN
722 irow = iatom
723 icol = jatom
724 ELSE
725 irow = jatom
726 icol = iatom
727 END IF
728 ELSE
729 irow = iatom
730 icol = jatom
731 END IF
732
733 NULLIFY (sint(1)%block)
734 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
735 row=irow, col=icol, block=sint(1)%block, found=found)
736 cpassert(found)
737 trans = do_symmetric .AND. (iatom > jatom)
738
739 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
740 dab = sqrt(rab2)
741
742 DO iset = 1, nseta
743
744 ncoa = npgfa(iset)*ncoset(la_max(iset))
745 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
746 sgfa = first_sgfa(1, iset)
747
748 DO jset = 1, nsetb
749
750 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
751
752!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
753!$ hash = MOD(hash1 + hash2, nlock) + 1
754
755 ncob = npgfb(jset)*ncoset(lb_max(jset))
756 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
757 sgfb = first_sgfb(1, jset)
758
759 ! calculate integrals and derivatives
760 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
761 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
762 rab, sab=oint(:, :, 1))
763 ! Contraction
764 IF (.NOT. ldocart) THEN
765 CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
766 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=trans)
767 ELSE
768 CALL contraction(oint(:, :, 1), owork, ca=ccona(:, sgfa:), na=n1, ma=nsgfa(iset), &
769 cb=cconb(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=trans)
770 END IF
771!$OMP CRITICAL(blockadd)
772 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(1)%block, &
773 sgfa, sgfb, trans=trans)
774!$OMP END CRITICAL(blockadd)
775
776 END DO
777 END DO
778
779 END DO
780 DEALLOCATE (oint, owork)
781 DEALLOCATE (sint)
782!$OMP DO
783!$ DO lock_num = 1, nlock
784!$ call omp_destroy_lock(locks(lock_num))
785!$ END DO
786!$OMP END DO
787
788!$OMP SINGLE
789!$ DEALLOCATE (locks)
790!$OMP END SINGLE NOWAIT
791
792!$OMP END PARALLEL
793
794 CALL dbcsr_finalize(matrix_s(1)%matrix)
795 CALL dbcsr_filter(matrix_s(1)%matrix, dft_control%qs_control%eps_filter_matrix)
796
797 CALL timestop(handle)
798
799 END SUBROUTINE build_overlap_matrix_simple
800
801! **************************************************************************************************
802!> \brief Calculation of the force contribution from an overlap matrix
803!> over Cartesian Gaussian functions.
804!> \param ks_env the QS environment
805!> \param force holds the calcuated force Tr(P dS/dR)
806!> \param basis_type_a basis set to be used for bra in <a|b>
807!> \param basis_type_b basis set to be used for ket in <a|b>
808!> \param sab_nl pair list (must be consistent with basis sets!)
809!> \param matrix_p density matrix for force calculation
810!> \param matrixkp_p ...
811!> \date 11.03.2002
812!> \par History
813!> Enlarged functionality of this routine. Now overlap matrices based
814!> on different basis sets can be calculated, taking into account also
815!> mixed overlaps NOTE: the pointer to the overlap matrix must now be
816!> put into its corresponding env outside of this routine
817!> [Manuel Guidon]
818!> Generalized for derivatives and force calculations [JHU]
819!> This special version is only for forces [07.07.2014, JGH]
820!> \author MK/JGH
821!> \version 1.0
822! **************************************************************************************************
823 SUBROUTINE build_overlap_force(ks_env, force, basis_type_a, basis_type_b, &
824 sab_nl, matrix_p, matrixkp_p)
825
826 TYPE(qs_ks_env_type), POINTER :: ks_env
827 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
828 CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b
829 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
830 POINTER :: sab_nl
831 TYPE(dbcsr_type), OPTIONAL :: matrix_p
832 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL :: matrixkp_p
833
834 CHARACTER(len=*), PARAMETER :: routinen = 'build_overlap_force'
835
836 INTEGER :: handle, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, n1, n2, &
837 natom, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
838 INTEGER, DIMENSION(3) :: cell
839 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
840 npgfb, nsgfa, nsgfb
841 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
842 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
843 LOGICAL :: do_symmetric, dokp, found, trans, &
844 use_cell_mapping, use_virial
845 REAL(kind=dp) :: dab, f0, ff, rab2
846 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, sab
847 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: drab
848 REAL(kind=dp), DIMENSION(3) :: force_a, rab
849 REAL(kind=dp), DIMENSION(3, 3) :: virial_thread
850 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
851 REAL(kind=dp), DIMENSION(3, SIZE(force, 2)) :: force_thread
852 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
853 zeta, zetb
854 TYPE(dft_control_type), POINTER :: dft_control
855 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
856 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
857 TYPE(kpoint_type), POINTER :: kpoints
858 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
859 TYPE(virial_type), POINTER :: virial
860
861 CALL timeset(routinen, handle)
862
863 NULLIFY (qs_kind_set)
864 CALL get_ks_env(ks_env=ks_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
865 nimg = dft_control%nimages
866
867 ! test for matrices (kpoints or standard gamma point)
868 IF (PRESENT(matrix_p)) THEN
869 dokp = .false.
870 use_cell_mapping = .false.
871 ELSEIF (PRESENT(matrixkp_p)) THEN
872 dokp = .true.
873 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
874 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
875 use_cell_mapping = (SIZE(cell_to_index) > 1)
876 ELSE
877 cpabort("")
878 END IF
879
880 nkind = SIZE(qs_kind_set)
881 nder = 1
882
883 ! check for symmetry
884 cpassert(SIZE(sab_nl) > 0)
885 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
886
887 CALL get_ks_env(ks_env=ks_env, virial=virial)
888 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
889 virial_thread = 0.0_dp
890
891 ! set up basis sets
892 ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
893 CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set)
894 CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set)
895 ldsab = get_memory_usage(qs_kind_set, basis_type_a, basis_type_b)
896
897 natom = SIZE(force, 2)
898 force_thread = 0.0_dp
899
900!$OMP PARALLEL DEFAULT(NONE) &
901!$OMP SHARED (ldsab, do_symmetric, ncoset, nder, use_virial, force, virial, ndod, sab_nl, cell_to_index, &
902!$OMP matrix_p, basis_set_list_a, basis_set_list_b, dokp, use_cell_mapping, nimg, matrixkp_p) &
903!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, ic, &
904!$OMP basis_set_a, basis_set_b, sab, pab, drab, &
905!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
906!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, &
907!$OMP zetb, scon_a, scon_b, irow, icol, f0, ff, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
908!$OMP slot, cell) &
909!$OMP REDUCTION (+ : virial_thread, force_thread )
910
911 ! *** Allocate work storage ***
912 ALLOCATE (sab(ldsab, ldsab), pab(ldsab, ldsab))
913 ALLOCATE (drab(ldsab, ldsab, 3))
914
915 ! Loop over neighbor list
916!$OMP DO SCHEDULE(GUIDED)
917 DO slot = 1, sab_nl(1)%nl_size
918 ikind = sab_nl(1)%nlist_task(slot)%ikind
919 jkind = sab_nl(1)%nlist_task(slot)%jkind
920 iatom = sab_nl(1)%nlist_task(slot)%iatom
921 jatom = sab_nl(1)%nlist_task(slot)%jatom
922 cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
923 rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
924
925 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
926 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
927 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
928 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
929 ! basis ikind
930 first_sgfa => basis_set_a%first_sgf
931 la_max => basis_set_a%lmax
932 la_min => basis_set_a%lmin
933 npgfa => basis_set_a%npgf
934 nseta = basis_set_a%nset
935 nsgfa => basis_set_a%nsgf_set
936 rpgfa => basis_set_a%pgf_radius
937 set_radius_a => basis_set_a%set_radius
938 scon_a => basis_set_a%scon
939 zeta => basis_set_a%zet
940 ! basis jkind
941 first_sgfb => basis_set_b%first_sgf
942 lb_max => basis_set_b%lmax
943 lb_min => basis_set_b%lmin
944 npgfb => basis_set_b%npgf
945 nsetb = basis_set_b%nset
946 nsgfb => basis_set_b%nsgf_set
947 rpgfb => basis_set_b%pgf_radius
948 set_radius_b => basis_set_b%set_radius
949 scon_b => basis_set_b%scon
950 zetb => basis_set_b%zet
951
952 IF (use_cell_mapping) THEN
953 ic = cell_to_index(cell(1), cell(2), cell(3))
954 IF (ic < 1 .OR. ic > nimg) cycle
955 ELSE
956 ic = 1
957 END IF
958
959 IF (do_symmetric) THEN
960 IF (iatom <= jatom) THEN
961 irow = iatom
962 icol = jatom
963 ELSE
964 irow = jatom
965 icol = iatom
966 END IF
967 f0 = 2.0_dp
968 IF (iatom == jatom) f0 = 1.0_dp
969 ff = 2.0_dp
970 ELSE
971 irow = iatom
972 icol = jatom
973 f0 = 1.0_dp
974 ff = 1.0_dp
975 END IF
976 NULLIFY (p_block)
977 IF (dokp) THEN
978 CALL dbcsr_get_block_p(matrix=matrixkp_p(ic)%matrix, row=irow, col=icol, &
979 block=p_block, found=found)
980 ELSE
981 CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
982 block=p_block, found=found)
983 END IF
984
985 trans = do_symmetric .AND. (iatom > jatom)
986
987 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
988 dab = sqrt(rab2)
989
990 DO iset = 1, nseta
991
992 ncoa = npgfa(iset)*ncoset(la_max(iset))
993 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
994 sgfa = first_sgfa(1, iset)
995
996 DO jset = 1, nsetb
997
998 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
999
1000 ncob = npgfb(jset)*ncoset(lb_max(jset))
1001 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
1002 sgfb = first_sgfb(1, jset)
1003
1004 IF (ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
1005 ! Decontract P matrix block
1006 sab = 0.0_dp
1007 CALL block_add("OUT", sab, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
1008 CALL decontraction(sab, pab, scon_a(:, sgfa:), n1, nsgfa(iset), scon_b(:, sgfb:), n2, &
1009 nsgfb(jset), trans=trans)
1010 ! calculate integrals and derivatives
1011 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1012 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
1013 rab, dab=drab)
1014 CALL force_trace(force_a, drab, pab, n1, n2, 3)
1015 force_thread(1:3, iatom) = force_thread(1:3, iatom) - ff*force_a(1:3)
1016 force_thread(1:3, jatom) = force_thread(1:3, jatom) + ff*force_a(1:3)
1017 IF (use_virial) THEN
1018 CALL virial_pair_force(virial_thread, -f0, force_a, rab)
1019 END IF
1020 END IF
1021
1022 END DO
1023 END DO
1024
1025 END DO
1026 DEALLOCATE (sab, pab, drab)
1027!$OMP END PARALLEL
1028
1029!$OMP PARALLEL DEFAULT(NONE) &
1030!$OMP SHARED(force,force_thread,natom)
1031!$OMP WORKSHARE
1032 force(1:3, 1:natom) = force(1:3, 1:natom) + force_thread(1:3, 1:natom)
1033!$OMP END WORKSHARE
1034!$OMP END PARALLEL
1035 IF (use_virial) THEN
1036 virial%pv_virial = virial%pv_virial + virial_thread
1037 virial%pv_overlap = virial%pv_overlap + virial_thread
1038 END IF
1039
1040 DEALLOCATE (basis_set_list_a, basis_set_list_b)
1041
1042 CALL timestop(handle)
1043
1044 END SUBROUTINE build_overlap_force
1045
1046! **************************************************************************************************
1047!> \brief Setup the structure of a sparse matrix based on the overlap
1048!> neighbor list
1049!> \param ks_env The QS environment
1050!> \param matrix_s Matrices to be constructed
1051!> \param matrix_name Matrix base name
1052!> \param basis_set_list_a Basis set used for <a|
1053!> \param basis_set_list_b Basis set used for |b>
1054!> \param sab_nl Overlap neighbor list
1055!> \param symmetric Is symmetry used in the neighbor list?
1056!> \param lcart ...
1057! **************************************************************************************************
1058 SUBROUTINE create_sab_matrix_1d(ks_env, matrix_s, matrix_name, &
1059 basis_set_list_a, basis_set_list_b, sab_nl, symmetric, lcart)
1060
1061 TYPE(qs_ks_env_type), POINTER :: ks_env
1062 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1063 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
1064 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
1065 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1066 POINTER :: sab_nl
1067 LOGICAL, INTENT(IN) :: symmetric
1068 LOGICAL, INTENT(IN), OPTIONAL :: lcart
1069
1070 CHARACTER(LEN=12) :: cgfsym
1071 CHARACTER(LEN=32) :: symmetry_string
1072 CHARACTER(LEN=default_string_length) :: mname, name
1073 INTEGER :: i, maxs, natom
1074 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
1075 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1076 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1077 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1078
1079 LOGICAL:: my_lcart
1080
1081 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, &
1082 qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist)
1083
1084 natom = SIZE(particle_set)
1085 my_lcart = .false.
1086 IF (PRESENT(lcart)) my_lcart = lcart
1087
1088 IF (PRESENT(matrix_name)) THEN
1089 mname = matrix_name
1090 ELSE
1091 mname = "DUMMY"
1092 END IF
1093
1094 maxs = SIZE(matrix_s)
1095
1096 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
1097
1098 IF (.NOT. my_lcart) THEN
1099 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1100 basis=basis_set_list_a)
1101 CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, &
1102 basis=basis_set_list_b)
1103 ELSE
1104 CALL get_particle_set(particle_set, qs_kind_set, ncgf=row_blk_sizes, &
1105 basis=basis_set_list_a)
1106 CALL get_particle_set(particle_set, qs_kind_set, ncgf=col_blk_sizes, &
1107 basis=basis_set_list_b)
1108 END IF
1109
1110 ! prepare for allocation
1111 IF (symmetric) THEN
1112 symmetry_string = dbcsr_type_symmetric
1113 ELSE
1114 symmetry_string = dbcsr_type_no_symmetry
1115 END IF
1116
1117 DO i = 1, maxs
1118 IF (symmetric) THEN
1119 IF (ndod(i) == 1) THEN
1120 ! odd derivatives are anti-symmetric
1121 symmetry_string = dbcsr_type_antisymmetric
1122 ELSE
1123 symmetry_string = dbcsr_type_symmetric
1124 END IF
1125 ELSE
1126 symmetry_string = dbcsr_type_no_symmetry
1127 END IF
1128 cgfsym = cgf_symbol(1, indco(1:3, i))
1129 IF (i == 1) THEN
1130 name = mname
1131 ELSE
1132 name = trim(cgfsym(4:))//" DERIVATIVE OF THE "//trim(mname)// &
1133 " W.R.T. THE NUCLEAR COORDINATES"
1134 END IF
1135 CALL compress(name)
1136 CALL uppercase(name)
1137 ALLOCATE (matrix_s(i)%matrix)
1138 CALL dbcsr_create(matrix=matrix_s(i)%matrix, &
1139 name=trim(name), &
1140 dist=dbcsr_dist, matrix_type=symmetry_string, &
1141 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
1142 CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i)%matrix, sab_nl)
1143 END DO
1144
1145 DEALLOCATE (row_blk_sizes, col_blk_sizes)
1146
1147 END SUBROUTINE create_sab_matrix_1d
1148
1149! **************************************************************************************************
1150!> \brief Setup the structure of a sparse matrix based on the overlap
1151!> neighbor list, 2d version
1152!> \param ks_env The QS environment
1153!> \param matrix_s Matrices to be constructed
1154!> \param matrix_name Matrix base name
1155!> \param basis_set_list_a Basis set used for <a|
1156!> \param basis_set_list_b Basis set used for |b>
1157!> \param sab_nl Overlap neighbor list
1158!> \param symmetric Is symmetry used in the neighbor list?
1159!> \param lcart ...
1160! **************************************************************************************************
1161 SUBROUTINE create_sab_matrix_2d(ks_env, matrix_s, matrix_name, &
1162 basis_set_list_a, basis_set_list_b, sab_nl, symmetric, lcart)
1163
1164 TYPE(qs_ks_env_type), POINTER :: ks_env
1165 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1166 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
1167 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
1168 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1169 POINTER :: sab_nl
1170 LOGICAL, INTENT(IN) :: symmetric
1171 LOGICAL, INTENT(in), OPTIONAL :: lcart
1172
1173 CHARACTER(LEN=12) :: cgfsym
1174 CHARACTER(LEN=32) :: symmetry_string
1175 CHARACTER(LEN=default_string_length) :: mname, name
1176 INTEGER :: i1, i2, natom
1177 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
1178 LOGICAL :: my_lcart
1179 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1180 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1181 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1182
1183 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, &
1184 qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist)
1185
1186 natom = SIZE(particle_set)
1187 my_lcart = .false.
1188 IF (PRESENT(lcart)) my_lcart = lcart
1189
1190 IF (PRESENT(matrix_name)) THEN
1191 mname = matrix_name
1192 ELSE
1193 mname = "DUMMY"
1194 END IF
1195
1196 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
1197
1198 IF (.NOT. my_lcart) THEN
1199 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1200 basis=basis_set_list_a)
1201 CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, &
1202 basis=basis_set_list_b)
1203 ELSE
1204 CALL get_particle_set(particle_set, qs_kind_set, ncgf=row_blk_sizes, &
1205 basis=basis_set_list_a)
1206 CALL get_particle_set(particle_set, qs_kind_set, ncgf=col_blk_sizes, &
1207 basis=basis_set_list_b)
1208 END IF
1209
1210 ! prepare for allocation
1211 IF (symmetric) THEN
1212 symmetry_string = dbcsr_type_symmetric
1213 ELSE
1214 symmetry_string = dbcsr_type_no_symmetry
1215 END IF
1216
1217 DO i2 = 1, SIZE(matrix_s, 2)
1218 DO i1 = 1, SIZE(matrix_s, 1)
1219 IF (symmetric) THEN
1220 IF (ndod(i1) == 1) THEN
1221 ! odd derivatives are anti-symmetric
1222 symmetry_string = dbcsr_type_antisymmetric
1223 ELSE
1224 symmetry_string = dbcsr_type_symmetric
1225 END IF
1226 ELSE
1227 symmetry_string = dbcsr_type_no_symmetry
1228 END IF
1229 cgfsym = cgf_symbol(1, indco(1:3, i1))
1230 IF (i1 == 1) THEN
1231 name = mname
1232 ELSE
1233 name = trim(cgfsym(4:))//" DERIVATIVE OF THE "//trim(mname)// &
1234 " W.R.T. THE NUCLEAR COORDINATES"
1235 END IF
1236 CALL compress(name)
1237 CALL uppercase(name)
1238 ALLOCATE (matrix_s(i1, i2)%matrix)
1239 CALL dbcsr_create(matrix=matrix_s(i1, i2)%matrix, &
1240 name=trim(name), &
1241 dist=dbcsr_dist, matrix_type=symmetry_string, &
1242 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
1243 CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i1, i2)%matrix, sab_nl)
1244 END DO
1245 END DO
1246
1247 DEALLOCATE (row_blk_sizes, col_blk_sizes)
1248
1249 END SUBROUTINE create_sab_matrix_2d
1250
1251END 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, ccon)
...
Definition basis_set_types.F:643
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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:525
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, ncgf)
Get the components of a particle set.
Definition particle_methods.F:121
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_matrix_simple
subroutine, public build_overlap_matrix_simple(ks_env, matrix_s, basis_set_list_a, basis_set_list_b, sab_nl, lcart)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:575
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:825
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:744
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:171
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