(git:4dc5213)
Loading...
Searching...
No Matches
qs_operators_ao.F
Go to the documentation of this file.
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!> \par History
10!> created 07.2005
11!> \author MI (07.2005)
12! **************************************************************************************************
13MODULE qs_operators_ao
14 USE ai_angmom, ONLY: angmom
15 USE ai_moments, ONLY: contract_cossin,&
16 cossin
17 USE ai_overlap, ONLY: overlap_ab
18 USE basis_set_types, ONLY: gto_basis_set_p_type,&
19 gto_basis_set_type
20 USE block_p_types, ONLY: block_p_type
21 USE cell_types, ONLY: cell_type,&
22 pbc
23 USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
24 dbcsr_get_matrix_type,&
25 dbcsr_has_symmetry,&
26 dbcsr_p_type,&
27 dbcsr_set,&
28 dbcsr_type_antisymmetric
29 USE kinds, ONLY: dp
30 USE orbital_pointers, ONLY: init_orbital_pointers,&
31 ncoset
32 USE particle_types, ONLY: particle_type
33 USE qs_environment_types, ONLY: get_qs_env,&
34 qs_environment_type
35 USE qs_kind_types, ONLY: get_qs_kind,&
36 get_qs_kind_set,&
37 qs_kind_type
38 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
39 neighbor_list_iterate,&
40 neighbor_list_iterator_create,&
41 neighbor_list_iterator_p_type,&
42 neighbor_list_iterator_release,&
43 neighbor_list_set_p_type
44#include "./base/base_uses.f90"
45
46 IMPLICIT NONE
47 PRIVATE
48
49 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_operators_ao'
50
51! *** Public subroutines ***
52
53 PUBLIC :: build_exp_ikr_matrix
54 PUBLIC :: build_lin_mom_matrix, build_ang_mom_matrix
55
56CONTAINS
57
58! **************************************************************************************************
59!> \brief Build real and imaginary AO matrices for exp(i*k*r).
60!> \param qs_env ...
61!> \param op_sm_set op_sm_set(1, i) is the cosine matrix and op_sm_set(2, i) is the sine matrix
62!> for kvec(:, i)
63!> \param kvec Cartesian wave vectors, one column for each output matrix pair
64!> \param sab_orb_external optional neighbor list used instead of the default orbital list
65!> \param basis_type optional basis-set name
66!> \param force_periodic temporarily use all three periodic directions
67!> \param cell_external optional cell used for coordinate folding and periodicity
68!> \author CP2K developers
69!> \note The cosine and sine matrices must be preallocated with matching symmetric
70!> DBCSR topology. The matrices are overwritten from zero and only canonical
71!> atom-pair blocks are written.
72! **************************************************************************************************
73 SUBROUTINE build_exp_ikr_matrix(qs_env, op_sm_set, kvec, sab_orb_external, basis_type, force_periodic, cell_external)
74
75 TYPE(qs_environment_type), POINTER :: qs_env
76 TYPE(dbcsr_p_type), DIMENSION(:, :) :: op_sm_set
77 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: kvec
78 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
79 OPTIONAL, POINTER :: sab_orb_external
80 CHARACTER(LEN=*), OPTIONAL :: basis_type
81 LOGICAL, OPTIONAL :: force_periodic
82 TYPE(cell_type), OPTIONAL, POINTER :: cell_external
83
84 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_exp_ikr_matrix'
85
86 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
87 ldsa, ldsb, ldwork, ncoa, ncob, nkind, nkvec, nseta, nsetb, reim, sgfa, sgfb
88 INTEGER, DIMENSION(3) :: perd0
89 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
90 npgfb, nsgfa, nsgfb
91 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
92 LOGICAL :: found, my_force_periodic, new_atom_b
93 REAL(kind=dp) :: dab
94 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
95 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
96 REAL(kind=dp), DIMENSION(:, :), POINTER :: cosab, rpgfa, rpgfb, sinab, sphi_a, &
97 sphi_b, work, zeta, zetb
98 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: op_cossin
99 TYPE(cell_type), POINTER :: cell
100 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
101 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
102 TYPE(neighbor_list_iterator_p_type), &
103 DIMENSION(:), POINTER :: nl_iterator
104 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
105 POINTER :: sab_orb
106 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
107 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
108 TYPE(qs_kind_type), POINTER :: qs_kind
109
110 CALL timeset(routinen, handle)
111
112 cpassert(SIZE(kvec, 1) == 3)
113 nkvec = SIZE(kvec, 2)
114 cpassert(nkvec > 0)
115 cpassert(SIZE(op_sm_set, 1) == 2)
116 cpassert(SIZE(op_sm_set, 2) == nkvec)
117 DO i = 1, nkvec
118 DO reim = 1, 2
119 cpassert(ASSOCIATED(op_sm_set(reim, i)%matrix))
120 cpassert(dbcsr_has_symmetry(op_sm_set(reim, i)%matrix))
121 CALL dbcsr_set(op_sm_set(reim, i)%matrix, 0.0_dp)
122 END DO
123 END DO
124
125 NULLIFY (qs_kind, qs_kind_set, particle_set, sab_orb, cell)
126 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
127 particle_set=particle_set, cell=cell, sab_orb=sab_orb)
128 cpassert(ASSOCIATED(cell))
129 cpassert(ASSOCIATED(qs_kind_set))
130 cpassert(ASSOCIATED(particle_set))
131 IF (PRESENT(cell_external)) THEN
132 cpassert(ASSOCIATED(cell_external))
133 cell => cell_external
134 END IF
135
136 IF (PRESENT(sab_orb_external)) THEN
137 cpassert(ASSOCIATED(sab_orb_external))
138 sab_orb => sab_orb_external
139 END IF
140 cpassert(ASSOCIATED(sab_orb))
141
142 my_force_periodic = .false.
143 IF (PRESENT(force_periodic)) my_force_periodic = force_periodic
144 IF (my_force_periodic) THEN
145 perd0(:) = cell%perd(:)
146 cell%perd(:) = 1
147 END IF
148
149 nkind = SIZE(qs_kind_set)
150 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
151 ALLOCATE (cosab(ldwork, ldwork), sinab(ldwork, ldwork), work(ldwork, ldwork))
152
153 ALLOCATE (op_cossin(2, nkvec))
154
155 ALLOCATE (basis_set_list(nkind))
156 DO ikind = 1, nkind
157 qs_kind => qs_kind_set(ikind)
158 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
159 IF (ASSOCIATED(basis_set_a)) THEN
160 basis_set_list(ikind)%gto_basis_set => basis_set_a
161 ELSE
162 NULLIFY (basis_set_list(ikind)%gto_basis_set)
163 END IF
164 END DO
165
166 last_jatom = 0
167 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
168 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
169 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
170 iatom=iatom, jatom=jatom, r=rab)
171 basis_set_a => basis_set_list(ikind)%gto_basis_set
172 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
173 basis_set_b => basis_set_list(jkind)%gto_basis_set
174 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
175
176 ra(:) = pbc(particle_set(iatom)%r(:), cell)
177 rb(:) = ra(:) + rab(:)
178
179 first_sgfa => basis_set_a%first_sgf
180 la_max => basis_set_a%lmax
181 la_min => basis_set_a%lmin
182 npgfa => basis_set_a%npgf
183 nsgfa => basis_set_a%nsgf_set
184 rpgfa => basis_set_a%pgf_radius
185 set_radius_a => basis_set_a%set_radius
186 sphi_a => basis_set_a%sphi
187 zeta => basis_set_a%zet
188
189 first_sgfb => basis_set_b%first_sgf
190 lb_max => basis_set_b%lmax
191 lb_min => basis_set_b%lmin
192 npgfb => basis_set_b%npgf
193 nsgfb => basis_set_b%nsgf_set
194 rpgfb => basis_set_b%pgf_radius
195 set_radius_b => basis_set_b%set_radius
196 sphi_b => basis_set_b%sphi
197 zetb => basis_set_b%zet
198
199 nseta = basis_set_a%nset
200 nsetb = basis_set_b%nset
201 ldsa = SIZE(sphi_a, 1)
202 ldsb = SIZE(sphi_b, 1)
203 IF (inode == 1) last_jatom = 0
204
205 IF (jatom /= last_jatom) THEN
206 new_atom_b = .true.
207 last_jatom = jatom
208 ELSE
209 new_atom_b = .false.
210 END IF
211
212 IF (new_atom_b) THEN
213 IF (iatom <= jatom) THEN
214 irow = iatom
215 icol = jatom
216 ELSE
217 irow = jatom
218 icol = iatom
219 END IF
220
221 DO i = 1, nkvec
222 DO reim = 1, 2
223 NULLIFY (op_cossin(reim, i)%block)
224 CALL dbcsr_get_block_p(matrix=op_sm_set(reim, i)%matrix, &
225 row=irow, col=icol, block=op_cossin(reim, i)%block, found=found)
226 END DO
227 IF (ASSOCIATED(op_cossin(1, i)%block) .NEQV. ASSOCIATED(op_cossin(2, i)%block)) THEN
228 cpabort("cosine and sine blocks should have the same topology")
229 END IF
230 END DO
231 END IF
232
233 dab = norm2(rab)
234 DO iset = 1, nseta
235 ncoa = npgfa(iset)*ncoset(la_max(iset))
236 sgfa = first_sgfa(1, iset)
237 DO jset = 1, nsetb
238 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
239 ncob = npgfb(jset)*ncoset(lb_max(jset))
240 sgfb = first_sgfb(1, jset)
241
242 DO i = 1, nkvec
243 IF (.NOT. ASSOCIATED(op_cossin(1, i)%block)) cycle
244 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), &
245 la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), &
246 rpgfb(:, jset), lb_min(jset), ra, rb, kvec(:, i), cosab, sinab)
247 CALL contract_cossin(op_cossin(1, i)%block, op_cossin(2, i)%block, &
248 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
249 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
250 cosab, sinab, ldwork, work, ldwork)
251 END DO
252 END DO
253 END DO
254 END DO
255 CALL neighbor_list_iterator_release(nl_iterator)
256
257 DEALLOCATE (op_cossin, cosab, sinab, work, basis_set_list)
258
259 IF (my_force_periodic) cell%perd(:) = perd0(:)
260
261 CALL timestop(handle)
262 END SUBROUTINE build_exp_ikr_matrix
263
264! **************************************************************************************************
265!> \brief Calculation of the linear momentum matrix <mu|∂|nu> over
266!> Cartesian Gaussian functions.
267!> \param qs_env ...
268!> \param matrix ...
269!> \param minimum_image take into account only the first neighbors in the lists
270!> \date 27.02.2009
271!> \author VW
272!> \version 1.0
273! **************************************************************************************************
274 SUBROUTINE build_lin_mom_matrix(qs_env, matrix, minimum_image)
275
276 TYPE(qs_environment_type), POINTER :: qs_env
277 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
278 LOGICAL, INTENT(IN), OPTIONAL :: minimum_image
279
280 CHARACTER(len=*), PARAMETER :: routinen = 'build_lin_mom_matrix'
281
282 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
283 ldai, maxco, maxlgto, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
284 INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
285 nsgfb
286 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
287 LOGICAL :: do_symmetric, found, my_minimum_image, &
288 new_atom_b
289 REAL(kind=dp) :: alpha, dab, lxo2, lyo2, lzo2
290 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), TARGET :: rr_work
291 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
292 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
293 TARGET :: intab
294 REAL(kind=dp), DIMENSION(3) :: rab
295 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
296 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
297 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
298 TYPE(cell_type), POINTER :: cell
299 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
300 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
301 TYPE(neighbor_list_iterator_p_type), &
302 DIMENSION(:), POINTER :: nl_iterator
303 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
304 POINTER :: sab_nl
305 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
306 TYPE(qs_kind_type), POINTER :: qs_kind
307
308 CALL timeset(routinen, handle)
309
310 NULLIFY (cell, sab_nl, qs_kind_set)
311
312 CALL get_qs_env(qs_env=qs_env, &
313 qs_kind_set=qs_kind_set, &
314 cell=cell)
315
316 nkind = SIZE(qs_kind_set)
317
318 my_minimum_image = .false.
319 IF (PRESENT(minimum_image)) THEN
320 my_minimum_image = minimum_image
321 lxo2 = norm2(cell%hmat(:, 1))/2.0_dp
322 lyo2 = norm2(cell%hmat(:, 2))/2.0_dp
323 lzo2 = norm2(cell%hmat(:, 3))/2.0_dp
324 END IF
325
326 ! Take into account the symmetry of the input matrix
327 do_symmetric = dbcsr_has_symmetry(matrix(1)%matrix)
328 IF (do_symmetric) THEN
329 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
330 ELSE
331 CALL get_qs_env(qs_env=qs_env, sab_all=sab_nl)
332 END IF
333! *** Allocate work storage ***
334
335 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
336 maxco=maxco, &
337 maxlgto=maxlgto, &
338 maxsgf=maxsgf)
339
340 ldai = ncoset(maxlgto + 1)
341 CALL init_orbital_pointers(ldai)
342
343 ALLOCATE (rr_work(ldai*ldai*3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
344 rr_work(:) = 0.0_dp
345 intab(:, :, :) = 0.0_dp
346 work(:, :) = 0.0_dp
347
348 ALLOCATE (basis_set_list(nkind))
349 DO ikind = 1, nkind
350 qs_kind => qs_kind_set(ikind)
351 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
352 IF (ASSOCIATED(basis_set_a)) THEN
353 basis_set_list(ikind)%gto_basis_set => basis_set_a
354 ELSE
355 NULLIFY (basis_set_list(ikind)%gto_basis_set)
356 END IF
357 END DO
358 CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
359 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
360 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
361 iatom=iatom, jatom=jatom, r=rab)
362 basis_set_a => basis_set_list(ikind)%gto_basis_set
363 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
364 basis_set_b => basis_set_list(jkind)%gto_basis_set
365 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
366 ! basis ikind
367 first_sgfa => basis_set_a%first_sgf
368 la_max => basis_set_a%lmax
369 npgfa => basis_set_a%npgf
370 nseta = basis_set_a%nset
371 nsgfa => basis_set_a%nsgf_set
372 rpgfa => basis_set_a%pgf_radius
373 set_radius_a => basis_set_a%set_radius
374 sphi_a => basis_set_a%sphi
375 zeta => basis_set_a%zet
376 ! basis jkind
377 first_sgfb => basis_set_b%first_sgf
378 lb_max => basis_set_b%lmax
379 npgfb => basis_set_b%npgf
380 nsetb = basis_set_b%nset
381 nsgfb => basis_set_b%nsgf_set
382 rpgfb => basis_set_b%pgf_radius
383 set_radius_b => basis_set_b%set_radius
384 sphi_b => basis_set_b%sphi
385 zetb => basis_set_b%zet
386
387 IF (inode == 1) last_jatom = 0
388
389 IF (my_minimum_image) THEN
390 IF (any(abs(rab(:)) > [lxo2, lyo2, lzo2])) cycle
391 END IF
392
393 IF (jatom /= last_jatom) THEN
394 new_atom_b = .true.
395 last_jatom = jatom
396 ELSE
397 new_atom_b = .false.
398 END IF
399
400 IF (new_atom_b) THEN
401 alpha = 1.0_dp
402 IF (do_symmetric) THEN
403 IF (iatom <= jatom) THEN
404 irow = iatom
405 icol = jatom
406 ELSE
407 irow = jatom
408 icol = iatom
409 IF (dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_antisymmetric) THEN
410 alpha = -1.0_dp
411 END IF
412 END IF
413 ELSE
414 irow = iatom
415 icol = jatom
416 END IF
417
418 DO i = 1, 3
419 NULLIFY (integral(i)%block)
420 CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
421 row=irow, col=icol, block=integral(i)%block, found=found)
422 cpassert(ASSOCIATED(integral(i)%block))
423 END DO
424 END IF
425
426 dab = norm2(rab)
427
428 DO iset = 1, nseta
429
430 ncoa = npgfa(iset)*ncoset(la_max(iset))
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 ncob = npgfb(jset)*ncoset(lb_max(jset))
438 sgfb = first_sgfb(1, jset)
439
440 ! *** Calculate the primitive fermi contact integrals ***
441
442 ! Keep the full primitive layout used by the contraction below.
443 CALL overlap_ab(la_max(iset), 0, npgfa(iset), &
444 rpgfa(:, iset), zeta(:, iset), &
445 lb_max(jset), 0, npgfb(jset), &
446 rpgfb(:, jset), zetb(:, jset), &
447 rab, dab=intab, rr_work=rr_work)
448
449 ! *** Contraction step ***
450
451 DO i = 1, 3
452
453 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
454 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
455 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
456 0.0_dp, work(1, 1), SIZE(work, 1))
457
458 IF (do_symmetric) THEN
459 IF (iatom <= jatom) THEN
460
461 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
462 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
463 work(1, 1), SIZE(work, 1), &
464 1.0_dp, integral(i)%block(sgfa, sgfb), &
465 SIZE(integral(i)%block, 1))
466
467 ELSE
468
469 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
470 alpha, work(1, 1), SIZE(work, 1), &
471 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
472 1.0_dp, integral(i)%block(sgfb, sgfa), &
473 SIZE(integral(i)%block, 1))
474
475 END IF
476 ELSE
477 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
478 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
479 work(1, 1), SIZE(work, 1), &
480 1.0_dp, integral(i)%block(sgfa, sgfb), &
481 SIZE(integral(i)%block, 1))
482 END IF
483
484 END DO
485
486 END DO
487
488 END DO
489
490 END DO
491 CALL neighbor_list_iterator_release(nl_iterator)
492
493 ! *** Release work storage ***
494
495 DEALLOCATE (intab, rr_work, work, integral, basis_set_list)
496
497 CALL timestop(handle)
498
499 END SUBROUTINE build_lin_mom_matrix
500
501! **************************************************************************************************
502!> \brief Calculation of the angular momentum matrix over
503!> Cartesian Gaussian functions.
504!> \param qs_env ...
505!> \param matrix ...
506!> \param rc ...
507!> \date 27.02.2009
508!> \author VW
509!> \version 1.0
510! **************************************************************************************************
511
512 SUBROUTINE build_ang_mom_matrix(qs_env, matrix, rc)
513
514 TYPE(qs_environment_type), POINTER :: qs_env
515 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
516 REAL(dp), DIMENSION(:), INTENT(IN) :: rc
517
518 CHARACTER(len=*), PARAMETER :: routinen = 'build_ang_mom_matrix'
519
520 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
521 maxco, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
522 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
523 npgfb, nsgfa, nsgfb
524 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
525 LOGICAL :: found, new_atom_b
526 REAL(kind=dp) :: dab
527 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
528 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab
529 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rbc
530 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
531 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
532 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
533 TYPE(cell_type), POINTER :: cell
534 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
535 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
536 TYPE(neighbor_list_iterator_p_type), &
537 DIMENSION(:), POINTER :: nl_iterator
538 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
539 POINTER :: sab_all
540 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
541 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
542 TYPE(qs_kind_type), POINTER :: qs_kind
543
544 CALL timeset(routinen, handle)
545
546 CALL get_qs_env(qs_env=qs_env, &
547 qs_kind_set=qs_kind_set, &
548 particle_set=particle_set, &
549 sab_all=sab_all, &
550 cell=cell)
551
552 nkind = SIZE(qs_kind_set)
553
554! *** Allocate work storage ***
555
556 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
557 maxco=maxco, &
558 maxsgf=maxsgf)
559
560 ALLOCATE (intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
561 intab(:, :, :) = 0.0_dp
562 work(:, :) = 0.0_dp
563
564 ALLOCATE (basis_set_list(nkind))
565 DO ikind = 1, nkind
566 qs_kind => qs_kind_set(ikind)
567 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
568 IF (ASSOCIATED(basis_set_a)) THEN
569 basis_set_list(ikind)%gto_basis_set => basis_set_a
570 ELSE
571 NULLIFY (basis_set_list(ikind)%gto_basis_set)
572 END IF
573 END DO
574 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
575 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
576 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
577 iatom=iatom, jatom=jatom, r=rab)
578 basis_set_a => basis_set_list(ikind)%gto_basis_set
579 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
580 basis_set_b => basis_set_list(jkind)%gto_basis_set
581 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
582 ra = pbc(particle_set(iatom)%r, cell)
583 ! basis ikind
584 first_sgfa => basis_set_a%first_sgf
585 la_max => basis_set_a%lmax
586 la_min => basis_set_a%lmin
587 npgfa => basis_set_a%npgf
588 nseta = basis_set_a%nset
589 nsgfa => basis_set_a%nsgf_set
590 rpgfa => basis_set_a%pgf_radius
591 set_radius_a => basis_set_a%set_radius
592 sphi_a => basis_set_a%sphi
593 zeta => basis_set_a%zet
594 ! basis jkind
595 first_sgfb => basis_set_b%first_sgf
596 lb_max => basis_set_b%lmax
597 lb_min => basis_set_b%lmin
598 npgfb => basis_set_b%npgf
599 nsetb = basis_set_b%nset
600 nsgfb => basis_set_b%nsgf_set
601 rpgfb => basis_set_b%pgf_radius
602 set_radius_b => basis_set_b%set_radius
603 sphi_b => basis_set_b%sphi
604 zetb => basis_set_b%zet
605
606 IF (inode == 1) last_jatom = 0
607
608 IF (jatom /= last_jatom) THEN
609 new_atom_b = .true.
610 last_jatom = jatom
611 ELSE
612 new_atom_b = .false.
613 END IF
614
615 IF (new_atom_b) THEN
616 !IF (iatom <= jatom) THEN
617 irow = iatom
618 icol = jatom
619 !ELSE
620 ! irow = jatom
621 ! icol = iatom
622 !END IF
623
624 DO i = 1, 3
625 NULLIFY (integral(i)%block)
626 CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
627 row=irow, col=icol, block=integral(i)%block, found=found)
628 cpassert(found .AND. ASSOCIATED(integral(i)%block))
629 END DO
630 END IF
631
632 dab = norm2(rab)
633
634 DO iset = 1, nseta
635
636 ncoa = npgfa(iset)*ncoset(la_max(iset))
637 sgfa = first_sgfa(1, iset)
638
639 DO jset = 1, nsetb
640
641 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
642
643 !IF(PRESENT(wancen)) THEN
644 ! rc = wancen
645 rac = pbc(rc, ra, cell)
646 rbc = rac + rab
647 !ELSE
648 ! rc(1:3) = rb(1:3)
649 ! rac(1:3) = -rab(1:3)
650 ! rbc(1:3) = 0.0_dp
651 !ENDIF
652
653 ncob = npgfb(jset)*ncoset(lb_max(jset))
654 sgfb = first_sgfb(1, jset)
655
656 ! *** Calculate the primitive angular momentum integrals ***
657
658 CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
659 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), rac, rbc, intab)
660
661 ! *** Contraction step ***
662
663 DO i = 1, 3
664
665 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
666 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
667 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
668 0.0_dp, work(1, 1), SIZE(work, 1))
669
670 !IF (iatom <= jatom) THEN
671
672 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
673 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
674 work(1, 1), SIZE(work, 1), &
675 1.0_dp, integral(i)%block(sgfa, sgfb), &
676 SIZE(integral(i)%block, 1))
677
678 !ELSE
679 !
680 ! CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncoa,&
681 ! -1.0_dp,work(1,1),SIZE(work,1),&
682 ! sphi_a(1,sgfa),SIZE(sphi_a,1),&
683 ! 1.0_dp,integral(i)%block(sgfb,sgfa),&
684 ! SIZE(integral(i)%block,1))
685 !
686 !ENDIF
687
688 END DO
689
690 END DO
691
692 END DO
693
694 END DO
695 CALL neighbor_list_iterator_release(nl_iterator)
696
697 ! *** Release work storage ***
698
699 DEALLOCATE (intab, work, integral, basis_set_list)
700
701! *** Print the spin orbit matrix, if requested ***
702
703 !IF (BTEST(cp_print_key_should_output(logger%iter_info,&
704 ! qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM"),cp_p_file)) THEN
705 ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM",&
706 ! extension=".Log")
707 ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw)
708 ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw)
709 ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw)
710 ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,&
711 ! "DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM")
712 !END IF
713
714 CALL timestop(handle)
715
716 END SUBROUTINE build_ang_mom_matrix
717
718! **************************************************************************************************
719END MODULE qs_operators_ao
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cell_types::pbc
Definition cell_types.F:104
ai_angmom
Calculation of the angular momentum integrals over Cartesian Gaussian-type functions.
Definition ai_angmom.F:17
ai_angmom::angmom
subroutine, public angmom(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, rac, rbc, angab)
...
Definition ai_angmom.F:52
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition ai_moments.F:17
ai_moments::contract_cossin
subroutine, public contract_cossin(cos_block, sin_block, iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, cosab, sinab, ldab, work, ldwork)
...
Definition ai_moments.F:77
ai_moments::cossin
subroutine, public cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
...
Definition ai_moments.F:155
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab, rr_work)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:273
basis_set_types
Definition basis_set_types.F:15
block_p_types
collect pointers to a block of reals
Definition block_p_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_has_symmetry
logical function, public dbcsr_has_symmetry(matrix)
...
Definition cp_dbcsr_api.F:928
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:861
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition orbital_pointers.F:253
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:539
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:550
qs_operators_ao
Definition qs_operators_ao.F:13
qs_operators_ao::build_ang_mom_matrix
subroutine, public build_ang_mom_matrix(qs_env, matrix, rc)
Calculation of the angular momentum matrix over Cartesian Gaussian functions.
Definition qs_operators_ao.F:513
qs_operators_ao::build_lin_mom_matrix
subroutine, public build_lin_mom_matrix(qs_env, matrix, minimum_image)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
Definition qs_operators_ao.F:275
qs_operators_ao::build_exp_ikr_matrix
subroutine, public build_exp_ikr_matrix(qs_env, op_sm_set, kvec, sab_orb_external, basis_type, force_periodic, cell_external)
Build real and imaginary AO matrices for exp(i*k*r).
Definition qs_operators_ao.F:74
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
block_p_types::block_p_type
Definition block_p_types.F:23
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67