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 35520 : 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 35520 : 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 38208 : matrix_p, matrixkp_p, ext_kpoints, natom)
147 :
148 35520 : 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 35520 : 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 35520 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
193 : INTEGER, DIMENSION(3) :: cell
194 35520 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
195 35520 : npgfb, nsgfa, nsgfb
196 35520 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
197 35520 : 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 35520 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork, pmat
202 35520 : 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 71040 : REAL(KIND=dp), DIMENSION(3, natom) :: force_thread
206 35520 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
207 35520 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
208 35520 : zeta, zetb
209 35520 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
210 35520 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint
211 : TYPE(dft_control_type), POINTER :: dft_control
212 35520 : 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 35520 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
216 35520 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
217 : TYPE(virial_type), POINTER :: virial
218 :
219 : !$ INTEGER(kind=omp_lock_kind), &
220 35520 : !$ ALLOCATABLE, DIMENSION(:) :: locks
221 : !$ INTEGER :: lock_num, hash, hash1, hash2
222 : !$ INTEGER(KIND=int_8) :: iatom8
223 : !$ INTEGER, PARAMETER :: nlock = 501
224 :
225 35520 : CALL timeset(routineN, handle)
226 :
227 35520 : 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 35520 : dft_control=dft_control)
232 :
233 : ! test for matrices (kpoints or standard gamma point)
234 35520 : IF (PRESENT(matrix_s)) THEN
235 15083 : dokp = .FALSE.
236 15083 : use_cell_mapping = .FALSE.
237 15083 : nimg = dft_control%nimages
238 20437 : ELSEIF (PRESENT(matrixkp_s)) THEN
239 20437 : dokp = .TRUE.
240 20437 : IF (PRESENT(ext_kpoints)) THEN
241 0 : IF (ASSOCIATED(ext_kpoints)) THEN
242 0 : CALL get_kpoint_info(kpoint=ext_kpoints, cell_to_index=cell_to_index)
243 0 : use_cell_mapping = (SIZE(cell_to_index) > 1)
244 0 : nimg = SIZE(ext_kpoints%index_to_cell, 2)
245 : ELSE
246 0 : use_cell_mapping = .FALSE.
247 0 : nimg = 1
248 : END IF
249 : ELSE
250 20437 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
251 20437 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
252 81748 : use_cell_mapping = (SIZE(cell_to_index) > 1)
253 20437 : nimg = dft_control%nimages
254 : END IF
255 : ELSE
256 0 : CPABORT("")
257 : END IF
258 :
259 35520 : nkind = SIZE(qs_kind_set)
260 :
261 35520 : IF (PRESENT(calculate_forces)) THEN
262 21269 : do_forces = calculate_forces
263 : ELSE
264 : do_forces = .FALSE.
265 : END IF
266 :
267 35520 : IF (PRESENT(nderivative)) THEN
268 21743 : nder = nderivative
269 : ELSE
270 : nder = 0
271 : END IF
272 35520 : maxder = ncoset(nder)
273 :
274 : ! check for symmetry
275 35520 : CPASSERT(SIZE(sab_nl) > 0)
276 35520 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
277 35520 : IF (do_symmetric) THEN
278 34356 : CPASSERT(basis_type_a == basis_type_b)
279 : END IF
280 :
281 : ! set up basis set lists
282 277952 : ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
283 35520 : CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set)
284 35520 : CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set)
285 :
286 35520 : IF (dokp) THEN
287 20437 : 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 20437 : sab_nl, do_symmetric)
290 : ELSE
291 15083 : 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 17771 : sab_nl, do_symmetric)
294 : END IF
295 35520 : maxs = maxder
296 :
297 35520 : use_virial = .FALSE.
298 35520 : IF (do_forces) THEN
299 10041 : CALL get_ks_env(ks_env=ks_env, force=force, virial=virial)
300 10041 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
301 : END IF
302 :
303 669352 : force_thread = 0.0_dp
304 35520 : pv_thread = 0.0_dp
305 :
306 35520 : ldsab = get_memory_usage(qs_kind_set, basis_type_a, basis_type_b)
307 35520 : IF (do_forces) THEN
308 : ! we need density matrix for forces
309 10041 : IF (dokp) THEN
310 6503 : CPASSERT(PRESENT(matrixkp_p))
311 : ELSE
312 3538 : CPASSERT(PRESENT(matrix_p))
313 : END IF
314 10041 : nder = MAX(nder, 1)
315 : END IF
316 35520 : 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 35520 : !$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 35520 : IF (do_forces) THEN
522 10041 : 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 35241 : atom_a = atom_of_kind(iatom)
526 35241 : ikind = kind_of(iatom)
527 140964 : force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) + force_thread(:, iatom)
528 : END DO
529 : !$OMP END DO
530 : END IF
531 10041 : IF (do_forces .AND. use_virial) THEN
532 10972 : virial%pv_overlap = virial%pv_overlap + pv_thread
533 10972 : virial%pv_virial = virial%pv_virial + pv_thread
534 : END IF
535 :
536 35520 : IF (dokp) THEN
537 46502 : DO i = 1, maxs
538 127545 : DO ic = 1, nimg
539 81043 : CALL dbcsr_finalize(matrixkp_s(i, ic)%matrix)
540 : CALL dbcsr_filter(matrixkp_s(i, ic)%matrix, &
541 107108 : dft_control%qs_control%eps_filter_matrix)
542 : END DO
543 : END DO
544 : ELSE
545 45166 : DO i = 1, maxs
546 30083 : CALL dbcsr_finalize(matrix_s(i)%matrix)
547 : CALL dbcsr_filter(matrix_s(i)%matrix, &
548 45166 : dft_control%qs_control%eps_filter_matrix)
549 : END DO
550 : END IF
551 :
552 : ! *** Release work storage ***
553 35520 : DEALLOCATE (basis_set_list_a, basis_set_list_b)
554 :
555 35520 : CALL timestop(handle)
556 :
557 106560 : END SUBROUTINE build_overlap_matrix_low
558 :
559 : ! **************************************************************************************************
560 : !> \brief Calculation of the overlap matrix over Cartesian Gaussian functions.
561 : !> \param ks_env the QS env
562 : !> \param matrix_s The overlap matrix to be calculated
563 : !> \param basis_set_list_a basis set list to be used for bra in <a|b>
564 : !> \param basis_set_list_b basis set list to be used for ket in <a|b>
565 : !> \param sab_nl pair list (must be consistent with basis sets!)
566 : !> \date 11.03.2016
567 : !> \par History
568 : !> Simplified version of build_overlap_matrix
569 : !> \author MK
570 : !> \version 1.0
571 : ! **************************************************************************************************
572 162 : SUBROUTINE build_overlap_matrix_simple(ks_env, matrix_s, &
573 : basis_set_list_a, basis_set_list_b, sab_nl)
574 :
575 : TYPE(qs_ks_env_type), POINTER :: ks_env
576 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
577 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
578 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
579 : POINTER :: sab_nl
580 :
581 : CHARACTER(len=*), PARAMETER :: routineN = 'build_overlap_matrix_simple'
582 :
583 : INTEGER :: handle, iatom, icol, ikind, irow, iset, &
584 : jatom, jkind, jset, ldsab, m1, m2, n1, &
585 : n2, natom, ncoa, ncob, nkind, nseta, &
586 : nsetb, sgfa, sgfb, slot
587 162 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
588 162 : npgfb, nsgfa, nsgfb
589 162 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
590 : LOGICAL :: do_symmetric, found, trans
591 : REAL(KIND=dp) :: dab, rab2
592 162 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
593 162 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint
594 : REAL(KIND=dp), DIMENSION(3) :: rab
595 162 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
596 162 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
597 162 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
598 162 : TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: sint
599 : TYPE(dft_control_type), POINTER :: dft_control
600 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
601 162 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
602 :
603 : !$ INTEGER(kind=omp_lock_kind), &
604 162 : !$ ALLOCATABLE, DIMENSION(:) :: locks
605 : !$ INTEGER :: lock_num, hash, hash1, hash2
606 : !$ INTEGER(KIND=int_8) :: iatom8
607 : !$ INTEGER, PARAMETER :: nlock = 501
608 :
609 162 : NULLIFY (dft_control)
610 :
611 162 : CALL timeset(routineN, handle)
612 :
613 162 : NULLIFY (atomic_kind_set)
614 : CALL get_ks_env(ks_env, &
615 : atomic_kind_set=atomic_kind_set, &
616 : natom=natom, &
617 : qs_kind_set=qs_kind_set, &
618 162 : dft_control=dft_control)
619 :
620 : ! check for symmetry
621 162 : CPASSERT(SIZE(sab_nl) > 0)
622 162 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
623 :
624 162 : nkind = SIZE(qs_kind_set)
625 :
626 162 : CALL dbcsr_allocate_matrix_set(matrix_s, 1)
627 : CALL create_sab_matrix(ks_env, matrix_s, "Matrix", basis_set_list_a, basis_set_list_b, &
628 162 : sab_nl, do_symmetric)
629 :
630 162 : ldsab = 0
631 474 : DO ikind = 1, nkind
632 312 : basis_set_a => basis_set_list_a(ikind)%gto_basis_set
633 312 : CALL get_gto_basis_set(gto_basis_set=basis_set_a, maxco=m1, nsgf=m2)
634 312 : ldsab = MAX(m1, m2, ldsab)
635 312 : basis_set_b => basis_set_list_b(ikind)%gto_basis_set
636 312 : CALL get_gto_basis_set(gto_basis_set=basis_set_b, maxco=m1, nsgf=m2)
637 474 : ldsab = MAX(m1, m2, ldsab)
638 : END DO
639 :
640 : !$OMP PARALLEL DEFAULT(NONE) &
641 : !$OMP SHARED (ldsab,sab_nl,do_symmetric,ncoset,natom,&
642 : !$OMP matrix_s,basis_set_list_a,basis_set_list_b,locks) &
643 : !$OMP PRIVATE (oint,owork,sint,ikind,jkind,iatom,jatom,rab,basis_set_a,basis_set_b,&
644 : !$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, &
645 : !$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, dab, &
646 : !$OMP zetb, scon_a, scon_b, irow, icol, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
647 162 : !$OMP slot, lock_num, hash, hash1, hash2, iatom8 )
648 :
649 : !$OMP SINGLE
650 : !$ ALLOCATE (locks(nlock))
651 : !$OMP END SINGLE
652 :
653 : !$OMP DO
654 : !$ DO lock_num = 1, nlock
655 : !$ call omp_init_lock(locks(lock_num))
656 : !$ END DO
657 : !$OMP END DO
658 :
659 : ALLOCATE (oint(ldsab, ldsab, 1), owork(ldsab, ldsab))
660 : ALLOCATE (sint(1))
661 : NULLIFY (sint(1)%block)
662 :
663 : !$OMP DO SCHEDULE(GUIDED)
664 : DO slot = 1, sab_nl(1)%nl_size
665 : ikind = sab_nl(1)%nlist_task(slot)%ikind
666 : jkind = sab_nl(1)%nlist_task(slot)%jkind
667 : iatom = sab_nl(1)%nlist_task(slot)%iatom
668 : jatom = sab_nl(1)%nlist_task(slot)%jatom
669 : rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
670 : !$ iatom8 = (iatom - 1)*natom + jatom
671 : !$ hash1 = INT(MOD(iatom8, INT(nlock, int_8)) + 1)
672 :
673 : basis_set_a => basis_set_list_a(ikind)%gto_basis_set
674 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
675 : basis_set_b => basis_set_list_b(jkind)%gto_basis_set
676 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
677 : ! basis ikind
678 : first_sgfa => basis_set_a%first_sgf
679 : la_max => basis_set_a%lmax
680 : la_min => basis_set_a%lmin
681 : npgfa => basis_set_a%npgf
682 : nseta = basis_set_a%nset
683 : nsgfa => basis_set_a%nsgf_set
684 : rpgfa => basis_set_a%pgf_radius
685 : set_radius_a => basis_set_a%set_radius
686 : scon_a => basis_set_a%scon
687 : zeta => basis_set_a%zet
688 : ! basis jkind
689 : first_sgfb => basis_set_b%first_sgf
690 : lb_max => basis_set_b%lmax
691 : lb_min => basis_set_b%lmin
692 : npgfb => basis_set_b%npgf
693 : nsetb = basis_set_b%nset
694 : nsgfb => basis_set_b%nsgf_set
695 : rpgfb => basis_set_b%pgf_radius
696 : set_radius_b => basis_set_b%set_radius
697 : scon_b => basis_set_b%scon
698 : zetb => basis_set_b%zet
699 :
700 : IF (do_symmetric) THEN
701 : IF (iatom <= jatom) THEN
702 : irow = iatom
703 : icol = jatom
704 : ELSE
705 : irow = jatom
706 : icol = iatom
707 : END IF
708 : ELSE
709 : irow = iatom
710 : icol = jatom
711 : END IF
712 :
713 : NULLIFY (sint(1)%block)
714 : CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
715 : row=irow, col=icol, BLOCK=sint(1)%block, found=found)
716 : CPASSERT(found)
717 : trans = do_symmetric .AND. (iatom > jatom)
718 :
719 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
720 : dab = SQRT(rab2)
721 :
722 : DO iset = 1, nseta
723 :
724 : ncoa = npgfa(iset)*ncoset(la_max(iset))
725 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
726 : sgfa = first_sgfa(1, iset)
727 :
728 : DO jset = 1, nsetb
729 :
730 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
731 :
732 : !$ hash2 = MOD((iset - 1)*nsetb + jset, nlock) + 1
733 : !$ hash = MOD(hash1 + hash2, nlock) + 1
734 :
735 : ncob = npgfb(jset)*ncoset(lb_max(jset))
736 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
737 : sgfb = first_sgfb(1, jset)
738 :
739 : ! calculate integrals and derivatives
740 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
741 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
742 : rab, sab=oint(:, :, 1))
743 : ! Contraction
744 : CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
745 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=trans)
746 : !$OMP CRITICAL(blockadd)
747 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(1)%block, &
748 : sgfa, sgfb, trans=trans)
749 : !$OMP END CRITICAL(blockadd)
750 :
751 : END DO
752 : END DO
753 :
754 : END DO
755 : DEALLOCATE (oint, owork)
756 : DEALLOCATE (sint)
757 : !$OMP DO
758 : !$ DO lock_num = 1, nlock
759 : !$ call omp_destroy_lock(locks(lock_num))
760 : !$ END DO
761 : !$OMP END DO
762 :
763 : !$OMP SINGLE
764 : !$ DEALLOCATE (locks)
765 : !$OMP END SINGLE NOWAIT
766 :
767 : !$OMP END PARALLEL
768 :
769 162 : CALL dbcsr_finalize(matrix_s(1)%matrix)
770 162 : CALL dbcsr_filter(matrix_s(1)%matrix, dft_control%qs_control%eps_filter_matrix)
771 :
772 162 : CALL timestop(handle)
773 :
774 324 : END SUBROUTINE build_overlap_matrix_simple
775 :
776 : ! **************************************************************************************************
777 : !> \brief Calculation of the force contribution from an overlap matrix
778 : !> over Cartesian Gaussian functions.
779 : !> \param ks_env the QS environment
780 : !> \param force holds the calcuated force Tr(P dS/dR)
781 : !> \param basis_type_a basis set to be used for bra in <a|b>
782 : !> \param basis_type_b basis set to be used for ket in <a|b>
783 : !> \param sab_nl pair list (must be consistent with basis sets!)
784 : !> \param matrix_p density matrix for force calculation
785 : !> \param matrixkp_p ...
786 : !> \date 11.03.2002
787 : !> \par History
788 : !> Enlarged functionality of this routine. Now overlap matrices based
789 : !> on different basis sets can be calculated, taking into account also
790 : !> mixed overlaps NOTE: the pointer to the overlap matrix must now be
791 : !> put into its corresponding env outside of this routine
792 : !> [Manuel Guidon]
793 : !> Generalized for derivatives and force calculations [JHU]
794 : !> This special version is only for forces [07.07.2014, JGH]
795 : !> \author MK/JGH
796 : !> \version 1.0
797 : ! **************************************************************************************************
798 2464 : SUBROUTINE build_overlap_force(ks_env, force, basis_type_a, basis_type_b, &
799 2464 : sab_nl, matrix_p, matrixkp_p)
800 :
801 : TYPE(qs_ks_env_type), POINTER :: ks_env
802 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force
803 : CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b
804 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
805 : POINTER :: sab_nl
806 : TYPE(dbcsr_type), OPTIONAL :: matrix_p
807 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL :: matrixkp_p
808 :
809 : CHARACTER(len=*), PARAMETER :: routineN = 'build_overlap_force'
810 :
811 : INTEGER :: handle, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, n1, n2, &
812 : natom, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
813 : INTEGER, DIMENSION(3) :: cell
814 2464 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
815 2464 : npgfb, nsgfa, nsgfb
816 2464 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
817 2464 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
818 : LOGICAL :: do_symmetric, dokp, found, trans, &
819 : use_cell_mapping, use_virial
820 : REAL(KIND=dp) :: dab, f0, ff, rab2
821 2464 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, sab
822 2464 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: drab
823 : REAL(KIND=dp), DIMENSION(3) :: force_a, rab
824 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_thread
825 2464 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
826 4928 : REAL(KIND=dp), DIMENSION(3, SIZE(force, 2)) :: force_thread
827 2464 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
828 2464 : zeta, zetb
829 : TYPE(dft_control_type), POINTER :: dft_control
830 2464 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
831 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
832 : TYPE(kpoint_type), POINTER :: kpoints
833 2464 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
834 : TYPE(virial_type), POINTER :: virial
835 :
836 2464 : CALL timeset(routineN, handle)
837 :
838 2464 : NULLIFY (qs_kind_set)
839 2464 : CALL get_ks_env(ks_env=ks_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
840 2464 : nimg = dft_control%nimages
841 :
842 : ! test for matrices (kpoints or standard gamma point)
843 2464 : IF (PRESENT(matrix_p)) THEN
844 : dokp = .FALSE.
845 : use_cell_mapping = .FALSE.
846 64 : ELSEIF (PRESENT(matrixkp_p)) THEN
847 64 : dokp = .TRUE.
848 64 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
849 64 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
850 256 : use_cell_mapping = (SIZE(cell_to_index) > 1)
851 : ELSE
852 0 : CPABORT("")
853 : END IF
854 :
855 2464 : nkind = SIZE(qs_kind_set)
856 2464 : nder = 1
857 :
858 : ! check for symmetry
859 2464 : CPASSERT(SIZE(sab_nl) > 0)
860 2464 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
861 :
862 2464 : CALL get_ks_env(ks_env=ks_env, virial=virial)
863 2464 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
864 2464 : virial_thread = 0.0_dp
865 :
866 : ! set up basis sets
867 18824 : ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
868 2464 : CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set)
869 2464 : CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set)
870 2464 : ldsab = get_memory_usage(qs_kind_set, basis_type_a, basis_type_b)
871 :
872 2464 : natom = SIZE(force, 2)
873 30800 : force_thread = 0.0_dp
874 :
875 : !$OMP PARALLEL DEFAULT(NONE) &
876 : !$OMP SHARED (ldsab, do_symmetric, ncoset, nder, use_virial, force, virial, ndod, sab_nl, cell_to_index, &
877 : !$OMP matrix_p, basis_set_list_a, basis_set_list_b, dokp, use_cell_mapping, nimg, matrixkp_p) &
878 : !$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, ic, &
879 : !$OMP basis_set_a, basis_set_b, sab, pab, drab, &
880 : !$OMP first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a, ncoa, ncob, force_a, &
881 : !$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb, p_block, dab, &
882 : !$OMP zetb, scon_a, scon_b, irow, icol, f0, ff, found, trans, rab2, n1, n2, sgfa, sgfb, iset, jset, &
883 : !$OMP slot, cell) &
884 2464 : !$OMP REDUCTION (+ : virial_thread, force_thread )
885 :
886 : ! *** Allocate work storage ***
887 : ALLOCATE (sab(ldsab, ldsab), pab(ldsab, ldsab))
888 : ALLOCATE (drab(ldsab, ldsab, 3))
889 :
890 : ! Loop over neighbor list
891 : !$OMP DO SCHEDULE(GUIDED)
892 : DO slot = 1, sab_nl(1)%nl_size
893 : ikind = sab_nl(1)%nlist_task(slot)%ikind
894 : jkind = sab_nl(1)%nlist_task(slot)%jkind
895 : iatom = sab_nl(1)%nlist_task(slot)%iatom
896 : jatom = sab_nl(1)%nlist_task(slot)%jatom
897 : cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
898 : rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
899 :
900 : basis_set_a => basis_set_list_a(ikind)%gto_basis_set
901 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
902 : basis_set_b => basis_set_list_b(jkind)%gto_basis_set
903 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
904 : ! basis ikind
905 : first_sgfa => basis_set_a%first_sgf
906 : la_max => basis_set_a%lmax
907 : la_min => basis_set_a%lmin
908 : npgfa => basis_set_a%npgf
909 : nseta = basis_set_a%nset
910 : nsgfa => basis_set_a%nsgf_set
911 : rpgfa => basis_set_a%pgf_radius
912 : set_radius_a => basis_set_a%set_radius
913 : scon_a => basis_set_a%scon
914 : zeta => basis_set_a%zet
915 : ! basis jkind
916 : first_sgfb => basis_set_b%first_sgf
917 : lb_max => basis_set_b%lmax
918 : lb_min => basis_set_b%lmin
919 : npgfb => basis_set_b%npgf
920 : nsetb = basis_set_b%nset
921 : nsgfb => basis_set_b%nsgf_set
922 : rpgfb => basis_set_b%pgf_radius
923 : set_radius_b => basis_set_b%set_radius
924 : scon_b => basis_set_b%scon
925 : zetb => basis_set_b%zet
926 :
927 : IF (use_cell_mapping) THEN
928 : ic = cell_to_index(cell(1), cell(2), cell(3))
929 : IF (ic < 1 .OR. ic > nimg) CYCLE
930 : ELSE
931 : ic = 1
932 : END IF
933 :
934 : IF (do_symmetric) THEN
935 : IF (iatom <= jatom) THEN
936 : irow = iatom
937 : icol = jatom
938 : ELSE
939 : irow = jatom
940 : icol = iatom
941 : END IF
942 : f0 = 2.0_dp
943 : IF (iatom == jatom) f0 = 1.0_dp
944 : ff = 2.0_dp
945 : ELSE
946 : irow = iatom
947 : icol = jatom
948 : f0 = 1.0_dp
949 : ff = 1.0_dp
950 : END IF
951 : NULLIFY (p_block)
952 : IF (dokp) THEN
953 : CALL dbcsr_get_block_p(matrix=matrixkp_p(ic)%matrix, row=irow, col=icol, &
954 : block=p_block, found=found)
955 : ELSE
956 : CALL dbcsr_get_block_p(matrix=matrix_p, row=irow, col=icol, &
957 : block=p_block, found=found)
958 : END IF
959 :
960 : trans = do_symmetric .AND. (iatom > jatom)
961 :
962 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
963 : dab = SQRT(rab2)
964 :
965 : DO iset = 1, nseta
966 :
967 : ncoa = npgfa(iset)*ncoset(la_max(iset))
968 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
969 : sgfa = first_sgfa(1, iset)
970 :
971 : DO jset = 1, nsetb
972 :
973 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
974 :
975 : ncob = npgfb(jset)*ncoset(lb_max(jset))
976 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
977 : sgfb = first_sgfb(1, jset)
978 :
979 : IF (ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial)) THEN
980 : ! Decontract P matrix block
981 : sab = 0.0_dp
982 : CALL block_add("OUT", sab, nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
983 : CALL decontraction(sab, pab, scon_a(:, sgfa:), n1, nsgfa(iset), scon_b(:, sgfb:), n2, &
984 : nsgfb(jset), trans=trans)
985 : ! calculate integrals and derivatives
986 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
987 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
988 : rab, dab=drab)
989 : CALL force_trace(force_a, drab, pab, n1, n2, 3)
990 : force_thread(1:3, iatom) = force_thread(1:3, iatom) - ff*force_a(1:3)
991 : force_thread(1:3, jatom) = force_thread(1:3, jatom) + ff*force_a(1:3)
992 : IF (use_virial) THEN
993 : CALL virial_pair_force(virial_thread, -f0, force_a, rab)
994 : END IF
995 : END IF
996 :
997 : END DO
998 : END DO
999 :
1000 : END DO
1001 : DEALLOCATE (sab, pab, drab)
1002 : !$OMP END PARALLEL
1003 :
1004 : !$OMP PARALLEL DEFAULT(NONE) &
1005 2464 : !$OMP SHARED(force,force_thread,natom)
1006 : !$OMP WORKSHARE
1007 : force(1:3, 1:natom) = force(1:3, 1:natom) + force_thread(1:3, 1:natom)
1008 : !$OMP END WORKSHARE
1009 : !$OMP END PARALLEL
1010 2464 : IF (use_virial) THEN
1011 3432 : virial%pv_virial = virial%pv_virial + virial_thread
1012 3432 : virial%pv_overlap = virial%pv_overlap + virial_thread
1013 : END IF
1014 :
1015 2464 : DEALLOCATE (basis_set_list_a, basis_set_list_b)
1016 :
1017 2464 : CALL timestop(handle)
1018 :
1019 7392 : END SUBROUTINE build_overlap_force
1020 :
1021 : ! **************************************************************************************************
1022 : !> \brief Setup the structure of a sparse matrix based on the overlap
1023 : !> neighbor list
1024 : !> \param ks_env The QS environment
1025 : !> \param matrix_s Matrices to be constructed
1026 : !> \param matrix_name Matrix base name
1027 : !> \param basis_set_list_a Basis set used for <a|
1028 : !> \param basis_set_list_b Basis set used for |b>
1029 : !> \param sab_nl Overlap neighbor list
1030 : !> \param symmetric Is symmetry used in the neighbor list?
1031 : ! **************************************************************************************************
1032 15331 : SUBROUTINE create_sab_matrix_1d(ks_env, matrix_s, matrix_name, &
1033 : basis_set_list_a, basis_set_list_b, sab_nl, symmetric)
1034 :
1035 : TYPE(qs_ks_env_type), POINTER :: ks_env
1036 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1037 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
1038 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
1039 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1040 : POINTER :: sab_nl
1041 : LOGICAL, INTENT(IN) :: symmetric
1042 :
1043 : CHARACTER(LEN=12) :: cgfsym
1044 : CHARACTER(LEN=32) :: symmetry_string
1045 : CHARACTER(LEN=default_string_length) :: mname, name
1046 : INTEGER :: i, maxs, natom
1047 15331 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
1048 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1049 15331 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1050 15331 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1051 :
1052 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, &
1053 15331 : qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist)
1054 :
1055 15331 : natom = SIZE(particle_set)
1056 :
1057 15331 : IF (PRESENT(matrix_name)) THEN
1058 12643 : mname = matrix_name
1059 : ELSE
1060 2688 : mname = "DUMMY"
1061 : END IF
1062 :
1063 15331 : maxs = SIZE(matrix_s)
1064 :
1065 61324 : ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
1066 :
1067 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1068 15331 : basis=basis_set_list_a)
1069 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, &
1070 15331 : basis=basis_set_list_b)
1071 :
1072 : ! prepare for allocation
1073 15331 : IF (symmetric) THEN
1074 15039 : symmetry_string = dbcsr_type_symmetric
1075 : ELSE
1076 292 : symmetry_string = dbcsr_type_no_symmetry
1077 : END IF
1078 :
1079 45734 : DO i = 1, maxs
1080 30403 : IF (symmetric) THEN
1081 30111 : IF (ndod(i) == 1) THEN
1082 : ! odd derivatives are anti-symmetric
1083 13488 : symmetry_string = dbcsr_type_antisymmetric
1084 : ELSE
1085 16623 : symmetry_string = dbcsr_type_symmetric
1086 : END IF
1087 : ELSE
1088 292 : symmetry_string = dbcsr_type_no_symmetry
1089 : END IF
1090 30403 : cgfsym = cgf_symbol(1, indco(1:3, i))
1091 30403 : IF (i == 1) THEN
1092 15331 : name = mname
1093 : ELSE
1094 : name = TRIM(cgfsym(4:))//" DERIVATIVE OF THE "//TRIM(mname)// &
1095 15072 : " W.R.T. THE NUCLEAR COORDINATES"
1096 : END IF
1097 30403 : CALL compress(name)
1098 30403 : CALL uppercase(name)
1099 30403 : ALLOCATE (matrix_s(i)%matrix)
1100 : CALL dbcsr_create(matrix=matrix_s(i)%matrix, &
1101 : name=TRIM(name), &
1102 : dist=dbcsr_dist, matrix_type=symmetry_string, &
1103 30403 : row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
1104 45734 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i)%matrix, sab_nl)
1105 : END DO
1106 :
1107 15331 : DEALLOCATE (row_blk_sizes, col_blk_sizes)
1108 :
1109 15331 : END SUBROUTINE create_sab_matrix_1d
1110 :
1111 : ! **************************************************************************************************
1112 : !> \brief Setup the structure of a sparse matrix based on the overlap
1113 : !> neighbor list, 2d version
1114 : !> \param ks_env The QS environment
1115 : !> \param matrix_s Matrices to be constructed
1116 : !> \param matrix_name Matrix base name
1117 : !> \param basis_set_list_a Basis set used for <a|
1118 : !> \param basis_set_list_b Basis set used for |b>
1119 : !> \param sab_nl Overlap neighbor list
1120 : !> \param symmetric Is symmetry used in the neighbor list?
1121 : ! **************************************************************************************************
1122 46104 : SUBROUTINE create_sab_matrix_2d(ks_env, matrix_s, matrix_name, &
1123 : basis_set_list_a, basis_set_list_b, sab_nl, symmetric)
1124 :
1125 : TYPE(qs_ks_env_type), POINTER :: ks_env
1126 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1127 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: matrix_name
1128 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
1129 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1130 : POINTER :: sab_nl
1131 : LOGICAL, INTENT(IN) :: symmetric
1132 :
1133 : CHARACTER(LEN=12) :: cgfsym
1134 : CHARACTER(LEN=32) :: symmetry_string
1135 : CHARACTER(LEN=default_string_length) :: mname, name
1136 : INTEGER :: i1, i2, natom
1137 46104 : INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes
1138 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1139 46104 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1140 46104 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1141 :
1142 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, &
1143 46104 : qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist)
1144 :
1145 46104 : natom = SIZE(particle_set)
1146 :
1147 46104 : IF (PRESENT(matrix_name)) THEN
1148 45346 : mname = matrix_name
1149 : ELSE
1150 758 : mname = "DUMMY"
1151 : END IF
1152 :
1153 184416 : ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
1154 :
1155 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1156 46104 : basis=basis_set_list_a)
1157 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes, &
1158 46104 : basis=basis_set_list_b)
1159 :
1160 : ! prepare for allocation
1161 46104 : IF (symmetric) THEN
1162 45174 : symmetry_string = dbcsr_type_symmetric
1163 : ELSE
1164 930 : symmetry_string = dbcsr_type_no_symmetry
1165 : END IF
1166 :
1167 267980 : DO i2 = 1, SIZE(matrix_s, 2)
1168 528304 : DO i1 = 1, SIZE(matrix_s, 1)
1169 260324 : IF (symmetric) THEN
1170 256854 : IF (ndod(i1) == 1) THEN
1171 : ! odd derivatives are anti-symmetric
1172 38448 : symmetry_string = dbcsr_type_antisymmetric
1173 : ELSE
1174 218406 : symmetry_string = dbcsr_type_symmetric
1175 : END IF
1176 : ELSE
1177 3470 : symmetry_string = dbcsr_type_no_symmetry
1178 : END IF
1179 260324 : cgfsym = cgf_symbol(1, indco(1:3, i1))
1180 260324 : IF (i1 == 1) THEN
1181 221876 : name = mname
1182 : ELSE
1183 : name = TRIM(cgfsym(4:))//" DERIVATIVE OF THE "//TRIM(mname)// &
1184 38448 : " W.R.T. THE NUCLEAR COORDINATES"
1185 : END IF
1186 260324 : CALL compress(name)
1187 260324 : CALL uppercase(name)
1188 260324 : ALLOCATE (matrix_s(i1, i2)%matrix)
1189 : CALL dbcsr_create(matrix=matrix_s(i1, i2)%matrix, &
1190 : name=TRIM(name), &
1191 : dist=dbcsr_dist, matrix_type=symmetry_string, &
1192 260324 : row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes)
1193 482200 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_s(i1, i2)%matrix, sab_nl)
1194 : END DO
1195 : END DO
1196 :
1197 46104 : DEALLOCATE (row_blk_sizes, col_blk_sizes)
1198 :
1199 46104 : END SUBROUTINE create_sab_matrix_2d
1200 :
1201 : END MODULE qs_overlap
|
