(git:e7e05ae)
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-2024 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 ! **************************************************************************************************
13 MODULE 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_log_handling, ONLY: cp_get_default_logger,&
24  cp_logger_type
25  USE dbcsr_api, ONLY: dbcsr_get_block_p,&
26  dbcsr_get_matrix_type,&
27  dbcsr_has_symmetry,&
28  dbcsr_p_type,&
29  dbcsr_type_antisymmetric,&
30  dbcsr_type_no_symmetry
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 
62 CONTAINS
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 
1642 END MODULE qs_operators_ao
1643 
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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:195
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_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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
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, 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_r3d_rs_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_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
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)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
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