(git:ed6f26b)
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-2025 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_moments, ONLY: diff_momop,&
15 diffop,&
16 moment
17 USE ai_os_rr, ONLY: os_rr_ovlp
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_type_antisymmetric,&
28 dbcsr_type_no_symmetry
29 USE cp_log_handling, ONLY: cp_get_default_logger,&
30 cp_logger_type
31 USE kinds, ONLY: default_string_length,&
32 dp
33 USE mathconstants, ONLY: pi
34 USE message_passing, ONLY: mp_para_env_type
35 USE orbital_pointers, ONLY: coset,&
36 init_orbital_pointers,&
37 ncoset
38 USE particle_types, ONLY: particle_type
39 USE qs_environment_types, ONLY: get_qs_env,&
40 qs_environment_type
41 USE qs_kind_types, ONLY: get_qs_kind,&
42 get_qs_kind_set,&
43 qs_kind_type
44 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
45 neighbor_list_iterate,&
46 neighbor_list_iterator_create,&
47 neighbor_list_iterator_p_type,&
48 neighbor_list_iterator_release,&
49 neighbor_list_set_p_type
50#include "./base/base_uses.f90"
51
52 IMPLICIT NONE
53 PRIVATE
54
55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_operators_ao'
56
57! *** Public subroutines ***
58
59 PUBLIC :: p_xyz_ao, rrc_xyz_ao, rrc_xyz_der_ao
60 PUBLIC :: build_lin_mom_matrix, build_ang_mom_matrix
61
62CONTAINS
63
64! **************************************************************************************************
65!> \brief Calculation of the linear momentum matrix <mu|∂|nu> over
66!> Cartesian Gaussian functions.
67!> \param qs_env ...
68!> \param matrix ...
69!> \date 27.02.2009
70!> \author VW
71!> \version 1.0
72! **************************************************************************************************
73 SUBROUTINE build_lin_mom_matrix(qs_env, matrix)
74
75 TYPE(qs_environment_type), POINTER :: qs_env
76 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
77
78 CHARACTER(len=*), PARAMETER :: routinen = 'build_lin_mom_matrix'
79
80 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
81 ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, nseta, nsetb, &
82 sgfa, sgfb
83 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
84 npgfb, nsgfa, nsgfb
85 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
86 LOGICAL :: do_symmetric, found, new_atom_b
87 REAL(kind=dp) :: dab, rab2
88 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
89 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab, rr_work
90 REAL(kind=dp), DIMENSION(3) :: rab
91 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
92 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
93 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
94 TYPE(cell_type), POINTER :: cell
95 TYPE(cp_logger_type), POINTER :: logger
96 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
97 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
98 TYPE(mp_para_env_type), POINTER :: para_env
99 TYPE(neighbor_list_iterator_p_type), &
100 DIMENSION(:), POINTER :: nl_iterator
101 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
102 POINTER :: sab_nl
103 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
104 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
105 TYPE(qs_kind_type), POINTER :: qs_kind
106
107 CALL timeset(routinen, handle)
108
109 NULLIFY (cell, sab_nl, qs_kind_set, particle_set, para_env)
110 NULLIFY (logger)
111
112 logger => cp_get_default_logger()
113
114 CALL get_qs_env(qs_env=qs_env, &
115 qs_kind_set=qs_kind_set, &
116 particle_set=particle_set, &
117 neighbor_list_id=neighbor_list_id, &
118 para_env=para_env, &
119 cell=cell)
120
121 nkind = SIZE(qs_kind_set)
122 natom = SIZE(particle_set)
123
124 ! Take into account the symmetry of the input matrix
125 do_symmetric = dbcsr_has_symmetry(matrix(1)%matrix)
126 IF (do_symmetric) THEN
127 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_nl)
128 ELSE
129 CALL get_qs_env(qs_env=qs_env, sab_all=sab_nl)
130 END IF
131! *** Allocate work storage ***
132
133 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
134 maxco=maxco, &
135 maxlgto=maxlgto, &
136 maxsgf=maxsgf)
137
138 ldai = ncoset(maxlgto + 1)
139 CALL init_orbital_pointers(ldai)
140
141 ALLOCATE (rr_work(ldai, ldai, 3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
142 rr_work(:, :, :) = 0.0_dp
143 intab(:, :, :) = 0.0_dp
144 work(:, :) = 0.0_dp
145
146 ALLOCATE (basis_set_list(nkind))
147 DO ikind = 1, nkind
148 qs_kind => qs_kind_set(ikind)
149 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
150 IF (ASSOCIATED(basis_set_a)) THEN
151 basis_set_list(ikind)%gto_basis_set => basis_set_a
152 ELSE
153 NULLIFY (basis_set_list(ikind)%gto_basis_set)
154 END IF
155 END DO
156 CALL neighbor_list_iterator_create(nl_iterator, sab_nl)
157 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
158 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
159 iatom=iatom, jatom=jatom, r=rab)
160 basis_set_a => basis_set_list(ikind)%gto_basis_set
161 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
162 basis_set_b => basis_set_list(jkind)%gto_basis_set
163 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
164 ! basis ikind
165 first_sgfa => basis_set_a%first_sgf
166 la_max => basis_set_a%lmax
167 la_min => basis_set_a%lmin
168 npgfa => basis_set_a%npgf
169 nseta = basis_set_a%nset
170 nsgfa => basis_set_a%nsgf_set
171 rpgfa => basis_set_a%pgf_radius
172 set_radius_a => basis_set_a%set_radius
173 sphi_a => basis_set_a%sphi
174 zeta => basis_set_a%zet
175 ! basis jkind
176 first_sgfb => basis_set_b%first_sgf
177 lb_max => basis_set_b%lmax
178 lb_min => basis_set_b%lmin
179 npgfb => basis_set_b%npgf
180 nsetb = basis_set_b%nset
181 nsgfb => basis_set_b%nsgf_set
182 rpgfb => basis_set_b%pgf_radius
183 set_radius_b => basis_set_b%set_radius
184 sphi_b => basis_set_b%sphi
185 zetb => basis_set_b%zet
186
187 IF (inode == 1) last_jatom = 0
188 IF (jatom /= last_jatom) THEN
189 new_atom_b = .true.
190 last_jatom = jatom
191 ELSE
192 new_atom_b = .false.
193 END IF
194
195 IF (new_atom_b) THEN
196 IF (do_symmetric) THEN
197 IF (iatom <= jatom) THEN
198 irow = iatom
199 icol = jatom
200 ELSE
201 irow = jatom
202 icol = iatom
203 END IF
204 ELSE
205 irow = iatom
206 icol = jatom
207 END IF
208
209 DO i = 1, 3
210 NULLIFY (integral(i)%block)
211 CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
212 row=irow, col=icol, block=integral(i)%block, found=found)
213 cpassert(ASSOCIATED(integral(i)%block))
214 END DO
215 END IF
216
217 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
218 dab = sqrt(rab2)
219
220 DO iset = 1, nseta
221
222 ncoa = npgfa(iset)*ncoset(la_max(iset))
223 sgfa = first_sgfa(1, iset)
224
225 DO jset = 1, nsetb
226
227 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
228
229 ncob = npgfb(jset)*ncoset(lb_max(jset))
230 sgfb = first_sgfb(1, jset)
231
232 ! *** Calculate the primitive fermi contact integrals ***
233
234 CALL lin_mom(la_max(iset), la_min(iset), npgfa(iset), &
235 rpgfa(:, iset), zeta(:, iset), &
236 lb_max(jset), lb_min(jset), npgfb(jset), &
237 rpgfb(:, jset), zetb(:, jset), &
238 rab, intab, SIZE(rr_work, 1), rr_work)
239
240 ! *** Contraction step ***
241
242 DO i = 1, 3
243
244 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
245 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
246 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
247 0.0_dp, work(1, 1), SIZE(work, 1))
248
249 IF (do_symmetric) THEN
250 IF (iatom <= jatom) THEN
251
252 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
253 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
254 work(1, 1), SIZE(work, 1), &
255 1.0_dp, integral(i)%block(sgfa, sgfb), &
256 SIZE(integral(i)%block, 1))
257
258 ELSE
259
260 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
261 -1.0_dp, work(1, 1), SIZE(work, 1), &
262 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
263 1.0_dp, integral(i)%block(sgfb, sgfa), &
264 SIZE(integral(i)%block, 1))
265
266 END IF
267 ELSE
268 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
269 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
270 work(1, 1), SIZE(work, 1), &
271 1.0_dp, integral(i)%block(sgfa, sgfb), &
272 SIZE(integral(i)%block, 1))
273 END IF
274
275 END DO
276
277 END DO
278
279 END DO
280
281 END DO
282 CALL neighbor_list_iterator_release(nl_iterator)
283
284 ! *** Release work storage ***
285
286 DEALLOCATE (intab, work, integral, basis_set_list)
287
288! *** Print the spin orbit matrix, if requested ***
289
290 !IF (BTEST(cp_print_key_should_output(logger%iter_info,&
291 ! qs_env%input,"DFT%PRINT%AO_MATRICES/LINEAR_MOMENTUM"),cp_p_file)) THEN
292 ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/LINEA_MOMENTUM",&
293 ! extension=".Log")
294 ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw)
295 ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw)
296 ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw)
297 ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,&
298 ! "DFT%PRINT%AO_MATRICES/LINEAR_MOMENTUM")
299 !END IF
300
301 CALL timestop(handle)
302
303 END SUBROUTINE build_lin_mom_matrix
304
305! **************************************************************************************************
306!> \brief Calculation of the primitive paramagnetic spin orbit integrals over
307!> Cartesian Gaussian-type functions.
308!> \param la_max ...
309!> \param la_min ...
310!> \param npgfa ...
311!> \param rpgfa ...
312!> \param zeta ...
313!> \param lb_max ...
314!> \param lb_min ...
315!> \param npgfb ...
316!> \param rpgfb ...
317!> \param zetb ...
318!> \param rab ...
319!> \param intab ...
320!> \param ldrr ...
321!> \param rr ...
322!> \date 02.03.2009
323!> \author VW
324!> \version 1.0
325! **************************************************************************************************
326 SUBROUTINE lin_mom(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, &
327 rab, intab, ldrr, rr)
328 INTEGER, INTENT(IN) :: la_max, la_min, npgfa
329 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
330 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
331 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
332 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
333 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: intab
334 INTEGER, INTENT(IN) :: ldrr
335 REAL(dp), DIMENSION(0:ldrr-1, 0:ldrr-1, 3), &
336 INTENT(INOUT) :: rr
337
338 INTEGER :: ax, ay, az, bx, by, bz, coa, cob, i, &
339 ipgf, j, jpgf, la, lb, ma, mb, na, nb
340 REAL(dp) :: dab, dumx, dumy, dumz, f0, rab2, xhi, zet
341 REAL(dp), DIMENSION(3) :: rap, rbp
342
343! *** Calculate the distance of the centers a and c ***
344
345 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
346 dab = sqrt(rab2)
347
348 ! *** Loop over all pairs of primitive Gaussian-type functions ***
349
350 na = 0
351
352 DO ipgf = 1, npgfa
353
354 nb = 0
355
356 DO jpgf = 1, npgfb
357
358 ! *** Screening ***
359
360 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
361 DO j = nb + 1, nb + ncoset(lb_max)
362 DO i = na + 1, na + ncoset(la_max)
363 intab(i, j, 1) = 0.0_dp
364 intab(i, j, 2) = 0.0_dp
365 intab(i, j, 3) = 0.0_dp
366 END DO
367 END DO
368 nb = nb + ncoset(lb_max)
369 cycle
370 END IF
371
372 ! *** Calculate some prefactors ***
373 zet = zeta(ipgf) + zetb(jpgf)
374 xhi = zeta(ipgf)*zetb(jpgf)/zet
375 rap = zetb(jpgf)*rab/zet
376 rbp = -zeta(ipgf)*rab/zet
377
378 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
379
380 ! *** Calculate the recurrence relation ***
381
382 CALL os_rr_ovlp(rap, la_max + 1, rbp, lb_max, zet, ldrr, rr)
383
384 ! *** Calculate the primitive linear momentum integrals ***
385 DO lb = lb_min, lb_max
386 DO bx = 0, lb
387 DO by = 0, lb - bx
388 bz = lb - bx - by
389 cob = coset(bx, by, bz)
390 mb = nb + cob
391 DO la = la_min, la_max
392 DO ax = 0, la
393 DO ay = 0, la - ax
394 az = la - ax - ay
395 coa = coset(ax, ay, az)
396 ma = na + coa
397 !
398 !
399 ! (a|p_x|b) = 2*a*(a+1x|b) - N_x(a)*(a-1_x|b)
400 dumx = 2.0_dp*zeta(ipgf)*rr(ax + 1, bx, 1)
401 IF (ax .GT. 0) dumx = dumx - real(ax, dp)*rr(ax - 1, bx, 1)
402 intab(ma, mb, 1) = f0*dumx*rr(ay, by, 2)*rr(az, bz, 3)
403 !
404 ! (a|p_y|b)
405 dumy = 2.0_dp*zeta(ipgf)*rr(ay + 1, by, 2)
406 IF (ay .GT. 0) dumy = dumy - real(ay, dp)*rr(ay - 1, by, 2)
407 intab(ma, mb, 2) = f0*rr(ax, bx, 1)*dumy*rr(az, bz, 3)
408 !
409 ! (a|p_z|b)
410 dumz = 2.0_dp*zeta(ipgf)*rr(az + 1, bz, 3)
411 IF (az .GT. 0) dumz = dumz - real(az, dp)*rr(az - 1, bz, 3)
412 intab(ma, mb, 3) = f0*rr(ax, bx, 1)*rr(ay, by, 2)*dumz
413 !
414 END DO
415 END DO
416 END DO !la
417
418 END DO
419 END DO
420 END DO !lb
421
422 nb = nb + ncoset(lb_max)
423
424 END DO
425
426 na = na + ncoset(la_max)
427
428 END DO
429
430 END SUBROUTINE lin_mom
431
432! **************************************************************************************************
433!> \brief Calculation of the angular momentum matrix over
434!> Cartesian Gaussian functions.
435!> \param qs_env ...
436!> \param matrix ...
437!> \param rc ...
438!> \date 27.02.2009
439!> \author VW
440!> \version 1.0
441! **************************************************************************************************
442
443 SUBROUTINE build_ang_mom_matrix(qs_env, matrix, rc)
444
445 TYPE(qs_environment_type), POINTER :: qs_env
446 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix
447 REAL(dp), DIMENSION(:), INTENT(IN) :: rc
448
449 CHARACTER(len=*), PARAMETER :: routinen = 'build_ang_mom_matrix'
450
451 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
452 ldai, maxco, maxlgto, maxsgf, natom, ncoa, ncob, neighbor_list_id, nkind, nseta, nsetb, &
453 sgfa, sgfb
454 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
455 npgfb, nsgfa, nsgfb
456 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
457 LOGICAL :: found, new_atom_b
458 REAL(kind=dp) :: dab, rab2
459 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
460 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: intab, rr_work
461 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rbc
462 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
463 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
464 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: integral
465 TYPE(cell_type), POINTER :: cell
466 TYPE(cp_logger_type), POINTER :: logger
467 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
468 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
469 TYPE(mp_para_env_type), POINTER :: para_env
470 TYPE(neighbor_list_iterator_p_type), &
471 DIMENSION(:), POINTER :: nl_iterator
472 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
473 POINTER :: sab_all
474 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
475 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
476 TYPE(qs_kind_type), POINTER :: qs_kind
477
478 CALL timeset(routinen, handle)
479
480 NULLIFY (cell, sab_all, qs_kind_set, particle_set, para_env)
481 NULLIFY (logger)
482
483 logger => cp_get_default_logger()
484
485 CALL get_qs_env(qs_env=qs_env, &
486 qs_kind_set=qs_kind_set, &
487 particle_set=particle_set, &
488 neighbor_list_id=neighbor_list_id, &
489 para_env=para_env, &
490 sab_all=sab_all, &
491 cell=cell)
492
493 nkind = SIZE(qs_kind_set)
494 natom = SIZE(particle_set)
495
496! *** Allocate work storage ***
497
498 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
499 maxco=maxco, &
500 maxlgto=maxlgto, &
501 maxsgf=maxsgf)
502
503 ldai = ncoset(maxlgto + 1)
504 CALL init_orbital_pointers(ldai)
505
506 ALLOCATE (rr_work(ldai, ldai, 3), intab(maxco, maxco, 3), work(maxco, maxsgf), integral(3))
507 rr_work(:, :, :) = 0.0_dp
508 intab(:, :, :) = 0.0_dp
509 work(:, :) = 0.0_dp
510
511 ALLOCATE (basis_set_list(nkind))
512 DO ikind = 1, nkind
513 qs_kind => qs_kind_set(ikind)
514 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
515 IF (ASSOCIATED(basis_set_a)) THEN
516 basis_set_list(ikind)%gto_basis_set => basis_set_a
517 ELSE
518 NULLIFY (basis_set_list(ikind)%gto_basis_set)
519 END IF
520 END DO
521 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
522 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
523 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
524 iatom=iatom, jatom=jatom, r=rab)
525 basis_set_a => basis_set_list(ikind)%gto_basis_set
526 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
527 basis_set_b => basis_set_list(jkind)%gto_basis_set
528 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
529 ra = pbc(particle_set(iatom)%r, cell)
530 ! basis ikind
531 first_sgfa => basis_set_a%first_sgf
532 la_max => basis_set_a%lmax
533 la_min => basis_set_a%lmin
534 npgfa => basis_set_a%npgf
535 nseta = basis_set_a%nset
536 nsgfa => basis_set_a%nsgf_set
537 rpgfa => basis_set_a%pgf_radius
538 set_radius_a => basis_set_a%set_radius
539 sphi_a => basis_set_a%sphi
540 zeta => basis_set_a%zet
541 ! basis jkind
542 first_sgfb => basis_set_b%first_sgf
543 lb_max => basis_set_b%lmax
544 lb_min => basis_set_b%lmin
545 npgfb => basis_set_b%npgf
546 nsetb = basis_set_b%nset
547 nsgfb => basis_set_b%nsgf_set
548 rpgfb => basis_set_b%pgf_radius
549 set_radius_b => basis_set_b%set_radius
550 sphi_b => basis_set_b%sphi
551 zetb => basis_set_b%zet
552
553 IF (inode == 1) last_jatom = 0
554
555 IF (jatom /= last_jatom) THEN
556 new_atom_b = .true.
557 last_jatom = jatom
558 ELSE
559 new_atom_b = .false.
560 END IF
561
562 IF (new_atom_b) THEN
563 !IF (iatom <= jatom) THEN
564 irow = iatom
565 icol = jatom
566 !ELSE
567 ! irow = jatom
568 ! icol = iatom
569 !END IF
570
571 DO i = 1, 3
572 NULLIFY (integral(i)%block)
573 CALL dbcsr_get_block_p(matrix=matrix(i)%matrix, &
574 row=irow, col=icol, block=integral(i)%block, found=found)
575 cpassert(ASSOCIATED(integral(i)%block))
576 END DO
577 END IF
578
579 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
580 dab = sqrt(rab2)
581
582 DO iset = 1, nseta
583
584 ncoa = npgfa(iset)*ncoset(la_max(iset))
585 sgfa = first_sgfa(1, iset)
586
587 DO jset = 1, nsetb
588
589 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
590
591 !IF(PRESENT(wancen)) THEN
592 ! rc = wancen
593 rac = pbc(rc, ra, cell)
594 rbc = rac + rab
595 !ELSE
596 ! rc(1:3) = rb(1:3)
597 ! rac(1:3) = -rab(1:3)
598 ! rbc(1:3) = 0.0_dp
599 !ENDIF
600
601 ncob = npgfb(jset)*ncoset(lb_max(jset))
602 sgfb = first_sgfb(1, jset)
603
604 ! *** Calculate the primitive angular momentum integrals ***
605
606 CALL ang_mom(la_max(iset), la_min(iset), npgfa(iset), &
607 rpgfa(:, iset), zeta(:, iset), &
608 lb_max(jset), lb_min(jset), npgfb(jset), &
609 rpgfb(:, jset), zetb(:, jset), &
610 rab, rac, intab, SIZE(rr_work, 1), rr_work)
611
612 ! *** Contraction step ***
613
614 DO i = 1, 3
615
616 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
617 1.0_dp, intab(1, 1, i), SIZE(intab, 1), &
618 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
619 0.0_dp, work(1, 1), SIZE(work, 1))
620
621 !IF (iatom <= jatom) THEN
622
623 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
624 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
625 work(1, 1), SIZE(work, 1), &
626 1.0_dp, integral(i)%block(sgfa, sgfb), &
627 SIZE(integral(i)%block, 1))
628
629 !ELSE
630 !
631 ! CALL dgemm("T","N",nsgfb(jset),nsgfa(iset),ncoa,&
632 ! -1.0_dp,work(1,1),SIZE(work,1),&
633 ! sphi_a(1,sgfa),SIZE(sphi_a,1),&
634 ! 1.0_dp,integral(i)%block(sgfb,sgfa),&
635 ! SIZE(integral(i)%block,1))
636 !
637 !ENDIF
638
639 END DO
640
641 END DO
642
643 END DO
644
645 END DO
646 CALL neighbor_list_iterator_release(nl_iterator)
647
648 ! *** Release work storage ***
649
650 DEALLOCATE (intab, work, integral, basis_set_list)
651
652! *** Print the spin orbit matrix, if requested ***
653
654 !IF (BTEST(cp_print_key_should_output(logger%iter_info,&
655 ! qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM"),cp_p_file)) THEN
656 ! iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM",&
657 ! extension=".Log")
658 ! CALL cp_dbcsr_write_sparse_matrix(matrix(1)%matrix,4,6,qs_env,para_env,output_unit=iw)
659 ! CALL cp_dbcsr_write_sparse_matrix(matrix(2)%matrix,4,6,qs_env,para_env,output_unit=iw)
660 ! CALL cp_dbcsr_write_sparse_matrix(matrix(3)%matrix,4,6,qs_env,para_env,output_unit=iw)
661 ! CALL cp_print_key_finished_output(iw,logger,qs_env%input,&
662 ! "DFT%PRINT%AO_MATRICES/ANGULAR_MOMENTUM")
663 !END IF
664
665 CALL timestop(handle)
666
667 END SUBROUTINE build_ang_mom_matrix
668
669! **************************************************************************************************
670!> \brief Calculation of the primitive paramagnetic spin orbit integrals over
671!> Cartesian Gaussian-type functions.
672!> \param la_max ...
673!> \param la_min ...
674!> \param npgfa ...
675!> \param rpgfa ...
676!> \param zeta ...
677!> \param lb_max ...
678!> \param lb_min ...
679!> \param npgfb ...
680!> \param rpgfb ...
681!> \param zetb ...
682!> \param rab ...
683!> \param rac ...
684!> \param intab ...
685!> \param ldrr ...
686!> \param rr ...
687!> \date 02.03.2009
688!> \author VW
689!> \version 1.0
690! **************************************************************************************************
691 SUBROUTINE ang_mom(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, &
692 rab, rac, intab, ldrr, rr)
693 INTEGER, INTENT(IN) :: la_max, la_min, npgfa
694 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
695 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
696 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
697 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab, rac
698 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: intab
699 INTEGER, INTENT(IN) :: ldrr
700 REAL(dp), DIMENSION(0:ldrr-1, 0:ldrr-1, 3), &
701 INTENT(INOUT) :: rr
702
703 INTEGER :: ax, ay, az, bx, by, bz, coa, cob, i, &
704 ipgf, j, jpgf, la, lb, ma, mb, na, nb
705 REAL(dp) :: dab, dumx, dumy, dumz, f0, rab2, xhi, zet
706 REAL(dp), DIMENSION(3) :: rap, rbp
707
708! *** Calculate the distance of the centers a and c ***
709
710 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
711 dab = sqrt(rab2)
712
713 ! *** Loop over all pairs of primitive Gaussian-type functions ***
714
715 na = 0
716
717 DO ipgf = 1, npgfa
718
719 nb = 0
720
721 DO jpgf = 1, npgfb
722
723 ! *** Screening ***
724
725 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
726 DO j = nb + 1, nb + ncoset(lb_max)
727 DO i = na + 1, na + ncoset(la_max)
728 intab(i, j, 1) = 0.0_dp
729 intab(i, j, 2) = 0.0_dp
730 intab(i, j, 3) = 0.0_dp
731 END DO
732 END DO
733 nb = nb + ncoset(lb_max)
734 cycle
735 END IF
736
737 ! *** Calculate some prefactors ***
738 zet = zeta(ipgf) + zetb(jpgf)
739 xhi = zeta(ipgf)*zetb(jpgf)/zet
740 rap = zetb(jpgf)*rab/zet
741 rbp = -zeta(ipgf)*rab/zet
742
743 f0 = (pi/zet)**(1.5_dp)*exp(-xhi*rab2)
744
745 ! *** Calculate the recurrence relation ***
746
747 CALL os_rr_ovlp(rap, la_max + 1, rbp, lb_max, zet, ldrr, rr)
748
749 ! *** Calculate the primitive Fermi contact integrals ***
750
751 DO lb = lb_min, lb_max
752 DO bx = 0, lb
753 DO by = 0, lb - bx
754 bz = lb - bx - by
755 cob = coset(bx, by, bz)
756 mb = nb + cob
757 DO la = la_min, la_max
758 DO ax = 0, la
759 DO ay = 0, la - ax
760 az = la - ax - ay
761 coa = coset(ax, ay, az)
762 ma = na + coa
763 !
764 dumx = -2.0_dp*zeta(ipgf)*rr(ax + 1, bx, 1)
765 dumy = -2.0_dp*zeta(ipgf)*rr(ay + 1, by, 2)
766 dumz = -2.0_dp*zeta(ipgf)*rr(az + 1, bz, 3)
767 IF (ax .GT. 0) dumx = dumx + real(ax, dp)*rr(ax - 1, bx, 1)
768 IF (ay .GT. 0) dumy = dumy + real(ay, dp)*rr(ay - 1, by, 2)
769 IF (az .GT. 0) dumz = dumz + real(az, dp)*rr(az - 1, bz, 3)
770 !
771 ! (a|l_z|b)
772 intab(ma, mb, 1) = -f0*rr(ax, bx, 1)*( &
773 & (rr(ay + 1, by, 2) + rac(2)*rr(ay, by, 2))*dumz &
774 & - (rr(az + 1, bz, 3) + rac(3)*rr(az, bz, 3))*dumy)
775 !
776 ! (a|l_y|b)
777 intab(ma, mb, 2) = -f0*rr(ay, by, 2)*( &
778 & (rr(az + 1, bz, 3) + rac(3)*rr(az, bz, 3))*dumx &
779 & - (rr(ax + 1, bx, 1) + rac(1)*rr(ax, bx, 1))*dumz)
780 !
781 ! (a|l_z|b)
782 intab(ma, mb, 3) = -f0*rr(az, bz, 3)*( &
783 & (rr(ax + 1, bx, 1) + rac(1)*rr(ax, bx, 1))*dumy &
784 & - (rr(ay + 1, by, 2) + rac(2)*rr(ay, by, 2))*dumx)
785 !
786 END DO
787 END DO
788 END DO !la
789
790 END DO
791 END DO
792 END DO !lb
793
794 nb = nb + ncoset(lb_max)
795
796 END DO
797
798 na = na + ncoset(la_max)
799
800 END DO
801
802 END SUBROUTINE ang_mom
803
804! **************************************************************************************************
805!> \brief Calculation of the components of the dipole operator in the velocity form
806!> The elements of the sparse matrices are the integrals in the
807!> basis functions
808!> \param op matrix representation of the p operator
809!> calculated in terms of the contracted basis functions
810!> \param qs_env environment for the lists and the basis sets
811!> \param minimum_image take into account only the first neighbors in the lists
812!> \par History
813!> 06.2005 created [MI]
814!> \author MI
815! **************************************************************************************************
816
817 SUBROUTINE p_xyz_ao(op, qs_env, minimum_image)
818
819 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op
820 TYPE(qs_environment_type), POINTER :: qs_env
821 LOGICAL, INTENT(IN), OPTIONAL :: minimum_image
822
823 CHARACTER(len=*), PARAMETER :: routinen = 'p_xyz_ao'
824
825 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
826 ldab, ldsa, ldsb, ldwork, maxl, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
827 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
828 npgfb, nsgfa, nsgfb
829 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
830 LOGICAL :: found, my_minimum_image, new_atom_b
831 REAL(kind=dp) :: alpha, dab, lxo2, lyo2, lzo2, rab2
832 REAL(kind=dp), DIMENSION(3) :: ra, rab
833 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
834 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
835 zeta, zetb
836 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: difab
837 TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip
838 TYPE(cell_type), POINTER :: cell
839 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
840 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
841 TYPE(neighbor_list_iterator_p_type), &
842 DIMENSION(:), POINTER :: nl_iterator
843 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
844 POINTER :: sab_orb
845 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
846 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
847 TYPE(qs_kind_type), POINTER :: qs_kind
848
849 CALL timeset(routinen, handle)
850
851 NULLIFY (qs_kind, qs_kind_set)
852 NULLIFY (cell, particle_set)
853 NULLIFY (sab_orb)
854 NULLIFY (difab, op_dip, work)
855 NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb)
856 NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
857
858 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
859 cell=cell, particle_set=particle_set, &
860 sab_orb=sab_orb)
861
862 nkind = SIZE(qs_kind_set)
863
864 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
865 maxco=ldwork, maxlgto=maxl)
866
867 my_minimum_image = .false.
868 IF (PRESENT(minimum_image)) THEN
869 my_minimum_image = minimum_image
870 lxo2 = sqrt(sum(cell%hmat(:, 1)**2))/2.0_dp
871 lyo2 = sqrt(sum(cell%hmat(:, 2)**2))/2.0_dp
872 lzo2 = sqrt(sum(cell%hmat(:, 3)**2))/2.0_dp
873 END IF
874
875 ldab = ldwork
876
877 ALLOCATE (difab(ldab, ldab, 3))
878 difab(1:ldab, 1:ldab, 1:3) = 0.0_dp
879 ALLOCATE (work(ldwork, ldwork))
880 work(1:ldwork, 1:ldwork) = 0.0_dp
881 ALLOCATE (op_dip(3))
882
883 DO i = 1, 3
884 NULLIFY (op_dip(i)%block)
885 END DO
886
887 ALLOCATE (basis_set_list(nkind))
888 DO ikind = 1, nkind
889 qs_kind => qs_kind_set(ikind)
890 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
891 IF (ASSOCIATED(basis_set_a)) THEN
892 basis_set_list(ikind)%gto_basis_set => basis_set_a
893 ELSE
894 NULLIFY (basis_set_list(ikind)%gto_basis_set)
895 END IF
896 END DO
897 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
898 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
899 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
900 iatom=iatom, jatom=jatom, r=rab)
901 basis_set_a => basis_set_list(ikind)%gto_basis_set
902 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
903 basis_set_b => basis_set_list(jkind)%gto_basis_set
904 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
905 ra = pbc(particle_set(iatom)%r, cell)
906 ! basis ikind
907 first_sgfa => basis_set_a%first_sgf
908 la_max => basis_set_a%lmax
909 la_min => basis_set_a%lmin
910 npgfa => basis_set_a%npgf
911 nseta = basis_set_a%nset
912 nsgfa => basis_set_a%nsgf_set
913 rpgfa => basis_set_a%pgf_radius
914 set_radius_a => basis_set_a%set_radius
915 sphi_a => basis_set_a%sphi
916 zeta => basis_set_a%zet
917 ! basis jkind
918 first_sgfb => basis_set_b%first_sgf
919 lb_max => basis_set_b%lmax
920 lb_min => basis_set_b%lmin
921 npgfb => basis_set_b%npgf
922 nsetb = basis_set_b%nset
923 nsgfb => basis_set_b%nsgf_set
924 rpgfb => basis_set_b%pgf_radius
925 set_radius_b => basis_set_b%set_radius
926 sphi_b => basis_set_b%sphi
927 zetb => basis_set_b%zet
928
929 IF (inode == 1) THEN
930 last_jatom = 0
931 alpha = 1.0_dp
932 END IF
933 ldsa = SIZE(sphi_a, 1)
934 ldsb = SIZE(sphi_b, 1)
935
936 IF (my_minimum_image) THEN
937 IF (abs(rab(1)) > lxo2 .OR. abs(rab(2)) > lyo2 .OR. abs(rab(3)) > lzo2) cycle
938 END IF
939
940 IF (jatom /= last_jatom) THEN
941 new_atom_b = .true.
942 last_jatom = jatom
943 ELSE
944 new_atom_b = .false.
945 END IF
946
947 IF (new_atom_b) THEN
948 IF (iatom <= jatom) THEN
949 irow = iatom
950 icol = jatom
951 alpha = 1.0_dp
952 ELSE
953 irow = jatom
954 icol = iatom
955 IF (dbcsr_get_matrix_type(op(1)%matrix) == dbcsr_type_antisymmetric) THEN
956 !IF(op(1)%matrix%symmetry=="antisymmetric") THEN
957 alpha = -1.0_dp
958 END IF
959 END IF
960
961 DO i = 1, 3
962 NULLIFY (op_dip(i)%block)
963 CALL dbcsr_get_block_p(matrix=op(i)%matrix, &
964 row=irow, col=icol, block=op_dip(i)%block, found=found)
965 cpassert(ASSOCIATED(op_dip(i)%block))
966 END DO
967 END IF ! new_atom_b
968 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
969 dab = sqrt(rab2)
970
971 DO iset = 1, nseta
972
973 ncoa = npgfa(iset)*ncoset(la_max(iset))
974 sgfa = first_sgfa(1, iset)
975
976 DO jset = 1, nsetb
977
978 ncob = npgfb(jset)*ncoset(lb_max(jset))
979 sgfb = first_sgfb(1, jset)
980
981 IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
982
983! *** Calculate the primitive overlap integrals ***
984 CALL diffop(la_max(iset), npgfa(iset), zeta(:, iset), &
985 rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
986 zetb(:, jset), rpgfb(:, jset), lb_min(jset), rab, difab)
987
988! *** Contraction ***
989 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
990 alpha, difab(1, 1, 1), ldab, sphi_b(1, sgfb), ldsb, &
991 0.0_dp, work(1, 1), ldwork)
992 IF (iatom <= jatom) THEN
993 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
994 1.0_dp, sphi_a(1, sgfa), ldsa, &
995 work(1, 1), ldwork, &
996 1.0_dp, op_dip(1)%block(sgfa, sgfb), &
997 SIZE(op_dip(1)%block, 1))
998
999 ELSE
1000 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1001 1.0_dp, work(1, 1), ldwork, &
1002 sphi_a(1, sgfa), ldsa, &
1003 1.0_dp, op_dip(1)%block(sgfb, sgfa), &
1004 SIZE(op_dip(1)%block, 1))
1005
1006 END IF
1007
1008! *** Contraction ***
1009 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1010 alpha, difab(1, 1, 2), ldab, sphi_b(1, sgfb), ldsb, &
1011 0.0_dp, work(1, 1), ldwork)
1012 IF (iatom <= jatom) THEN
1013 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1014 1.0_dp, sphi_a(1, sgfa), ldsa, &
1015 work(1, 1), ldwork, &
1016 1.0_dp, op_dip(2)%block(sgfa, sgfb), &
1017 SIZE(op_dip(2)%block, 1))
1018 ELSE
1019 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1020 1.0_dp, work(1, 1), ldwork, &
1021 sphi_a(1, sgfa), ldsa, &
1022 1.0_dp, op_dip(2)%block(sgfb, sgfa), &
1023 SIZE(op_dip(2)%block, 1))
1024 END IF
1025
1026! *** Contraction ***
1027 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1028 alpha, difab(1, 1, 3), ldab, sphi_b(1, sgfb), ldsb, &
1029 0.0_dp, work(1, 1), ldwork)
1030 IF (iatom <= jatom) THEN
1031 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1032 1.0_dp, sphi_a(1, sgfa), ldsa, &
1033 work(1, 1), ldwork, &
1034 1.0_dp, op_dip(3)%block(sgfa, sgfb), &
1035 SIZE(op_dip(3)%block, 1))
1036 ELSE
1037 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1038 1.0_dp, work(1, 1), ldwork, &
1039 sphi_a(1, sgfa), ldsa, &
1040 1.0_dp, op_dip(3)%block(sgfb, sgfa), &
1041 SIZE(op_dip(3)%block, 1))
1042 END IF
1043 END IF ! >= dab
1044
1045 END DO ! jset
1046
1047 END DO ! iset
1048
1049 END DO
1050 CALL neighbor_list_iterator_release(nl_iterator)
1051
1052 DO i = 1, 3
1053 NULLIFY (op_dip(i)%block)
1054 END DO
1055 DEALLOCATE (op_dip)
1056
1057 DEALLOCATE (difab, work, basis_set_list)
1058
1059 CALL timestop(handle)
1060
1061 END SUBROUTINE p_xyz_ao
1062
1063! **************************************************************************************************
1064!> \brief Calculation of the components of the dipole operator in the length form
1065!> by taking the relative position operator r-Rc, with respect a reference point Rc
1066!> Probably it does not work for PBC, or maybe yes if the wfn are
1067!> sufficiently localized
1068!> The elements of the sparse matrices are the integrals in the
1069!> basis functions
1070!> \param op matrix representation of the p operator
1071!> calculated in terms of the contracted basis functions
1072!> \param qs_env environment for the lists and the basis sets
1073!> \param rc reference vector position
1074!> \param order maximum order of the momentum, for the dipole order = 1, order = -2 for quad only
1075!> \param minimum_image take into account only the first neighbors in the lists
1076!> \param soft ...
1077!> \par History
1078!> 03.2006 created [MI]
1079!> 06.2019 added quarupole only option (A.Bussy)
1080!> \author MI
1081! **************************************************************************************************
1082
1083 SUBROUTINE rrc_xyz_ao(op, qs_env, rc, order, minimum_image, soft)
1084
1085 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op
1086 TYPE(qs_environment_type), POINTER :: qs_env
1087 REAL(dp) :: rc(3)
1088 INTEGER, INTENT(IN) :: order
1089 LOGICAL, INTENT(IN), OPTIONAL :: minimum_image, soft
1090
1091 CHARACTER(len=*), PARAMETER :: routinen = 'rRc_xyz_ao'
1092
1093 CHARACTER(LEN=default_string_length) :: basis_type
1094 INTEGER :: handle, iatom, icol, ikind, imom, inode, irow, iset, jatom, jkind, jset, &
1095 last_jatom, ldab, ldsa, ldsb, ldwork, m_dim, maxl, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
1096 sgfb, smom
1097 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, npgfa, npgfb, &
1098 nsgfa, nsgfb
1099 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1100 LOGICAL :: found, my_minimum_image, my_soft, &
1101 new_atom_b
1102 REAL(kind=dp) :: dab, lxo2, lyo2, lzo2, rab2
1103 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
1104 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1105 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
1106 zeta, zetb
1107 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab
1108 TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip
1109 TYPE(cell_type), POINTER :: cell
1110 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1111 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1112 TYPE(neighbor_list_iterator_p_type), &
1113 DIMENSION(:), POINTER :: nl_iterator
1114 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1115 POINTER :: sab_orb
1116 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1117 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1118 TYPE(qs_kind_type), POINTER :: qs_kind
1119
1120 CALL timeset(routinen, handle)
1121
1122 NULLIFY (qs_kind, qs_kind_set)
1123 NULLIFY (cell, particle_set)
1124 NULLIFY (sab_orb)
1125 NULLIFY (mab, op_dip, work)
1126 NULLIFY (la_max, la_min, lb_max, npgfa, npgfb, nsgfa, nsgfb)
1127 NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
1128
1129 my_soft = .false.
1130 IF (PRESENT(soft)) my_soft = soft
1131 IF (my_soft) THEN
1132 basis_type = "ORB_SOFT"
1133 ELSE
1134 basis_type = "ORB"
1135 END IF
1136
1137 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
1138 cell=cell, particle_set=particle_set, sab_orb=sab_orb)
1139
1140 nkind = SIZE(qs_kind_set)
1141
1142 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1143 maxco=ldwork, maxlgto=maxl)
1144
1145 my_minimum_image = .false.
1146 IF (PRESENT(minimum_image)) THEN
1147 my_minimum_image = minimum_image
1148 lxo2 = sqrt(sum(cell%hmat(:, 1)**2))/2.0_dp
1149 lyo2 = sqrt(sum(cell%hmat(:, 2)**2))/2.0_dp
1150 lzo2 = sqrt(sum(cell%hmat(:, 3)**2))/2.0_dp
1151 END IF
1152
1153 ldab = ldwork
1154
1155 smom = 1
1156 IF (order == -2) smom = 4
1157 m_dim = ncoset(abs(order)) - 1
1158 cpassert(m_dim <= SIZE(op, 1))
1159
1160 ALLOCATE (mab(ldab, ldab, 1:m_dim))
1161 mab(1:ldab, 1:ldab, 1:m_dim) = 0.0_dp
1162 ALLOCATE (work(ldwork, ldwork))
1163 work(1:ldwork, 1:ldwork) = 0.0_dp
1164 ALLOCATE (op_dip(smom:m_dim))
1165
1166 DO imom = smom, m_dim
1167 NULLIFY (op_dip(imom)%block)
1168 END DO
1169
1170 ALLOCATE (basis_set_list(nkind))
1171 DO ikind = 1, nkind
1172 qs_kind => qs_kind_set(ikind)
1173 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1174 IF (ASSOCIATED(basis_set_a)) THEN
1175 basis_set_list(ikind)%gto_basis_set => basis_set_a
1176 ELSE
1177 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1178 END IF
1179 END DO
1180 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1181 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1182 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1183 iatom=iatom, jatom=jatom, r=rab)
1184 basis_set_a => basis_set_list(ikind)%gto_basis_set
1185 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1186 basis_set_b => basis_set_list(jkind)%gto_basis_set
1187 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1188 ra = pbc(particle_set(iatom)%r, cell)
1189 ! basis ikind
1190 first_sgfa => basis_set_a%first_sgf
1191 la_max => basis_set_a%lmax
1192 la_min => basis_set_a%lmin
1193 npgfa => basis_set_a%npgf
1194 nseta = basis_set_a%nset
1195 nsgfa => basis_set_a%nsgf_set
1196 rpgfa => basis_set_a%pgf_radius
1197 set_radius_a => basis_set_a%set_radius
1198 sphi_a => basis_set_a%sphi
1199 zeta => basis_set_a%zet
1200 ! basis jkind
1201 first_sgfb => basis_set_b%first_sgf
1202 lb_max => basis_set_b%lmax
1203 npgfb => basis_set_b%npgf
1204 nsetb = basis_set_b%nset
1205 nsgfb => basis_set_b%nsgf_set
1206 rpgfb => basis_set_b%pgf_radius
1207 set_radius_b => basis_set_b%set_radius
1208 sphi_b => basis_set_b%sphi
1209 zetb => basis_set_b%zet
1210
1211 ldsa = SIZE(sphi_a, 1)
1212 ldsb = SIZE(sphi_b, 1)
1213 IF (inode == 1) last_jatom = 0
1214
1215 IF (my_minimum_image) THEN
1216 IF (abs(rab(1)) > lxo2 .OR. abs(rab(2)) > lyo2 .OR. abs(rab(3)) > lzo2) cycle
1217 END IF
1218
1219 rb = rab + ra
1220
1221 IF (jatom /= last_jatom) THEN
1222 new_atom_b = .true.
1223 last_jatom = jatom
1224 ELSE
1225 new_atom_b = .false.
1226 END IF
1227
1228 IF (new_atom_b) THEN
1229 IF (iatom <= jatom) THEN
1230 irow = iatom
1231 icol = jatom
1232 ELSE
1233 irow = jatom
1234 icol = iatom
1235 END IF
1236
1237 DO imom = smom, m_dim
1238 NULLIFY (op_dip(imom)%block)
1239 CALL dbcsr_get_block_p(matrix=op(imom)%matrix, &
1240 row=irow, col=icol, block=op_dip(imom)%block, found=found)
1241 cpassert(ASSOCIATED(op_dip(imom)%block))
1242 END DO ! imom
1243 END IF ! new_atom_b
1244
1245 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1246 dab = sqrt(rab2)
1247
1248 DO iset = 1, nseta
1249
1250 ncoa = npgfa(iset)*ncoset(la_max(iset))
1251 sgfa = first_sgfa(1, iset)
1252
1253 DO jset = 1, nsetb
1254
1255 ncob = npgfb(jset)*ncoset(lb_max(jset))
1256 sgfb = first_sgfb(1, jset)
1257
1258 IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
1259
1260 rac = pbc(rc, ra, cell)
1261 rbc = pbc(rc, rb, cell)
1262
1263! *** Calculate the primitive overlap integrals ***
1264 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
1265 rpgfa(:, iset), la_min(iset), &
1266 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), &
1267 abs(order), rac, rbc, mab)
1268
1269 DO imom = smom, m_dim
1270! *** Contraction ***
1271 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1272 1.0_dp, mab(1, 1, imom), ldab, sphi_b(1, sgfb), ldsb, &
1273 0.0_dp, work(1, 1), ldwork)
1274 IF (iatom <= jatom) THEN
1275 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1276 1.0_dp, sphi_a(1, sgfa), ldsa, &
1277 work(1, 1), ldwork, &
1278 1.0_dp, op_dip(imom)%block(sgfa, sgfb), &
1279 SIZE(op_dip(imom)%block, 1))
1280 ELSE
1281 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1282 1.0_dp, work(1, 1), ldwork, &
1283 sphi_a(1, sgfa), ldsa, &
1284 1.0_dp, op_dip(imom)%block(sgfb, sgfa), &
1285 SIZE(op_dip(imom)%block, 1))
1286 END IF
1287
1288 END DO ! imom
1289 END IF ! >= dab
1290
1291 END DO ! jset
1292
1293 END DO ! iset
1294
1295 END DO
1296 CALL neighbor_list_iterator_release(nl_iterator)
1297
1298 DO imom = smom, m_dim
1299 NULLIFY (op_dip(imom)%block)
1300 END DO
1301 DEALLOCATE (op_dip)
1302
1303 DEALLOCATE (mab, work, basis_set_list)
1304
1305 CALL timestop(handle)
1306
1307 END SUBROUTINE rrc_xyz_ao
1308
1309! **************************************************************************************************
1310!> \brief Calculation of the multipole operators integrals
1311!> and of its derivatives of the type
1312!> [\mu | op | d(\nu)/dR(\nu)]-[d(\mu)/dR(\mu)| op | \nu]
1313!> by taking the relative position operator r-Rc, with respect a reference point Rc
1314!> The derivative are with respect to the primitive position,
1315!> The multipole operator is symmetric and if it does not depend on R(\mu) or R(\nu)
1316!> therefore [\mu | op | d(\nu)/dR(\nu)] = -[d(\mu)/dR(\mu)| op | \nu]
1317!> [\mu|op|d(\nu)/dR]-[d(\mu)/dR|op|\nu]=2[\mu|op|d(\nu)/dR]
1318!> When it is not the case a correction term is needed
1319!>
1320!> The momentum operator [\mu|M|\nu] is symmetric, the number of components is
1321!> determined by the order: 3 for order 1 (x,y,x), 9 for order 2(xx,xy,xz,yy,yz,zz)
1322!> The derivative of the type [\mu | op | d(\nu)/dR_i(\nu)], where
1323!> i indicates the cartesian direction, is antisymmetric only when
1324!> the no component M =(r_i) or (r_i r_j) is in the same cartesian
1325!> direction of the derivative, indeed
1326!> d([\mu|M|\nu])/dr_i = [d(\mu)/dr_i|M|\nu] + [\mu|M|d(\nu)/dr_i] + [\mu |d(M)/dr_i|\nu]
1327!> d([\mu|M|\nu])/dr_i = -[d(\mu)/dR_i(\mu)|M|\nu] -[\mu|M|d(\nu)/dR_i(\nu)] + [\mu |d(M)/dr_i|\nu]
1328!> Therefore we cannot use an antisymmetric matrix
1329!>
1330!> The same holds for the derivative with respect to the electronic position r
1331!> taking into account that [\mu|op|d(\nu)/dR] = -[\mu|op|d(\nu)/dr]
1332!> \param op matrix representation of the p operator
1333!> calculated in terms of the contracted basis functions
1334!> \param op_der ...
1335!> \param qs_env environment for the lists and the basis sets
1336!> \param rc reference vector position
1337!> \param order maximum order of the momentum, for the dipole order = 1
1338!> \param minimum_image take into account only the first neighbors in the lists
1339!> \param soft ...
1340!> \par History
1341!> 03.2006 created [MI]
1342!> \author MI
1343!> \note
1344!> Probably it does not work for PBC, or maybe yes if the wfn are
1345!> sufficiently localized
1346!> The elements of the sparse matrices are the integrals in the
1347!> basis functions
1348! **************************************************************************************************
1349 SUBROUTINE rrc_xyz_der_ao(op, op_der, qs_env, rc, order, minimum_image, soft)
1350
1351 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op
1352 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_der
1353 TYPE(qs_environment_type), POINTER :: qs_env
1354 REAL(dp) :: rc(3)
1355 INTEGER, INTENT(IN) :: order
1356 LOGICAL, INTENT(IN), OPTIONAL :: minimum_image, soft
1357
1358 CHARACTER(len=*), PARAMETER :: routinen = 'rRc_xyz_der_ao'
1359
1360 CHARACTER(LEN=default_string_length) :: basis_type
1361 INTEGER :: handle, i, iatom, icol, idir, ikind, imom, inode, ipgf, irow, iset, j, jatom, &
1362 jkind, jpgf, jset, last_jatom, lda_min, ldab, ldb_min, ldsa, ldsb, ldwork, m_dim, maxl, &
1363 na, nb, ncoa, ncob, nda, ndb, nkind, nseta, nsetb, sgfa, sgfb
1364 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
1365 npgfb, nsgfa, nsgfb
1366 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
1367 LOGICAL :: my_minimum_image, my_soft, new_atom_b, &
1368 op_der_found, op_found
1369 REAL(kind=dp) :: alpha, alpha_der, dab, lxo2, lyo2, lzo2, &
1370 rab2
1371 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
1372 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
1373 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, work, &
1374 zeta, zetb
1375 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab, mab_tmp
1376 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: difmab
1377 TYPE(block_p_type), DIMENSION(:), POINTER :: op_dip
1378 TYPE(block_p_type), DIMENSION(:, :), POINTER :: op_dip_der
1379 TYPE(cell_type), POINTER :: cell
1380 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1381 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1382 TYPE(neighbor_list_iterator_p_type), &
1383 DIMENSION(:), POINTER :: nl_iterator
1384 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1385 POINTER :: sab_all
1386 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1387 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1388 TYPE(qs_kind_type), POINTER :: qs_kind
1389
1390 CALL timeset(routinen, handle)
1391
1392 cpassert(ASSOCIATED(op))
1393 cpassert(ASSOCIATED(op_der))
1394 !IF(.NOT.op_sm_der(1,1)%matrix%symmetry=="none") THEN
1395 IF (.NOT. dbcsr_get_matrix_type(op_der(1, 1)%matrix) .EQ. dbcsr_type_no_symmetry) THEN
1396 cpabort("")
1397 END IF
1398
1399 NULLIFY (qs_kind, qs_kind_set)
1400 NULLIFY (cell, particle_set)
1401 NULLIFY (sab_all)
1402 NULLIFY (difmab, mab, mab_tmp)
1403 NULLIFY (op_dip, op_dip_der, work)
1404 NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb)
1405 NULLIFY (set_radius_a, set_radius_b, rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb)
1406
1407 my_soft = .false.
1408 IF (PRESENT(soft)) my_soft = soft
1409 IF (my_soft) THEN
1410 basis_type = "ORB_SOFT"
1411 ELSE
1412 basis_type = "ORB"
1413 END IF
1414
1415 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
1416 cell=cell, particle_set=particle_set, &
1417 sab_all=sab_all)
1418
1419 nkind = SIZE(qs_kind_set)
1420
1421 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1422 maxco=ldwork, maxlgto=maxl)
1423
1424 my_minimum_image = .false.
1425 IF (PRESENT(minimum_image)) THEN
1426 my_minimum_image = minimum_image
1427 lxo2 = sqrt(sum(cell%hmat(:, 1)**2))/2.0_dp
1428 lyo2 = sqrt(sum(cell%hmat(:, 2)**2))/2.0_dp
1429 lzo2 = sqrt(sum(cell%hmat(:, 3)**2))/2.0_dp
1430 END IF
1431
1432 ldab = ldwork
1433
1434 m_dim = ncoset(order) - 1
1435 cpassert(m_dim <= SIZE(op, 1))
1436
1437 ALLOCATE (mab(ldab, ldab, m_dim))
1438 mab(1:ldab, 1:ldab, 1:m_dim) = 0.0_dp
1439 ALLOCATE (difmab(ldab, ldab, m_dim, 3))
1440 difmab(1:ldab, 1:ldab, 1:m_dim, 1:3) = 0.0_dp
1441
1442 ALLOCATE (work(ldwork, ldwork))
1443 work(1:ldwork, 1:ldwork) = 0.0_dp
1444 ALLOCATE (op_dip(m_dim))
1445 ALLOCATE (op_dip_der(m_dim, 3))
1446
1447 DO imom = 1, m_dim
1448 NULLIFY (op_dip(imom)%block)
1449 DO i = 1, 3
1450 NULLIFY (op_dip_der(imom, i)%block)
1451 END DO
1452 END DO
1453
1454 ALLOCATE (basis_set_list(nkind))
1455 DO ikind = 1, nkind
1456 qs_kind => qs_kind_set(ikind)
1457 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1458 IF (ASSOCIATED(basis_set_a)) THEN
1459 basis_set_list(ikind)%gto_basis_set => basis_set_a
1460 ELSE
1461 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1462 END IF
1463 END DO
1464 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
1465 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1466 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1467 iatom=iatom, jatom=jatom, r=rab)
1468 basis_set_a => basis_set_list(ikind)%gto_basis_set
1469 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1470 basis_set_b => basis_set_list(jkind)%gto_basis_set
1471 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1472 ra = pbc(particle_set(iatom)%r, cell)
1473 ! basis ikind
1474 first_sgfa => basis_set_a%first_sgf
1475 la_max => basis_set_a%lmax
1476 la_min => basis_set_a%lmin
1477 npgfa => basis_set_a%npgf
1478 nseta = basis_set_a%nset
1479 nsgfa => basis_set_a%nsgf_set
1480 rpgfa => basis_set_a%pgf_radius
1481 set_radius_a => basis_set_a%set_radius
1482 sphi_a => basis_set_a%sphi
1483 zeta => basis_set_a%zet
1484 ! basis jkind
1485 first_sgfb => basis_set_b%first_sgf
1486 lb_max => basis_set_b%lmax
1487 lb_min => basis_set_b%lmin
1488 npgfb => basis_set_b%npgf
1489 nsetb = basis_set_b%nset
1490 nsgfb => basis_set_b%nsgf_set
1491 rpgfb => basis_set_b%pgf_radius
1492 set_radius_b => basis_set_b%set_radius
1493 sphi_b => basis_set_b%sphi
1494 zetb => basis_set_b%zet
1495
1496 ldsa = SIZE(sphi_a, 1)
1497 IF (ldsa .EQ. 0) cycle
1498 ldsb = SIZE(sphi_b, 1)
1499 IF (ldsb .EQ. 0) cycle
1500 IF (inode == 1) last_jatom = 0
1501
1502 IF (my_minimum_image) THEN
1503 IF (abs(rab(1)) > lxo2 .OR. abs(rab(2)) > lyo2 .OR. abs(rab(3)) > lzo2) cycle
1504 END IF
1505
1506 rb = rab + ra
1507
1508 IF (jatom /= last_jatom) THEN
1509 new_atom_b = .true.
1510 last_jatom = jatom
1511 ELSE
1512 new_atom_b = .false.
1513 END IF
1514
1515 IF (new_atom_b) THEN
1516 irow = iatom
1517 icol = jatom
1518 alpha_der = 2.0_dp
1519
1520 DO imom = 1, m_dim
1521 NULLIFY (op_dip(imom)%block)
1522 CALL dbcsr_get_block_p(matrix=op(imom)%matrix, &
1523 row=irow, col=icol, &
1524 block=op_dip(imom)%block, &
1525 found=op_found)
1526 cpassert(op_found .AND. ASSOCIATED(op_dip(imom)%block))
1527 DO idir = 1, 3
1528 NULLIFY (op_dip_der(imom, idir)%block)
1529 CALL dbcsr_get_block_p(matrix=op_der(imom, idir)%matrix, &
1530 row=irow, col=icol, &
1531 block=op_dip_der(imom, idir)%block, &
1532 found=op_der_found)
1533 cpassert(op_der_found .AND. ASSOCIATED(op_dip_der(imom, idir)%block))
1534 END DO ! idir
1535 END DO ! imom
1536 END IF ! new_atom_b
1537
1538 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1539 dab = sqrt(rab2)
1540
1541 DO iset = 1, nseta
1542
1543 ncoa = npgfa(iset)*ncoset(la_max(iset))
1544 sgfa = first_sgfa(1, iset)
1545
1546 DO jset = 1, nsetb
1547
1548 ncob = npgfb(jset)*ncoset(lb_max(jset))
1549 sgfb = first_sgfb(1, jset)
1550
1551 IF (set_radius_a(iset) + set_radius_b(jset) >= dab) THEN
1552
1553 rac = pbc(rc, ra, cell)
1554 rbc = rac + rab
1555! rac = pbc(rc,ra,cell)
1556! rbc = pbc(rc,rb,cell)
1557
1558 ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
1559 npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(order) - 1))
1560
1561 lda_min = max(0, la_min(iset) - 1)
1562 ldb_min = max(0, lb_min(jset) - 1)
1563! *** Calculate the primitive overlap integrals ***
1564 CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
1565 rpgfa(:, iset), lda_min, &
1566 lb_max(jset) + 1, npgfb(jset), zetb(:, jset), rpgfb(:, jset), &
1567 order, rac, rbc, mab_tmp)
1568
1569! *** Calculate the derivatives
1570 CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
1571 rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
1572 zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
1573 difmab, mab_ext=mab_tmp)
1574
1575! Contract and copy in the sparse matrix
1576 mab = 0.0_dp
1577 DO imom = 1, m_dim
1578 na = 0
1579 nda = 0
1580 DO ipgf = 1, npgfa(iset)
1581 nb = 0
1582 ndb = 0
1583 DO jpgf = 1, npgfb(jset)
1584 DO j = 1, ncoset(lb_max(jset))
1585 DO i = 1, ncoset(la_max(iset))
1586 mab(i + na, j + nb, imom) = mab_tmp(i + nda, j + ndb, imom)
1587 END DO ! i
1588 END DO ! j
1589 nb = nb + ncoset(lb_max(jset))
1590 ndb = ndb + ncoset(lb_max(jset) + 1)
1591 END DO ! jpgf
1592 na = na + ncoset(la_max(iset))
1593 nda = nda + ncoset(la_max(iset) + 1)
1594 END DO ! ipgf
1595
1596! *** Contraction ***
1597 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1598 1.0_dp, mab(1, 1, imom), ldab, sphi_b(1, sgfb), ldsb, &
1599 0.0_dp, work(1, 1), ldwork)
1600 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1601 1.0_dp, sphi_a(1, sgfa), ldsa, &
1602 work(1, 1), ldwork, &
1603 1.0_dp, op_dip(imom)%block(sgfa, sgfb), &
1604 SIZE(op_dip(imom)%block, 1))
1605
1606 alpha = -1.0_dp !-alpha_der
1607 DO idir = 1, 3
1608 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1609 alpha, difmab(1, 1, imom, idir), ldab, sphi_b(1, sgfb), ldsb, &
1610 0.0_dp, work(1, 1), ldwork)
1611 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1612 1.0_dp, sphi_a(1, sgfa), ldsa, &
1613 work(1, 1), ldwork, &
1614 1.0_dp, op_dip_der(imom, idir)%block(sgfa, sgfb), &
1615 SIZE(op_dip_der(imom, idir)%block, 1))
1616
1617 END DO ! idir
1618
1619 END DO ! imom
1620
1621 DEALLOCATE (mab_tmp)
1622 END IF ! >= dab
1623
1624 END DO ! jset
1625
1626 END DO ! iset
1627
1628 END DO
1629 CALL neighbor_list_iterator_release(nl_iterator)
1630
1631 DO i = 1, 3
1632 NULLIFY (op_dip(i)%block)
1633 END DO
1634 DEALLOCATE (op_dip, op_dip_der)
1635
1636 DEALLOCATE (mab, difmab, work, basis_set_list)
1637
1638 CALL timestop(handle)
1639
1640 END SUBROUTINE rrc_xyz_der_ao
1641
1642END MODULE qs_operators_ao
1643
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:206
cell_types::pbc
Definition cell_types.F:92
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition ai_moments.F:17
ai_moments::diff_momop
subroutine, public diff_momop(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, mab_ext)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the position of the pr...
Definition ai_moments.F:1583
ai_moments::diffop
subroutine, public diffop(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, rab, difab)
This returns the derivative of the overlap integrals [a|b], with respect to the position of the primi...
Definition ai_moments.F:1518
ai_moments::moment
subroutine, public moment(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab)
...
Definition ai_moments.F:986
ai_os_rr
Definition ai_os_rr.F:8
ai_os_rr::os_rr_ovlp
subroutine, public os_rr_ovlp(rap, la_max, rbp, lb_max, zet, ldrr, rr)
Calculation of the basic Obara-Saika recurrence relation.
Definition ai_os_rr.F:39
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:914
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:847
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:679
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
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
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
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, 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, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, 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, 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, init_u_ramping_each_scf, reltmat, ghost, 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:454
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, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:880
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:549
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:444
qs_operators_ao::rrc_xyz_der_ao
subroutine, public rrc_xyz_der_ao(op, op_der, qs_env, rc, order, minimum_image, soft)
Calculation of the multipole operators integrals and of its derivatives of the type [\mu | op | d(\nu...
Definition qs_operators_ao.F:1350
qs_operators_ao::p_xyz_ao
subroutine, public p_xyz_ao(op, qs_env, minimum_image)
Calculation of the components of the dipole operator in the velocity form The elements of the sparse ...
Definition qs_operators_ao.F:818
qs_operators_ao::build_lin_mom_matrix
subroutine, public build_lin_mom_matrix(qs_env, matrix)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
Definition qs_operators_ao.F:74
qs_operators_ao::rrc_xyz_ao
subroutine, public rrc_xyz_ao(op, qs_env, rc, order, minimum_image, soft)
Calculation of the components of the dipole operator in the length form by taking the relative positi...
Definition qs_operators_ao.F:1084
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:55
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
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