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