(git:b279b6b)
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 ! **************************************************************************************************
19 MODULE 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 
90 CONTAINS
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 
1187 END MODULE qs_overlap
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.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