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