(git:374b731)
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-2024 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
33 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
34 USE dbcsr_api, ONLY: &
35 dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_get_block_p, &
36 dbcsr_p_type, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
37 dbcsr_type_symmetric
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!> \date 11.03.2002
106!> \par History
107!> Enlarged functionality of this routine. Now overlap matrices based
108!> on different basis sets can be calculated, taking into account also
109!> mixed overlaps NOTE: the pointer to the overlap matrix must now be
110!> put into its corresponding env outside of this routine
111!> [Manuel Guidon]
112!> Generalized for derivatives and force calculations [JHU]
113!> Kpoints, returns overlap matrices in real space index form
114!> \author MK
115!> \version 1.0
116! **************************************************************************************************
117 SUBROUTINE build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, &
118 nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, &
119 matrix_p, matrixkp_p)
120
121 TYPE(qs_ks_env_type), POINTER :: ks_env
122 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
123 POINTER :: matrix_s
124 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
125 POINTER :: matrixkp_s
126 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
127 INTEGER, INTENT(IN), OPTIONAL :: nderivative
128 CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b
129 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
130 POINTER :: sab_nl
131 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
132 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p
133 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
134 POINTER :: matrixkp_p
135
136 INTEGER :: natom
137
138 mark_used(int_8)
139
140 CALL get_ks_env(ks_env, natom=natom)
141
142 CALL build_overlap_matrix_low(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, &
143 basis_type_a, basis_type_b, sab_nl, calculate_forces, &
144 matrix_p, matrixkp_p, natom)
145
146 END SUBROUTINE build_overlap_matrix
147
148! **************************************************************************************************
149!> \brief Implementation of build_overlap_matrix with the additional natom argument passed in to
150!> allow for explicitly shaped force_thread array which is needed for OMP REDUCTION.
151!> \param ks_env ...
152!> \param matrix_s ...
153!> \param matrixkp_s ...
154!> \param matrix_name ...
155!> \param nderivative ...
156!> \param basis_type_a ...
157!> \param basis_type_b ...
158!> \param sab_nl ...
159!> \param calculate_forces ...
160!> \param matrix_p ...
161!> \param matrixkp_p ...
162!> \param natom ...
163! **************************************************************************************************
164 SUBROUTINE build_overlap_matrix_low(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, &
165 basis_type_a, basis_type_b, sab_nl, calculate_forces, &
166 matrix_p, matrixkp_p, natom)
167
168 TYPE(qs_ks_env_type), POINTER :: ks_env
169 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
170 POINTER :: matrix_s
171 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
172 POINTER :: matrixkp_s
173 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
174 INTEGER, INTENT(IN), OPTIONAL :: nderivative
175 CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b
176 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
177 POINTER :: sab_nl
178 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
179 TYPE(dbcsr_type), OPTIONAL, POINTER :: matrix_p
180 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
181 POINTER :: matrixkp_p
182 INTEGER, INTENT(IN) :: natom
183
184 CHARACTER(len=*), PARAMETER :: routinen = 'build_overlap_matrix_low'
185
186 INTEGER :: atom_a, handle, i, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, &
187 maxder, maxs, n1, n2, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
188 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
189 INTEGER, DIMENSION(3) :: cell
190 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
191 npgfb, nsgfa, nsgfb
192 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
193 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
194 LOGICAL :: do_forces, do_symmetric, dokp, found, &
195 trans, use_cell_mapping, use_virial
196 REAL(kind=dp) :: dab, f, f0, ff, rab2
197 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: owork, pmat
198 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint
199 REAL(kind=dp), DIMENSION(3) :: force_a, rab
200 REAL(kind=dp), DIMENSION(3, 3) :: pv_thread
201 REAL(kind=dp), DIMENSION(3, natom) :: force_thread
202 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
203 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
204 zeta, zetb
205 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
206 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint
207 TYPE(dft_control_type), POINTER :: dft_control
208 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
209 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
210 TYPE(kpoint_type), POINTER :: kpoints
211 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
212 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
213 TYPE(virial_type), POINTER :: virial
214
215!$ INTEGER(kind=omp_lock_kind), &
216!$ ALLOCATABLE, DIMENSION(:) :: locks
217!$ INTEGER :: lock_num, hash, hash1, hash2
218!$ INTEGER(KIND=int_8) :: iatom8
219!$ INTEGER, PARAMETER :: nlock = 501
220
221 CALL timeset(routinen, handle)
222
223 ! test for matrices (kpoints or standard gamma point)
224 IF (PRESENT(matrix_s)) THEN
225 dokp = .false.
226 use_cell_mapping = .false.
227 ELSEIF (PRESENT(matrixkp_s)) THEN
228 dokp = .true.
229 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
230 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
231 use_cell_mapping = (SIZE(cell_to_index) > 1)
232 ELSE
233 cpabort("")
234 END IF
235
236 NULLIFY (atomic_kind_set)
237 CALL get_ks_env(ks_env, &
238 atomic_kind_set=atomic_kind_set, &
239 qs_kind_set=qs_kind_set, &
240 dft_control=dft_control)
241
242 nimg = dft_control%nimages
243 nkind = SIZE(qs_kind_set)
244
245 IF (PRESENT(calculate_forces)) THEN
246 do_forces = calculate_forces
247 ELSE
248 do_forces = .false.
249 END IF
250
251 IF (PRESENT(nderivative)) THEN
252 nder = nderivative
253 ELSE
254 nder = 0
255 END IF
256 maxder = ncoset(nder)
257
258 ! check for symmetry
259 cpassert(SIZE(sab_nl) > 0)
260 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
261 IF (do_symmetric) THEN
262 cpassert(basis_type_a == basis_type_b)
263 END IF
264
265 ! set up basis set lists
266 ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
267 CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set)
268 CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set)
269
270 IF (dokp) THEN
271 CALL dbcsr_allocate_matrix_set(matrixkp_s, maxder, nimg)
272 CALL create_sab_matrix(ks_env, matrixkp_s, matrix_name, basis_set_list_a, basis_set_list_b, &
273 sab_nl, do_symmetric)
274 ELSE
275 CALL dbcsr_allocate_matrix_set(matrix_s, maxder)
276 CALL create_sab_matrix(ks_env, matrix_s, matrix_name, basis_set_list_a, basis_set_list_b, &
277 sab_nl, do_symmetric)
278 END IF
279 maxs = maxder
280
281 use_virial = .false.
282 IF (do_forces) THEN
283 CALL get_ks_env(ks_env=ks_env, force=force, virial=virial)
284 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
285 END IF
286
287 force_thread = 0.0_dp
288 pv_thread = 0.0_dp
289
290 ldsab = get_memory_usage(qs_kind_set, basis_type_a, basis_type_b)
291 IF (do_forces) THEN
292 ! we need density matrix for forces
293 IF (dokp) THEN
294 cpassert(PRESENT(matrixkp_p))
295 ELSE
296 cpassert(PRESENT(matrix_p))
297 END IF
298 nder = max(nder, 1)
299 END IF
300 maxder = ncoset(nder)
301
302!$OMP PARALLEL DEFAULT(NONE) &
303!$OMP SHARED (do_forces, ldsab, maxder, use_cell_mapping, do_symmetric, maxs, dokp, &
304!$OMP ncoset, nder, use_virial, ndod, sab_nl, nimg,&
305!$OMP matrix_s, matrix_p,basis_set_list_a, basis_set_list_b, cell_to_index, &
306!$OMP matrixkp_s, matrixkp_p, locks, natom) &
307!$OMP PRIVATE (oint, owork, pmat, sint, ikind, jkind, iatom, jatom, rab, cell, &
308!$OMP basis_set_a, basis_set_b, &
309!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
310!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, f, &
311!$OMP zetb, scon_a, scon_b, ic, irow, icol, f0, ff, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
312!$OMP hash, hash1, hash2, iatom8, slot ) &
313!$OMP REDUCTION (+ : pv_thread, force_thread )
314
315!$OMP SINGLE
316!$ ALLOCATE (locks(nlock))
317!$OMP END SINGLE
318
319!$OMP DO
320!$ DO lock_num = 1, nlock
321!$ call omp_init_lock(locks(lock_num))
322!$ END DO
323!$OMP END DO
324
325 NULLIFY (p_block)
326 ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
327 IF (do_forces) ALLOCATE (pmat(ldsab, ldsab))
328 ALLOCATE (sint(maxs))
329 DO i = 1, maxs
330 NULLIFY (sint(i)%block)
331 END DO
332
333!$OMP DO SCHEDULE(GUIDED)
334 DO slot = 1, sab_nl(1)%nl_size
335
336 ikind = sab_nl(1)%nlist_task(slot)%ikind
337 jkind = sab_nl(1)%nlist_task(slot)%jkind
338 iatom = sab_nl(1)%nlist_task(slot)%iatom
339 jatom = sab_nl(1)%nlist_task(slot)%jatom
340 cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
341 rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
342
343 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
344 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
345 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
346 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
347
348!$ iatom8 = INT(iatom - 1, int_8)*INT(natom, int_8) + INT(jatom, int_8)
349!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
350
351 ! basis ikind
352 first_sgfa => basis_set_a%first_sgf
353 la_max => basis_set_a%lmax
354 la_min => basis_set_a%lmin
355 npgfa => basis_set_a%npgf
356 nseta = basis_set_a%nset
357 nsgfa => basis_set_a%nsgf_set
358 rpgfa => basis_set_a%pgf_radius
359 set_radius_a => basis_set_a%set_radius
360 scon_a => basis_set_a%scon
361 zeta => basis_set_a%zet
362 ! basis jkind
363 first_sgfb => basis_set_b%first_sgf
364 lb_max => basis_set_b%lmax
365 lb_min => basis_set_b%lmin
366 npgfb => basis_set_b%npgf
367 nsetb = basis_set_b%nset
368 nsgfb => basis_set_b%nsgf_set
369 rpgfb => basis_set_b%pgf_radius
370 set_radius_b => basis_set_b%set_radius
371 scon_b => basis_set_b%scon
372 zetb => basis_set_b%zet
373
374 IF (use_cell_mapping) THEN
375 ic = cell_to_index(cell(1), cell(2), cell(3))
376 IF (ic < 1 .OR. ic > nimg) cycle
377 ELSE
378 ic = 1
379 END IF
380
381 IF (do_symmetric) THEN
382 IF (iatom <= jatom) THEN
383 irow = iatom
384 icol = jatom
385 ELSE
386 irow = jatom
387 icol = iatom
388 END IF
389 f0 = 2.0_dp
390 ff = 2.0_dp
391 IF (iatom == jatom) f0 = 1.0_dp
392 ELSE
393 irow = iatom
394 icol = jatom
395 f0 = 1.0_dp
396 ff = 1.0_dp
397 END IF
398 DO i = 1, maxs
399 NULLIFY (sint(i)%block)
400 IF (dokp) THEN
401 CALL dbcsr_get_block_p(matrix=matrixkp_s(i, ic)%matrix, &
402 row=irow, col=icol, block=sint(i)%block, found=found)
403 cpassert(found)
404 ELSE
405 CALL dbcsr_get_block_p(matrix=matrix_s(i)%matrix, &
406 row=irow, col=icol, block=sint(i)%block, found=found)
407 cpassert(found)
408 END IF
409 END DO
410 IF (do_forces) THEN
411 NULLIFY (p_block)
412 IF (dokp) THEN
413 CALL dbcsr_get_block_p(matrix=matrixkp_p(1, ic)%matrix, &
414 row=irow, col=icol, block=p_block, found=found)
415 cpassert(found)
416 ELSE
417 CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
418 block=p_block, found=found)
419 cpassert(found)
420 END IF
421 END IF
422 trans = do_symmetric .AND. (iatom > jatom)
423
424 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
425 dab = sqrt(rab2)
426
427 DO iset = 1, nseta
428
429 ncoa = npgfa(iset)*ncoset(la_max(iset))
430 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
431 sgfa = first_sgfa(1, iset)
432
433 DO jset = 1, nsetb
434
435 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
436
437!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
438!$ hash = MOD(hash1 + hash2, nlock) + 1
439
440 ncob = npgfb(jset)*ncoset(lb_max(jset))
441 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
442 sgfb = first_sgfb(1, jset)
443
444 ! calculate integrals and derivatives
445 SELECT CASE (nder)
446 CASE (0)
447 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
448 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
449 rab, sab=oint(:, :, 1))
450 CASE (1)
451 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
452 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
453 rab, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
454 CASE (2)
455 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
456 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
457 rab, sab=oint(:, :, 1), dab=oint(:, :, 2:4), ddab=oint(:, :, 5:10))
458 CASE DEFAULT
459 cpabort("")
460 END SELECT
461 IF (do_forces .AND. ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
462 ! Decontract P matrix block
463 owork = 0.0_dp
464 CALL block_add("OUT", owork, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
465 CALL decontraction(owork, pmat, scon_a(:, sgfa:), n1, nsgfa(iset), scon_b(:, sgfb:), n2, &
466 nsgfb(jset), trans=trans)
467 CALL force_trace(force_a, oint(:, :, 2:4), pmat, n1, n2, 3)
468 force_thread(:, iatom) = force_thread(:, iatom) - ff*force_a(:)
469 force_thread(:, jatom) = force_thread(:, jatom) + ff*force_a(:)
470 IF (use_virial) THEN
471 CALL virial_pair_force(pv_thread, -f0, force_a, rab)
472 END IF
473 END IF
474 ! Contraction
475 DO i = 1, maxs
476 f = 1.0_dp
477 IF (ndod(i) == 1 .AND. trans) f = -1.0_dp
478 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
479 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=f, trans=trans)
480!$ CALL omp_set_lock(locks(hash))
481 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(i)%block, &
482 sgfa, sgfb, trans=trans)
483!$ CALL omp_unset_lock(locks(hash))
484 END DO
485
486 END DO
487 END DO
488
489 END DO
490 IF (do_forces) DEALLOCATE (pmat)
491 DEALLOCATE (oint, owork)
492 DEALLOCATE (sint)
493!$OMP DO
494!$ DO lock_num = 1, nlock
495!$ call omp_destroy_lock(locks(lock_num))
496!$ END DO
497!$OMP END DO
498
499!$OMP SINGLE
500!$ DEALLOCATE (locks)
501!$OMP END SINGLE NOWAIT
502
503!$OMP END PARALLEL
504
505 IF (do_forces) THEN
506 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
507!$OMP DO
508 DO iatom = 1, natom
509 atom_a = atom_of_kind(iatom)
510 ikind = kind_of(iatom)
511 force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) + force_thread(:, iatom)
512 END DO
513!$OMP END DO
514 END IF
515 IF (do_forces .AND. use_virial) THEN
516 virial%pv_overlap = virial%pv_overlap + pv_thread
517 virial%pv_virial = virial%pv_virial + pv_thread
518 END IF
519
520 IF (dokp) THEN
521 DO i = 1, maxs
522 DO ic = 1, nimg
523 CALL dbcsr_finalize(matrixkp_s(i, ic)%matrix)
524 CALL dbcsr_filter(matrixkp_s(i, ic)%matrix, &
525 dft_control%qs_control%eps_filter_matrix)
526 END DO
527 END DO
528 ELSE
529 DO i = 1, maxs
530 CALL dbcsr_finalize(matrix_s(i)%matrix)
531 CALL dbcsr_filter(matrix_s(i)%matrix, &
532 dft_control%qs_control%eps_filter_matrix)
533 END DO
534 END IF
535
536 ! *** Release work storage ***
537 DEALLOCATE (basis_set_list_a, basis_set_list_b)
538
539 CALL timestop(handle)
540
541 END SUBROUTINE build_overlap_matrix_low
542
543! **************************************************************************************************
544!> \brief Calculation of the overlap matrix over Cartesian Gaussian functions.
545!> \param ks_env the QS env
546!> \param matrix_s The overlap matrix to be calculated
547!> \param basis_set_list_a basis set list to be used for bra in <a|b>
548!> \param basis_set_list_b basis set list to be used for ket in <a|b>
549!> \param sab_nl pair list (must be consistent with basis sets!)
550!> \date 11.03.2016
551!> \par History
552!> Simplified version of build_overlap_matrix
553!> \author MK
554!> \version 1.0
555! **************************************************************************************************
556 SUBROUTINE build_overlap_matrix_simple(ks_env, matrix_s, &
557 basis_set_list_a, basis_set_list_b, sab_nl)
558
559 TYPE(qs_ks_env_type), POINTER :: ks_env
560 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
561 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
562 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
563 POINTER :: sab_nl
564
565 CHARACTER(len=*), PARAMETER :: routinen = 'build_overlap_matrix_simple'
566
567 INTEGER :: handle, iatom, icol, ikind, irow, iset, &
568 jatom, jkind, jset, ldsab, m1, m2, n1, &
569 n2, natom, ncoa, ncob, nkind, nseta, &
570 nsetb, sgfa, sgfb, slot
571 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
572 npgfb, nsgfa, nsgfb
573 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
574 LOGICAL :: do_symmetric, found, trans
575 REAL(kind=dp) :: dab, rab2
576 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
577 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint
578 REAL(kind=dp), DIMENSION(3) :: rab
579 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
580 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
581 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
582 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint
583 TYPE(dft_control_type), POINTER :: dft_control
584 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
585 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
586
587!$ INTEGER(kind=omp_lock_kind), &
588!$ ALLOCATABLE, DIMENSION(:) :: locks
589!$ INTEGER :: lock_num, hash, hash1, hash2
590!$ INTEGER(KIND=int_8) :: iatom8
591!$ INTEGER, PARAMETER :: nlock = 501
592
593 NULLIFY (dft_control)
594
595 CALL timeset(routinen, handle)
596
597 NULLIFY (atomic_kind_set)
598 CALL get_ks_env(ks_env, &
599 atomic_kind_set=atomic_kind_set, &
600 natom=natom, &
601 qs_kind_set=qs_kind_set, &
602 dft_control=dft_control)
603
604 ! check for symmetry
605 cpassert(SIZE(sab_nl) > 0)
606 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
607
608 nkind = SIZE(qs_kind_set)
609
610 CALL dbcsr_allocate_matrix_set(matrix_s, 1)
611 CALL create_sab_matrix(ks_env, matrix_s, "Matrix", basis_set_list_a, basis_set_list_b, &
612 sab_nl, do_symmetric)
613
614 ldsab = 0
615 DO ikind = 1, nkind
616 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
617 CALL get_gto_basis_set(gto_basis_set=basis_set_a, maxco=m1, nsgf=m2)
618 ldsab = max(m1, m2, ldsab)
619 basis_set_b => basis_set_list_b(ikind)%gto_basis_set
620 CALL get_gto_basis_set(gto_basis_set=basis_set_b, maxco=m1, nsgf=m2)
621 ldsab = max(m1, m2, ldsab)
622 END DO
623
624!$OMP PARALLEL DEFAULT(NONE) &
625!$OMP SHARED (ldsab,sab_nl,do_symmetric,ncoset,natom,&
626!$OMP matrix_s,basis_set_list_a,basis_set_list_b,locks) &
627!$OMP PRIVATE (oint,owork,sint,ikind,jkind,iatom,jatom,rab,basis_set_a,basis_set_b,&
628!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, &
629!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, dab, &
630!$OMP zetb, scon_a, scon_b, irow, icol, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
631!$OMP slot, lock_num, hash, hash1, hash2, iatom8 )
632
633!$OMP SINGLE
634!$ ALLOCATE (locks(nlock))
635!$OMP END SINGLE
636
637!$OMP DO
638!$ DO lock_num = 1, nlock
639!$ call omp_init_lock(locks(lock_num))
640!$ END DO
641!$OMP END DO
642
643 ALLOCATE (oint(ldsab, ldsab, 1), owork(ldsab, ldsab))
644 ALLOCATE (sint(1))
645 NULLIFY (sint(1)%block)
646
647!$OMP DO SCHEDULE(GUIDED)
648 DO slot = 1, sab_nl(1)%nl_size
649 ikind = sab_nl(1)%nlist_task(slot)%ikind
650 jkind = sab_nl(1)%nlist_task(slot)%jkind
651 iatom = sab_nl(1)%nlist_task(slot)%iatom
652 jatom = sab_nl(1)%nlist_task(slot)%jatom
653 rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
654!$ iatom8 = (iatom - 1)*natom + jatom
655!$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
656
657 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
658 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
659 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
660 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
661 ! basis ikind
662 first_sgfa => basis_set_a%first_sgf
663 la_max => basis_set_a%lmax
664 la_min => basis_set_a%lmin
665 npgfa => basis_set_a%npgf
666 nseta = basis_set_a%nset
667 nsgfa => basis_set_a%nsgf_set
668 rpgfa => basis_set_a%pgf_radius
669 set_radius_a => basis_set_a%set_radius
670 scon_a => basis_set_a%scon
671 zeta => basis_set_a%zet
672 ! basis jkind
673 first_sgfb => basis_set_b%first_sgf
674 lb_max => basis_set_b%lmax
675 lb_min => basis_set_b%lmin
676 npgfb => basis_set_b%npgf
677 nsetb = basis_set_b%nset
678 nsgfb => basis_set_b%nsgf_set
679 rpgfb => basis_set_b%pgf_radius
680 set_radius_b => basis_set_b%set_radius
681 scon_b => basis_set_b%scon
682 zetb => basis_set_b%zet
683
684 IF (do_symmetric) THEN
685 IF (iatom <= jatom) THEN
686 irow = iatom
687 icol = jatom
688 ELSE
689 irow = jatom
690 icol = iatom
691 END IF
692 ELSE
693 irow = iatom
694 icol = jatom
695 END IF
696
697 NULLIFY (sint(1)%block)
698 CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
699 row=irow, col=icol, block=sint(1)%block, found=found)
700 cpassert(found)
701 trans = do_symmetric .AND. (iatom > jatom)
702
703 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
704 dab = sqrt(rab2)
705
706 DO iset = 1, nseta
707
708 ncoa = npgfa(iset)*ncoset(la_max(iset))
709 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
710 sgfa = first_sgfa(1, iset)
711
712 DO jset = 1, nsetb
713
714 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
715
716!$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
717!$ hash = MOD(hash1 + hash2, nlock) + 1
718
719 ncob = npgfb(jset)*ncoset(lb_max(jset))
720 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
721 sgfb = first_sgfb(1, jset)
722
723 ! calculate integrals and derivatives
724 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
725 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
726 rab, sab=oint(:, :, 1))
727 ! Contraction
728 CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
729 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=trans)
730!$OMP CRITICAL(blockadd)
731 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(1)%block, &
732 sgfa, sgfb, trans=trans)
733!$OMP END CRITICAL(blockadd)
734
735 END DO
736 END DO
737
738 END DO
739 DEALLOCATE (oint, owork)
740 DEALLOCATE (sint)
741!$OMP DO
742!$ DO lock_num = 1, nlock
743!$ call omp_destroy_lock(locks(lock_num))
744!$ END DO
745!$OMP END DO
746
747!$OMP SINGLE
748!$ DEALLOCATE (locks)
749!$OMP END SINGLE NOWAIT
750
751!$OMP END PARALLEL
752
753 CALL dbcsr_finalize(matrix_s(1)%matrix)
754 CALL dbcsr_filter(matrix_s(1)%matrix, dft_control%qs_control%eps_filter_matrix)
755
756 CALL timestop(handle)
757
758 END SUBROUTINE build_overlap_matrix_simple
759
760! **************************************************************************************************
761!> \brief Calculation of the force contribution from an overlap matrix
762!> over Cartesian Gaussian functions.
763!> \param ks_env the QS environment
764!> \param force holds the calcuated force Tr(P dS/dR)
765!> \param basis_type_a basis set to be used for bra in <a|b>
766!> \param basis_type_b basis set to be used for ket in <a|b>
767!> \param sab_nl pair list (must be consistent with basis sets!)
768!> \param matrix_p density matrix for force calculation
769!> \param matrixkp_p ...
770!> \date 11.03.2002
771!> \par History
772!> Enlarged functionality of this routine. Now overlap matrices based
773!> on different basis sets can be calculated, taking into account also
774!> mixed overlaps NOTE: the pointer to the overlap matrix must now be
775!> put into its corresponding env outside of this routine
776!> [Manuel Guidon]
777!> Generalized for derivatives and force calculations [JHU]
778!> This special version is only for forces [07.07.2014, JGH]
779!> \author MK/JGH
780!> \version 1.0
781! **************************************************************************************************
782 SUBROUTINE build_overlap_force(ks_env, force, basis_type_a, basis_type_b, &
783 sab_nl, matrix_p, matrixkp_p)
784
785 TYPE(qs_ks_env_type), POINTER :: ks_env
786 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force
787 CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b
788 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
789 POINTER :: sab_nl
790 TYPE(dbcsr_type), OPTIONAL :: matrix_p
791 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL :: matrixkp_p
792
793 CHARACTER(len=*), PARAMETER :: routinen = 'build_overlap_force'
794
795 INTEGER :: handle, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, n1, n2, &
796 natom, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
797 INTEGER, DIMENSION(3) :: cell
798 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
799 npgfb, nsgfa, nsgfb
800 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
801 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
802 LOGICAL :: do_symmetric, dokp, found, trans, &
803 use_cell_mapping, use_virial
804 REAL(kind=dp) :: dab, f0, ff, rab2
805 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, sab
806 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: drab
807 REAL(kind=dp), DIMENSION(3) :: force_a, rab
808 REAL(kind=dp), DIMENSION(3, 3) :: virial_thread
809 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
810 REAL(kind=dp), DIMENSION(3, SIZE(force, 2)) :: force_thread
811 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
812 zeta, zetb
813 TYPE(dft_control_type), POINTER :: dft_control
814 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
815 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
816 TYPE(kpoint_type), POINTER :: kpoints
817 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
818 TYPE(virial_type), POINTER :: virial
819
820 CALL timeset(routinen, handle)
821
822 NULLIFY (qs_kind_set)
823 CALL get_ks_env(ks_env=ks_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
824 nimg = dft_control%nimages
825
826 ! test for matrices (kpoints or standard gamma point)
827 IF (PRESENT(matrix_p)) THEN
828 dokp = .false.
829 use_cell_mapping = .false.
830 ELSEIF (PRESENT(matrixkp_p)) THEN
831 dokp = .true.
832 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
833 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
834 use_cell_mapping = (SIZE(cell_to_index) > 1)
835 ELSE
836 cpabort("")
837 END IF
838
839 nkind = SIZE(qs_kind_set)
840 nder = 1
841
842 ! check for symmetry
843 cpassert(SIZE(sab_nl) > 0)
844 CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
845
846 CALL get_ks_env(ks_env=ks_env, virial=virial)
847 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
848 virial_thread = 0.0_dp
849
850 ! set up basis sets
851 ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
852 CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set)
853 CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set)
854 ldsab = get_memory_usage(qs_kind_set, basis_type_a, basis_type_b)
855
856 natom = SIZE(force, 2)
857 force_thread = 0.0_dp
858
859!$OMP PARALLEL DEFAULT(NONE) &
860!$OMP SHARED (ldsab, do_symmetric, ncoset, nder, use_virial, force, virial, ndod, sab_nl, cell_to_index, &
861!$OMP matrix_p, basis_set_list_a, basis_set_list_b, dokp, use_cell_mapping, nimg, matrixkp_p) &
862!$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, ic, &
863!$OMP basis_set_a, basis_set_b, sab, pab, drab, &
864!$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
865!$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, &
866!$OMP zetb, scon_a, scon_b, irow, icol, f0, ff, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
867!$OMP slot, cell) &
868!$OMP REDUCTION (+ : virial_thread, force_thread )
869
870 ! *** Allocate work storage ***
871 ALLOCATE (sab(ldsab, ldsab), pab(ldsab, ldsab))
872 ALLOCATE (drab(ldsab, ldsab, 3))
873
874 ! Loop over neighbor list
875!$OMP DO SCHEDULE(GUIDED)
876 DO slot = 1, sab_nl(1)%nl_size
877 ikind = sab_nl(1)%nlist_task(slot)%ikind
878 jkind = sab_nl(1)%nlist_task(slot)%jkind
879 iatom = sab_nl(1)%nlist_task(slot)%iatom
880 jatom = sab_nl(1)%nlist_task(slot)%jatom
881 cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
882 rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
883
884 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
885 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
886 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
887 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
888 ! basis ikind
889 first_sgfa => basis_set_a%first_sgf
890 la_max => basis_set_a%lmax
891 la_min => basis_set_a%lmin
892 npgfa => basis_set_a%npgf
893 nseta = basis_set_a%nset
894 nsgfa => basis_set_a%nsgf_set
895 rpgfa => basis_set_a%pgf_radius
896 set_radius_a => basis_set_a%set_radius
897 scon_a => basis_set_a%scon
898 zeta => basis_set_a%zet
899 ! basis jkind
900 first_sgfb => basis_set_b%first_sgf
901 lb_max => basis_set_b%lmax
902 lb_min => basis_set_b%lmin
903 npgfb => basis_set_b%npgf
904 nsetb = basis_set_b%nset
905 nsgfb => basis_set_b%nsgf_set
906 rpgfb => basis_set_b%pgf_radius
907 set_radius_b => basis_set_b%set_radius
908 scon_b => basis_set_b%scon
909 zetb => basis_set_b%zet
910
911 IF (use_cell_mapping) THEN
912 ic = cell_to_index(cell(1), cell(2), cell(3))
913 IF (ic < 1 .OR. ic > nimg) cycle
914 ELSE
915 ic = 1
916 END IF
917
918 IF (do_symmetric) THEN
919 IF (iatom <= jatom) THEN
920 irow = iatom
921 icol = jatom
922 ELSE
923 irow = jatom
924 icol = iatom
925 END IF
926 f0 = 2.0_dp
927 IF (iatom == jatom) f0 = 1.0_dp
928 ff = 2.0_dp
929 ELSE
930 irow = iatom
931 icol = jatom
932 f0 = 1.0_dp
933 ff = 1.0_dp
934 END IF
935 NULLIFY (p_block)
936 IF (dokp) THEN
937 CALL dbcsr_get_block_p(matrix=matrixkp_p(ic)%matrix, row=irow, col=icol, &
938 block=p_block, found=found)
939 ELSE
940 CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
941 block=p_block, found=found)
942 END IF
943
944 trans = do_symmetric .AND. (iatom > jatom)
945
946 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
947 dab = sqrt(rab2)
948
949 DO iset = 1, nseta
950
951 ncoa = npgfa(iset)*ncoset(la_max(iset))
952 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
953 sgfa = first_sgfa(1, iset)
954
955 DO jset = 1, nsetb
956
957 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
958
959 ncob = npgfb(jset)*ncoset(lb_max(jset))
960 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
961 sgfb = first_sgfb(1, jset)
962
963 IF (ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
964 ! Decontract P matrix block
965 sab = 0.0_dp
966 CALL block_add("OUT", sab, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
967 CALL decontraction(sab, pab, scon_a(:, sgfa:), n1, nsgfa(iset), scon_b(:, sgfb:), n2, &
968 nsgfb(jset), trans=trans)
969 ! calculate integrals and derivatives
970 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
971 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
972 rab, dab=drab)
973 CALL force_trace(force_a, drab, pab, n1, n2, 3)
974 force_thread(1:3, iatom) = force_thread(1:3, iatom) - ff*force_a(1:3)
975 force_thread(1:3, jatom) = force_thread(1:3, jatom) + ff*force_a(1:3)
976 IF (use_virial) THEN
977 CALL virial_pair_force(virial_thread, -f0, force_a, rab)
978 END IF
979 END IF
980
981 END DO
982 END DO
983
984 END DO
985 DEALLOCATE (sab, pab, drab)
986!$OMP END PARALLEL
987
988!$OMP PARALLEL DEFAULT(NONE) &
989!$OMP SHARED(force,force_thread,natom)
990!$OMP WORKSHARE
991 force(1:3, 1:natom) = force(1:3, 1:natom) + force_thread(1:3, 1:natom)
992!$OMP END WORKSHARE
993!$OMP END PARALLEL
994 IF (use_virial) THEN
995 virial%pv_virial = virial%pv_virial + virial_thread
996 virial%pv_overlap = virial%pv_overlap + virial_thread
997 END IF
998
999 DEALLOCATE (basis_set_list_a, basis_set_list_b)
1000
1001 CALL timestop(handle)
1002
1003 END SUBROUTINE build_overlap_force
1004
1005! **************************************************************************************************
1006!> \brief Setup the structure of a sparse matrix based on the overlap
1007!> neighbor list
1008!> \param ks_env The QS environment
1009!> \param matrix_s Matrices to be constructed
1010!> \param matrix_name Matrix base name
1011!> \param basis_set_list_a Basis set used for <a|
1012!> \param basis_set_list_b Basis set used for |b>
1013!> \param sab_nl Overlap neighbor list
1014!> \param symmetric Is symmetry used in the neighbor list?
1015! **************************************************************************************************
1016 SUBROUTINE create_sab_matrix_1d(ks_env, matrix_s, matrix_name, &
1017 basis_set_list_a, basis_set_list_b, sab_nl, symmetric)
1018
1019 TYPE(qs_ks_env_type), POINTER :: ks_env
1020 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1021 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
1022 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
1023 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1024 POINTER :: sab_nl
1025 LOGICAL, INTENT(IN) :: symmetric
1026
1027 CHARACTER(LEN=12) :: cgfsym
1028 CHARACTER(LEN=32) :: symmetry_string
1029 CHARACTER(LEN=default_string_length) :: mname, name
1030 INTEGER :: i, maxs, natom
1031 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
1032 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1033 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1034 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1035
1036 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, &
1037 qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist)
1038
1039 natom = SIZE(particle_set)
1040
1041 IF (PRESENT(matrix_name)) THEN
1042 mname = matrix_name
1043 ELSE
1044 mname = "DUMMY"
1045 END IF
1046
1047 maxs = SIZE(matrix_s)
1048
1049 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
1050
1051 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1052 basis=basis_set_list_a)
1053 CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, &
1054 basis=basis_set_list_b)
1055
1056 ! prepare for allocation
1057 IF (symmetric) THEN
1058 symmetry_string = dbcsr_type_symmetric
1059 ELSE
1060 symmetry_string = dbcsr_type_no_symmetry
1061 END IF
1062
1063 DO i = 1, maxs
1064 IF (symmetric) THEN
1065 IF (ndod(i) == 1) THEN
1066 ! odd derivatives are anti-symmetric
1067 symmetry_string = dbcsr_type_antisymmetric
1068 ELSE
1069 symmetry_string = dbcsr_type_symmetric
1070 END IF
1071 ELSE
1072 symmetry_string = dbcsr_type_no_symmetry
1073 END IF
1074 cgfsym = cgf_symbol(1, indco(1:3, i))
1075 IF (i == 1) THEN
1076 name = mname
1077 ELSE
1078 name = trim(cgfsym(4:))//" DERIVATIVE OF THE "//trim(mname)// &
1079 " W.R.T. THE NUCLEAR COORDINATES"
1080 END IF
1081 CALL compress(name)
1082 CALL uppercase(name)
1083 ALLOCATE (matrix_s(i)%matrix)
1084 CALL dbcsr_create(matrix=matrix_s(i)%matrix, &
1085 name=trim(name), &
1086 dist=dbcsr_dist, matrix_type=symmetry_string, &
1087 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, &
1088 nze=0)
1089 CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i)%matrix, sab_nl)
1090 END DO
1091
1092 DEALLOCATE (row_blk_sizes, col_blk_sizes)
1093
1094 END SUBROUTINE create_sab_matrix_1d
1095
1096! **************************************************************************************************
1097!> \brief Setup the structure of a sparse matrix based on the overlap
1098!> neighbor list, 2d version
1099!> \param ks_env The QS environment
1100!> \param matrix_s Matrices to be constructed
1101!> \param matrix_name Matrix base name
1102!> \param basis_set_list_a Basis set used for <a|
1103!> \param basis_set_list_b Basis set used for |b>
1104!> \param sab_nl Overlap neighbor list
1105!> \param symmetric Is symmetry used in the neighbor list?
1106! **************************************************************************************************
1107 SUBROUTINE create_sab_matrix_2d(ks_env, matrix_s, matrix_name, &
1108 basis_set_list_a, basis_set_list_b, sab_nl, symmetric)
1109
1110 TYPE(qs_ks_env_type), POINTER :: ks_env
1111 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1112 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
1113 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
1114 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1115 POINTER :: sab_nl
1116 LOGICAL, INTENT(IN) :: symmetric
1117
1118 CHARACTER(LEN=12) :: cgfsym
1119 CHARACTER(LEN=32) :: symmetry_string
1120 CHARACTER(LEN=default_string_length) :: mname, name
1121 INTEGER :: i1, i2, natom
1122 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
1123 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1124 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1125 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1126
1127 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, &
1128 qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist)
1129
1130 natom = SIZE(particle_set)
1131
1132 IF (PRESENT(matrix_name)) THEN
1133 mname = matrix_name
1134 ELSE
1135 mname = "DUMMY"
1136 END IF
1137
1138 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
1139
1140 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1141 basis=basis_set_list_a)
1142 CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, &
1143 basis=basis_set_list_b)
1144
1145 ! prepare for allocation
1146 IF (symmetric) THEN
1147 symmetry_string = dbcsr_type_symmetric
1148 ELSE
1149 symmetry_string = dbcsr_type_no_symmetry
1150 END IF
1151
1152 DO i2 = 1, SIZE(matrix_s, 2)
1153 DO i1 = 1, SIZE(matrix_s, 1)
1154 IF (symmetric) THEN
1155 IF (ndod(i1) == 1) THEN
1156 ! odd derivatives are anti-symmetric
1157 symmetry_string = dbcsr_type_antisymmetric
1158 ELSE
1159 symmetry_string = dbcsr_type_symmetric
1160 END IF
1161 ELSE
1162 symmetry_string = dbcsr_type_no_symmetry
1163 END IF
1164 cgfsym = cgf_symbol(1, indco(1:3, i1))
1165 IF (i1 == 1) THEN
1166 name = mname
1167 ELSE
1168 name = trim(cgfsym(4:))//" DERIVATIVE OF THE "//trim(mname)// &
1169 " W.R.T. THE NUCLEAR COORDINATES"
1170 END IF
1171 CALL compress(name)
1172 CALL uppercase(name)
1173 ALLOCATE (matrix_s(i1, i2)%matrix)
1174 CALL dbcsr_create(matrix=matrix_s(i1, i2)%matrix, &
1175 name=trim(name), &
1176 dist=dbcsr_dist, matrix_type=symmetry_string, &
1177 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, &
1178 nze=0)
1179 CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i1, i2)%matrix, sab_nl)
1180 END DO
1181 END DO
1182
1183 DEALLOCATE (row_blk_sizes, col_blk_sizes)
1184
1185 END SUBROUTINE create_sab_matrix_2d
1186
1187END 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_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:81
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:680
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)
...
Definition basis_set_types.F:615
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_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:333
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:98
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, 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, 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_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, 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:330
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:993
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
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)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:120
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:784
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:558
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:3190
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3361
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:139
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:93
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
block_p_types::block_p_type
Definition block_p_types.F:23
cp_control_types::dft_control_type
Definition cp_control_types.F:559
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:138
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:171
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:135
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