(git:31b3539)
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_expikr,&
16 expikr
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 COMPLEX(KIND=dp), ALLOCATABLE :: expab(:, :)
87 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
88 ldsa, ldsb, ldwork, ncoa, ncob, nkind, nkvec, nseta, nsetb, reim, sgfa, sgfb
89 INTEGER, DIMENSION(3) :: perd0
90 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
91 npgfb, nsgfa, nsgfb
92 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
93 LOGICAL :: found, my_force_periodic, new_atom_b
94 REAL(kind=dp) :: dab
95 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
96 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
97 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
98 zeta, zetb
99 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: op_cossin
100 TYPE(cell_type), POINTER :: cell
101 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
102 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
103 TYPE(neighbor_list_iterator_p_type), &
104 DIMENSION(:), POINTER :: nl_iterator
105 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
106 POINTER :: sab_orb
107 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
109 TYPE(qs_kind_type), POINTER :: qs_kind
110
111 CALL timeset(routinen, handle)
112
113 cpassert(SIZE(kvec, 1) == 3)
114 nkvec = SIZE(kvec, 2)
115 cpassert(nkvec > 0)
116 cpassert(SIZE(op_sm_set, 1) == 2)
117 cpassert(SIZE(op_sm_set, 2) == nkvec)
118 DO i = 1, nkvec
119 DO reim = 1, 2
120 cpassert(ASSOCIATED(op_sm_set(reim, i)%matrix))
121 cpassert(dbcsr_has_symmetry(op_sm_set(reim, i)%matrix))
122 CALL dbcsr_set(op_sm_set(reim, i)%matrix, 0.0_dp)
123 END DO
124 END DO
125
126 NULLIFY (qs_kind, qs_kind_set, particle_set, sab_orb, cell)
127 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
128 particle_set=particle_set, cell=cell, sab_orb=sab_orb)
129 cpassert(ASSOCIATED(cell))
130 cpassert(ASSOCIATED(qs_kind_set))
131 cpassert(ASSOCIATED(particle_set))
132 IF (PRESENT(cell_external)) THEN
133 cpassert(ASSOCIATED(cell_external))
134 cell => cell_external
135 END IF
136
137 IF (PRESENT(sab_orb_external)) THEN
138 cpassert(ASSOCIATED(sab_orb_external))
139 sab_orb => sab_orb_external
140 END IF
141 cpassert(ASSOCIATED(sab_orb))
142
143 my_force_periodic = .false.
144 IF (PRESENT(force_periodic)) my_force_periodic = force_periodic
145 perd0(:) = -1
146 IF (my_force_periodic) THEN
147 perd0(:) = cell%perd(:)
148 cell%perd(:) = 1
149 END IF
150
151 nkind = SIZE(qs_kind_set)
152 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
153 ALLOCATE (expab(ldwork, ldwork), work(ldwork, ldwork))
154
155 ALLOCATE (op_cossin(2, nkvec))
156
157 ALLOCATE (basis_set_list(nkind))
158 DO ikind = 1, nkind
159 qs_kind => qs_kind_set(ikind)
160 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
161 IF (ASSOCIATED(basis_set_a)) THEN
162 basis_set_list(ikind)%gto_basis_set => basis_set_a
163 ELSE
164 NULLIFY (basis_set_list(ikind)%gto_basis_set)
165 END IF
166 END DO
167
168 last_jatom = 0
169 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
170 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
171 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
172 iatom=iatom, jatom=jatom, r=rab)
173 basis_set_a => basis_set_list(ikind)%gto_basis_set
174 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
175 basis_set_b => basis_set_list(jkind)%gto_basis_set
176 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
177
178 ra(:) = pbc(particle_set(iatom)%r(:), cell)
179 rb(:) = ra(:) + rab(:)
180
181 first_sgfa => basis_set_a%first_sgf
182 la_max => basis_set_a%lmax
183 la_min => basis_set_a%lmin
184 npgfa => basis_set_a%npgf
185 nsgfa => basis_set_a%nsgf_set
186 rpgfa => basis_set_a%pgf_radius
187 set_radius_a => basis_set_a%set_radius
188 sphi_a => basis_set_a%sphi
189 zeta => basis_set_a%zet
190
191 first_sgfb => basis_set_b%first_sgf
192 lb_max => basis_set_b%lmax
193 lb_min => basis_set_b%lmin
194 npgfb => basis_set_b%npgf
195 nsgfb => basis_set_b%nsgf_set
196 rpgfb => basis_set_b%pgf_radius
197 set_radius_b => basis_set_b%set_radius
198 sphi_b => basis_set_b%sphi
199 zetb => basis_set_b%zet
200
201 nseta = basis_set_a%nset
202 nsetb = basis_set_b%nset
203 ldsa = SIZE(sphi_a, 1)
204 ldsb = SIZE(sphi_b, 1)
205 IF (inode == 1) last_jatom = 0
206
207 IF (jatom /= last_jatom) THEN
208 new_atom_b = .true.
209 last_jatom = jatom
210 ELSE
211 new_atom_b = .false.
212 END IF
213
214 IF (new_atom_b) THEN
215 IF (iatom <= jatom) THEN
216 irow = iatom
217 icol = jatom
218 ELSE
219 irow = jatom
220 icol = iatom
221 END IF
222
223 DO i = 1, nkvec
224 DO reim = 1, 2
225 NULLIFY (op_cossin(reim, i)%block)
226 CALL dbcsr_get_block_p(matrix=op_sm_set(reim, i)%matrix, &
227 row=irow, col=icol, block=op_cossin(reim, i)%block, found=found)
228 END DO
229 IF (ASSOCIATED(op_cossin(1, i)%block) .NEQV. ASSOCIATED(op_cossin(2, i)%block)) THEN
230 cpabort("cosine and sine blocks should have the same topology")
231 END IF
232 END DO
233 END IF
234
235 dab = norm2(rab)
236 DO iset = 1, nseta
237 ncoa = npgfa(iset)*ncoset(la_max(iset))
238 sgfa = first_sgfa(1, iset)
239 DO jset = 1, nsetb
240 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
241 ncob = npgfb(jset)*ncoset(lb_max(jset))
242 sgfb = first_sgfb(1, jset)
243
244 DO i = 1, nkvec
245 IF (.NOT. ASSOCIATED(op_cossin(1, i)%block)) cycle
246 CALL expikr(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), &
247 la_min(iset), lb_max(jset), npgfb(jset), zetb(:, jset), &
248 rpgfb(:, jset), lb_min(jset), ra, rb, kvec(:, i), expab)
249 CALL contract_expikr(op_cossin(1, i)%block, op_cossin(2, i)%block, &
250 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
251 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
252 expab, work, ldwork)
253 END DO
254 END DO
255 END DO
256 END DO
257 CALL neighbor_list_iterator_release(nl_iterator)
258
259 DEALLOCATE (op_cossin, expab, work, basis_set_list)
260
261 IF (my_force_periodic) cell%perd(:) = perd0(:)
262
263 CALL timestop(handle)
264 END SUBROUTINE build_exp_ikr_matrix
265
266! **************************************************************************************************
267!> \brief Calculation of the linear momentum matrix <mu|∂|nu> over
268!> Cartesian Gaussian functions.
269!> \param qs_env ...
270!> \param matrix ...
271!> \param minimum_image take into account only the first neighbors in the lists
272!> \date 27.02.2009
273!> \author VW
274!> \version 1.0
275! **************************************************************************************************
276 SUBROUTINE build_lin_mom_matrix(qs_env, matrix, minimum_image)
277
278 TYPE(qs_environment_type), POINTER :: qs_env
279 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
280 LOGICAL, INTENT(IN), OPTIONAL :: minimum_image
281
282 CHARACTER(len=*), PARAMETER :: routinen = 'build_lin_mom_matrix'
283
284 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
285 ldai, maxco, maxlgto, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
286 INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
287 nsgfb
288 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
289 LOGICAL :: do_symmetric, found, my_minimum_image, &
290 new_atom_b
291 REAL(kind=dp) :: alpha, dab, lxo2, lyo2, lzo2
292 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), TARGET :: rr_work
293 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
294 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
295 TARGET :: intab
296 REAL(kind=dp), DIMENSION(3) :: rab
297 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
298 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
299 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
300 TYPE(cell_type), POINTER :: cell
301 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
302 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
303 TYPE(neighbor_list_iterator_p_type), &
304 DIMENSION(:), POINTER :: nl_iterator
305 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
306 POINTER :: sab_nl
307 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
308 TYPE(qs_kind_type), POINTER :: qs_kind
309
310 CALL timeset(routinen, handle)
311
312 NULLIFY (cell, sab_nl, qs_kind_set)
313
314 CALL get_qs_env(qs_env=qs_env, &
315 qs_kind_set=qs_kind_set, &
316 cell=cell)
317
318 nkind = SIZE(qs_kind_set)
319
320 my_minimum_image = .false.
321 IF (PRESENT(minimum_image)) THEN
322 my_minimum_image = minimum_image
323 lxo2 = norm2(cell%hmat(:, 1))/2.0_dp
324 lyo2 = norm2(cell%hmat(:, 2))/2.0_dp
325 lzo2 = norm2(cell%hmat(:, 3))/2.0_dp
326 END IF
327
328 ! Take into account the symmetry of the input matrix
329 do_symmetric = dbcsr_has_symmetry(matrix(1)%matrix)
330 IF (do_symmetric) THEN
331 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
332 ELSE
333 CALL get_qs_env(qs_env=qs_env, sab_all=sab_nl)
334 END IF
335! *** Allocate work storage ***
336
337 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
338 maxco=maxco, &
339 maxlgto=maxlgto, &
340 maxsgf=maxsgf)
341
342 ldai = ncoset(maxlgto + 1)
343 CALL init_orbital_pointers(ldai)
344
345 ALLOCATE (rr_work(ldai*ldai*3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
346 rr_work(:) = 0.0_dp
347 intab(:, :, :) = 0.0_dp
348 work(:, :) = 0.0_dp
349
350 ALLOCATE (basis_set_list(nkind))
351 DO ikind = 1, nkind
352 qs_kind => qs_kind_set(ikind)
353 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
354 IF (ASSOCIATED(basis_set_a)) THEN
355 basis_set_list(ikind)%gto_basis_set => basis_set_a
356 ELSE
357 NULLIFY (basis_set_list(ikind)%gto_basis_set)
358 END IF
359 END DO
360 CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
361 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
362 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
363 iatom=iatom, jatom=jatom, r=rab)
364 basis_set_a => basis_set_list(ikind)%gto_basis_set
365 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
366 basis_set_b => basis_set_list(jkind)%gto_basis_set
367 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
368 ! basis ikind
369 first_sgfa => basis_set_a%first_sgf
370 la_max => basis_set_a%lmax
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 sphi_a => basis_set_a%sphi
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 npgfb => basis_set_b%npgf
382 nsetb = basis_set_b%nset
383 nsgfb => basis_set_b%nsgf_set
384 rpgfb => basis_set_b%pgf_radius
385 set_radius_b => basis_set_b%set_radius
386 sphi_b => basis_set_b%sphi
387 zetb => basis_set_b%zet
388
389 IF (inode == 1) last_jatom = 0
390
391 IF (my_minimum_image) THEN
392 IF (any(abs(rab(:)) > [lxo2, lyo2, lzo2])) cycle
393 END IF
394
395 IF (jatom /= last_jatom) THEN
396 new_atom_b = .true.
397 last_jatom = jatom
398 ELSE
399 new_atom_b = .false.
400 END IF
401
402 IF (new_atom_b) THEN
403 alpha = 1.0_dp
404 IF (do_symmetric) THEN
405 IF (iatom <= jatom) THEN
406 irow = iatom
407 icol = jatom
408 ELSE
409 irow = jatom
410 icol = iatom
411 IF (dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_antisymmetric) THEN
412 alpha = -1.0_dp
413 END IF
414 END IF
415 ELSE
416 irow = iatom
417 icol = jatom
418 END IF
419
420 DO i = 1, 3
421 NULLIFY (integral(i)%block)
422 CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
423 row=irow, col=icol, block=integral(i)%block, found=found)
424 cpassert(ASSOCIATED(integral(i)%block))
425 END DO
426 END IF
427
428 dab = norm2(rab)
429
430 DO iset = 1, nseta
431
432 ncoa = npgfa(iset)*ncoset(la_max(iset))
433 sgfa = first_sgfa(1, iset)
434
435 DO jset = 1, nsetb
436
437 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
438
439 ncob = npgfb(jset)*ncoset(lb_max(jset))
440 sgfb = first_sgfb(1, jset)
441
442 ! *** Calculate the primitive fermi contact integrals ***
443
444 ! Keep the full primitive layout used by the contraction below.
445 CALL overlap_ab(la_max(iset), 0, npgfa(iset), &
446 rpgfa(:, iset), zeta(:, iset), &
447 lb_max(jset), 0, npgfb(jset), &
448 rpgfb(:, jset), zetb(:, jset), &
449 rab, dab=intab, rr_work=rr_work)
450
451 ! *** Contraction step ***
452
453 DO i = 1, 3
454
455 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
456 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
457 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
458 0.0_dp, work(1, 1), SIZE(work, 1))
459
460 IF (do_symmetric) THEN
461 IF (iatom <= jatom) THEN
462
463 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
464 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
465 work(1, 1), SIZE(work, 1), &
466 1.0_dp, integral(i)%block(sgfa, sgfb), &
467 SIZE(integral(i)%block, 1))
468
469 ELSE
470
471 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
472 alpha, work(1, 1), SIZE(work, 1), &
473 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
474 1.0_dp, integral(i)%block(sgfb, sgfa), &
475 SIZE(integral(i)%block, 1))
476
477 END IF
478 ELSE
479 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
480 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
481 work(1, 1), SIZE(work, 1), &
482 1.0_dp, integral(i)%block(sgfa, sgfb), &
483 SIZE(integral(i)%block, 1))
484 END IF
485
486 END DO
487
488 END DO
489
490 END DO
491
492 END DO
493 CALL neighbor_list_iterator_release(nl_iterator)
494
495 ! *** Release work storage ***
496
497 DEALLOCATE (intab, rr_work, work, integral, basis_set_list)
498
499 CALL timestop(handle)
500
501 END SUBROUTINE build_lin_mom_matrix
502
503! **************************************************************************************************
504!> \brief Calculation of the angular momentum matrix over
505!> Cartesian Gaussian functions.
506!> \param qs_env ...
507!> \param matrix ...
508!> \param rc ...
509!> \date 27.02.2009
510!> \author VW
511!> \version 1.0
512! **************************************************************************************************
513
514 SUBROUTINE build_ang_mom_matrix(qs_env, matrix, rc)
515
516 TYPE(qs_environment_type), POINTER :: qs_env
517 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
518 REAL(dp), DIMENSION(:), INTENT(IN) :: rc
519
520 CHARACTER(len=*), PARAMETER :: routinen = 'build_ang_mom_matrix'
521
522 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
523 maxco, maxsgf, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
524 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
525 npgfb, nsgfa, nsgfb
526 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
527 LOGICAL :: found, new_atom_b
528 REAL(kind=dp) :: dab
529 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
530 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab
531 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rbc
532 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
533 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
534 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
535 TYPE(cell_type), POINTER :: cell
536 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
537 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
538 TYPE(neighbor_list_iterator_p_type), &
539 DIMENSION(:), POINTER :: nl_iterator
540 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
541 POINTER :: sab_all
542 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
543 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
544 TYPE(qs_kind_type), POINTER :: qs_kind
545
546 CALL timeset(routinen, handle)
547
548 CALL get_qs_env(qs_env=qs_env, &
549 qs_kind_set=qs_kind_set, &
550 particle_set=particle_set, &
551 sab_all=sab_all, &
552 cell=cell)
553
554 nkind = SIZE(qs_kind_set)
555
556! *** Allocate work storage ***
557
558 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
559 maxco=maxco, &
560 maxsgf=maxsgf)
561
562 ALLOCATE (intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
563 intab(:, :, :) = 0.0_dp
564 work(:, :) = 0.0_dp
565
566 ALLOCATE (basis_set_list(nkind))
567 DO ikind = 1, nkind
568 qs_kind => qs_kind_set(ikind)
569 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
570 IF (ASSOCIATED(basis_set_a)) THEN
571 basis_set_list(ikind)%gto_basis_set => basis_set_a
572 ELSE
573 NULLIFY (basis_set_list(ikind)%gto_basis_set)
574 END IF
575 END DO
576 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
577 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
578 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
579 iatom=iatom, jatom=jatom, r=rab)
580 basis_set_a => basis_set_list(ikind)%gto_basis_set
581 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
582 basis_set_b => basis_set_list(jkind)%gto_basis_set
583 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
584 ra = pbc(particle_set(iatom)%r, cell)
585 ! basis ikind
586 first_sgfa => basis_set_a%first_sgf
587 la_max => basis_set_a%lmax
588 la_min => basis_set_a%lmin
589 npgfa => basis_set_a%npgf
590 nseta = basis_set_a%nset
591 nsgfa => basis_set_a%nsgf_set
592 rpgfa => basis_set_a%pgf_radius
593 set_radius_a => basis_set_a%set_radius
594 sphi_a => basis_set_a%sphi
595 zeta => basis_set_a%zet
596 ! basis jkind
597 first_sgfb => basis_set_b%first_sgf
598 lb_max => basis_set_b%lmax
599 lb_min => basis_set_b%lmin
600 npgfb => basis_set_b%npgf
601 nsetb = basis_set_b%nset
602 nsgfb => basis_set_b%nsgf_set
603 rpgfb => basis_set_b%pgf_radius
604 set_radius_b => basis_set_b%set_radius
605 sphi_b => basis_set_b%sphi
606 zetb => basis_set_b%zet
607
608 IF (inode == 1) last_jatom = 0
609
610 IF (jatom /= last_jatom) THEN
611 new_atom_b = .true.
612 last_jatom = jatom
613 ELSE
614 new_atom_b = .false.
615 END IF
616
617 IF (new_atom_b) THEN
618 !IF (iatom <= jatom) THEN
619 irow = iatom
620 icol = jatom
621 !ELSE
622 ! irow = jatom
623 ! icol = iatom
624 !END IF
625
626 DO i = 1, 3
627 NULLIFY (integral(i)%block)
628 CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
629 row=irow, col=icol, block=integral(i)%block, found=found)
630 cpassert(found .AND. ASSOCIATED(integral(i)%block))
631 END DO
632 END IF
633
634 dab = norm2(rab)
635
636 DO iset = 1, nseta
637
638 ncoa = npgfa(iset)*ncoset(la_max(iset))
639 sgfa = first_sgfa(1, iset)
640
641 DO jset = 1, nsetb
642
643 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
644
645 !IF(PRESENT(wancen)) THEN
646 ! rc = wancen
647 rac = pbc(rc, ra, cell)
648 rbc = rac + rab
649 !ELSE
650 ! rc(1:3) = rb(1:3)
651 ! rac(1:3) = -rab(1:3)
652 ! rbc(1:3) = 0.0_dp
653 !ENDIF
654
655 ncob = npgfb(jset)*ncoset(lb_max(jset))
656 sgfb = first_sgfb(1, jset)
657
658 ! *** Calculate the primitive angular momentum integrals ***
659
660 CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
661 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), rac, rbc, intab)
662
663 ! *** Contraction step ***
664
665 DO i = 1, 3
666
667 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
668 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
669 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
670 0.0_dp, work(1, 1), SIZE(work, 1))
671
672 !IF (iatom <= jatom) THEN
673
674 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
675 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
676 work(1, 1), SIZE(work, 1), &
677 1.0_dp, integral(i)%block(sgfa, sgfb), &
678 SIZE(integral(i)%block, 1))
679
680 !ELSE
681 !
682 ! CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncoa,&
683 ! -1.0_dp,work(1,1),SIZE(work,1),&
684 ! sphi_a(1,sgfa),SIZE(sphi_a,1),&
685 ! 1.0_dp,integral(i)%block(sgfb,sgfa),&
686 ! SIZE(integral(i)%block,1))
687 !
688 !ENDIF
689
690 END DO
691
692 END DO
693
694 END DO
695
696 END DO
697 CALL neighbor_list_iterator_release(nl_iterator)
698
699 ! *** Release work storage ***
700
701 DEALLOCATE (intab, work, integral, basis_set_list)
702
703! *** Print the spin orbit matrix, if requested ***
704
705 !IF (BTEST(cp_print_key_should_output(logger%iter_info,&
706 ! qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM"),cp_p_file)) THEN
707 ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM",&
708 ! extension=".Log")
709 ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw)
710 ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw)
711 ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw)
712 ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,&
713 ! "DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM")
714 !END IF
715
716 CALL timestop(handle)
717
718 END SUBROUTINE build_ang_mom_matrix
719
720! **************************************************************************************************
721END 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:105
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)
Real angular-momentum kernel (r-C) cross grad, equal to i times L.
Definition ai_angmom.F:54
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition ai_moments.F:17
ai_moments::expikr
subroutine, public expikr(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, expab, dexpab, dexpba, momentab)
Evaluate primitive exp(i k.r) integrals through the complex overlap recurrence.
Definition ai_moments.F:163
ai_moments::contract_expikr
subroutine, public contract_expikr(cos_block, sin_block, iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, expab, work, ldwork, ordered)
Contract complex primitive integrals into real cosine and sine AO blocks.
Definition ai_moments.F:80
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:249
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:544
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:515
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:277
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:223
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